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

Comments on Example Code

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

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(33);
    DenseVector<Array<Complex> >        b(3);

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

    cout << "inv(T) = " << T << 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-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)