# Matrix Views

## Matrix Views Referencing Matrix Parts

If you have setup a general matrix then you can use matrix views to reference matrix parts. Matrix views act like references in C++. Just like you can have constant references and non-constant references in C++ you can have constant matrix views and non-constant matrix views.

### Example Code

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

using namespace flens;
using namespace std;

int
main()
{
typedef FullStorage<double, ColMajor>  FS;
typedef GeMatrix<FS>                   GeMatrix;

typedef GeMatrix::ConstView            GeMatrixConstView;
typedef GeMatrix::View                 GeMatrixView;

GeMatrix  A(3,4);
A = 1,  2,  3,  4,
5,  6,  7,  8,
9101112;

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

const Underscore<GeMatrix::IndexType>   _;

GeMatrixView       B = A(_(2,3),_(2,4));

//
//  error:
//  GeMatrixView       B;
//

GeMatrixConstView  C = A(_,_(2,4));

B = 111222333,
444555666;

cout << "Changed B:" << endl;
cout << "A = " << A << endl;
cout << "B = " << B << endl;
cout << "C = " << C << endl;

B(1,1) = 999;

cout << "Changed B(1,1):" << endl;
cout << "A = " << A << endl;
cout << "B = " << B << endl;
cout << "C = " << C << endl;

A(_,_(1,2)) = A(_,_(3,4));

cout << "After copying A(:,1:2) = A(:,3:4)" << endl;
cout << "A = " << A << endl;

//
//  error:
//  C(1,1) = 3;
//  C = 1, 2,
//      3, 4,
//      5, 6;
}

We start with a typedef for a GeMatrix with elements of type double that are stored column-wise in memory.

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

Then we define matrix view types for referencing matrix parts.

typedef GeMatrix::ConstView            GeMatrixConstView;
typedef GeMatrix::View                 GeMatrixView;

We setup some matrix.

GeMatrix  A(3,4);
A = 1,  2,  3,  4,
5,  6,  7,  8,
9101112;

We also define a “range-operator” to ease the selection of matrix parts. for this purpose we simply define an object called “_” of type Underscore<I> where I denotes the type of index used in the matrix.

const Underscore<GeMatrix::IndexType>   _;

Now we let B reference the matrix part A(2:3,2:4). The range operator lets you specify row/column ranges by _(from, to) similar to the Matlab notation from:to.

GeMatrixView       B = A(_(2,3),_(2,4));

Note that you can not define something like an uninitialized matrix view. Just like you can not define uninitialized references in C++. I.e. this would give you an error

//
//  error:
//  GeMatrixView       B;
//

Next we define a matrix view that provides read-only access to A(:,2,4). So by the way, notice that you can also select all rows if you do not pass arguments to the range operator “_”. Analogously you can select all columns of course.

GeMatrixConstView  C = A(_,_(2,4));

We overwrite B and by that A

B = 111222333,
444555666;

Note that B acts like a “regular” matrix and indices also start at One.

B(1,1) = 999;

Of course you also can use unnamed matrix views:

A(_,_(1,2)) = A(_,_(3,4));

And obviously it should not be possible to write to C. Indeed this would give you a compile-time error:

//
//  error:
//  C(1,1) = 3;
//  C = 1, 2,
//      3, 4,
//      5, 6;
}

### Compile and Run

