@Override
    public void solve(DMatrixRMaj B, DMatrixRMaj X) {
        if( B.numRows != numRows )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        X.reshape(numCols,B.numCols);

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                a[i] = B.data[i*BnumCols + colB];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = n*numRows + n + 1;

                double ub = a[n];
                // U^T*b
                for( int i = n+1; i < numRows; i++ , indexU++ ) {
                    ub += dataQR[indexU]*a[i];
                }

                // gamma*U^T*b
                ub *= gammas[n];

                a[n] -= ub;
                indexU = n*numRows + n + 1;
                for( int i = n+1; i < numRows; i++ , indexU++) {
                    a[i] -= dataQR[indexU]*ub;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_DDRM.solveU(U.data,a,numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.data[i*X.numCols+colB] = a[i];
            }
        }
    }
@ Override public void solve ( DMatrixRMaj B , DMatrixRMaj X ) { if ( B . numRows != numRows ) throw new IllegalArgumentException ( "Unexpected dimensions for X: X rows = " + X . numRows + " expected = " + numCols ) ; X . reshape ( numCols , B . numCols ) ; U = decomposer . getR ( U , true ) ; final double gammas [ ] = decomposer . getGammas ( ) ; final double dataQR [ ] = QR . data ; final int BnumCols = B . numCols ; for ( int colB = 0 ; colB < BnumCols ; colB ++ ) { for ( int i = 0 ; i < numRows ; i ++ ) { a [ i ] = B . data [ i * BnumCols + colB ] ; } for ( int n = 0 ; n < numCols ; n ++ ) { int indexU = n * numRows + n + 1 ; double ub = a [ n ] ; for ( int i = n + 1 ; i < numRows ; i ++ , indexU ++ ) { ub += dataQR [ indexU ] * a [ i ] ; } ub *= gammas [ n ] ; a [ n ] -= ub ; indexU = n * numRows + n + 1 ; for ( int i = n + 1 ; i < numRows ; i ++ , indexU ++ ) { a [ i ] -= dataQR [ indexU ] * ub ; } } TriangularSolver_DDRM . solveU ( U . data , a , numCols ) ; for ( int i = 0 ; i < numCols ; i ++ ) { X . data [ i * X . numCols + colB ] = a [ i ] ; } } }
Solves for X using the QR decomposition.

@param B A matrix that is n by m.  Not modified.
@param X An n by m matrix where the solution is written to.  Modified.