MatrixBase.hpp

Go to the documentation of this file.

/*
 * ------------------------------------------------------------------------------------------------------------
 * SPDX-License-Identifier: LGPL-2.1-only
 *
 * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
 * Copyright (c) 2018-2024 TotalEnergies
 * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
 * Copyright (c) 2023-2024 Chevron
 * Copyright (c) 2019-     GEOS/GEOSX Contributors
 * All rights reserved
 *
 * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
 * ------------------------------------------------------------------------------------------------------------
 */
 
#ifndef GEOS_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_
#define GEOS_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_
 
#include "linearAlgebra/common/common.hpp"
#include "linearAlgebra/common/LinearOperator.hpp"
#include "LvArray/src/output.hpp"
 
namespace geos
{
 
class DofManager;
 
enum class RowSumType
{
  SumValues,
  SumAbsValues,
  SumSqrValues,
  MaxAbsValues
};
 
enum class MatrixPatternOp
{
  Equal,    // Caller guarantees patterns of arguments are exactly the same
  Subset,   // Caller guarantees pattern of second argument is a subset of the first
  Restrict, // Restrict pattern of second argument, ignoring any entries that don't exist in first
  Extend    // Extend pattern of the first argument with entries from the second
};
 
template< typename MATRIX, typename VECTOR >
class MatrixBase : public virtual LinearOperator< VECTOR >
{
protected:
 
  using Base = LinearOperator< VECTOR >;
 
  using Matrix = MATRIX;
 
  using Vector = VECTOR;
 
  using Base::numGlobalRows;
  using Base::numGlobalCols;
  using Base::numGlobalNonzeros;
  using Base::numLocalRows;
  using Base::numLocalCols;
  using Base::numLocalNonzeros;
 
 
  inline bool closed() const { return m_closed; }
 
  inline bool assembled() const { return m_assembled; }
 
  inline bool ready() const { return closed() && assembled(); }
 
  inline bool modifiable() const { return !closed() && assembled(); }
 
  inline bool insertable() const { return !closed() && !assembled(); }
 
  virtual bool created() const = 0;
 
 
 
  void setDofManager( DofManager const * const dofManager )
  {
    m_dofManager = dofManager;
  }
 
  DofManager const * dofManager() const
  {
    return m_dofManager;
  }
 
 
 
  virtual void
  createWithLocalSize( localIndex const localSize,
                       localIndex const maxEntriesPerRow,
                       MPI_Comm const & comm )
  {
    createWithLocalSize( localSize, localSize, maxEntriesPerRow, comm );
  }
 
  virtual void
  createWithGlobalSize( globalIndex const globalSize,
                        localIndex const maxEntriesPerRow,
                        MPI_Comm const & comm )
  {
    createWithGlobalSize( globalSize, globalSize, maxEntriesPerRow, comm );
  }
 
  virtual void
  createWithLocalSize( localIndex const localRows,
                       localIndex const localCols,
                       localIndex const maxEntriesPerRow,
                       MPI_Comm const & comm ) = 0;
 
  virtual void
  createWithGlobalSize( globalIndex const globalRows,
                        globalIndex const globalCols,
                        localIndex const maxEntriesPerRow,
                        MPI_Comm const & comm ) = 0;
 
  virtual void create( CRSMatrixView< real64 const, globalIndex const > const & localMatrix,
                       localIndex const numLocalColumns,
                       MPI_Comm const & comm )
  {
    localMatrix.move( hostMemorySpace, false );
 
    localIndex maxEntriesPerRow = 0;
    for( localIndex i = 0; i < localMatrix.numRows(); ++i )
    {
      maxEntriesPerRow = std::max( maxEntriesPerRow, localMatrix.numNonZeros( i ) );
    }
 
    createWithLocalSize( localMatrix.numRows(),
                         numLocalColumns,
                         maxEntriesPerRow,
                         comm );
 
    globalIndex const rankOffset = ilower();
 
    open();
    for( localIndex localRow = 0; localRow < localMatrix.numRows(); ++localRow )
    {
      insert( localRow + rankOffset, localMatrix.getColumns( localRow ), localMatrix.getEntries( localRow ) );
    }
    close();
  }
 
 
 
  virtual void open() = 0;
 
  virtual void close() = 0;
 
  virtual void reset()
  {
    m_assembled = false;
    m_closed = true;
  }
 
 
 
  virtual void set( real64 const value ) = 0;
 
  virtual void zero() = 0;
 
 
 
  virtual void add( globalIndex const rowIndex,
                    globalIndex const colIndex,
                    real64 const value ) = 0;
 