$shell> cd flens/examples$shell> g++ -Wall -std=c++11 -I../.. tut01-page03-example01.cc
$shell> ./a.out A = 1 2 3 4 5 6 7 8 9 10 11 12 Changed B: A = 1 2 3 4 5 111 222 333 9 444 555 666 B = 111 222 333 444 555 666 C = 2 3 4 111 222 333 444 555 666 Changed B(1,1): A = 1 2 3 4 5 999 222 333 9 444 555 666 B = 999 222 333 444 555 666 C = 2 3 4 999 222 333 444 555 666 After copying A(:,1:2) = A(:,3:4) A = 3 4 3 4 222 333 222 333 555 666 555 666  ## Creating a Matrix View from a C-Array If you have a plain C-array you can create a matrix view that references this C-array. That means you wrap the C-array into a FLENS matrix without copying any data. This trick is essential for interfacing FLENS with other libraries. ### Example Code #include <flens/flens.cxx> #include <iostream> using namespace flens; using namespace std; int main() { double data[4*4] = { 1, 2, 3, 4, 5, 6, 7, 8, 9101112, 13141516}; typedef FullStorageView<double, ColMajor> FSView; typedef GeMatrix<FSView> GeMatrixView; GeMatrixView A = FSView(44, data, 4); cout << "A = " << A << endl; GeMatrixView B = FSView(32, data+54); cout << "B = " << B << endl; typedef GeMatrix<FullStorage<double, ColMajor> > GeMatrixNoView; GeMatrixNoView M = GeMatrixView(FSView(44, data, 4)); // GeMatrixNoView M = FSView(4, 4, data, 4); // error: This would not compile because // GeMatrixNoView::Engine is of type "FullStorage<...>" // FSView is of type "FullStorageView<...>" M(2,2) = -666; B(1,2) = 666; cout << "now: A = " << A << endl; cout << "Data dump of C-array: " << endl; for (int i=0; i<16; ++i) { cout << data[i] << " "; } cout << endl; return 0; } ### Comments on Example Code So this is our plain C-array double data[4*4] = { 1, 2, 3, 4, 5, 6, 7, 8, 9101112, 13141516}; We define one typedef for a storage scheme that only references data and another typedef for a corresponding general matrix type. The trick is that you can construct GeMatrix objects from full storage objects. So what we later do is building a storage object that just keeps internally a pointer to the C-array. Note: The type GeMatrixView is actually identicall with the one from the previous example (there we defined it via GeMatrix::View). typedef FullStorageView<double, ColMajor> FSView; typedef GeMatrix<FSView> GeMatrixView; We finally create a matrix view that references the C-array as follows: • We wrap the C-array in an anonymous full storage view object • and construct with it a general matrix view. The syntax for the full storage view is FSView(numRows, numCols, data-pointer, leading dimension). GeMatrixView A = FSView(44, data, 4); With a little pointer arithmetic we also can create views for a sub-matrix directly (note that data+5points to the element of A(2,3). GeMatrixView B = FSView(32, data+54); Note that in the following the data of the C-array gets copied into matrix M. This is because GeMatrixNoView has a non-view storage scheme. typedef GeMatrix<FullStorage<double, ColMajor> > GeMatrixNoView; GeMatrixNoView M = GeMatrixView(FSView(44, data, 4)); Note that you only can construct a GeMatrix from a storage object if the type of the storage object is identical with GeMatrix::Engine (Engine is the storage scheme used by GeMatrix). Hence, the following would not compile: // GeMatrixNoView M = FSView(4, 4, data, 4); // error: This would not compile because // GeMatrixNoView::Engine is of type "FullStorage<...>" // FSView is of type "FullStorageView<...>" Let us demonstrate what is a view and what not: • We change M(2,2), • we change B(1,2), • we output A. Only the change of B(1,2) affects elements of A as both matrices share data. M(2,2) = -666; B(1,2) = 666; cout << "now: A = " << A << endl; cout << "Data dump of C-array: " << endl; for (int i=0; i<16; ++i) { cout << data[i] << " "; } cout << endl; ### Compile and Run $shell> cd flens/examples
$shell> g++ -Wall -std=c++11 -I../.. tut01-page03-example02.cc$shell> ./a.out
A =
1             5             9            13
2             6            10            14
3             7            11            15
4             8            12            16
B =
6            10
7            11
8            12
now: A =
1             5             9            13
2             6           666            14
3             7            11            15
4             8            12            16
Data dump of C-array:
1  2  3  4  5  6  7  8  9  666  11  12  13  14  15  16


## More on Matrix View Types

You saw that there are different ways to define a matrix view types

• Get the matrix view types through typedefs in GeMatrix, e.g. by GeMatrix<MA>::View or GeMatrix<MA>::ConstView were MA denotes some full storage scheme.

typedef FullStorage<double, ColMajor>  MA;
typedef GeMatrix<MA>             NoView;

typedef GeMatrix<MA>::NoView     AlsoNoView;
typedef GeMatrix<MA>::View       View;
typedef GeMatrix<MA>::ConstView  ConstView;

Note that NoView and AlsoNoView are identical.

• You can achieve the same by directly defining different storage types for GeMatrix:

typedef FullStorage<double, ColMajor>           MA;
typedef FullStorageView<double, ColMajor>       MAView;
typedef CosntFullStorageView<double, ColMajor>  MAConstView;

