Content

Solving Systems of Linear Equations

In this example we will solve a system of linear equations \(Ax=b\). For this purpose we first compute the \(LU\) factorization of \(A\) with lapack::trf and then use the triangular solver lapack::trs to finish the job.

Example Code

#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

int
main()
{
    GeMatrix<FullStorage<double> >   A(44);
    DenseVector<Array<double> >      b(4);
    DenseVector<Array<int> >         piv(4);

    A =  2,   3,  -1,   0,
        -6,  -5,   0,   2,
         2,  -5,   6,  -6,
         4,   6,   2,  -3;

    b = 20,
       -33,
       -43,
        49;

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

    lapack::trf(A, piv);

    lapack::trs(NoTrans, A, piv, b);
    cout << "X = " << b << endl;
}

Comments on Example Code

Compute the \(PLU = A\) factorization with lapack::trf

    lapack::trf(A, piv);

Solve the system of linear equation \(PLUx = b\) using lapack::trs. Vector \(b\) gets overwritten with vector \(x\).

    lapack::trs(NoTrans, A, piv, b);
    cout << "X = " << b << endl;
}

Compile

$shell> cd flens/examples                                                       
$shell> g++ -std=c++11 -Wall -I../.. -o lapack-getrs lapack-getrs.cc            

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-getrs                                                          
A = 
            2             3            -1             0 
           -6            -5             0             2 
            2            -5             6            -6 
            4             6             2            -3 
b = 
           20            -33            -43             49 
X = 
            1              9              9              9 

Example with Complex Numbers

Example Code

#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

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;

    lapack::trf(A, piv);

    lapack::trs(NoTrans, A, piv, b);
    cout << "X = " << b << endl;
}

Comments on Example Code

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;

Compute the \(PLU = A\) factorization with lapack::trf

    lapack::trf(A, piv);

Solve the system of linear equation \(PLUx = b\) using lapack::trs. Vector \(b\) gets overwritten with vector \(x\).

    lapack::trs(NoTrans, A, piv, b);
    cout << "X = " << b << endl;
}

Compile

$shell> cd flens/examples                                                       
$shell> g++ -std=c++11 -Wall -I../.. -o lapack-complex-getrs lapack-complex-getrs.cc                 

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-getrs                                                  
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) 
X = 
           (16.5967,25.0647)           (-9.15476,-5.78498)             (32.7744,14.2709)              (39.2151,24.624)