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>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
GeMatrix<FullStorage<double> >   A(33);
DenseVector<Array<double> >      b(3);

A = 210,
021,
002;

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

Setup the raw data.

A = 210,
021,
002;

$$T$$ is an upper triangular matrix referencing the upper triangular part of matrix $$A$$.

auto T = A.upper();
cout << "T = " << T << endl;

Computes the inverse $$T^{-1}$$ of matrix $$T$$.

int info = lapack::tri(T);

If info is not zero then matrix $$T$$ was singular. To be more precise, in this case one diagonal element was exactly zero.

if (info!=0) {
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>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
typedef complex<double>             Complex;
const Complex                       I(0,1);

GeMatrix<FullStorage<Complex> >     A(33);
DenseVector<Array<Complex> >        b(3);

A = 2, I,    0,
021.+I,
00,    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;
}

Setup the raw data.

A = 2, I,    0,
021.+I,
00,    3;

$$T$$ is an upper triangular matrix referencing the upper triangular part of matrix $$A$$.

auto T = A.upper();
cout << "T = " << T << endl;

Computes the inverse $$T^{-1}$$ of matrix $$T$$.

int info = lapack::tri(T);

If info is not zero then matrix $$T$$ was singular. To be more precise, in this case one diagonal element was exactly zero.

if (info!=0) {
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++ -DUSE_CXXLAPACK -framework vecLib -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)