Content

Solving Systems of Linear Equations

In this example we solve a system of linear equations \(Ax = b\) were the coefficient matrix is symmetric positiv definite and banded. We first compute the Cholesky factorization \(S = L L^T\) with lapack::pbtrf which is the FLENS port of LAPACK's dpbtrf. We then solve \(Lu=b\) and then \(L^T x = u\) using lapack::pbtrs which is the port of dpbtrs.

Note that we might rename lapack::pbtrf to lapack::potrf and lapack::pbtrs to lapack::potrs.

Example Code

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

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{
    typedef SbMatrix<BandStorage<T> >           SymmetricBandMatrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef SymmetricBandMatrix::IndexType      IndexType;

    const IndexType n = 4;

    SymmetricBandMatrix     A(n, Upper, 1);

    Vector                  b(n);

    A.general().diag( 1) = -1;
    A.general().diag( 0) =  2;

    b =  1,
         2,
         3,
         4;

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

    lapack::pbtrf(A);

    cout << "L^T = " << A.general().upper() << endl;

    lapack::pbtrs(A, b);

    cout << "x = " << b << endl;
}

Comments on Example Code

Define some convenient typedefs for the matrix/vector types of our system of linear equations.

    typedef SbMatrix<BandStorage<T> >           SymmetricBandMatrix;
    typedef DenseVector<Array<T> >              Vector;

We also need an extra vector type for the pivots. The type of the pivots is taken for the system matrix.

    typedef SymmetricBandMatrix::IndexType      IndexType;

Set up the baby problem ...

    const IndexType n = 4;

We allocate a symmetric band matrix with one off-diagonal.

    SymmetricBandMatrix     A(n, Upper, 1);

We initialize the positiv definite \(n \times n\) tridiagonal matrix.

    A.general().diag( 1) = -1;
    A.general().diag( 0) =  2;

We factorize \(A\) with lapack::pbtrf

    lapack::pbtrf(A);

In the upper part of \(A\) now the triangular factor \(L^T\) is stored.

    cout << "L^T = " << A.general().upper() << endl;

Solve it with lapack::pbtrs

    lapack::pbtrs(A, b);

Compile

Note that we need to link against an external LAPACK implementation:

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

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-pbtrs                                                         
A = 
   2.000000  -1.000000                      
  -1.000000   2.000000  -1.000000           
             -1.000000   2.000000  -1.000000
                        -1.000000   2.000000
b = 
     1.000000       2.000000       3.000000       4.000000 
L^T = 
   1.414214   -0.707107                         
               1.224745   -0.816497             
                           1.154701   -0.866025 
                                       1.118034 
x = 
     4.000000       7.000000       8.000000       6.000000