# 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

#include <iostream>
#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> g++ -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