Home     Back     Up     Next Browse Files

Content

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

flens/examples/tut01-page03-example01.cc
#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,
        9, 10, 11, 12;

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

    const Underscore<GeMatrix::IndexType>   _;

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

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

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

    B = 111, 222, 333,
        444, 555, 666;

    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;
}

Comments on Example Code

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,
        9, 10, 11, 12;

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 = 111, 222, 333,
        444, 555, 666;

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

flens/examples/tut01-page03-example02.cc
#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,
                         9, 10, 11, 12,
                        13, 14, 15, 16};

    typedef FullStorageView<double, ColMajor>  FSView;
    typedef GeMatrix<FSView>                   GeMatrixView;

    GeMatrixView  A = FSView(4, 4, data, 4);

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

    GeMatrixView  B = FSView(3, 2, data+5, 4);

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

    typedef GeMatrix<FullStorage<double, ColMajor> >  GeMatrixNoView;
    GeMatrixNoView  M = GeMatrixView(FSView(4, 4, 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,
                         9, 10, 11, 12,
                        13, 14, 15, 16};

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:

The syntax for the full storage view is FSView(numRows, numCols, data-pointer, leading dimension).

    GeMatrixView  A = FSView(4, 4, 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(3, 2, data+5, 4);

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(4, 4, 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:

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

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.

flens/examples/tut01-page03-example03.cc
#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,
        9, 10, 11, 12;

    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)));
}

Comments on the Code

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,
        9, 10, 11, 12;

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

flens/examples/tut01-page03-example04.cc
#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,
        9, 10, 11, 12;

    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,
        9, 10, 11, 12;

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.

Home     Back     Up     Next