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

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