# Solving Systems of Linear Equations

In this example we solve a system of linear equations $$Ax = b$$ were the coefficient matrix is general, i.e. not necessarily symmetric. We solve the system with lapack::sv which is the FLENS port of LAPACK's dgesv.

## Example Code

#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{
typedef GeMatrix<FullStorage<T> >           Matrix;
typedef DenseVector<Array<T> >              Vector;

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

const IndexType n = 4;

Matrix         A(n,n);
Vector         b(n);
IndexVector    piv(n);

A =  2,   3,  -1,   0,
-6,  -5,   0,   2,
2,  -5,   6,  -6,
4,   6,   2,  -3;

b = 20,
-33,
-43,
49;

cerr << "A = " << A << endl;
cerr << "b = " << b << endl;

lapack::sv(A, piv, b);

cerr << "x = " << b << endl;
}

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

Define some convenient typedefs for the matrix/vector types of our system of linear equations.

typedef GeMatrix<FullStorage<T> >           Matrix;
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 Matrix::IndexType                   IndexType;
typedef DenseVector<Array<IndexType> >      IndexVector;

Set up the baby problem ...

const IndexType n = 4;

And solve it with lapack::sv

lapack::sv(A, piv, b);

## Compile

$shell> cd flens/examples$shell> g++ -std=c++11 -Wall -I../.. -o lapack-gesv lapack-gesv.cc


## Run

$shell> cd flens/examples$shell> ./lapack-gesv
A =
2             3            -1             0
-6            -5             0             2
2            -5             6            -6
4             6             2            -3
b =
20            -33            -43             49
x =
1              9              9              9


## Example with Complex Numbers

### Example Code

#include <complex>
#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef complex<double>   T;

int
main()
{
typedef GeMatrix<FullStorage<T> >           Matrix;
typedef DenseVector<Array<T> >              Vector;
typedef Matrix::IndexType                   IndexType;
typedef DenseVector<Array<IndexType> >      IndexVector;

const IndexType n = 3;

Matrix         A(n,n);
Vector         b(n);
IndexVector    piv(n);

A = T(10), T( 1,-1), T(  2,20),
T(01), T( 12), T(-108),
T(0,-1), T(-11), T( 406);

b = T( 10),
T(-11),
T(-2,-1);

cerr << "A = " << A << endl;
cerr << "b = " << b << endl;

lapack::sv(A, piv, b);

cerr << "A = " << A << endl;
cerr << "piv = " << piv << endl;
cerr << "x = " << b << endl;
}

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

Elements are now defined to be complex

typedef complex<double>   T;

Like before we define convenient matrix/vector types ...

typedef GeMatrix<FullStorage<T> >           Matrix;
typedef DenseVector<Array<T> >              Vector;
typedef Matrix::IndexType                   IndexType;
typedef DenseVector<Array<IndexType> >      IndexVector;

Then we setup another toy problem

const IndexType n = 3;

And solve it with lapack::sv

lapack::sv(A, piv, b);

### Compile

$shell> cd flens/examples$shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-gesv lapack-complex-gesv.cc


### Run

$shell> cd flens/examples$shell> ./lapack-complex-gesv
A =
(1,0)                       (1,-1)                       (2,20)
(0,1)                        (1,2)                      (-10,8)
(0,-1)                       (-1,1)                       (40,6)
b =
(1,0)                        (-1,1)                       (-2,-1)
A =
(1,0)                       (1,-1)                       (2,20)
(0,-1)                        (0,2)                       (20,8)
(0,1)                      (0.5,0)                        (0,2)
piv =
1              3              3
x =
(0,-1)                         (0,1)                         (0,0)