7. LvArray::SparsityPattern
and LvArray::CRSMatrix
¶
LvArray::SparsityPattern
represents just the sparsity pattern of a matrix while the LvArray::CRSMatrix
represents a sparse matrix. Both use a slightly modified version of the compressed row storage format. The modifications to the standard format are as follows:
- The columns of each row are sorted.
- Each row can have a capacity in addition to a size which means that the rows aren’t necessarily adjacent in memory.
7.1. Template arguments¶
The LvArray::SparsityPattern
has three template arguments
COL_TYPE
: The integral type used to enumerate the columns of the matrix.INDEX_TYPE
: An integral type used in index calculations, the suggested type isstd::ptrdiff_t
.BUFFER_TYPE
: A template template parameter specifying the buffer type used for allocation and de-allocation, theLvArray::SparsityPattern
contains aBUFFER_TYPE< COL_TYPE >
along with twoBUFFER_TYPE< INDEX_TYPE >
.
The LvArray::CRSMatrix
adds an additional template argument T
which is the type of the entries in the matrix. It also has an addition member of type BUFER_TYPE< T >
.
7.2. Usage¶
The LvArray::SparsityPattern
is just a LvArray::ArrayOfSets
by a different name. The only functional difference is that it doesn’t support inserting or removing rows from the matrix (inserting or removing inner sets). Both the LvArray::SparsityPattern
and LvArray::CRSMatrix
support insertNonZero
and insertNonZeros
for inserting entries into a row as well as removeNonZero
and removeNonZeros
for removing entries from a row. LvArray::CRSMatrix
also supports various addToRow
methods which will add to existing entries in a specific row.
It is worth noting that neither LvArray::SparsityPattern
nor LvArray::CRSMatrix
have an operator()
or operator[]
. Instead they both support getColumns
which returns a LvArray::ArraySlice< COL_TYPE const, 1, 0, INDEX_TYPE >
with the columns of the row and LvArray::CRSMatrix
supports getEntries
which returns a LvArray::ArraySlice< T, 1, 0, INDEX_TYPE >
with the entries of the row.
TEST( CRSMatrix, examples )
{
// Create a matrix with three rows and three columns.
LvArray::CRSMatrix< double, int, std::ptrdiff_t, LvArray::MallocBuffer > matrix( 2, 3 );
EXPECT_EQ( matrix.numRows(), 2 );
EXPECT_EQ( matrix.numColumns(), 3 );
// Insert two entries into the first row.
int const row0Columns[ 2 ] = { 0, 2 };
double const row0Values[ 2 ] = { 4, 3 };
matrix.insertNonZeros( 0, row0Columns, row0Values, 2 );
EXPECT_EQ( matrix.numNonZeros( 0 ), 2 );
// Insert three entries into the second row.
int const row1Columns[ 3 ] = { 0, 1, 2 };
double const row1Values[ 3 ] = { 55, -1, 4 };
matrix.insertNonZeros( 1, row1Columns, row1Values, 3 );
EXPECT_EQ( matrix.numNonZeros( 1 ), 3 );
// The entire matrix has five non zero entries.
EXPECT_EQ( matrix.numNonZeros(), 5 );
// Row 0 does not have an entry for column 1.
EXPECT_TRUE( matrix.empty( 0, 1 ) );
// Row 1 does have an entry for column 1.
EXPECT_FALSE( matrix.empty( 1, 1 ) );
LvArray::ArraySlice< int const, 1, 0, std::ptrdiff_t > columns = matrix.getColumns( 0 );
LvArray::ArraySlice< double, 1, 0, std::ptrdiff_t > entries = matrix.getEntries( 0 );
// Check the entries of the matrix.
EXPECT_EQ( columns.size(), 2 );
EXPECT_EQ( columns[ 0 ], row0Columns[ 0 ] );
EXPECT_EQ( entries[ 0 ], row0Values[ 0 ] );
EXPECT_EQ( columns[ 1 ], row0Columns[ 1 ] );
entries[ 1 ] += 10;
EXPECT_EQ( entries[ 1 ], 10 + row0Values[ 1 ] );
}
[Source: examples/exampleSparsityPatternAndCRSMatrix.cpp]
Third party packages such as MKL or cuSparse expect the rows of the CRS matrix to be adjacent in memory. Specifically they don’t accept a pointer that gives the size of each row they only accept an offsets pointer and calculate the sizes from that. To make an LvArray::CRSMatrix
conform to this layout you can call compress
which will set the capacity of each row equal to its size making the rows adjacent in memory.
TEST( CRSMatrix, compress )
{
// Create a matrix with two rows and four columns where each row has capacity 3.
LvArray::CRSMatrix< double, int, std::ptrdiff_t, LvArray::MallocBuffer > matrix( 2, 4, 3 );
// Insert two entries into the first row.
int const row0Columns[ 2 ] = { 0, 2 };
double const row0Values[ 2 ] = { 4, 3 };
matrix.insertNonZeros( 0, row0Columns, row0Values, 2 );
EXPECT_EQ( matrix.numNonZeros( 0 ), 2 );
EXPECT_EQ( matrix.nonZeroCapacity( 0 ), 3 );
// Insert two entries into the second row.
int const row1Columns[ 3 ] = { 0, 1 };
double const row1Values[ 3 ] = { 55, -1 };
matrix.insertNonZeros( 1, row1Columns, row1Values, 2 );
EXPECT_EQ( matrix.numNonZeros( 1 ), 2 );
EXPECT_EQ( matrix.nonZeroCapacity( 1 ), 3 );
// The rows are not adjacent in memory.
EXPECT_NE( &matrix.getColumns( 0 )[ 0 ] + matrix.numNonZeros( 0 ), &matrix.getColumns( 1 )[ 0 ] );
EXPECT_NE( &matrix.getEntries( 0 )[ 0 ] + matrix.numNonZeros( 0 ), &matrix.getEntries( 1 )[ 0 ] );
// After compression the rows are adjacent in memroy.
matrix.compress();
EXPECT_EQ( &matrix.getColumns( 0 )[ 0 ] + matrix.numNonZeros( 0 ), &matrix.getColumns( 1 )[ 0 ] );
EXPECT_EQ( &matrix.getEntries( 0 )[ 0 ] + matrix.numNonZeros( 0 ), &matrix.getEntries( 1 )[ 0 ] );
EXPECT_EQ( matrix.numNonZeros( 0 ), matrix.nonZeroCapacity( 0 ) );
EXPECT_EQ( matrix.numNonZeros( 1 ), matrix.nonZeroCapacity( 1 ) );
// The entries in the matrix are unchanged.
EXPECT_EQ( matrix.getColumns( 0 )[ 0 ], row0Columns[ 0 ] );
EXPECT_EQ( matrix.getEntries( 0 )[ 0 ], row0Values[ 0 ] );
EXPECT_EQ( matrix.getColumns( 0 )[ 1 ], row0Columns[ 1 ] );
EXPECT_EQ( matrix.getEntries( 0 )[ 1 ], row0Values[ 1 ] );
EXPECT_EQ( matrix.getColumns( 1 )[ 0 ], row1Columns[ 0 ] );
EXPECT_EQ( matrix.getEntries( 1 )[ 0 ], row1Values[ 0 ] );
EXPECT_EQ( matrix.getColumns( 1 )[ 1 ], row1Columns[ 1 ] );
EXPECT_EQ( matrix.getEntries( 1 )[ 1 ], row1Values[ 1 ] );
}
[Source: examples/exampleSparsityPatternAndCRSMatrix.cpp]
LvArray::CRSMatrix
also has an assimilate
method which takes an r-values reference to a LvArray::SparsityPattern
and converts it into a LvArray::CRSMatrix
. It takes a RAJA execution policy as a template parameter.
TEST( CRSMatrix, assimilate )
{
// Create a sparsity pattern with two rows and four columns where each row has capacity 3.
LvArray::SparsityPattern< int, std::ptrdiff_t, LvArray::MallocBuffer > sparsity( 2, 4, 3 );
// Insert two entries into the first row.
int const row0Columns[ 2 ] = { 0, 2 };
sparsity.insertNonZeros( 0, row0Columns, row0Columns + 2 );
// Insert two entries into the second row.
int const row1Columns[ 3 ] = { 0, 1 };
sparsity.insertNonZeros( 1, row1Columns, row1Columns + 2 );
// Create a matrix with the sparsity pattern
LvArray::CRSMatrix< double, int, std::ptrdiff_t, LvArray::MallocBuffer > matrix;
matrix.assimilate< RAJA::loop_exec >( std::move( sparsity ) );
// The sparsity is empty after being assimilated.
EXPECT_EQ( sparsity.numRows(), 0 );
EXPECT_EQ( sparsity.numColumns(), 0 );
// The matrix has the same shape as the sparsity pattern.
EXPECT_EQ( matrix.numRows(), 2 );
EXPECT_EQ( matrix.numColumns(), 4 );
EXPECT_EQ( matrix.numNonZeros( 0 ), 2 );
EXPECT_EQ( matrix.nonZeroCapacity( 0 ), 3 );
EXPECT_EQ( matrix.numNonZeros( 1 ), 2 );
EXPECT_EQ( matrix.nonZeroCapacity( 1 ), 3 );
// The entries in the matrix are zero initialized.
EXPECT_EQ( matrix.getColumns( 0 )[ 0 ], row0Columns[ 0 ] );
EXPECT_EQ( matrix.getEntries( 0 )[ 0 ], 0 );
EXPECT_EQ( matrix.getColumns( 0 )[ 1 ], row0Columns[ 1 ] );
EXPECT_EQ( matrix.getEntries( 0 )[ 1 ], 0 );
EXPECT_EQ( matrix.getColumns( 1 )[ 0 ], row1Columns[ 0 ] );
EXPECT_EQ( matrix.getEntries( 1 )[ 0 ], 0 );
EXPECT_EQ( matrix.getColumns( 1 )[ 1 ], row1Columns[ 1 ] );
EXPECT_EQ( matrix.getEntries( 1 )[ 1 ], 0 );
}
[Source: examples/exampleSparsityPatternAndCRSMatrix.cpp]
7.3. LvArray::CRSMatrixView
¶
The LvArray::SparsityPatternView
and LvArray::CRSMatrixView
are the view counterparts to LvArray::SparsityPattern
and LvArray::CRSMatrix
. They both have the same template arguments as their counterparts. The LvArray::SparsityPatternView
behaves exactly like an LvArray::ArrayOfSetsView
. From a LvArray::CRSMatrix< T, COL_TYPE, INDEX_TYPE, BUFFER_TYPE >
there are four view types you can create
toView()
returns aLvArray::CRSMatrixView< T, COL_TYPE, INDEX_TYPE const, BUFFER_TYPE >
which can modify the entries as well as insert and remove columns from the rows as long as the size of each row doesn’t exceed its capacity.toViewConstSizes()
returns aLvArray::CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
which can modify existing values but cannot add or remove columns from the rows.toViewConst()
returns aLvArray::CRSMatrixView< T const, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
which provides read only access to both the columns and values of the rows.toSparsityPatternView()
returns aLvArray::SparsityPatternView< COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
which provides read only access to the columns of the rows.
TEST( CRSMatrix, views )
{
// Create a 100x100 tri-diagonal sparsity pattern.
LvArray::SparsityPattern< int, std::ptrdiff_t, LvArray::MallocBuffer > sparsity( 100, 100, 3 );
// Since enough space has been preallocated in each row we can insert into the sparsity
// pattern in parallel.
LvArray::SparsityPatternView< int,
std::ptrdiff_t const,
LvArray::MallocBuffer > const sparsityView = sparsity.toView();
RAJA::forall< RAJA::omp_parallel_for_exec >(
RAJA::TypedRangeSegment< std::ptrdiff_t >( 0, sparsityView.numRows() ),
[sparsityView] ( std::ptrdiff_t const row )
{
int const columns[ 3 ] = { int( row - 1 ), int( row ), int( row + 1 ) };
int const begin = row == 0 ? 1 : 0;
int const end = row == sparsityView.numRows() - 1 ? 2 : 3;
sparsityView.insertNonZeros( row, columns + begin, columns + end );
}
);
// Create a matrix from the sparsity pattern
LvArray::CRSMatrix< double, int, std::ptrdiff_t, LvArray::MallocBuffer > matrix;
matrix.assimilate< RAJA::omp_parallel_for_exec >( std::move( sparsity ) );
// Assemble into the matrix.
LvArray::CRSMatrixView< double,
int const,
std::ptrdiff_t const,
LvArray::MallocBuffer > const matrixView = matrix.toViewConstSizes();
RAJA::forall< RAJA::omp_parallel_for_exec >(
RAJA::TypedRangeSegment< std::ptrdiff_t >( 0, matrixView.numRows() ),
[matrixView] ( std::ptrdiff_t const row )
{
// Some silly assembly where each diagonal entry ( r, r ) has a contribution to the row above ( r-1, r )
// The current row (r, r-1), (r, r), (r, r+1) and the row below (r+1, r).
int const columns[ 3 ] = { int( row - 1 ), int( row ), int( row + 1 ) };
double const contribution[ 3 ] = { double( row ), double( row ), double( row ) };
// Contribution to the row above
if( row > 0 )
{
matrixView.addToRow< RAJA::builtin_atomic >( row - 1, columns + 1, contribution, 1 );
}
// Contribution to the current row
int const begin = row == 0 ? 1 : 0;
int const end = row == matrixView.numRows() - 1 ? 2 : 3;
matrixView.addToRow< RAJA::builtin_atomic >( row, columns + begin, contribution, end - begin );
// Contribution to the row below
if( row < matrixView.numRows() - 1 )
{
matrixView.addToRow< RAJA::builtin_atomic >( row + 1, columns + 1, contribution, 1 );
}
}
);
// Check every row except for the first and the last.
for( std::ptrdiff_t row = 1; row < matrix.numRows() - 1; ++row )
{
EXPECT_EQ( matrix.numNonZeros( row ), 3 );
LvArray::ArraySlice< int const, 1, 0, std::ptrdiff_t > const rowColumns = matrix.getColumns( row );
LvArray::ArraySlice< double const, 1, 0, std::ptrdiff_t > const rowEntries = matrix.getEntries( row );
for( std::ptrdiff_t i = 0; i < matrix.numNonZeros( row ); ++i )
{
// The first first entry is the sum of the contributions of the row above and the current row.
EXPECT_EQ( rowColumns[ 0 ], row - 1 );
EXPECT_EQ( rowEntries[ 0 ], row + row - 1 );
// The second entry is the from the current row alone.
EXPECT_EQ( rowColumns[ 1 ], row );
EXPECT_EQ( rowEntries[ 1 ], row );
// The third entry is the sum of the contributions of the row below and the current row.
EXPECT_EQ( rowColumns[ 2 ], row + 1 );
EXPECT_EQ( rowEntries[ 2 ], row + row + 1 );
}
}
// Check the first row.
EXPECT_EQ( matrix.numNonZeros( 0 ), 2 );
EXPECT_EQ( matrix.getColumns( 0 )[ 0 ], 0 );
EXPECT_EQ( matrix.getEntries( 0 )[ 0 ], 0 );
EXPECT_EQ( matrix.getColumns( 0 )[ 1 ], 1 );
EXPECT_EQ( matrix.getEntries( 0 )[ 1 ], 1 );
// Check the last row.
EXPECT_EQ( matrix.numNonZeros( matrix.numRows() - 1 ), 2 );
EXPECT_EQ( matrix.getColumns( matrix.numRows() - 1 )[ 0 ], matrix.numRows() - 2 );
EXPECT_EQ( matrix.getEntries( matrix.numRows() - 1 )[ 0 ], matrix.numRows() - 1 + matrix.numRows() - 2 );
EXPECT_EQ( matrix.getColumns( matrix.numRows() - 1 )[ 1 ], matrix.numRows() - 1 );
EXPECT_EQ( matrix.getEntries( matrix.numRows() - 1 )[ 1 ], matrix.numRows() - 1 );
}
[Source: examples/exampleSparsityPatternAndCRSMatrix.cpp]
7.4. Usage with LvArray::ChaiBuffer
¶
The three types of LvArray::CRSMatrixView
obtainable from an LvArray::CRSMatrixs
all act differently when moved to a new memory space.
LvArray::CRSMatrixsView< T, COL_TYPE, INDEX_TYPE const, LvArray::ChaiBuffer >
, obtained by callingtoView()
. When it is moved to a new space the values, columns and sizes are all touched. The offsets are not touched.LvArray::CRSMatrixsView< T, COL_TYPE const, INDEX_TYPE const, LvArray::ChaiBuffer >
, obtained by callingtoViewConstSizes()
. When it is moved to a new space the values are touched but the columns, sizes and offsets aren’t.LvArray::CRSMatrixsView< T const, COL_TYPE const, INDEX_TYPE const, LvArray::ChaiBuffer >
, obtained by callingtoViewConst()
. None of the buffers are touched in the new space.
Calling the explicit move
method with the touch parameter set to true
on a view type has the behavior described above. However calling move( MemorySpace::host )
on an LvArray::CRSMatrix
or LvArray::SparsityPattern
will also touch the offsets (if moving to the GPU the offsets aren’t touched). This is the only way to touch the offsets so if an LvArray::CRSMatrix
was previously on the device then it must be explicitly moved and touched on the host before any modification to the offsets can safely take place.
7.5. Usage with LVARRAY_BOUNDS_CHECK
¶
When LVARRAY_BOUNDS_CHECK
is defined access all row and column access is checked. Methods which expect a sorted unique set of columns check that the columns are indeed sorted and unique. In addition if addToRow
checks that all the given columns are present in the row.
7.6. Guidelines¶
As with all the LvArray
containers it is important to pass around the most restrictive form. A function should only accept a LvArray::CRSMatrix
if it needs to resize the matrix or might bust the capacity of a row. If a function only needs to be able to modify existing entries it should accept a LvArray::CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
. If a function only needs to examine the sparsity pattern of the matrix it should accept a LvArray::SparsityPatternView< COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
.
Like the LvArray::ArrayOfArrays
and LvArray::ArrayOfSets
when constructing an LvArray::SparsityPattern
or LvArray::CRSMatrix
it is important to preallocate space for each row in order to achieve decent performance.
A common pattern with sparse matrices is that the sparsity pattern need only be assembled once but the matrix is used multiple times with different values. For example when using the finite element method on an unstructured mesh you can generate the sparsity pattern once at the beginning of the simulation and then each time step you repopulate the entries of the matrix. When this is the case it is usually best to do the sparsity generation with a LvArray::SparsityPattern
and then assimilate that into a LvArray::CRSMatrix
. Once you have a matrix with the proper sparsity pattern create a LvArray::CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
via toViewConstSizes()
and you can then assemble into the matrix in parallel with the addToRow
methods. Finally before beginning the next iteration you can zero out the entries in the matrix by calling setValues
.