# Solving a Symmetric Positive Definite System of Linear Equations

In this example we solve a system of linear equations $$Sx = b$$ were the coefficient matrix is symmetric and positiv definite. We first compute the Cholesky factorization $$S = L L^T$$ with lapack::potrf which is the FLENS port of LAPACK's dpotrf. We then solve $$Lu=b$$ and then $$L^T x = u$$ using lapack::potrs which is the port of dpotrs.

## Example Code

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

using namespace flens;
using namespace std;

int
main()
{
GeMatrix<FullStorage<double> >   A(33);
DenseVector<Array<double> >      b(3);

A = 210,
021,
002;

b = 1,
2,
3;

auto S = A.upper().symmetric();

cout << "S = " << S << endl;
cout << "b = " << b << endl;

int info = lapack::potrf(S);

if (info!=0) {
cerr << "S is singular" << endl;
return info;
}

lapack::potrs(S, b);
cout << "inv(S)*b = " << b << endl;

return 0;
}

Setup the raw data.

A = 210,
021,
002;

$$S$$ is a symmetric matrix view constructed from the upper triangular part of $$A$$. Note that $$S$$ only references data from $$A$$.

auto S = A.upper().symmetric();

Computes the Cholesky factorization $$S = L L^T$$ where $$L$$ is upper triangular. Matrix $$S$$ (i.e. the lower triangular part of $$A$$) gets overwritten with $$L^T$$.

int info = lapack::potrf(S);

If info is not zero the factorization could not be computed as the matrix was singular.

if (info!=0) {
cerr << "S is singular" << endl;
return info;
}

Compute $$x = S^{-1}b$$. Vector $$b$$ gets overwritten with the solution vector $$x$$.

lapack::potrs(S, b);
cout << "inv(S)*b = " << b << endl;

## Compile

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


## Run

$shell> cd flens/examples$shell> ./lapack-potrs
S =
2             1             0
1             2             1
0             1             2
b =
1              2              3
inv(S)*b =
0.5              0            1.5


## Example with Complex Numbers

### Example Code

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

using namespace flens;
using namespace std;

int
main()
{
typedef complex<double>             Complex;
const Complex                       I(0,1);

GeMatrix<FullStorage<Complex> >     A(33);
DenseVector<Array<Complex> >        b(3);

A = 2, I,    0,
021.+I,
00,    3;

b = 1,
2,
3;

auto H = A.upper().hermitian();

cout << "H = " << H << endl;
cout << "b = " << b << endl;

int info = lapack::potrf(H);

if (info!=0) {
cerr << "H is singular" << endl;
return info;
}

lapack::potrs(H, b);
cout << "inv(H)*b = " << b << endl;

return 0;
}

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

using namespace flens;
using namespace std;

int
main()
{
typedef complex<double>             Complex;
const Complex                       I(0,1);

GeMatrix<FullStorage<Complex> >     A(33);
DenseVector<Array<Complex> >        b(3);

Setup the raw data.

A = 2, I,    0,
021.+I,
00,    3;

b = 1,
2,
3;

$$H$$ is a symmetric matrix view constructed from the upper triangular part of $$A$$. Note that $$H$$ only references data from $$A$$.

auto H = A.upper().hermitian();

cout << "H = " << H << endl;
cout << "b = " << b << endl;

Computes the Cholesky factorization $$H = L L^T$$ where $$L$$ is upper triangular. Matrix $$H$$ (i.e. the lower triangular part of $$A$$) gets overwritten with $$L^T$$.

int info = lapack::potrf(H);

If info is not zero the factorization could not be computed as the matrix was singular.

if (info!=0) {
cerr << "H is singular" << endl;
return info;
}

Compute $$x = H^{-1}b$$. Vector $$b$$ gets overwritten with the solution vector $$x$$.

lapack::potrs(H, b);
cout << "inv(H)*b = " << b << endl;

return 0;
}

### Compile

Note that an external LAPACK implementation is needed for the complex variant of lapack::potrs

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


### Run

$shell> cd flens/examples$shell> ./lapack-complex-potrs
H =
(2,0)                       (0,1)                       (0,0)
(0,-1)                       (2,0)                       (1,1)
(0,-0)                      (1,-1)                       (3,0)
b =
(1,0)                         (2,0)                         (3,0)
inv(H)*b =
(0.2,-0.6)                    (1.2,-0.6)                     (0.8,0.6)