CRSMatrix.hpp

Go to the documentation of this file.

/*
 * Copyright (c) 2021, Lawrence Livermore National Security, LLC and LvArray contributors.
 * All rights reserved.
 * See the LICENSE file for details.
 * SPDX-License-Identifier: (BSD-3-Clause)
 */

#pragma once

#include "CRSMatrixView.hpp"
#include "arrayManipulation.hpp"

namespace LvArray
{

template< typename COL_TYPE, typename INDEX_TYPE, template< typename > class BUFFER_TYPE >
class SparsityPattern;

template< typename T,
          typename COL_TYPE,
          typename INDEX_TYPE,
          template< typename > class BUFFER_TYPE >
class CRSMatrix : protected CRSMatrixView< T, COL_TYPE, INDEX_TYPE, BUFFER_TYPE >
{

  using ParentClass = CRSMatrixView< T, COL_TYPE, INDEX_TYPE, BUFFER_TYPE >;

public:

  using typename ParentClass::EntryType;
  using typename ParentClass::ColType;
  using typename ParentClass::IndexType;


  CRSMatrix( INDEX_TYPE const nrows=0,
             INDEX_TYPE const ncols=0,
             INDEX_TYPE const initialRowCapacity=0 ):
    ParentClass( true )
  {
    resize( nrows, ncols, initialRowCapacity );
    setName( "" );
  }

  inline
  CRSMatrix( CRSMatrix const & src ):
    ParentClass( true )
  { *this = src; }

  inline
  CRSMatrix( CRSMatrix && ) = default;

  ~CRSMatrix()
  { ParentClass::free( this->m_entries ); }



  inline
  CRSMatrix & operator=( CRSMatrix const & src )
  {
    this->m_numCols = src.m_numCols;
    ParentClass::setEqualTo( src.m_numArrays,
                             src.m_offsets[ src.m_numArrays ],
                             src.m_offsets,
                             src.m_sizes,
                             src.m_values,
                             typename ParentClass::template PairOfBuffers< T >( this->m_entries, src.m_entries ) );
    return *this;
  }

  inline
  CRSMatrix & operator=( CRSMatrix && src )
  {
    ParentClass::free( this->m_entries );
    ParentClass::operator=( std::move( src ) );
    return *this;
  }

  template< typename POLICY >
  inline
  void assimilate( SparsityPattern< COL_TYPE, INDEX_TYPE, BUFFER_TYPE > && src )
  {
    // Destroy the current entries.
    if( !std::is_trivially_destructible< T >::value )
    {
      CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE > const view = toViewConstSizes();
      RAJA::forall< POLICY >( RAJA::TypedRangeSegment< INDEX_TYPE >( 0, numRows() ),
                              [view] LVARRAY_HOST_DEVICE ( INDEX_TYPE const row )
        {
          INDEX_TYPE const nnz = view.numNonZeros( row );
          T * const entries = view.getEntries( row );
          arrayManipulation::destroy( entries, nnz );
        } );
    }

    // Reallocate to the appropriate length
    bufferManipulation::reserve( this->m_entries, 0, MemorySpace::host, src.nonZeroCapacity() );

    ParentClass::assimilate( reinterpret_cast< SparsityPatternView< COL_TYPE, INDEX_TYPE, BUFFER_TYPE > && >( src ) );

    CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE > const view = toViewConstSizes();
    RAJA::forall< POLICY >( RAJA::TypedRangeSegment< INDEX_TYPE >( 0, numRows() ),
                            [view] LVARRAY_HOST_DEVICE ( INDEX_TYPE const row )
      {
        T * const rowEntries = view.getEntries( row );
        for( INDEX_TYPE i = 0; i < view.numNonZeros( row ); ++i )
        {
          new ( rowEntries + i ) T();
        }
      } );

    setName( "" );
  }

