# 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;
}

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