typedef GeMatrix<MA>           NoView;
typedef GeMatrix<MAView>       View;
typedef GeMatrix<MAConstView>  ConstView;
• As you can see: (Regular) matrices, matrix views and const matrix views are just GeMatrix types with different stroage schemes:

• FullStorage is a storage scheme that actually allocates memory in its constructor and releases it in its destructor.

• FullStorageView does not allocate or release memory but besides that acts like FullStorage and

• ConstFullStorage acts like FullStorageView but only provides read-only methods and operators.

• In an assignment

A = B;

(were A and B are general matrices) data gets always copied (except A and B are identical objects). More precise, the memory referenced by B gets copied to the memory referenced by A. Again, this is analogously to references in C++.

• If you create a matrix view from an object A of type GeMatrix<MA> it will either be of type GeMatrix<MA>::View or GeMatrix<MA>::NoView:

• If A is in const context any matrix view created from it will be of type GeMatrix<MA>::ConstView.

• If A is in non-const context any matrix view created from it will be of type GeMatrix<MA>::View.

Attention: If a const matrix view is in non-const context and you create a view from it it would be a non-const matrix view. Obviously this would not make sense and give you a compile time error. We give an example for this below and also show how to prevent such a problem.

### Using the auto keyword

The auto keyword can make it easier to work with the matrix views. This is because it can save typing. However, you should always know what view type gets actually represented through an auto type.

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

using namespace flens;
using namespace std;

template <typename MA>
void
dummy(const GeMatrix<MA> &A)
{

typedef typename GeMatrix<MA>::IndexType  IndexType;
const Underscore<IndexType>  _;

const IndexType n = A.numCols();
const IndexType k = n/2;

const auto A1 = A(_,_(1,k));
const auto A2 = A(_,_(k+1,n));

cout << "A1 = " << A1 << endl;
cout << "A2 = " << A2 << endl;
}

int
main()
{
typedef FullStorage<double, ColMajor>  FS;
typedef GeMatrix<FS>                   GeMatrix;
const Underscore<GeMatrix::IndexType>  _;

GeMatrix  A(3,4);
A = 1,  2,  3,  4,
5,  6,  7,  8,
9101112;

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

dummy(A);

dummy(A(_(1,3),_));

auto B = A(_(1,3),_);
dummy(B);

const auto C = A(_(1,3),_);
dummy(C(_,_(1,4)));
}

This function receives a general matrix and MA can be of type FullStorage, FullStorageView or ConstFullStorageView. Note that the function receives the matrix as a const reference.

template <typename MA>
void
dummy(const GeMatrix<MA> &A)
{

We setup a range operator for this matrix type ...

typedef typename GeMatrix<MA>::IndexType  IndexType;
const Underscore<IndexType>  _;

... and get some numbers we will use to split the matrix vertically.

const IndexType n = A.numCols();
const IndexType k = n/2;

As A is in this function in const context the created matrix views will always be of type GeMatrix<ConstFullStorageView<..> >. And this is also the type represented by auto. So A1 and A2 will be of type const GeMatrix<ConstFullStorageView<..> >.

const auto A1 = A(_,_(1,k));
const auto A2 = A(_,_(k+1,n));

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 FullStorage<double, ColMajor>  FS;
typedef GeMatrix<FS>                   GeMatrix;
const Underscore<GeMatrix::IndexType>  _;

Then we setup some matrix.

GeMatrix  A(3,4);
A = 1,  2,  3,  4,
5,  6,  7,  8,
9101112;

Now we call dummy with different view types.

Call it with matrix of type GeMatrix<FullStorage<..> >

dummy(A);

Call it with matrix of type GeMatrix<FullStorageView<..> >

dummy(A(_(1,3),_));

Call it again with matrix of type GeMatrix<FullStorageView<..> >

auto B = A(_(1,3),_);
dummy(B);

Call it with matrix of type GeMatrix<ConstFullStorageView<..> >. Here C is of type const GeMatrix<FullStorageView>, i.e. a (non-const) matrix view in const context.

const auto C = A(_(1,3),_);
dummy(C(_,_(1,4)));
}

### Compile and Run

