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 <iostream>
using namespace flens;
using namespace std;
int
main()
{
GeMatrix<FullStorage<double> > A(3,3);
DenseVector<Array<int> > piv(3);
A = 1, 2, 3,
4, 5, 6,
7, 8, 1;
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.
4, 5, 6,
7, 8, 1;
Compute the \(LU\) factoriation of \(A\).
If info is not zero then the matrix was singular. More precise, the upper triangular matrix \(U\) computed by trf is exactly singular.
cerr << "A is singular" << endl;
return info;
}
Compute \(A^{-1}\) from the \(LU\) factorization of \(A\).
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(4, 4);
DenseVector<Array<Complex> > b(4);
DenseVector<Array<int> > piv(4);
A = 2, 3, -I, 0,
-6, -5, 0, 2.+I,
2.*I, -5, 6, -6,
4, 6.*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
#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(4, 4);
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>.
-6, -5, 0, 2.+I,
2.*I, -5, 6, -6,
4, 6.*I, 2, -3;
b = 20,
-33.*I,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
Compute the \(LU\) factoriation of \(A\).
If info is not zero then the matrix was singular. More precise, the upper triangular matrix \(U\) computed by trf is exactly singular.
cerr << "A is singular" << endl;
return info;
}
Compute \(A^{-1}\) from the \(LU\) factorization of \(A\).
cout << "inv(A) = " << A << endl;
Compute \(x = A^{-1}b\).
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)