# 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;
}

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

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

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;

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

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

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;

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

return 0;
}

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