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

flens/examples/lapack-gesv.cc

#include <cxxstd/iostream.h>

#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

flens/examples/lapack-complex-gesv.cc

#include <cxxstd/complex.h>

#include <cxxstd/iostream.h>

#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(1, 0), T( 1,-1), T(  2,20),

        T(0, 1), T( 1, 2), T(-10, 8),

        T(0,-1), T(-1, 1), T( 40, 6);

    b = T( 1, 0),

        T(-1, 1),

        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)