  template< typename POLICY >
  void resizeFromRowCapacities( INDEX_TYPE const nRows, INDEX_TYPE const nCols, INDEX_TYPE const * const rowCapacities )
  {
    LVARRAY_ERROR_IF( !arrayManipulation::isPositive( nCols ), "nCols must be positive." );
    LVARRAY_ERROR_IF( nCols - 1 > std::numeric_limits< COL_TYPE >::max(),
                      "COL_TYPE must be able to hold the range of columns: [0, " << nCols - 1 << "]." );

    this->m_numCols = nCols;
    ParentClass::template resizeFromCapacities< POLICY >( nRows, rowCapacities, this->m_entries );
  }



  constexpr inline
  CRSMatrixView< T, COL_TYPE, INDEX_TYPE const, BUFFER_TYPE >
  toView() const &
  { return ParentClass::toView(); }

  constexpr inline
  CRSMatrixView< T, COL_TYPE, INDEX_TYPE const, BUFFER_TYPE >
  toView() const && = delete;

  LVARRAY_HOST_DEVICE constexpr inline
  CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
  toViewConstSizes() const &
  { return ParentClass::toViewConstSizes(); }

  LVARRAY_HOST_DEVICE constexpr inline
  CRSMatrixView< T, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
  toViewConstSizes() const && = delete;

  LVARRAY_HOST_DEVICE constexpr inline
  CRSMatrixView< T const, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
  toViewConst() const &
  { return ParentClass::toViewConst(); }

  LVARRAY_HOST_DEVICE constexpr inline
  CRSMatrixView< T const, COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
  toViewConst() const && = delete;

  using ParentClass::toSparsityPatternView;

  LVARRAY_HOST_DEVICE constexpr inline
  SparsityPatternView< COL_TYPE const, INDEX_TYPE const, BUFFER_TYPE >
  toSparsityPatternView() const && = delete;



  using ParentClass::numRows;
  using ParentClass::numColumns;
  using ParentClass::numNonZeros;
  using ParentClass::nonZeroCapacity;
  using ParentClass::empty;



  using ParentClass::getOffsets;
  using ParentClass::getColumns;
  using ParentClass::getEntries;



  inline
  void reserveNonZeros( INDEX_TYPE const nnz )
  { ParentClass::reserveValues( nnz, this->m_entries ); }

  inline
  void reserveNonZeros( INDEX_TYPE const row, INDEX_TYPE const nnz )
  {
    if( nonZeroCapacity( row ) >= nnz ) return;
    setRowCapacity( row, nnz );
  }

  inline
  void setRowCapacity( INDEX_TYPE const row, INDEX_TYPE newCapacity )
  {
    if( newCapacity > numColumns() ) newCapacity = numColumns();
    ParentClass::setCapacityOfArray( row, newCapacity, this->m_entries );
  }

  void resize( INDEX_TYPE const nRows, INDEX_TYPE const nCols, INDEX_TYPE const initialRowCapacity )
  { ParentClass::resize( nRows, nCols, initialRowCapacity, this->m_entries ); }

  inline
  void compress()
  { ParentClass::compress( this->m_entries ); }



  inline
  bool insertNonZero( INDEX_TYPE const row, COL_TYPE const col, T const & entry )
  { return ParentClass::insertIntoSetImpl( row, col, CallBacks( *this, row, &entry ) ); }

  inline
  INDEX_TYPE insertNonZeros( INDEX_TYPE const row,
                             COL_TYPE const * const cols,
                             T const * const entriesToInsert,
                             INDEX_TYPE const ncols )
  { return ParentClass::insertIntoSetImpl( row, cols, cols + ncols, CallBacks( *this, row, entriesToInsert ) ); }

  using ParentClass::removeNonZero;
  using ParentClass::removeNonZeros;



  template< typename POLICY >
  inline
  void setValues( T const & value ) const
  { ParentClass::template setValues< POLICY >( value ); }

  using ParentClass::zero;

