Content |
Compute the Inverse of a Triangular Matrix
In this example we compute the inverse \(T^{-1}\) of a triangular matrix \(T\). The inverse \(T^{-1}\) gets computed by lapack::tri which is the FLENS port of LAPACK's dtrtri.
Example Code
#include <flens/flens.cxx>
using namespace flens;
using namespace std;
int
main()
{
GeMatrix<FullStorage<double> > A(3, 3);
DenseVector<Array<double> > b(3);
A = 2, 1, 0,
0, 2, 1,
0, 0, 2;
auto T = A.upper();
cout << "T = " << T << endl;
int info = lapack::tri(T);
if (info!=0) {
cerr << "T is singular" << endl;
return info;
}
cout << "inv(T) = " << T << endl;
return 0;
}
Comments on Example Code
Setup the raw data.
0, 2, 1,
0, 0, 2;
\(T\) is an upper triangular matrix referencing the upper triangular part of matrix \(A\).
cout << "T = " << T << endl;
Computes the inverse \(T^{-1}\) of matrix \(T\).
If info is not zero then matrix \(T\) was singular. To be more precise, in this case one diagonal element was exactly zero.
cerr << "T is singular" << endl;
return info;
}
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-trtri lapack-trtri.cc
Run
$shell> cd flens/examples $shell> ./lapack-trtri T = 2 1 0 0 2 1 0 0 2 inv(T) = 0.5 -0.25 0.125 0 0.5 -0.25 0 0 0.5
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);
DenseVector<Array<Complex> > b(3);
A = 2, I, 0,
0, 2, 1.+I,
0, 0, 3;
auto T = A.upper();
cout << "T = " << T << endl;
int info = lapack::tri(T);
if (info!=0) {
cerr << "T is singular" << endl;
return info;
}
cout << "inv(T) = " << T << endl;
return 0;
}
Comments on Example Code
Setup the raw data.
0, 2, 1.+I,
0, 0, 3;
\(T\) is an upper triangular matrix referencing the upper triangular part of matrix \(A\).
cout << "T = " << T << endl;
Computes the inverse \(T^{-1}\) of matrix \(T\).
If info is not zero then matrix \(T\) was singular. To be more precise, in this case one diagonal element was exactly zero.
cerr << "T is singular" << endl;
return info;
}
Compile
Note that an external LAPACK implementation is needed for the complex variant of lapack::tri
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-trtri lapack-complex-trtri.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-trtri T = (2,0) (0,1) (0,0) (0,0) (2,0) (1,1) (0,0) (0,0) (3,0) inv(T) = (0.5,0) (0,-0.25) (-0.0833333,0.0833333) (0,0) (0.5,0) (-0.166667,-0.166667) (0,0) (0,0) (0.333333,0)