  virtual void set( globalIndex const rowIndex,
                    globalIndex const colIndex,
                    real64 const value ) = 0;
 
  virtual void insert( globalIndex const rowIndex,
                       globalIndex const colIndex,
                       real64 const value ) = 0;
 
  virtual void add( globalIndex const rowIndex,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const size ) = 0;
 
  virtual void set( globalIndex const rowIndex,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const size ) = 0;
 
  virtual void insert( globalIndex const rowIndex,
                       globalIndex const * colIndices,
                       real64 const * values,
                       localIndex const size ) = 0;
 
  virtual void add( globalIndex const rowIndex,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice1d< real64 const > const & values ) = 0;
 
  virtual void set( globalIndex const rowIndex,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice1d< real64 const > const & values ) = 0;
 
  virtual void insert( globalIndex const rowIndex,
                       arraySlice1d< globalIndex const > const & colIndices,
                       arraySlice1d< real64 const > const & values ) = 0;
 
  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice2d< real64 const > const & values ) = 0;
 
  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
                    arraySlice1d< globalIndex const > const & colIndices,
                    arraySlice2d< real64 const > const & values ) = 0;
 
  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
                       arraySlice1d< globalIndex const > const & colIndices,
                       arraySlice2d< real64 const > const & values ) = 0;
 
  virtual void add( globalIndex const * rowIndices,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const numRows,
                    localIndex const numCols ) = 0;
 
  virtual void set( globalIndex const * rowIndices,
                    globalIndex const * colIndices,
                    real64 const * values,
                    localIndex const numRows,
                    localIndex const numCols ) = 0;
 
  virtual void insert( globalIndex const * rowIndices,
                       globalIndex const * colIndices,
                       real64 const * values,
                       localIndex const numRows,
                       localIndex const numCols ) = 0;
 
  virtual void insert( arrayView1d< globalIndex const > const & rowIndices,
                       arrayView1d< globalIndex const > const & colIndices,
                       arrayView1d< real64 const > const & values ) = 0;
 
 
 
  virtual void apply( Vector const & src, Vector & dst ) const override = 0;
 
  virtual void applyTranspose( Vector const & src, Vector & dst ) const = 0;
 
  virtual void residual( Vector const & x, Vector const & b, Vector & r ) const override
  {
    if( &b != &r )
    {
      r.copy( b );
    }
    gemv( -1.0, x, 1.0, r );
  }
 
  virtual void multiply( Matrix const & src,
                         Matrix & dst ) const = 0;
 
  virtual void leftMultiplyTranspose( Matrix const & src,
                                      Matrix & dst ) const = 0;
 
  virtual void rightMultiplyTranspose( Matrix const & src,
                                       Matrix & dst ) const = 0;
 
  virtual void multiplyRAP( Matrix const & R,
                            Matrix const & P,
                            Matrix & dst ) const
  {
    GEOS_LAI_ASSERT( ready() );
    GEOS_LAI_ASSERT( R.ready() );
    GEOS_LAI_ASSERT( P.ready() );
    GEOS_LAI_ASSERT_EQ( numGlobalRows(), R.numGlobalCols() );
    GEOS_LAI_ASSERT_EQ( numGlobalCols(), P.numGlobalRows() );
 
    Matrix AP;
    multiply( P, AP );
    R.multiply( AP, dst );
  }
 
  virtual void multiplyP1tAP2( Matrix const & P1,
                               Matrix const & P2,
                               Matrix & dst ) const
  {
    GEOS_LAI_ASSERT( ready() );
    GEOS_LAI_ASSERT( P1.ready() );
    GEOS_LAI_ASSERT( P2.ready() );
    GEOS_LAI_ASSERT_EQ( numLocalRows(), P1.numLocalRows() );
    GEOS_LAI_ASSERT_EQ( numLocalCols(), P2.numLocalRows() );
 
    Matrix AP2;
    multiply( P2, AP2 );
    P1.leftMultiplyTranspose( AP2, dst );
  }
 
  virtual void multiplyPtAP( Matrix const & P,
                             Matrix & dst ) const
  {
    multiplyP1tAP2( P, P, dst );
  }
 
  virtual void multiplyRARt( Matrix const & R,
                             Matrix & dst ) const
  {
    Matrix P;
    R.transpose( P );
    multiplyPtAP( P, dst );
  }
 
  virtual void gemv( real64 const alpha,
                     Vector const & x,
                     real64 const beta,
                     Vector & y,
                     bool useTranspose = false ) const = 0;
 
  virtual void separateComponentFilter( Matrix & dst,
                                        integer const dofPerPoint ) const = 0;
 
 
 
  virtual void scale( real64 const scalingFactor ) = 0;
 
  virtual void leftScale( Vector const & vec ) = 0;
 
  virtual void rightScale( Vector const & vec ) = 0;
 
  virtual void leftRightScale( Vector const & vecLeft,
                               Vector const & vecRight ) = 0;
 
  virtual void computeScalingVector( Vector & scaling ) const = 0;
 
  virtual void rescaleRows( arrayView1d< globalIndex const > const & rowIndices,
                            RowSumType const rowSumType ) = 0;
 
  virtual void transpose( Matrix & dst ) const = 0;
 
  virtual real64 clearRow( globalIndex const row,
                           bool const keepDiag = false,
                           real64 const diagValue = 0.0 ) = 0;
 
  virtual void addEntries( Matrix const & src,
                           MatrixPatternOp const op,
                           real64 const scale ) = 0;
 
  virtual void addDiagonal( Vector const & src,
                            real64 const scale ) = 0;
 
  virtual void clampEntries( real64 const lo,
                             real64 const hi,
                             bool const excludeDiag ) = 0;
 
 
 
  virtual localIndex maxRowLengthLocal() const = 0;
 
  virtual localIndex maxRowLength() const
  {
    return MpiWrapper::max( maxRowLengthLocal(), this->comm() );
  }
 
  virtual localIndex rowLength( globalIndex const globalRowIndex ) const = 0;
 
  virtual void getRowLengths( arrayView1d< localIndex > const & lengths ) const = 0;
 
  virtual void getRowLocalLengths( arrayView1d< localIndex > const & lengths ) const = 0;
 
  virtual void getRowCopy( globalIndex const globalRow,
                           arraySlice1d< globalIndex > const & colIndices,
                           arraySlice1d< real64 > const & values ) const = 0;
 
  virtual void extractDiagonal( Vector & dst ) const = 0;
 
  virtual void extract( CRSMatrixView< real64, globalIndex > const & localMat ) const = 0;
 
  virtual void extract( CRSMatrixView< real64, globalIndex const > const & localMat ) const = 0;
 
  CRSMatrix< real64, globalIndex > extract() const
  {
    array1d< localIndex > rowLengths( numLocalRows() );
    getRowLengths( rowLengths.toView() );
    rowLengths.move( LvArray::MemorySpace::host );
    CRSMatrix< real64, globalIndex > localMat;
    localMat.resizeFromRowCapacities< parallelHostPolicy >( numLocalRows(), numGlobalCols(), rowLengths.data() );
    extract( localMat.toView() );
    return localMat;
  }
 
  virtual void extractLocal( CRSMatrixView< real64, localIndex > const & localMat ) const = 0;
 
  CRSMatrix< real64, localIndex > extractLocal() const
  {
    array1d< localIndex > rowLengths( numLocalRows() );
    getRowLocalLengths( rowLengths.toView() );
    CRSMatrix< real64, localIndex > localMat;
    localMat.resizeFromRowCapacities< parallelHostPolicy >( numLocalRows(), numGlobalCols(), rowLengths.data() );
    extractLocal( localMat.toView() );
    return localMat;
  }
 
  virtual void getRowSums( Vector & dst, RowSumType const rowSumType ) const = 0;
 
  virtual globalIndex ilower() const = 0;
 
  virtual globalIndex iupper() const = 0;
 
  virtual globalIndex jlower() const = 0;
 
  virtual globalIndex jupper() const = 0;
 
  virtual real64 normInf() const = 0;
 
  virtual real64 norm1() const = 0;
 
  virtual real64 normFrobenius() const = 0;
 
  virtual real64 normMax() const = 0;
 
  virtual real64 normMax( arrayView1d< globalIndex const > const & rowIndices ) const = 0;
 
  virtual localIndex getLocalRowID( globalIndex const index ) const = 0;
 
  virtual globalIndex getGlobalRowID( localIndex const index ) const = 0;
 
 
 
  virtual void print( std::ostream & os = std::cout ) const = 0;
 
  virtual void write( string const & filename,
                      LAIOutputFormat const format = LAIOutputFormat::MATRIX_MARKET ) const = 0;
 
 
  friend std::ostream & operator<<( std::ostream & os, Matrix const & matrix )
  {
    matrix.print( os );
    return os;
  }
 
  bool m_closed = true;
 
  bool m_assembled = false;
 
  DofManager const * m_dofManager{};
 
};
 
} // namespace geos
 
#endif //GEOS_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_