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:

  1. The columns of each row are sorted.
  2. 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

  1. COL_TYPE: The integral type used to enumerate the columns of the matrix.
  2. INDEX_TYPE: An integral type used in index calculations, the suggested type is std::ptrdiff_t.
  3. BUFFER_TYPE: A template template parameter specifying the buffer type used for allocation and de-allocation, the LvArray::SparsityPattern contains a BUFFER_TYPE< COL_TYPE > along with two BUFFER_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 a LvArray::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 a LvArray::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 a LvArray::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 a LvArray::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 calling toView(). 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 calling toViewConstSizes(). 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 calling toViewConst(). 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.