Content

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. We solve the system with lapack::posv which is the FLENS port of LAPACK's dposv.

Example Code

#include <flens/flens.cxx>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
    GeMatrix<FullStorage<double> >   A(33);
    DenseVector<Array<double> >      b(3);

    A = 200,
        120,
        012;

    b = 1,
        2,
        3;

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

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

    int info = lapack::posv(S, b);

    if (info!=0) {
        cerr << "S is singular" << endl;
        return info;
    }

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

    return 0;
}

Comments on Example Code

Setup the raw data.

    A = 200,
        120,
        012;

\(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();
    int info = lapack::posv(S, b);

If info is not zero the factorization could not be computed as the matrix was singular.

    if (info!=0) {
        cerr << "S is singular" << endl;
        return info;
    }

Compile

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

Run

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

Example with Complex Numbers

Example Code

#include <flens/flens.cxx>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
    typedef complex<double>             Complex;
    const Complex                       I(0,1);

    GeMatrix<FullStorage<Complex> >     A(33);
    DenseVector<Array<Complex> >        b(3);

    A = 2,    00,
        I,    20,
        01.+I, 3;

    b =    1,
        2.+I,
           3;

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

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

    int info = lapack::posv(H, b);

    if (info!=0) {
        cerr << "H is singular" << endl;
        return info;
    }

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

    return 0;
}

Comments on Example Code

Setup the raw data.

    A = 2,    00,
        I,    20,
        01.+I, 3;

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

    auto H = A.lower().hermitian();
    int info = lapack::posv(H, b);

If info is not zero the factorization could not be computed as the matrix was singular.

    if (info!=0) {
        cerr << "H is singular" << endl;
        return info;
    }

Compile

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

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-posv                                                   
H = 
                       (2,0)                      (0,-1)                      (0,-0)
                       (0,1)                       (2,0)                      (1,-1)
                       (0,0)                       (1,1)                       (3,0)
b = 
                       (1,0)                         (2,1)                         (3,0) 
inv(H)*b = 
                  (-0.4,0.6)                     (1.2,1.8)                      (1.2,-1)