  template< typename AtomicPolicy >
  inline
  void addToRow( INDEX_TYPE const row,
                 COL_TYPE const * const LVARRAY_RESTRICT cols,
                 T const * const LVARRAY_RESTRICT vals,
                 INDEX_TYPE const nCols ) const
  { ParentClass::template addToRow< AtomicPolicy >( row, cols, vals, nCols ); }

  template< typename AtomicPolicy >
  inline
  void addToRowBinarySearch( INDEX_TYPE const row,
                             COL_TYPE const * const LVARRAY_RESTRICT cols,
                             T const * const LVARRAY_RESTRICT vals,
                             INDEX_TYPE const nCols ) const
  { ParentClass::template addToRowBinarySearch< AtomicPolicy >( row, cols, vals, nCols ); }

  template< typename AtomicPolicy >
  inline
  void addToRowLinearSearch( INDEX_TYPE const row,
                             COL_TYPE const * const LVARRAY_RESTRICT cols,
                             T const * const LVARRAY_RESTRICT vals,
                             INDEX_TYPE const nCols ) const
  { ParentClass::template addToRowLinearSearch< AtomicPolicy >( row, cols, vals, nCols ); }

  template< typename AtomicPolicy >
  inline
  void addToRowBinarySearchUnsorted( INDEX_TYPE const row,
                                     COL_TYPE const * const LVARRAY_RESTRICT cols,
                                     T const * const LVARRAY_RESTRICT vals,
                                     INDEX_TYPE const nCols ) const
  { ParentClass::template addToRowBinarySearchUnsorted< AtomicPolicy >( row, cols, vals, nCols ); }



  void move( MemorySpace const space, bool const touch=true ) const
  { return ParentClass::move( space, touch ); }


  void setName( std::string const & name )
  { ParentClass::template setName< decltype( *this ) >( name ); }

private:

  void dynamicallyGrowRow( INDEX_TYPE const row, INDEX_TYPE const newNNZ )
  { setRowCapacity( row, 2 * newNNZ ); }

  class CallBacks
  {
public:

    CallBacks( CRSMatrix & crsM,
               INDEX_TYPE const row, T const * const entriesToInsert ):
      m_crsM( crsM ),
      m_row( row ),
      m_rowNNZ( crsM.numNonZeros( row ) ),
      m_rowCapacity( crsM.nonZeroCapacity( row ) ),
      m_entries( nullptr ),
      m_entriesToInsert( entriesToInsert )
    {}

    inline
    COL_TYPE * incrementSize( COL_TYPE * const curPtr, INDEX_TYPE const nToAdd )
    {
      LVARRAY_UNUSED_VARIABLE( curPtr );
      if( m_rowNNZ + nToAdd > m_rowCapacity )
      {
        m_crsM.dynamicallyGrowRow( m_row, m_rowNNZ + nToAdd );
      }

      m_entries = m_crsM.getEntries( m_row );
      return m_crsM.getSetValues( m_row );
    }

    inline
    void insert( INDEX_TYPE const pos ) const
    { arrayManipulation::emplace( m_entries, m_rowNNZ, pos, m_entriesToInsert[0] ); }

    inline
    void set( INDEX_TYPE const pos, INDEX_TYPE const colPos ) const
    { new (&m_entries[pos]) T( m_entriesToInsert[colPos] ); }

    inline
    void insert( INDEX_TYPE const nLeftToInsert,
                 INDEX_TYPE const colPos,
                 INDEX_TYPE const pos,
                 INDEX_TYPE const prevPos ) const
    {
      arrayManipulation::shiftUp( m_entries, prevPos, pos, nLeftToInsert );
      new (&m_entries[pos + nLeftToInsert - 1]) T( m_entriesToInsert[colPos] );
    }

private:
    CRSMatrix & m_crsM;

    INDEX_TYPE const m_row;

    INDEX_TYPE const m_rowNNZ;

    INDEX_TYPE const m_rowCapacity;

    T * m_entries;

    T const * const m_entriesToInsert;
  };
};

} /* namespace LvArray */