Solving a Symmetric Positive Definite System of Linear Equations

In this example we solve a system of linear equations \(Sx = b\) were the coefficient matrix is symmetric and positiv definite. Basically we do the same as in example lapack-posv. However, in this example we explicitly use lapack::trs which is the FLENS port of LAPACK's dtrtrs.

Example Code

flens/examples/lapack-trtrs.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    GeMatrix<FullStorage<double> >   A(3, 3);

    DenseVector<Array<double> >      b(3);

    A = 2, 0, 0,

        1, 2, 0,

        0, 1, 2;

    b = 1,

        2,

        3;

    auto S = A.lower().symmetric();

    cout << "S = " << S << endl;

    cout << "b = " << b << endl;

    int info = lapack::potrf(S);

    if (info!=0) {

        cerr << "S is singular" << endl;

        return info;

    }

    lapack::trs(NoTrans, A.lower(), b);

    lapack::trs(Trans, A.lower(), b);

    cout << "x = inv(S)*b = " << b << endl;

    return 0;

}

Comments on Example Code

Setup the raw data.

    A = 2, 0, 0,

        1, 2, 0,

        0, 1, 2;

\(S\) is a symmetric matrix view constructed from the lower triangular part of \(A\). Note that \(S\) only references data from \(A\).

auto S = A.lower().symmetric();

Computes the Cholesky factorization \(S = L L^T\) where \(L\) is lower triangular. Matrix \(S\) (i.e. the lower triangular part of \(A\)) gets overwritten with \(L\).

int info = lapack::potrf(S);

Solve \(L u = b\)

lapack::trs(NoTrans, A.lower(), b);

Solve \(L^T x = u\)

lapack::trs(Trans, A.lower(), b);

Compile

$shell> cd flens/examples                                                      
$shell> g++ -std=c++11 -Wall -I../.. -o lapack-trtrs lapack-trtrs.cc

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-trtrs                                                         
S = 
            2             1             0 
            1             2             1 
            0             1             2 
b = 
            1              2              3 
x = inv(S)*b = 
          0.5              0            1.5

Example with Complex Numbers

Example Code

flens/examples/lapack-complex-trtrs.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef complex<double>             Complex;

    const Complex                       I(0,1);

    GeMatrix<FullStorage<Complex> >     A(3, 3);

    DenseVector<Array<Complex> >        b(3);

    A = 2, I,    0,

        0, 2, 1.+I,

        0, 0,    3;

    b = 1,

        2,

        3;

    auto H = A.upper().hermitian();

    cout << "H = " << H << endl;

    cout << "b = " << b << endl;

    int info = lapack::potrf(H);

    if (info!=0) {

        cerr << "H is singular" << endl;

        return info;

    }

    lapack::trs(NoTrans, A.lower(), b);

    lapack::trs(Trans, A.lower(), b);

    cout << "x = inv(H)*b = " << b << endl;

    return 0;

}

Comments on Example Code

Setup the raw data.

    A = 2, I,    0,

        0, 2, 1.+I,

        0, 0,    3;

\(H\) is a symmetric matrix view constructed from the lower triangular part of \(A\). Note that \(H\) only references data from \(A\).

auto H = A.upper().hermitian();

Computes the Cholesky factorization \(H = L L^T\) where \(L\) is lower triangular. Matrix \(H\) (i.e. the lower triangular part of \(A\)) gets overwritten with \(L\).

int info = lapack::potrf(H);

Solve \(L u = b\)

lapack::trs(NoTrans, A.lower(), b);

Solve \(L^T x = u\)

lapack::trs(Trans, A.lower(), b);

Compile

Note that an external LAPACK implementation is needed for the complex variant of lapack::trs

$shell> cd flens/examples                                                       
$shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-trtrs lapack-complex-trtrs.cc

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-trtrs                                                  
H = 
                       (2,0)                       (0,1)                       (0,0)
                      (0,-1)                       (2,0)                       (1,1)
                      (0,-0)                      (1,-1)                       (3,0)
b = 
                       (1,0)                         (2,0)                         (3,0) 
x = inv(H)*b = 
                     (0.5,0)                   (1.33333,0)                       (1.8,0)