Content

Full Storage Schemes
'GeMatrix': General Matrix with Full Storage
Simple Example for using GeMatrix
Simple Example for using GeMatrix Views
Think in C, Write in C++
More on Selecting Rectangular Parts from GeMatrix
Interface in Detail: Creating Views Referencing Matrix Parts
Matrix-View Gotchas
Using auto when dealing with GeMatrix Views
Why you need 3 GeMatrix Types
Matrix Assignments are neither a Deep or Shallow copies!

Full Storage Schemes and General Matrices

The term general matrix refers to a matrix that is not necessarily square or symmetric. Further, full storage denotes that for a \(m \times n\) matrix all its \(m \cdot n\) elements are stored somewhere in memory.

In FLENS a matrix type always includes an underlying storage type. The general idea of this is well known in the LAPACK/BLAS world. However, in the traditional Fortran LAPACK/BLAS this is merley done by conventions. It is completely up to the programmer how to interpretate some data. When you pass a double array to dgetrs its treated as a general matrix, if you pass it to dtrtrs as a triangular matrix. On the on hand it is flexible on the other hand it can be error prone.

In FLENS/C++ you get language support through strong typing. That means errors can be detected by the compiler. You will see on later slides that at the same time we do not loose flexibility.

Full Storage Schemes

FLENS defines three template classes for full storage schemes:

Template Class FullStorage

Constructors of FullStorage allocate memory and the destructor deallocates it. The first template parameter specifies the element type. The second parameter specifies whether elements get stored row-wise (RowMajor) or column-wise (ColMahor).

It is guaranteed that elements are stored contiguously in memory, i.e. without any gaps.

Examples:

FullStorage<double, ColMajor>	Full storage scheme that allocates memory when instantiated. Elements are stored contiguously and column-wise.
FullStorage<double, RowMajor>	Full storage scheme that allocates memory when instantiated. Elements are stored contiguously and row-wise.
FullStorage<double>	The second template parameter of `FullStorage` defaults to `ColMajor`. Hence this is equivalent to `FullStorage<double, ColMajor>`.

Creating an instance of type FullStorage involves dynamic memory allocation. So you do not want to create an instance inside a loop or function.

Template Class FullStorageView

A storage view does not allocate or release any memory. It is just used to reference elements that were originally allocated by a still existing and living FullStorage instance.

An instance of type FullStorageView just contains a pointer to raw data and some integers for matrix dimensions and element strides. If you are familiar with BLAS/LAPACK That is essentially what you have to pass to a BLAS/LAPACK function: dimensions, pointer and leading dimension.

Creating an instance of type FullStorageView is cheap as it just initializes a few integer values of a struct/class. The destructor does nothing.
Template Class ConstFullStorageView

Instances of ConstFullStorageView are used for const-referencing elements from an existing and living FullStorage instance.

The creation/destruction of an ConstFullStorageView instance is cheap.

