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 solve the system with lapack::pbsv which is the FLENS interface for LAPACK's dpbsv.

Note that we might rename lapack::pbsv to lapack::posv.

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::pbsv(A, b);

    cout << "(L\\R) = " << A << endl;

    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;

And solve it with lapack::pbsv

    lapack::pbsv(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-pbsv lapack-pbsv.cc                                      

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-pbsv                                                          
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\R) = 
   1.414214  -0.707107                      
  -0.707107   1.224745  -0.816497           
             -0.816497   1.154701  -0.866025
                        -0.866025   1.118034
x = 
     4.000000       7.000000       8.000000       6.000000