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 general and banded. We solve the system with lapack::sv which is the FLENS interface for LAPACK's dgbsv.

If you compute the \(LU\) factorization of a band matrix with \(k_l\) sub-diagonals and \(k_u\) super-diagonals then additional \(k_l\) super-diagonals are needed. I.e. the factorization requires \(k_l\) sub-diagonals and \(k_u+k_l\) super-diagonals.

For this reason FLENS-LAPACK requires that this additional storage is allocated. So a tridiagonal matrix (\(k_l=k_u=1\)) must be stored in a band matrix with two super-diagonals:

\[ A = \begin{pmatrix} a_{1,1} & a_{1,2} & * & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & * & 0 \\ 0 & a_{3,2} & a_{3,3} & a_{3,4} & * \\ 0 & 0 & a_{4,3} & a_{4,4} & a_{4,5} \\ 0 & 0 & 0 & a_{5,4} & a_{5,5} \\ \end{pmatrix} \]

After calling lapack::sv the matrix is overwritten with the factorization:

\[ A_{LU} = \begin{pmatrix} u_{1,1} & u_{1,2} & u_{1,3} & 0 & 0 \\ m_{2,1} & u_{2,2} & u_{2,3} & u_{2,4} & 0 \\ 0 & m_{3,2} & u_{3,3} & u_{3,4} & u_{3,5} \\ 0 & 0 & m_{4,3} & u_{4,4} & u_{4,5} \\ 0 & 0 & 0 & m_{5,4} & u_{5,5} \\ \end{pmatrix} \\ \]

Example Code

flens/examples/lapack-gbsv.cc
#include <cxxstd/iostream.h>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{
    typedef GbMatrix<BandStorage<T> >           BandMatrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef BandMatrix::IndexType               IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    BandMatrix     A(n,n,1,2);

    Vector         b(n);
    IndexVector    piv(n);

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

    b =  4,
        -7,
         0,
         1;

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

    lapack::sv(A, piv, 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 GbMatrix<BandStorage<T> >           BandMatrix;
    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 BandMatrix::IndexType               IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

Set up the baby problem ...

    const IndexType n = 4;

We allocate a band matrix with one sub-diagonal and two superdiagonals. Note, if you compute the LU factorization of a tridiagonal matrix then an extra superdiagonal is needed!

    BandMatrix     A(n,n,1,2);

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

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

And solve it with lapack::sv

    lapack::sv(A, piv, 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-gbsv lapack-gbsv.cc

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-gbsv                                                          
A = 
   2.000000  -1.000000   0.000000           
  -3.000000   2.000000  -1.000000   0.000000
             -3.000000   2.000000  -1.000000
                        -3.000000   2.000000
b = 
     4.000000      -7.000000       0.000000       1.000000 
(L\R) = 
  -3.000000   2.000000  -1.000000           
  -0.666667  -3.000000   2.000000  -1.000000
             -0.111111  -3.000000   2.000000
                         0.148148  -0.407407
x = 
     2.000000      -0.000000       1.000000       2.000000
Home