Element-wise Matrix/Vector Operations

For-Loops provides the easiest, most readable and most efficient way to do element-wise operations. However, in the case of matrices this requires you to order loops depending on row and column major storage. In FLENS you can use special index variables to iterate over a general or triangular matrix the best order.

Example for GeMatrix

flens/examples/gematrix-indexvariable.cc

#include <complex>

#include <iostream>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double> >             RealGeMatrix;

    RealGeMatrix   A(3,4);

    RealGeMatrix::IndexVariable i, j;

    A(i,j) = 10*i + j;

    RealGeMatrix   B(4,3);

    B(i,j) = Abs(i+j)*(i-j);

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

    cout << "B = " << B << endl;

    RealGeMatrix   C(3,4);

    C(i,j) = A(i,j) * Sin(M_PI*B(j,i)/4);

    cout << "C = " << C << endl;

}

Some notes

We allocate a matrix and define index variables.

    RealGeMatrix   A(3,4);

    RealGeMatrix::IndexVariable i, j;

In the following statement \(i\) and \(j\) will automatically iterate depending on the matrix size of \(A\). Obviously \(i\) and \(j\) are not some plain integer values. But you can do some simple computations with them.

A(i,j) = 10*i + j;

Another example which is using an Abs function that can operate on such index variables.

    RealGeMatrix   B(4,3);

    B(i,j) = Abs(i+j)*(i-j);

Now we also can do some simple element-wise operation. Note that the range of \(i\) and \(j\) is determined by the matrix size of \(C\)

    RealGeMatrix   C(3,4);

    C(i,j) = A(i,j) * Sin(M_PI*B(j,i)/4);

So check out the results

$shell> cd flens/examples                                                         
$shell> g++ -Wall -std=c++11 -I../.. gematrix-indexvariable.cc                    
$shell> ./a.out                                                                   
A = 
           11            12            13            14 
           21            22            23            24 
           31            32            33            34 
B = 
            0            -3            -8 
            3             0            -5 
            8             5             0 
           15            12             7 
C = 
            0       8.48528  -3.18408e-15      -9.89949 
     -14.8492             0      -16.2635   8.81746e-15 
  7.59281e-15       22.6274             0      -24.0416

Example for TrMatrix, SyMatrix, HeMatrix

flens/examples/trmatrix-indexvariable.cc

#include <complex>

#include <iostream>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double> >             RealGeMatrix;

    typedef TrMatrix<FullStorage<double> >             RealTrMatrix;

    typedef SyMatrix<FullStorage<double> >             RealSyMatrix;

    typedef HeMatrix<FullStorage<complex<double> > >   ComplexHeMatrix;

    typedef GeMatrix<FullStorage<complex<double> > >   ComplexGeMatrix;

    RealTrMatrix   A(3, Upper);

    RealTrMatrix::IndexVariable i, j;

    A(i,j) = 10*i+j;

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

    RealSyMatrix      B(3, Upper);

    ComplexHeMatrix   C(3, Lower);

    B.triangular()(i,j) = 2+i-j;

    C.triangular()(i,j) = Complex(i+j,j);

    cout << "B = " << B << endl;

    cout << "C = " << C << endl;

    ComplexGeMatrix   D(3,4);

    D.strictLower()(i,j) = Complex(i,0);

    D.diag(0)(i)         = Complex(666,i);

    D.strictUpper()(i,j) = Complex(i,666);

    cout << "D = " << D << endl;

}

Some notes

We allocate a triangular matrix, define index variables and use them. As the matrix is upper triangular \(i\) and \(j\) will only iterate over the upper triangular part.

    RealTrMatrix   A(3, Upper);

    RealTrMatrix::IndexVariable i, j;

    A(i,j) = 10*i+j;

For symmetric and hermitian matrices ...

    RealSyMatrix      B(3, Upper);

    ComplexHeMatrix   C(3, Lower);

... you simply initialize their dominant triangular part

    B.triangular()(i,j) = 2+i-j;

    C.triangular()(i,j) = Complex(i+j,j);

Finally we initialize different triangular parts of a general complex valued matrix

    ComplexGeMatrix   D(3,4);

    D.strictLower()(i,j) = Complex(i,0);

    D.diag(0)(i)         = Complex(666,i);

    D.strictUpper()(i,j) = Complex(i,666);

So check out the results

$shell> cd flens/examples                                                         
$shell> g++ -Wall -std=c++11 -I../.. trmatrix-indexvariable.cc                    
$shell> ./a.out                                                                   
A = 
           11            12            13 
            0            22            23 
            0             0            33 
B = 
            2             1             0 
            1             2             1 
            0             1             2 
C = 
                       (2,0)                      (3,-1)                      (4,-1)
                       (3,1)                       (4,0)                      (5,-2)
                       (4,1)                       (5,2)                       (6,0)
D = 
                     (666,1)                      (1,666)                      (1,666)                      (1,666) 
                       (1,0)                      (666,2)                      (2,666)                      (2,666) 
                       (2,0)                        (2,0)                      (666,3)                      (3,666)