Home Browse Files

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

flens/examples/lapack-pbsv.cc
#include <cxxstd/iostream.h>
#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. Note that currently the method diag(d) is only provided for matrices of type GbMatrix not SbMatrix. So we access the method via a GbMatrix view. If frequently needed we can provide the method also through the SbMatrix interface directly.

    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> clang++ -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
Home