Further template parameters of the storage classes allow to change the default index type (which is `long), the default index base (which is 1) and the memory allocator.

'GeMatrix': General Matrix with Full Storage

The GeMatrix class has only one template parameter that defines the underlying storage scheme. That means we give the storage scheme a mathematical meaning.

Examples:

GeMatrix<FullStorage<double, ColMajor> >

A general matrix where memory gets allocated when instanciated and released at the end of its scope.

GeMatrix<FullStorage<double, ColMajor> >::View

A general matrix that references a rectangular part of another GeMatrix instance. Memory neither gets allocated nor released at any point. Creation/destruction is cheap.

View is a publib typedef in GeMatrix and in this case defined as GeMatrix<FullStorageView<double, ColMajor> >

GeMatrix<FullStorage<double, ColMajor> >::ConstView

A general matrix that references a rectangular part of another GeMatrix instance. Memory neither gets allocated nor released at any point. Creation/destruction is cheap

Calling any non-const methods will cause a compile time error. So this is a read-only matrix view.

ConstView is a public typedef in GeMatrix and in this case defined as GeMatrix<ConstFullStorageView<double, ColMajor> >

Some Public Typedefs

IndexType	Like its name suggests the type used for indexing.
ElementType	The element type.
View, ConstView	Types for referencing a rectangular part of the matrix.
NoView	Matrix type with memory allocation. If you want to copy matrix parts instead of referencing use this type.
EngineView, EngineConstView	Storage schemes for referencing parts from another storage scheme. For `GeMatrix` these are typedefs for `FullStorageView<..>` and `ConstFullStorageView<..>`.

Some Methods and operations for Matrix Manipulation

GeMatrix(IndexType m, IndexType n);	Constructor for a \(m \times n\) matrix.
IndexType numRows() const;	Return the number of rows.
IndexType numCols() const;	Return the number of columns.
const ElementType & operator()(i, j) const;	Return a reference to element in row \(i\) and column \(j\).
ElementType & operator()(i, j);
Initializer operator=(const ElementType &value);	List initializer. Allows to initialize the matrix with a comma separated list of values. If the list contains only a single value it is used to fill the matrix.

By default indices start at 1, i.e. the index base is 1. Usually that is what you want in numerical linear algebra. However, we will see that it is fairly easy to change this if desired. In this case the following methods become useful.

IndexType firstRow() const;	Return the index of the first row.
IndexType lastRow() const;	Return the index of the last row.
IndexType firstCol() const;	Return the index of the first column.
IndexType lastCol() const;	Return the index of the last column.

Simple Example for using `GeMatrix`

In a first example we show:

How to allocate a general matrix with full storage.
How to explicitly initialize all elements with the list initializer.
How to fill the matrix with a single value.
How to access/manipulate a particular matrix entry.
How to retrieve matrix dimensions and index ranges.

Example Code

flens/examples/tut01-page01-example1.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double, ColMajor> >  DGeMatrix;

    DGeMatrix A(4,4);

    A = 1, 2, 3,  4,

        5, 6, 7,  8,

        9, 8, 7,  6,

        5, 4, 3, 20;

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

    cout << "Dim. of A: " << A.numRows() << " x " << A.numCols() << endl;

    cout << endl;

    cout << "Row indices: " << A.firstRow() << ".." << A.lastRow() << endl;

    cout << endl;

    cout << "Col indices: " << A.firstCol() << ".." << A.lastCol() << endl;

    cout << endl;

    A(3,2) = 42;

    cout << "changed element: A(3,2) = " << A(3,2) << endl;

    cout << endl;

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

    A = 42;

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

    return 0;

}

Comments on Example Code

We simply include everything of FLENS

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

Typedef for a general matrix with elements of type double. Internally elements are stored in column major order.

typedef GeMatrix<FullStorage<double, ColMajor> > DGeMatrix;

Matrix A gets dynamically allocated and then initialized.

    DGeMatrix A(4,4);

    A = 1, 2, 3,  4,

        5, 6, 7,  8,

        9, 8, 7,  6,

        5, 4, 3, 20;

Print the matrix content using output streams:

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

We print some information about matrix dimensions and index ranges. You will see that by default in FLENS indices start at 1 (like in Fortran):

    cout << "Dim. of A: " << A.numRows() << " x " << A.numCols() << endl;

    cout << endl;

    cout << "Row indices: " << A.firstRow() << ".." << A.lastRow() << endl;

    cout << endl;

    cout << "Col indices: " << A.firstCol() << ".." << A.lastCol() << endl;

    cout << endl;

Also for element access (write) we provide a Fortran-Style interface:

A(3,2) = 42;

The same for read access:

cout << "changed element: A(3,2) = " << A(3,2) << endl;

You also can fill the whole matrix with a new value:

A = 42;

Compile and Run

$shell> cd flens/examples                                                       
$shell> g++ -Wall -std=c++11 -I../.. tut01-page01-example1.cc                   
$shell> ./a.out                                                                 
A = 
            1             2             3             4 
            5             6             7             8 
            9             8             7             6 
            5             4             3            20 
Dim. of A: 4 x 4
Row indices: 1..4
Col indices: 1..4
changed element: A(3,2) = 42
A = 
            1             2             3             4 
            5             6             7             8 
            9            42             7             6 
            5             4             3            20 
A = 
           42            42            42            42 
           42            42            42            42 
           42            42            42            42 
           42            42            42            42

Simple Example for using `GeMatrix` Views

In a first example we show:

How to create GeMatrix views that reference rectangular parts of another GeMatrix instance.
There are two view variants, i.e. const-views and (non-const-)views.

Example Code

flens/examples/tut01-page01-example2.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double, ColMajor> >      DGeMatrix;

    typedef GeMatrix<FullStorageView<double, ColMajor> >  DGeMatrixView;

    const Underscore<DGeMatrix::IndexType>  _;

    DGeMatrix  A(3,4);

    A = 1,  2,  3,  4,

        5,  6,  7,  8,

        9, 10, 11, 12;

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

    DGeMatrixView B = A(_(2,3),_);

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

    B(1,2) = 42;

    cout << "Changed B(1,2)." << std::endl;

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

    DGeMatrix C = B;

    cout << "C is a copy of B." << std::endl;

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

    C(1,2) = 666;

    cout << "Changed C(1,2)." << std::endl;

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

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

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

}

Comments on Example Code

In the main function we again start with a typedef for a GeMatrix with elements of type double that are stored column-wise in memory. Also, we define a suitable range operator.

int

main()

{

    typedef GeMatrix<FullStorage<double, ColMajor> >      DGeMatrix;

    typedef GeMatrix<FullStorageView<double, ColMajor> >  DGeMatrixView;

We define an underscore operator for selecting index ranges. The template class Underscore requires the index type of DGeMatrix as template parameter.

const Underscore<DGeMatrix::IndexType> _;

Then we setup some matrix.

    DGeMatrix  A(3,4);

    A = 1,  2,  3,  4,

        5,  6,  7,  8,

        9, 10, 11, 12;

Create a view B that references from A the second and third row.

DGeMatrixView B = A(_(2,3),_);

We change an element of B. So we also change A.

B(1,2) = 42;

Create a copy C of B

DGeMatrix C = B;

We change an element of C. So this neither will change A nor B as C has its own memory.

C(1,2) = 666;

Compile and Run

$shell> cd flens/examples                                                       
$shell> g++ -Wall -std=c++11 -I../.. tut01-page01-example2.cc                   
$shell> ./a.out                                                                 
A = 
            1             2             3             4 
            5             6             7             8 
            9            10            11            12 
B = 
            5             6             7             8 
            9            10            11            12 
Changed B(1,2).
A = 
            1             2             3             4 
            5            42             7             8 
            9            10            11            12 
C is a copy of B.
C = 
            5            42             7             8 
            9            10            11            12 
Changed C(1,2).
A = 
            1             2             3             4 
            5            42             7             8 
            9            10            11            12 
B = 
            5            42             7             8 
            9            10            11            12 
C = 
            5           666             7             8 
            9            10            11            12

Think in C, Write in C++

If you haven't done already you really have to checkout Write in C :-)

The advantage of C (or Fortran) is its simplicity. So thinking in C means you know what is going on. Coding the idea in C++ allows to simplify and generalize this idea. If you would use structs in you C code using C++ features like constructors makes code less error prone. Features like templates add type safety to a macro-solution in C.

Just keep it simple. Don't use C++ features for the sake of using cool features. Start with a straight forward and simple idea in C and only use C++ features if they help to improve things. E.g. don't use iterators just because you think in C++ everybody uses iterators. Ask yourself if iterators make sense in the context of your application.

FLENS was developed in the spirit: How would you do it in C or Fortan? And How can it be generalized or simplified without sacrifices in C++? So a good way of understanding concepts in FLENS is boiling them down to simple C.

How would you do it in C?

So here we port back some of the essential parts of the previous example to C. This pretty much explains what gets done:

What does the constructor/destructor of GeMatrix: allocate/release memory
What does the constructor/destructor of GeMatrix views: Not much/Nothing
What happen when you create a view: Just setting a pointer to the right place
What happens when you access elements using the function operator: Some pointer arithmetic done correct and dereferencing

flens/examples/tut01-page01-example3.cc

#include <cstdio>

#include <cstdlib>

struct DGeMatrix

{

    double *data;

    long   numRows, numCols;

};

struct DGeMatrixView

{

    double *data;

    long   numRows, numCols;

    long   leadingDimension;

};

// Constructor for DGeMatrix

void

DGeMatrix_alloc(DGeMatrix *A, long m, long n)

{

    A->data    = (double *)malloc(m*n*sizeof(double));

    A->numRows = m;

    A->numCols = n;

}

// Destructor for DGeMatrix

void

DGeMatrix_free(DGeMatrix *A)

{

    free(A->data);

    A->data = 0;

}

// View from DGematrix

void

DGeMatrix_view(const DGeMatrix *A, DGeMatrixView *B,

               long fromRow, long toRow,

               long fromCol, long toCol)

{

    long ldA = A->numRows;

    B->data             = A->data + fromRow + fromCol*ldA;

    B->numRows          = toRow-fromRow+1;

    B->numCols          = toCol-fromCol+1;

    B->leadingDimension = ldA;

}

// Change entry in DGeMatrix

void

DGeMatrix_set(const DGeMatrix *A, long row, long col, double value)

{

    A->data[row+col*A->numRows] = value;

}

// Get entry from DGeMatrix

double

DGeMatrix_get(const DGeMatrix *A, long row, long col)

{

    return A->data[row+col*A->numRows];

}

// Change entry in DGeMatrixView

void

DGeMatrixView_set(const DGeMatrixView *A, long row, long col, double value)

{

    A->data[row+col*A->leadingDimension] = value;

}

// Get entry from DGeMatrixView

double

DGeMatrixView_get(const DGeMatrixView *A, long row, long col)

{

    return A->data[row+col*A->leadingDimension];

}

int

main()

{

    // DGeMatrix  A(3,4);

    DGeMatrix A;

    DGeMatrix_alloc(&A, 3, 4);

    // A = 1,  2,  3,  4,

    //     5,  6,  7,  8,

    //     9, 10, 11, 12;

    //

    // The list initializer always fills the matrix row wise.  Even if col wise

    // would be faster.  But anyway, the sole purpose of the list initializer is

    // for writing compact examples in the tutorial

    //

    long counter = 1;

    for (long j=0; j<A.numCols; ++j) {

        for (long i=0; i<A.numRows; ++i) {

            DGeMatrix_set(&A, i, j, counter);

            ++counter;

        }

    }

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

    printf("A =\n");

    for (long i=0; i<A.numRows; ++i) {

        for (long j=0; j<A.numCols; ++j) {

            printf("%8.3lf", DGeMatrix_get(&A, i, j));

        }

        printf("\n");

    }

    printf("\n");

    // DGeMatrixView B = A(_(2,3),_);

    DGeMatrixView B;

    DGeMatrix_view(&A, &B, 1, 2, 0, 3);

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

    printf("B =\n");

    for (long i=0; i<B.numRows; ++i) {

        for (long j=0; j<B.numCols; ++j) {

            printf("%8.3lf", DGeMatrixView_get(&B, i, j));

        }

        printf("\n");

    }

    printf("\n");

    // B(1,2) = 42;

    DGeMatrixView_set(&B, 0, 1, 42);

    // cout << "Changed B(1,2)." << std::endl;

    printf("Changed B(1,2).\n");

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

    printf("A =\n");

    for (long i=0; i<A.numRows; ++i) {

        for (long j=0; j<A.numCols; ++j) {

            printf("%8.3lf", DGeMatrix_get(&A, i, j));

        }

        printf("\n");

    }

    printf("\n");

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

    printf("B =\n");

    for (long i=0; i<B.numRows; ++i) {

        for (long j=0; j<B.numCols; ++j) {

            printf("%8.3lf", DGeMatrixView_get(&B, i, j));

        }

        printf("\n");

    }

    printf("\n");

    // Destructor for A

    DGeMatrix_free(&A);

}

Compile and Run

$shell> cd flens/examples                                                       
$shell> g++ -Wall -std=c++11 -I../.. tut01-page01-example3.cc                   
$shell> ./a.out                                                                 
A =
   1.000   4.000   7.000  10.000
   2.000   5.000   8.000  11.000
   3.000   6.000   9.000  12.000
B =
   2.000   5.000   8.000  11.000
   3.000   6.000   9.000  12.000
Changed B(1,2).
A =
   1.000   4.000   7.000  10.000
   2.000  42.000   8.000  11.000
   3.000   6.000   9.000  12.000
B =
   2.000  42.000   8.000  11.000
   3.000   6.000   9.000  12.000

More on Selecting Rectangular Parts from `GeMatrix`

From a GeMatrix instance you can select rectangular parts be specifying row and column ranges. This can have the following forms:

all
from:to
from:step:to

In FLENS this is realized through the template classes Underscore and Range. The template parameter specifies the index type using in the GeMatrix instance, which is usually long:

An object of type Underscore always means all. And giving such an object the name _ gives a compact notation.

Underscore<long> _;
Objects of type Range can be created by calling the function operator of Underscore:

Range<long> range1 _(1,5); // 1:5
Range<long> range2 _(1,2,5); // 1:2:5
If for some reason you want a view of a complete matrix you don't have to select ranges. Simply construct the view from the original matrix:

typedef GeMatrix<FullStorage<double, ColMajor> >      DGeMatrix;
typedef DGeMatrix::View                               DGeMatrixView;
typedef DGeMatrix::ConstView                          DGeMatrixConstView;

DGeMatrix                  A(5, 8);
DGeMatrixView              X = A;      // X is a view of A
const DGeMatrixConstView   Y = A;      // Y is a const-view of A

Note that if you use from:step:to ranges the resulting GeMatrix view will be neither ColMajor nor RowMajor. We call the storage order in this case Grid. To our knowledge only BLIS and ulmBLAS efficiently support linear algebra operations for this kind of views.

Interface in Detail: Creating Views Referencing Matrix Parts

Methods for selecting rectangular matrix parts:

const ConstView operator()(const Range<IndexType> &rows, const Range<IndexType> &cols) const;	Select a rectangualr area of the form `(fromRow:toRow, fromCol:toCol)` or `(fromRow:stepRow:toRow, fromCol:stepRow:toCol)`
View operator()(const Range<IndexType> &rows, const Range<IndexType> &cols);
template <typename RHS> const ConstView operator()(const GeMatrix<RHS> &A) const;	The area gets selected by the index ranges of matrix `A`. That means the range is `(A.firstRow():A.lastRow(), A.firstCol():A.lastCol())`
template <typename RHS> View operator()(const GeMatrix<RHS> &A);
const ConstView operator()(const Underscore<IndexType> &, const Range<IndexType> &cols) const;	All rows are selected. So the range has the form `(all, from:to)`.
View operator()(const Underscore<IndexType> &, const Range<IndexType> &cols);
const ConstView operator()(const Range<IndexType> &rows, const Underscore<IndexType> &) const;	All columns are selected. So the range has the form `(from:to, all)`.
View operator()(const Range<IndexType> &rows, const Underscore<IndexType> &);

Methods for selecting vector parts (the DenseVector class will be introduced on the next page):

const ConstVectorView operator()(IndexType row, const Underscore<IndexType> &) const;	Select `(row, all)`, i.e. the row with index `row`.
VectorView operator()(IndexType row, const Underscore<IndexType> &);
const ConstVectorView operator()(IndexType row, const Range<IndexType> &cols) const;	Select `(row, from:to)`, i.e. from the row with index `row` the elements in the column in range `from:to`.
VectorView operator()(IndexType row, const Range<IndexType> &cols);
const ConstVectorView operator()(const Underscore<IndexType> &, IndexType col) const;	Select `(all, col)`, i.e. the column with index `col`.
VectorView operator()(const Underscore<IndexType> &, IndexType col);
const ConstVectorView operator()(const Range<IndexType> &rows, IndexType col) const;	Select `(from:to, col)`, i.e. from the column with index `col` the elements in the row range `from:to`.
VectorView operator()(const Range<IndexType> &rows, IndexType col);
const ConstVectorView diag(IndexType d) const;	Select all elements on the `d`-the diagonal. The main diagonal is specified by , the first super diagonal by `1`, the first sub diagonal by `-1`, etc.
VectorView diag(IndexType d);
const ConstVectorView antiDiag(IndexType d) const;	Select all elements on the `d`-the anti-diagonal. The main anti-diagonal is specified by , the first super anti-diagonal by `1`, the first sub anti-diagonal by `-1`, etc.
VectorView antiDiag(IndexType d);

Matrix-View Gotchas

Short and simple version: Always declare a const-view as const.

A bit longer: As the underlying storage scheme of a const-view only implements const-methods calling a non-const method of GeMatrix always causes a compile time error. Using static asserts we try to make this error messages more readable.

Technically it would be possible to disable in such a cause all non-const methods of GeMatrix using the std::enable_if trait. However, declaring a const-view not as const is a flaw. And we don't want to encourage flaws.

Example Code

flens/examples/tut01-page01-example4.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double, ColMajor> >      DGeMatrix;

    typedef DGeMatrix::ConstView                          DGeMatrixConstView;

    const Underscore<DGeMatrix::IndexType>  _;

    DGeMatrix  A(3,4);

    A = 1,  2,  3,  4,

        5,  6,  7,  8,

        9, 10, 11, 12;

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

    DGeMatrixConstView  B = A(_(2,3),_);

    cout << "B(1,1) = " << B(1,1) << endl;

}

Comments on Example Code

We want to create a const-view. However, we do not declare B as const. Ideally the instantiation here would already cause a compile time error. You see next, why.

DGeMatrixConstView B = A(_(2,3),_);

GeMatrix implements two function operators for element access. For a const-view merely the read-only access must be allowed. But as B is not declared as const the read-write access gets called.

cout << "B(1,1) = " << B(1,1) << endl;

Compile and Run (Yes, this should fail!!)

$shell> cd flens/examples                                                       
$shell> g++ -Wall -std=c++11 -I../.. tut01-page01-example4.cc                   
In file included from tut01-page01-example4.cc:2:
In file included from ../../flens/flens.cxx:60:
In file included from ../../flens/flens.tcc:41:
In file included from ../../flens/matrixtypes/matrixtypes.tcc:38:
In file included from ../../flens/matrixtypes/general/impl/impl.tcc:40:
../../flens/matrixtypes/general/impl/gematrix.tcc:268:12: error: binding of reference to type 'ElementType' (aka 'double') to a value of type 'const ElementType' (aka 'const double') drops qualifiers
    return engine_(row, col);
           ^~~~~~~~~~~~~~~~~
tut01-page01-example4.cc:34:29: note: in instantiation of member function 'flens::GeMatrix, std::__1::allocator > >::operator()' requested here
    cout << "B(1,1) = " << B(1,1) << endl;
                            ^
1 error generated.

Using `auto` when dealing with `GeMatrix` Views

The auto feature in C++11 is a great way to simplify C++ code. But it is common sense that you actually understand this C++ feature. Otherwise don't use it. The C++ compiler knows the rules that determine the actual type represented by auto. And so should you! If you do, the code becomes more readable for you. If you don't, the code become obscure.

There is a nice talk by Scott Meyers Type Deduction and Why You Care. Have a look at it before you blindly use auto.

Fortunately deducing the type of auto becomes simple and transparent if you recall:

If A is of type GeMatrix<FullStorage<double, ColMajor> > then A(_(1,3),_(2,5)) returns
- the type GeMatrix<ConstFullStorageView<double, ColMajor> > if A is in const context and
- the type GeMatrix<FullStorageView<double, ColMajor> > otherwise.
If A is of type GeMatrix<FullStorageView<double, ColMajor> > then A(_(1,3),_(2,5)) returns
- the type GeMatrix<ConstFullStorageView<double, ColMajor> > if A is in const context and
- the type GeMatrix<FullStorageView<double, ColMajor> > otherwise.
If A is of type GeMatrix<ConstFullStorageView<double, ColMajor> > then A must be in const context and always returns the type GeMatrix<ConstFullStorageView<double, ColMajor> >.

Example

So consider the following example and deduce the type of auto in each case:

flens/examples/tut01-page01-example5.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

using namespace flens;

using namespace std;

int

main()

{

    typedef GeMatrix<FullStorage<double, ColMajor> >      DGeMatrix;

    typedef DGeMatrix::View                               DGeMatrixView;

    typedef DGeMatrix::ConstView                          DGeMatrixConstView;

    const Underscore<DGeMatrix::IndexType>  _;

    DGeMatrix  A(3,4);

    A = 1,  2,  3,  4,

        5,  6,  7,  8,

        9, 10, 11, 12;

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

    auto       B = A(_(2,3),_);       // auto: DGeMatrixView,

                                      // B:    DGeMatrixView

    const auto C = B(_,_(2,4));       // auto: DGeMatrixView,

                                      // C:    const DGeMatrixView

    const auto D = C;                 // auto: DGeMatrixView,

                                      // D:    const DGeMatrixView

    const auto E = C(_,_(1,2));       // auto: DGeMatrixConstView,

                                      // E: const DGeMatrixConstView

    // auto F = C(_,_(1,2));          // auto: DGeMatrixConstView

                                      // F:    DGeMatrixConstView

    // cout << "F(1,1) = " << F(1,1)  // error: static assertion failure as you

    //      << endl;                  //        call the non-const operator.

}

Why you need 3 `GeMatrix` Types

We have essentially 3 types of GeMatrix<FS> as FS can be FullStorage, FullStorageView or ConstFullStorageView. Technically we even have more, using different storage orders, element types, etc. lead to even more types. However, using generic programming all these types can be used in a unified ways. So FLENS users can think of a single matrix type. And a matrix view just acts like C++ references of buit-in types:

C++ Built-in Types	Analogous Role of FLENS Types
int	`GeMatrix<FullStorage<double, ColMajor> >`
int &	`GeMatrix<FullStorage<double, ColMajor> >::View`
const int &	`GeMatrix<FullStorage<double, ColMajor> >::ConstView`

Just as you technically have 3 types int, int & and const int & you have 3 analogous types for GeMatrix`. Now recall how these are realized in C++:

int i;	Variable `i` is a integer value. So actual the data you want to use.
int &j = i;	`j` is a reference/alias of `i`. Its internally a pointer to the address of `i`. The reference just does the dereferencing of the pointer automatically. Also, the internal pointer will always point to the same address. So in C that's a `int * const`.
const int &k = i;	A const reference is like a `const int * const` in C that automatically gets dereferenced. You can not change the content of `i` through `k`. And `k` internally always points to the address of `i`.

Matrix Assignments are neither a Deep or Shallow copies!

You have to distinguish between creating a new matrices, i.e.

    DGeMatrix  A(5,4);

    DGeMatrixConstView  B = A(_(2,3),_);

and assignments, i.e.

    B         = C;

    A(_(4,5)) = B;

In an assignment you have existing matrices on both sides of the equal sign. And you always copy the value of the right hand side to the left hand side. That's a simple rule as it is exactly what you expect.

There is no way in FLENS that a matrix instance becomes an alias of another matrix. If you want an alias you have to create one. After creation it remains an alias throughout its lifetime. Again, just like for references in C++.

Content

Full Storage Schemes and General Matrices

Full Storage Schemes

'GeMatrix': General Matrix with Full Storage

Some Public Typedefs

Some Methods and operations for Matrix Manipulation

Simple Example for using GeMatrix

Example Code

Comments on Example Code

Compile and Run

Simple Example for using GeMatrix Views

Example Code

Comments on Example Code

Compile and Run

Think in C, Write in C++

How would you do it in C?

Compile and Run

More on Selecting Rectangular Parts from GeMatrix

Interface in Detail: Creating Views Referencing Matrix Parts

Matrix-View Gotchas

Example Code

Comments on Example Code

Compile and Run (Yes, this should fail!!)

Using auto when dealing with GeMatrix Views

Example

Why you need 3 GeMatrix Types

Matrix Assignments are neither a Deep or Shallow copies!

Simple Example for using `GeMatrix`

Simple Example for using `GeMatrix` Views

More on Selecting Rectangular Parts from `GeMatrix`

Using `auto` when dealing with `GeMatrix` Views

Why you need 3 `GeMatrix` Types