Content

Solving Systems of Linear Equations

In this example we solve a system of linear equations \(Ax = b\) were the coefficient matrix is general, i.e. not necessarily symmetric. We solve the system with lapack::sv which is the FLENS port of LAPACK's dgesv.

Example Code

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

using namespace std;
using namespace flens;

typedef double   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 4;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

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

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

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    lapack::sv(A, piv, b);

    cerr << "x = " << b << endl;
}

Comments on Example Code

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

Define some convenient typedefs for the matrix/vector types of our system of linear equations.

    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;

We also need an extra vector type for the pivots. The type of the pivots is taken for the system matrix.

    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

Set up the baby problem ...

    const IndexType n = 4;

And solve it with lapack::sv

    lapack::sv(A, piv, b);

Compile

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

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-gesv                                                          
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 <complex>
#include <iostream>
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef complex<double>   T;

int
main()
{
    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;
    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

    const IndexType n = 3;

    Matrix         A(n,n);
    Vector         b(n);
    IndexVector    piv(n);

    A = T(10), T( 1,-1), T(  2,20),
        T(01), T( 12), T(-108),
        T(0,-1), T(-11), T( 406);

    b = T( 10),
        T(-11),
        T(-2,-1);

    cerr << "A = " << A << endl;
    cerr << "b = " << b << endl;

    lapack::sv(A, piv, b);

    cerr << "A = " << A << endl;
    cerr << "piv = " << piv << endl;
    cerr << "x = " << b << endl;
}

Comments on Example Code

With header flens.cxx all of FLENS gets included.

#include <flens/flens.cxx>

Elements are now defined to be complex

typedef complex<double>   T;

Like before we define convenient matrix/vector types ...

    typedef GeMatrix<FullStorage<T> >           Matrix;
    typedef DenseVector<Array<T> >              Vector;
    typedef Matrix::IndexType                   IndexType;
    typedef DenseVector<Array<IndexType> >      IndexVector;

Then we setup another toy problem

    const IndexType n = 3;

And solve it with lapack::sv

    lapack::sv(A, piv, b);

Compile

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

Run

$shell> cd flens/examples                                                      
$shell> ./lapack-complex-gesv                                                  
A = 
                       (1,0)                       (1,-1)                       (2,20) 
                       (0,1)                        (1,2)                      (-10,8) 
                      (0,-1)                       (-1,1)                       (40,6) 
b = 
                       (1,0)                        (-1,1)                       (-2,-1) 
A = 
                       (1,0)                       (1,-1)                       (2,20) 
                      (0,-1)                        (0,2)                       (20,8) 
                       (0,1)                      (0.5,0)                        (0,2) 
piv = 
            1              3              3 
x = 
                      (0,-1)                         (0,1)                         (0,0)