Compute the Inverse of a Symmetric Positive Definite Matrix

In this example we compute the inverse $$S^{-1}$$ of a symmetric positive definite matrix $$S$$. We first compute the Cholesky factoriation $$S=L L^T$$. The inverse $$S^{-1}$$ then gets computed by lapack::potri which is the FLENS port of LAPACK's dpotri.

Example Code

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

using namespace flens;
using namespace std;

int
main()
{
GeMatrix<FullStorage<double> >   A(33);

A = 200,
120,
012;

auto S = A.lower().symmetric();
cout << "S = " << S << endl;

int info = lapack::potrf(S);

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

lapack::potri(S);
cout << "inv(S) = " << S << endl;

return 0;
}

Comments on Example Code

Setup the raw data.

A = 200,
120,
012;

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

auto S = A.lower().symmetric();
cout << "S = " << S << endl;

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

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 $$S^{-1}$$ from the Cholesky factorization $$S = L L^T$$.

lapack::potri(S);
cout << "inv(S) = " << S << endl;

Compile

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


Run

$shell> cd flens/examples$shell> ./lapack-potri
S =
2             1             0
1             2             1
0             1             2
inv(S) =
0.75          -0.5          0.25
-0.5             1          -0.5
0.25          -0.5          0.75


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);

A = 2,    00,
I,    20,
01.+I, 2;

auto H = A.lower().hermitian();
cout << "H = " << H << endl;

int info = lapack::potrf(H);

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

cout << "H = " << H << endl;
lapack::potri(H);
cout << "inv(H) = " << H << endl;

return 0;
}

Comments on Example Code

Setup the raw data.

A = 2,    00,
I,    20,
01.+I, 2;

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

auto H = A.lower().hermitian();
cout << "H = " << H << endl;

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

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 $$H^{-1}$$ from the Cholesky factorization $$H = L L^T$$.

cout << "H = " << H << endl;
lapack::potri(H);
cout << "inv(H) = " << H << endl;

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-potri lapack-complex-potri.cc


Run

$shell> cd flens/examples$shell> ./lapack-complex-potri
H =
(2,0)                      (0,-1)                      (0,-0)
(0,1)                       (2,0)                      (1,-1)
(0,0)                       (1,1)                       (2,0)
H =
(1.41421,0)               (0,-0.707107)                      (0,-0)
(0,0.707107)                 (1.22474,0)        (0.816497,-0.816497)
(0,0)         (0.816497,0.816497)                (0.816497,0)
inv(H) =
(1,0)                       (0,1)                 (-0.5,-0.5)
(0,-1)                       (2,0)                      (-1,1)
(-0.5,0.5)                     (-1,-1)                     (1.5,0)