Computing the QR Factorization

In this example we compute the \(QR\) factorization and use it for solving a system of linear equations.

Consider the system of linear equations \(Ax=b\). We store the coefficient matrix \(A\) and the right-hand side \(b\) in a single matrix \(Ab = (A,b)\). We then compute the \(QR\) factorization of \(Ab\) with lapack::qrf such that \(QR = (A,b).\) Hence, we have \(R = (Q^T A, Q^T b)\) where \(R\) is an upper trapezoidal matrix and \(Q^T A\) upper triangular. We finally use the triangular solver blas::sm to obtain the solution of \((Q^T A) x = (Q^T b)\).

Function lapack::qrf is the FLENS port of LAPACK's geqrf and blas::sm the port of the BLAS routine trsm.

Example Code

flens/examples/lapack-geqrf.cc

#include <iostream>

#include <flens/flens.cxx>

using namespace std;

using namespace flens;

typedef double   T;

int

main()

{

    typedef GeMatrix<FullStorage<T> >   GeMatrix;

    typedef DenseVector<Array<T> >      DenseVector;

    typedef GeMatrix::IndexType   IndexType;

    const Underscore<IndexType>   _;

    const IndexType m = 4,

                    n = 5;

    GeMatrix  Ab(m, n);

    Ab =  2,   3,  -1,   0,   20,

         -6,  -5,   0,   2,  -33,

          2,  -5,   6,  -6,  -43,

          4,   6,   2,  -3,   49;

    cout << "Ab = " << Ab << endl;

    DenseVector tau(min(m,n));

    DenseVector work;

    lapack::qrf(Ab, tau, work);

    cout << "Ab = " << Ab << endl;

    cout << "tau = " << tau << endl;

    const auto A = Ab(_,_(1,m));

    auto       B = Ab(_,_(m+1,n));

    blas::sm(Left, NoTrans, 1, A.upper(), B);

    cout << "X = " << B << endl;

}

Comments on Example Code

#include <iostream>

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

using namespace std;

using namespace flens;

typedef double   T;

int

main()

{

Define some convenient typedefs for the matrix/vector types

    typedef GeMatrix<FullStorage<T> >   GeMatrix;

    typedef DenseVector<Array<T> >      DenseVector;

Define an underscore operator for convenient matrix slicing

    typedef GeMatrix::IndexType   IndexType;

    const Underscore<IndexType>   _;

    const IndexType m = 4,

                    n = 5;

    GeMatrix  Ab(m, n);

    Ab =  2,   3,  -1,   0,   20,

         -6,  -5,   0,   2,  -33,

          2,  -5,   6,  -6,  -43,

          4,   6,   2,  -3,   49;

    cout << "Ab = " << Ab << endl;

tau will contain the scalar factors of the elementary reflectors (see dgeqrf for details). Vector work is used as workspace and, if empty, gets resized by lapack::qrf to optimal size.

    DenseVector tau(min(m,n));

    DenseVector work;

Compute the \(QR\) factorization of \(A\) with lapack::qrf the FLENS port of LAPACK's dgeqrf.

    lapack::qrf(Ab, tau, work);

    cout << "Ab = " << Ab << endl;

    cout << "tau = " << tau << endl;

Solve the system of linear equation \(Ax =B\) using the triangular solver blas::sm which is the FLENS implementation of the BLAS routine trsm. Note that A.upper() creates a triangular view.

    const auto A = Ab(_,_(1,m));

    auto       B = Ab(_,_(m+1,n));

    blas::sm(Left, NoTrans, 1, A.upper(), B);

    cout << "X = " << B << endl;

}

Compile

$shell> cd flens/examples                                                      
$shell> clang++ -std=c++0x -Wall -I../.. -o lapack-geqrf lapack-geqrf.cc

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-geqrf                                                         
Ab = 
                   2                    3                   -1                    0                   20 
                  -6                   -5                    0                    2                  -33 
                   2                   -5                    6                   -6                  -43 
                   4                    6                    2                   -3                   49 
Ab = 
            -7.74597             -6.45497             -2.32379              4.64758             -44.9266 
           -0.615639             -7.30297               4.9295             -4.38178             -60.7972 
            0.205213            -0.854313             -3.36155              2.85583             -4.55148 
            0.410426             0.260891             0.453851             0.210352              1.89316 
tau = 
     1.2582       1.1124       1.6584            0 
X = 
                   1 
                   9 
                   9 
                   9