$shell> cd flens/examples$shell> g++ -Wall -std=c++11 -I../.. tut01-page03-example03.cc
$shell> ./a.out A = 1 2 3 4 5 6 7 8 9 10 11 12 A1 = 1 2 5 6 9 10 A2 = 3 4 7 8 11 12 A1 = 1 2 5 6 9 10 A2 = 3 4 7 8 11 12 A1 = 1 2 5 6 9 10 A2 = 3 4 7 8 11 12 A1 = 1 2 5 6 9 10 A2 = 3 4 7 8 11 12  ## More on auto and Const Context In a shot: Always make sure that a const matrix view is in const context. ### Example Code #include <flens/flens.cxx> #include <iostream> using namespace flens; using namespace std; int main() { typedef FullStorage<double, ColMajor> FS; typedef GeMatrix<FS> GeMatrix; const Underscore<GeMatrix::IndexType> _; GeMatrix A(3,4); A = 1, 2, 3, 4, 5, 6, 7, 8, 9101112; cout << "A = " << A << endl; const auto B = A(_(1,3),_); auto C = B(_(1,3),_); auto D = C(_(1,3),_); D(1,1) = 666; cout << "D = " << D << endl; } ### Comments on the 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 FullStorage<double, ColMajor> FS; typedef GeMatrix<FS> GeMatrix; const Underscore<GeMatrix::IndexType> _; Then we setup some matrix. GeMatrix A(3,4); A = 1, 2, 3, 4, 5, 6, 7, 8, 9101112; Matrix B is of type const GeMatrix::View. So in particular this matrix is in const context. const auto B = A(_(1,3),_); Matrix C is of type GeMatrix::ConstView but it is not in const context. auto C = B(_(1,3),_); Matrix D would be of type GeMatrix::View and would reference a part of a matrix of type GeMatrix::ConstView. This would violate const correctness. auto D = C(_(1,3),_); ### Compile and Analyze Compile Error $shell> cd flens/examples
\$shell> g++ -Wall -std=c++11 -I../.. tut01-page03-example04.cc
In file included from ../../flens/matrixtypes/general/impl/impl.tcc:40:0,
from ../../flens/matrixtypes/matrixtypes.tcc:38,
from ../../flens/flens.tcc:41,
from ../../flens/flens.cxx:53,
from tut01-page03-example04.cc:1:
../../flens/matrixtypes/general/impl/gematrix.tcc: In instantiation of 'flens::GeMatrix::View flens::GeMatrix::operator()(const flens::Range&, const flens::Underscore&) [with FS = flens::ConstFullStorageView, std::allocator >; flens::GeMatrix::View = flens::GeMatrix, std::allocator > >; typename FS::View = flens::FullStorageView, std::allocator >; typename FS::IndexType = int]':
tut01-page03-example04.cc:45:24:   required from here
../../flens/matrixtypes/general/impl/gematrix.tcc:631:53: error: could not convert 'flens::ConstFullStorageView::view(flens::ConstFullStorageView::IndexType, flens::ConstFullStorageView::IndexType, flens::ConstFullStorageView::IndexType, flens::ConstFullStorageView::IndexType, flens::ConstFullStorageView::IndexType, flens::ConstFullStorageView::IndexType) const [with T = double; cxxblas::StorageOrder Order = (cxxblas::StorageOrder)1u; I = flens::IndexOptions<>; A = std::allocator; flens::ConstFullStorageView::IndexType = int]((& rows)->flens::Range::firstIndex(), flens::GeMatrix::firstCol, std::allocator > >(), (& rows)->flens::Range::lastIndex(), flens::GeMatrix::lastCol, std::allocator > >(), 1, 1)' from 'const flens::ConstFullStorageView, std::allocator >' to 'flens::GeMatrix, std::allocator > >::View {aka flens::GeMatrix, std::allocator > >}'
../../flens/matrixtypes/general/impl/gematrix.tcc: In member function 'flens::GeMatrix::View flens::GeMatrix::operator()(const flens::Range&, const flens::Underscore&) [with FS = flens::ConstFullStorageView, std::allocator >; flens::GeMatrix::View = flens::GeMatrix, std::allocator > >; typename FS::View = flens::FullStorageView, std::allocator >; typename FS::IndexType = int]':
../../flens/matrixtypes/general/impl/gematrix.tcc:632:1: warning: control reaches end of non-void function [-Wreturn-type]


The error message basically says that a conversion from ConstFullStorageView to FullStorageView (aka GeMatrix::View) would be require. As this is not possible we get a compile time error.