getting matrix diagonal

Klaus Iglberger

changed status to resolved

Summary

The feature has been implemented, tested, optimized (including parallelization) and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 3.3.

Band Views

Bands provide views on a specific band of a dense or sparse matrix. As such, bands act as a reference to a specific band. This reference is valid and can be used in every way any other vector can be used as long as the matrix containing the band is not resized or entirely destroyed. The band also acts as an alias to the band elements: Changes made to the elements (e.g. modifying values, inserting or erasing elements) are immediately visible in the matrix and changes made via the matrix are immediately visible in the band.

The Band Template

The blaze::Band template represents a reference to a specific band of a dense or sparse matrix primitive. It can be included via the header file

#include <blaze/math/Band.h>

The type of the matrix is specified via template parameter:

template< typename MT >
class Band;

MT specifies the type of the matrix primitive. Band can be used with every matrix primitive, but does not work with any matrix expression type.

Setup of Bands

A reference to a dense or sparse band can be created very conveniently via the band() function. The band index must be in the range from [1-M..N-1], where M is the total number of rows and N is the total number of columns, and can be specified both at compile time or at runtime:

using DenseMatrixType = blaze::DynamicMatrix<double,blaze::rowMajor>;

DenseMatrixType A;
// ... Resizing and initialization

// Creating a reference to the 1st lower band of matrix A (compile time index)
blaze::Band<DenseMatrixType,-1L> band1 = band<-1L>( A );

// Creating a reference to the 2nd upper band of matrix A (runtime index)
blaze::Band<DenseMatrixType> band2 = band( A, 2L );

In addition, the diagonal() function provides a convenient shortcut for the setup of a view on the diagonal of a dense or sparse matrix. It has the same effect as calling the band() function with a compile time index of 0:

using DenseMatrixType = blaze::DynamicMatrix<double,blaze::rowMajor>;
using DiagonalType1   = blaze::Band<DenseMatrixType,0L>;
using DiagonalType2   = blaze::Diagonal<DenseMatrixType>;

DenseMatrixType A;
// ... Resizing and initialization

// Creating a reference to the diagonal of matrix A via the band() and diagonal() functions
DiagonalType1 diag1 = band<0L>( A );
DiagonalType2 diag2 = diagonal( A );

static_assert( blaze::IsSame<DiagonalType1,DiagonalType2>::value, "Non-identical types detected" );

The resulting reference can be treated as any other vector, i.e. it can be assigned to, it can be copied from, and it can be used in arithmetic operations. By default, bands are considered column vectors, but this setting can be changed via the defaultTransposeFlag switch. The reference can also be used on both sides of an assignment: The band can either be used as an alias to grant write access to a specific band of a matrix primitive on the left-hand side of an assignment or to grant read-access to a specific band of a matrix primitive or expression on the right-hand side of an assignment. The following example demonstrates this in detail:

using DenseVectorType  = blaze::DynamicVector<double,blaze::rowVector>;
using SparseVectorType = blaze::CompressedVector<double,blaze::rowVector>;
using DenseMatrixType  = blaze::DynamicMatrix<double,blaze::rowMajor>;
using SparseMatrixType = blaze::CompressedMatrix<double,blaze::rowMajor>;

DenseVectorType  x;
SparseVectorType y;
DenseMatrixType  A, B;
SparseMatrixType C, D;
// ... Resizing and initialization

// Setting the 2nd upper band of matrix A to x
blaze::Band<DenseMatrixType> band2 = band( A, 2L );
band2 = x;

// Setting the 3rd upper band of matrix B to y
band( B, 3L ) = y;

// Setting x to the 2nd lower band of the result of the matrix multiplication
x = band( A * B, -2L );

// Setting y to the 2nd upper band of the result of the sparse matrix multiplication
y = band( C * D, 2L );

The band() function can be used on any dense or sparse matrix, including expressions, as illustrated by the source code example. However, bands cannot be instantiated for expression types, but only for matrix primitives, respectively, i.e. for matrix types that offer write access.

Element Access

A dense or sparse band can be used like any other vector. For instance, the elements of a band can be directly accessed with the subscript operator:

using MatrixType = blaze::DynamicMatrix<double,blaze::rowMajor>;
MatrixType A;
// ... Resizing and initialization

// Creating a view on the 4th upper band of matrix A
blaze::Band<MatrixType> band4 = band( A, 4L );

// Setting the 1st element of the dense band, which corresponds
// to the 1st element in the 4th upper band of matrix A
band4[1] = 2.0;

Alternatively, the elements of a band can be traversed via iterators. Just as with vectors, in case of non-const band, begin() and end() return an Iterator, which allows a manipulation of the non-zero values, in case of constant band a ConstIterator is returned:

using MatrixType = blaze::DynamicMatrix<int,blaze::rowMajor>;
using BandType   = blaze::Band<MatrixType>;

MatrixType A( 128UL, 256UL );
// ... Resizing and initialization

// Creating a reference to the 5th upper band of matrix A
BandType band5 = band( A, 5L );

for( BandType::Iterator it=band5.begin(); it!=band5.end(); ++it ) {
   *it = ...;  // OK; Write access to the dense band value
   ... = *it;  // OK: Read access to the dense band value.
}

for( BandType::ConstIterator it=band5.begin(); it!=band5.end(); ++it ) {
   *it = ...;  // Compilation error: Assignment to the value via a ConstIterator is invalid.
   ... = *it;  // OK: Read access to the dense band value.
}

