Home Browse Files

Content

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

flens/examples/lapack-potri.cc
#include <flens/flens.cxx>
#include <iostream>

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.

    A = 2, 0, 0,
        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\).

    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

flens/examples/lapack-complex-potri.cc
#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(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

#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(3, 3);

Setup the raw data.

    A = 2,    0, 0,
        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\).

    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;

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