Content

Compute the Inverse of a General Matrix

In this example we compute the inverse \(A^{-1}\) of a general \(m \times n\) matrix \(A\). We first compute the \(LU\) factorization. The inverse \(A^{-1}\) then gets computed by lapack::tri which is the FLENS port of LAPACK's dgetri.

Example Code

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

using namespace flens;
using namespace std;

int
main()
{
    GeMatrix<FullStorage<double> >   A(3,3);
    DenseVector<Array<int> >         piv(3);

    A = 123,
        456,
        781;

    cout << "A = " << A << endl;

    int info = lapack::trf(A, piv);

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

    lapack::tri(A, piv);
    cout << "inv(A) = " << A << endl;

    return 0;
}

Comments on Example Code

Setup the toy example.

    A = 123,
        456,
        781;

Compute the \(LU\) factoriation of \(A\).

    int info = lapack::trf(A, piv);

If info is not zero then the matrix was singular. More precise, the upper triangular matrix \(U\) computed by trf is exactly singular.

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

Compute \(A^{-1}\) from the \(LU\) factorization of \(A\).

    lapack::tri(A, piv);
    cout << "inv(A) = " << A << endl;

Compile

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

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-getri                                                         
A = 
            1             2             3 
            4             5             6 
            7             8             1 
inv(A) = 
     -1.79167      0.916667        -0.125 
      1.58333     -0.833333          0.25 
       -0.125          0.25        -0.125 

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(44);
    DenseVector<Array<Complex> >        b(4);
    DenseVector<Array<int> >            piv(4);

    A =   2,    3,  -I,    0,
         -6,   -5,   02.+I,
       2.*I,   -5,   6,   -6,
          46.*I,   2,   -3;

    b = 20,
       -33.*I,
       -43,
        49;

    cout << "A = " << A << endl;
    cout << "b = " << b << endl;

    int info = lapack::trf(A, piv);

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

    lapack::tri(A, piv);
    cout << "inv(A) = " << A << endl;

    DenseVector<Array<Complex> >        x;
    x = A*b;

    cout << "x = inv(A)*b = " << x << endl;

    return 0;
}

Comments on Example Code

Setup matrix \(A\). Note that 2+I or 2*I would not work. That's because the STL does not define the operation int+complex<double> or int*complex<double>.

    A =   2,    3,  -I,    0,
         -6,   -5,   02.+I,
       2.*I,   -5,   6,   -6,
          46.*I,   2,   -3;

Compute the \(LU\) factoriation of \(A\).

    int info = lapack::trf(A, piv);

If info is not zero then the matrix was singular. More precise, the upper triangular matrix \(U\) computed by trf is exactly singular.

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

Compute \(A^{-1}\) from the \(LU\) factorization of \(A\).

    lapack::tri(A, piv);
    cout << "inv(A) = " << A << endl;

Compute \(x = A^{-1}b\).

    DenseVector<Array<Complex> >        x;
    x = A*b;

Compile

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

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-getri                                                  
A = 
                       (2,0)                        (3,0)                      (-0,-1)                        (0,0) 
                      (-6,0)                       (-5,0)                        (0,0)                        (2,1) 
                       (0,2)                       (-5,0)                        (6,0)                       (-6,0) 
                       (4,0)                        (0,6)                        (2,0)                       (-3,0) 
b = 
                      (20,0)                      (-0,-33)                       (-43,0)                        (49,0) 
inv(A) = 
         (0.244811,0.432046)        (-0.0798371,0.285037)      (-0.0677862,-0.0410065)        (-0.0126646,0.245425) 
       (-0.116046,0.0285651)      (-0.0920553,-0.0927248)        (0.0161794,0.0344789)       (-0.0628208,-0.161459) 
         (0.949788,0.858514)           (0.291899,0.43584)        (0.0214238,0.0870341)        (0.00647177,0.213792) 
         (0.902477,0.916313)            (0.2736,0.486499)        (-0.145057,0.0357063)         (-0.0229859,0.34412) 
x = inv(A)*b = 
           (16.5967,25.0647)           (-9.15476,-5.78498)             (32.7744,14.2709)              (39.2151,24.624)