Content |
Dense Vectors
Dense vectors have their elements stored linearly in memory separated from each other by a constant stride. If you create a dense vector that allocates its own memory then its stride will equal one. Strides different from one naturally arise when dealing with vector views. Vector views can reference vector slices or rows/columns of matrices.
First Steps
We start with the typical task: setting up a vector, initializing it, manipulating it and printing it.
Example Code
#include <iostream>
using namespace flens;
using namespace std;
int
main()
{
typedef DenseVector<Array<double> > DenseVector;
typedef DenseVector::IndexType IndexType;
DenseVector x(4);
x = 1, 2, 3, 4;
cout << "x.range() = " << x.range() << endl;
cout << "x.length() = " << x.length() << endl;
cout << "x = " << x << endl;
for (IndexType i=x.firstIndex(); i<=x.lastIndex(); ++i) {
x(i) = i*i;
}
const Underscore<IndexType> _;
DenseVector::View y = x(_(2,4));
y = 666;
DenseVector::NoView z = x(_(1,2));
z = 42;
cout << "x = " << x << endl;
cout << "y = " << y << endl;
cout << "z = " << z << endl;
}
Comments on Example Code
The storage scheme for DenseVector is Array, ArrayView or ConstArrayView. We setup a typedef for a “regular” dense vector that allocates its own memory.
typedef DenseVector::IndexType IndexType;
We define a dense vector of length 4 and initialize it.
x = 1, 2, 3, 4;
We retrieve some information like index range and vector length. You will see that the default index base is One.
cout << "x.length() = " << x.length() << endl;
We print the vector.
We iterate through the vector and reset its elements (such that \(x_i = i^2\)).
x(i) = i*i;
}
For selecting vector parts we define a range operator
We create a vector view y that references a part of x.
y = 666;
Here we create a new vector z that allocates its own memory and initializes it with a copy of a part of x. So z is no vector view.
z = 42;
Compile and Run
$shell> cd flens/examples $shell> clang++ -Wall -std=c++0x -I../.. tut01-page04-example01.cc $shell> ./a.out x.range() = [1,4] x.length() = 4 x = 1 2 3 4 x = 1 666 666 666 y = 666 666 666 z = 42 42
Vectors Referencing Matrix Rows, Columns and Diagonals
Referencing rows, columns or diagonals of a matrix is a common task in numerical algorithms. So here we show how to do this for general matrices where the corresponding vector views are dense vectors.
Example Code
#include <iostream>
using namespace flens;
using namespace std;
int
main()
{
typedef GeMatrix<FullStorage<double, ColMajor> > GeMatrix;
typedef DenseVector<Array<double> > DenseVector;
typedef DenseVector::IndexType IndexType;
const Underscore<IndexType> _;
GeMatrix A(4, 5);
A = 1, 2, 3, 4, 5,
6, 7, 8, 9, 10,
11, 12, 13, 14, 15,
16, 17, 18, 19, 20;
DenseVector::View y = A(2,_);
y = 666;
DenseVector::NoView z = A(3,_);
z = 42;
cout << "A = " << A << endl;
cout << "y = " << y << endl;
cout << "z = " << z << endl;
auto x = A(_(2,4),2);
x = 999;
auto d_1 = A.diag(-1);
auto d0 = A.diag(0);
auto d1 = A.diag(1);
auto d2 = A.diag(2);
cout << "d_1 = " << d_1 << endl;
cout << "d0 = " << d0 << endl;
cout << "d1 = " << d1 << endl;
cout << "d2 = " << d2 << endl;
d0 = -1;
d1 = -100, -200, -300, -400;
cout << "A = " << A << endl;
cout << "x = " << x << endl;
cout << "A.leadingDimension() = " << A.leadingDimension() << endl;
cout << "x.stride() = " << x.stride() << endl;
cout << "y.stride() = " << y.stride() << endl;
cout << "z.stride() = " << z.stride() << endl;
cout << "d_1.stride() = " << d_1.stride() << endl;
cout << "d0.stride() = " << d0.stride() << endl;
cout << "d1.stride() = " << d1.stride() << endl;
cout << "d2.stride() = " << d2.stride() << endl;
}
Comments on Example Code
We define typedefs for dense vectors and general matrices that allocate their own memory.
typedef DenseVector<Array<double> > DenseVector;
typedef DenseVector::IndexType IndexType;
For selecting vector parts we define a range operator
We setup some matrix.
A = 1, 2, 3, 4, 5,
6, 7, 8, 9, 10,
11, 12, 13, 14, 15,
16, 17, 18, 19, 20;
We create a vector view y that references the second row of A. Then we reset all elements of y to 666. This will also change elements in A.
y = 666;
We create a “regular” vector z and initialize it with the third row of A. Just to show that z is not a view we then reset all elements of z to 42. You will see that this will not change elements in A.
z = 42;
Of course we also can use the auto type instead of writing out the DenseVector::View type.
x = 999;
You also can reference the k-th off-diagonal of A. For \(k=0\) you get the main diagonal, for \(k>0\) the k-th super-diagonal and for \(k<0\) the corresponding sub-diagonal.
auto d0 = A.diag(0);
auto d1 = A.diag(1);
auto d2 = A.diag(2);
Next we set all elements on the diagonal to -1 and all elements on the first super-diagonal to -100, ..., -400:
d1 = -100, -200, -300, -400;
We finally print the strides of the vectors. Compare strides of vectors referencing rows, columns and diagonals.
cout << "x.stride() = " << x.stride() << endl;
cout << "y.stride() = " << y.stride() << endl;
cout << "z.stride() = " << z.stride() << endl;
cout << "d_1.stride() = " << d_1.stride() << endl;
cout << "d0.stride() = " << d0.stride() << endl;
cout << "d1.stride() = " << d1.stride() << endl;
cout << "d2.stride() = " << d2.stride() << endl;
}
Compile and Run
$shell> cd flens/examples $shell> clang++ -Wall -std=c++0x -I../.. tut01-page04-example02.cc $shell> ./a.out A = 1 2 3 4 5 666 666 666 666 666 11 12 13 14 15 16 17 18 19 20 y = 666 666 666 666 666 z = 42 42 42 42 42 d_1 = 666 999 18 d0 = 1 999 13 19 d1 = 2 666 14 20 d2 = 3 666 15 A = -1 -100 3 4 5 666 -1 -200 666 666 11 999 -1 -300 15 16 999 18 -1 -400 x = -1 999 999 A.leadingDimension() = 4 x.stride() = 1 y.stride() = 4 z.stride() = 1 d_1.stride() = 5 d0.stride() = 5 d1.stride() = 5 d2.stride() = 5