using MatrixType = blaze::CompressedMatrix<int,blaze::rowMajor>;
using BandType   = blaze::Band<MatrixType>;

MatrixType A( 128UL, 256UL );
// ... Resizing and initialization

// Creating a reference to the 5th band of matrix A
BandType band5 = band( A, 5L );

for( BandType::Iterator it=band5.begin(); it!=band5.end(); ++it ) {
   it->value() = ...;  // OK: Write access to the value of the non-zero element.
   ... = it->value();  // OK: Read access to the value of the non-zero element.
   it->index() = ...;  // Compilation error: The index of a non-zero element cannot be changed.
   ... = it->index();  // OK: Read access to the index of the sparse element.
}

for( BandType::ConstIterator it=band5.begin(); it!=band5.end(); ++it ) {
   it->value() = ...;  // Compilation error: Assignment to the value via a ConstIterator is invalid.
   ... = it->value();  // OK: Read access to the value of the non-zero element.
   it->index() = ...;  // Compilation error: The index of a non-zero element cannot be changed.
   ... = it->index();  // OK: Read access to the index of the sparse element.
}

Element Insertion

Inserting/accessing elements in a sparse band can be done by several alternative functions. The following example demonstrates all options:

using MatrixType = blaze::CompressedMatrix<double,blaze::rowMajor>;
MatrixType A( 10UL, 100UL );  // Non-initialized 10x100 matrix

using BandType = blaze::Band<MatrixType>;
BandType diag( band( A, 0L ) );  // Reference to the diagonal of A

// The subscript operator provides access to all possible elements of the sparse band,
// including the zero elements. In case the subscript operator is used to access an element
// that is currently not stored in the sparse band, the element is inserted into the band.
diag[42] = 2.0;

// The second operation for inserting elements is the set() function. In case the element
// is not contained in the band it is inserted into the band, if it is already contained in
// the band its value is modified.
diag.set( 45UL, -1.2 );

// An alternative for inserting elements into the band is the insert() function. However,
// it inserts the element only in case the element is not already contained in the band.
diag.insert( 50UL, 3.7 );

Common Operations

A band view can be used like any other column vector. This means that with only a few exceptions all basic vector operations and arithmetic operations can be used. For instance, the current number of band elements can be obtained via the size() function, the current capacity via the capacity() function, and the number of non-zero elements via the nonZeros() function. However, since bands are references to specific bands of a matrix, several operations are not possible, such as resizing and swapping. The following example shows this by means of a dense band view:

using MatrixType = blaze::DynamicMatrix<int,blaze::rowMajor>;
using BandType = blaze::Band<MatrixType>;

MatrixType A( 42UL, 42UL );
// ... Resizing and initialization

// Creating a reference to the 2nd upper band of matrix A
BandType band2 = band( A, 2L );

band2.size();          // Returns the number of elements in the band
band2.capacity();      // Returns the capacity of the band
band2.nonZeros();      // Returns the number of non-zero elements contained in the band

band2.resize( 84UL );  // Compilation error: Cannot resize a single band of a matrix

BandType band3 = band( A, 3L );
swap( band2, band3 );   // Compilation error: Swap operation not allowed

Arithmetic Operations

Both dense and sparse bands can be used in all arithmetic operations that any other dense or sparse vector can be used in. The following example gives an impression of the use of dense bands within arithmetic operations. All operations (addition, subtraction, multiplication, scaling, ...) can be performed on all possible combinations of dense and sparse bands with fitting element types:

blaze::DynamicVector<double,blaze::columnVector> a( 2UL, 2.0 ), b;
blaze::CompressedVector<double,blaze::columnVector> c( 2UL );
c[1] = 3.0;

using DenseMatrix = blaze::DynamicMatrix<double,blaze::rowMajor>;
DenseMatrix A( 4UL, 2UL );  // Non-initialized 4x2 matrix

using BandType = blaze::Band<DenseMatrix>;
BandType band1( band( A, 1L ) );  // Reference to the 1st upper band of A
BandType diag ( band( A, 0L ) );  // Reference to the diagonal of A

band1[0] = 0.0;      // Manual initialization of the 1st upper band of A
diag = 1.0;          // Homogeneous initialization of the diagonal of A
band( A, -1L ) = a;  // Dense vector initialization of the 1st lower band of A
band( A, -2L ) = c;  // Sparse vector initialization of the 2nd lower band of A

b = diag + a;               // Dense vector/dense vector addition
b = c + band( A, -1L );     // Sparse vector/dense vector addition
b = diag * band( A, -2L );  // Component-wise vector multiplication

band( A, -1L ) *= 2.0;     // In-place scaling of the 1st upper band
b = band( A, -1L ) * 2.0;  // Scaling of the 1st upper band
b = 2.0 * band( A, -1L );  // Scaling of the 1st upper band

band( A, -2L ) += a;              // Addition assignment
band( A, -2L ) -= c;              // Subtraction assignment
band( A, -2L ) *= band( A, 0L );  // Multiplication assignment

double scalar = trans( c ) * band( A, -1L );  // Scalar/dot/inner product between two vectors

A = band( A, -1L ) * trans( c );  // Outer product between two vectors

2017-08-28T14:44:52+00:00

Comments (4)