Content |
Compute the Inverse of a Symmetric Positive Definite Matrix
In this example we compute the inverse \(S^{-1}\) of a symmetric/hermitian 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/zpotri.
Example Code
#include <flens/flens.cxx>
using namespace flens;
using namespace std;
int
main()
{
GeMatrix<FullStorage<double> > A(3, 3);
A = 2, 0, 0,
1, 2, 0,
0, 1, 2;
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.
1, 2, 0,
0, 1, 2;
\(S\) is a symmetric matrix view constructed from the lower triangular part of \(A\). Note that \(S\) only references data from \(A\).
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\).
If info is not zero the factorization could not be computed as the matrix was singular.
cerr << "S is singular" << endl;
return info;
}
Compute \(S^{-1}\) from the Cholesky factorization \(S = L L^T\).
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>
using namespace flens;
using namespace std;
int
main()
{
typedef complex<double> Complex;
const Complex I(0,1);
GeMatrix<FullStorage<Complex> > A(3, 3);
A = 2, 0, 0,
I, 2, 0,
0, 1.+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.
I, 2, 0,
0, 1.+I, 2;
\(H\) is a symmetric matrix view constructed from the lower triangular part of \(A\). Note that \(H\) only references data from \(A\).
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\).
If info is not zero the factorization could not be computed as the matrix was singular.
cerr << "H is singular" << endl;
return info;
}
Compute \(H^{-1}\) from the Cholesky factorization \(H = L L^T\).
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++ -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)