GEOSX
MatrixBase.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University
7  * Copyright (c) 2018-2020 Total, S.A
8  * Copyright (c) 2019- GEOSX Contributors
9  * All rights reserved
10  *
11  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
12  * ------------------------------------------------------------------------------------------------------------
13  */
14 
19 #ifndef GEOSX_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_
20 #define GEOSX_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_
21 
22 #include "linearAlgebra/common.hpp"
24 //#include "LvArray/src/streamIO.hpp"
25 
26 namespace geosx
27 {
28 
48 template< typename MATRIX, typename VECTOR >
49 class MatrixBase : public virtual LinearOperator< VECTOR >
50 {
51 public:
52 
54  using Matrix = MATRIX;
55 
57  using Vector = VECTOR;
58 
59 protected:
60 
64 
70  : m_closed( true ),
71  m_assembled( false )
72  {}
73 
77  MatrixBase( MatrixBase const & ) = default;
78 
82  MatrixBase( MatrixBase && ) = default;
83 
88  MatrixBase & operator=( MatrixBase const & ) = default;
89 
94  MatrixBase & operator=( MatrixBase && ) = default;
95 
99  ~MatrixBase() = default;
100 
102 
106 
112  inline bool closed() const { return m_closed; }
113 
118  inline bool assembled() const { return m_assembled; }
119 
125  inline bool ready() const { return closed() && assembled(); }
126 
133  inline bool modifiable() const { return !closed() && assembled(); }
134 
140  inline bool insertable() const { return !closed() && !assembled(); }
141 
146  virtual bool created() const = 0;
147 
149 
153 
163  virtual void
164  createWithLocalSize( localIndex const localSize,
165  localIndex const maxEntriesPerRow,
166  MPI_Comm const & comm )
167  {
168  createWithLocalSize( localSize, localSize, maxEntriesPerRow, comm );
169  }
170 
181  virtual void
183  localIndex const maxEntriesPerRow,
184  MPI_Comm const & comm )
185  {
186  createWithGlobalSize( globalSize, globalSize, maxEntriesPerRow, comm );
187  }
188 
197  virtual void
198  createWithLocalSize( localIndex const localRows,
199  localIndex const localCols,
200  localIndex const maxEntriesPerRow,
201  MPI_Comm const & comm ) = 0;
202 
211  virtual void
212  createWithGlobalSize( globalIndex const globalRows,
213  globalIndex const globalCols,
214  localIndex const maxEntriesPerRow,
215  MPI_Comm const & comm ) = 0;
216 
225  virtual void create( CRSMatrixView< real64 const, globalIndex const > const & localMatrix,
226  MPI_Comm const & comm )
227  {
228  localMatrix.move( LvArray::MemorySpace::CPU, false );
229 
230  localIndex maxEntriesPerRow = 0;
231  for( localIndex i = 0; i < localMatrix.numRows(); ++i )
232  {
233  maxEntriesPerRow = std::max( maxEntriesPerRow, localMatrix.numNonZeros( i ) );
234  }
235 
236  createWithLocalSize( localMatrix.numRows(),
237  localMatrix.numRows(),
238  maxEntriesPerRow,
239  comm );
240 
241  globalIndex const rankOffset = ilower();
242 
243  open();
244  for( localIndex localRow = 0; localRow < localMatrix.numRows(); ++localRow )
245  {
246  insert( localRow + rankOffset, localMatrix.getColumns( localRow ), localMatrix.getEntries( localRow ) );
247  }
248  close();
249  }
250 
252 
256 
264  virtual void open() = 0;
265 
272  virtual void close() = 0;
273 
277  virtual void reset()
278  {
279  m_assembled = false;
280  m_closed = true;
281  }
282 
284 
288 
294  virtual void set( real64 const value ) = 0;
295 
299  virtual void zero() = 0;
300 
302 
313 
321  virtual void add( globalIndex const rowIndex,
322  globalIndex const colIndex,
323  real64 const value ) = 0;
324 
331  virtual void set( globalIndex const rowIndex,
332  globalIndex const colIndex,
333  real64 const value ) = 0;
334 
341  virtual void insert( globalIndex const rowIndex,
342  globalIndex const colIndex,
343  real64 const value ) = 0;
344 
352  virtual void add( globalIndex const rowIndex,
353  globalIndex const * colIndices,
354  real64 const * values,
355  localIndex const size ) = 0;
356 
364  virtual void set( globalIndex const rowIndex,
365  globalIndex const * colIndices,
366  real64 const * values,
367  localIndex const size ) = 0;
368 
376  virtual void insert( globalIndex const rowIndex,
377  globalIndex const * colIndices,
378  real64 const * values,
379  localIndex const size ) = 0;
380 
387  virtual void add( globalIndex const rowIndex,
388  arraySlice1d< globalIndex const > const & colIndices,
389  arraySlice1d< real64 const > const & values ) = 0;
390 
397  virtual void set( globalIndex const rowIndex,
398  arraySlice1d< globalIndex const > const & colIndices,
399  arraySlice1d< real64 const > const & values ) = 0;
400 
407  virtual void insert( globalIndex const rowIndex,
408  arraySlice1d< globalIndex const > const & colIndices,
409  arraySlice1d< real64 const > const & values ) = 0;
410 
417  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
418  arraySlice1d< globalIndex const > const & colIndices,
420 
427  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
428  arraySlice1d< globalIndex const > const & colIndices,
430 
437  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
438  arraySlice1d< globalIndex const > const & colIndices,
440 
447  virtual void add( arraySlice1d< globalIndex const > const & rowIndices,
448  arraySlice1d< globalIndex const > const & colIndices,
450 
457  virtual void set( arraySlice1d< globalIndex const > const & rowIndices,
458  arraySlice1d< globalIndex const > const & colIndices,
460 
467  virtual void insert( arraySlice1d< globalIndex const > const & rowIndices,
468  arraySlice1d< globalIndex const > const & colIndices,
470 
481  virtual void add( globalIndex const * rowIndices,
482  globalIndex const * colIndices,
483  real64 const * values,
484  localIndex const numRows,
485  localIndex const numCols ) = 0;
486 
497  virtual void set( globalIndex const * rowIndices,
498  globalIndex const * colIndices,
499  real64 const * values,
500  localIndex const numRows,
501  localIndex const numCols ) = 0;
502 
513  virtual void insert( globalIndex const * rowIndices,
514  globalIndex const * colIndices,
515  real64 const * values,
516  localIndex const numRows,
517  localIndex const numCols ) = 0;
518 
520 
524 
526  virtual void apply( Vector const & src, Vector & dst ) const override = 0;
527 
533  virtual void applyTranspose( Vector const & src, Vector & dst ) const = 0;
534 
546  virtual void residual( Vector const & x, Vector const & b, Vector & r ) const override
547  {
548  if( &b != &r )
549  {
550  r.copy( b );
551  }
552  gemv( -1.0, x, 1.0, r );
553  }
554 
568  virtual void multiply( Matrix const & src,
569  Matrix & dst ) const = 0;
570 
584  virtual void leftMultiplyTranspose( Matrix const & src,
585  Matrix & dst ) const = 0;
586 
600  virtual void rightMultiplyTranspose( Matrix const & src,
601  Matrix & dst ) const = 0;
602 
609  virtual void multiplyRAP( Matrix const & R,
610  Matrix const & P,
611  Matrix & dst ) const
612  {
613  GEOSX_LAI_ASSERT( ready() );
614  GEOSX_LAI_ASSERT( R.ready() );
615  GEOSX_LAI_ASSERT( P.ready() );
618 
619  Matrix AP;
620  multiply( P, AP );
621  R.multiply( AP, dst );
622  }
623 
629  virtual void multiplyPtAP( Matrix const & P,
630  Matrix & dst ) const
631  {
632  GEOSX_LAI_ASSERT( ready() );
633  GEOSX_LAI_ASSERT( P.ready() );
636 
637  Matrix AP;
638  multiply( P, AP );
639  P.leftMultiplyTranspose( AP, dst );
640  }
641 
655  virtual void gemv( real64 const alpha,
656  Vector const & x,
657  real64 const beta,
658  Vector & y,
659  bool useTranspose = false ) const = 0;
660 
665  virtual void scale( real64 const scalingFactor ) = 0;
666 
671  virtual void leftScale( Vector const & vec ) = 0;
672 
677  virtual void rightScale( Vector const & vec ) = 0;
678 
685  virtual void leftRightScale( Vector const & vecLeft,
686  Vector const & vecRight ) = 0;
687 
696  virtual void transpose( Matrix & dst ) const = 0;
697 
707  virtual real64 clearRow( globalIndex const row,
708  bool const keepDiag = false,
709  real64 const diagValue = 0.0 ) = 0;
710 
720  virtual void addEntries( Matrix const & src,
721  real64 const scale = 1.0,
722  bool const samePattern = true ) = 0;
723 
731  virtual void addDiagonal( Vector const & src ) = 0;
732 
734 
738 
746  virtual localIndex maxRowLength() const = 0;
747 
756  virtual localIndex localRowLength( localIndex localRowIndex ) const = 0;
757 
763  virtual localIndex globalRowLength( globalIndex const globalRowIndex ) const = 0;
764 
771  virtual void getRowCopy( globalIndex const globalRow,
772  arraySlice1d< globalIndex > const & colIndices,
773  arraySlice1d< real64 > const & values ) const = 0;
774 
780  virtual real64 getDiagValue( globalIndex globalRow ) const = 0;
781 
786  virtual void extractDiagonal( Vector & dst ) const = 0;
787 
792  virtual globalIndex numGlobalRows() const override = 0;
793 
798  virtual globalIndex numGlobalCols() const override = 0;
799 
804  virtual localIndex numLocalRows() const = 0;
805 
815  virtual localIndex numLocalCols() const = 0;
816 
821  virtual globalIndex ilower() const = 0;
822 
829  virtual globalIndex iupper() const = 0;
830 
840  virtual globalIndex jlower() const = 0;
841 
849  virtual globalIndex jupper() const = 0;
850 
855  virtual localIndex numLocalNonzeros() const = 0;
856 
861  virtual globalIndex numGlobalNonzeros() const = 0;
862 
867  virtual real64 normInf() const = 0;
868 
873  virtual real64 norm1() const = 0;
874 
879  virtual real64 normFrobenius() const = 0;
880 
886  virtual localIndex getLocalRowID( globalIndex const index ) const = 0;
887 
893  virtual globalIndex getGlobalRowID( localIndex const index ) const = 0;
894 
902  virtual MPI_Comm getComm() const = 0;
903 
905 
909 
915  virtual void print( std::ostream & os = std::cout ) const = 0;
916 
926  virtual void write( string const & filename,
927  LAIOutputFormat const format = LAIOutputFormat::MATRIX_MARKET ) const = 0;
928 
930 
937  friend std::ostream & operator<<( std::ostream & os, Matrix const & matrix )
938  {
939  matrix.print( os );
940  return os;
941  }
942 
944  bool m_closed;
945 
948 
949 };
950 
951 } // namespace geosx
952 
953 #endif //GEOSX_LINEARALGEBRA_INTERFACES_MATRIXBASE_HPP_
virtual void multiplyPtAP(Matrix const &P, Matrix &dst) const
Compute the triple product dst = P^T * this * P
Definition: MatrixBase.hpp:629
Wrapper class for Epetra&#39;s CrsMatrix.
MatrixBase & operator=(MatrixBase const &)=default
Copy assignment.
virtual void createWithLocalSize(localIndex const localSize, localIndex const maxEntriesPerRow, MPI_Comm const &comm)
Create a square matrix from local number of rows.
Definition: MatrixBase.hpp:164
bool assembled() const
Query matrix assembled status.
Definition: MatrixBase.hpp:118
virtual real64 normInf() const =0
Returns the infinity norm of the matrix.
virtual void transpose(Matrix &dst) const =0
Matrix transposition.
friend std::ostream & operator<<(std::ostream &os, Matrix const &matrix)
Stream insertion operator for all matrix types.
Definition: MatrixBase.hpp:937
virtual void leftMultiplyTranspose(Matrix const &src, Matrix &dst) const =0
Matrix/Matrix transpose multiplication.
bool m_assembled
Flag indicating whether the matrix (sparsity pattern) has been assembled.
Definition: MatrixBase.hpp:947
virtual globalIndex jupper() const =0
Returns index one past the last global col owned by that processor.
long long int globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:128
virtual real64 getDiagValue(globalIndex globalRow) const =0
get diagonal element value on a given row
This class serves to provide a sliced multidimensional interface to the family of LvArray classes...
Definition: ArraySlice.hpp:89
std::string format(int NDIM, INDEX_TYPE const *const dims)
This function returns a string that may be used as the "type" in a call to TV_ttf_add_row(). This will either be a single value or an array.
Definition: tv_helpers.hpp:37
virtual void getRowCopy(globalIndex const globalRow, arraySlice1d< globalIndex > const &colIndices, arraySlice1d< real64 > const &values) const =0
Returns a copy of the data in row globalRow.
void move(MemorySpace const space, bool const touch=true) const
Move this SparsityPattern to the given memory space and touch the values, sizes and offsets...
virtual localIndex globalRowLength(globalIndex const globalRowIndex) const =0
Get row length via global row index.
bool ready() const
Query matrix ready status.
Definition: MatrixBase.hpp:125
virtual bool created() const =0
Query matrix creation status.
virtual void close()=0
Assemble and compress the matrix.
virtual void addDiagonal(Vector const &src)=0
Add entries of a vector to the diagonal of this matrix.
virtual void multiply(Matrix const &src, Matrix &dst) const =0
Matrix/Matrix multiplication.
virtual globalIndex numGlobalNonzeros() const =0
Returns the total number of nonzeros in the matrix.
virtual void leftMultiplyTranspose(EpetraMatrix const &src, EpetraMatrix &dst) const override
Matrix/Matrix transpose multiplication.
constexpr INDEX_TYPE_NC numRows() const
virtual localIndex localRowLength(localIndex localRowIndex) const =0
Get row length via local row index.
virtual void applyTranspose(Vector const &src, Vector &dst) const =0
Apply transpose of the matrix to a vector.
double real64
64-bit floating point type.
Definition: DataTypes.hpp:136
virtual globalIndex iupper() const =0
Returns index one past the last global row owned by that processor.
virtual void leftScale(Vector const &vec)=0
Pre-multiplies (left) with diagonal matrix consisting of the values in vec.
virtual void print(std::ostream &os=std::cout) const =0
Print the matrix in Trilinos format to a stream.
MatrixBase()
Constructs a matrix in default state.
Definition: MatrixBase.hpp:69
constexpr ArraySlice< COL_TYPE const, 1, 0, INDEX_TYPE_NC > getColumns(INDEX_TYPE const row) const
virtual void add(globalIndex const rowIndex, globalIndex const colIndex, real64 const value)=0
Add to one element.
ArraySlice< T, 1, 0, INDEX_TYPE_NC > getEntries(INDEX_TYPE const row) const
virtual localIndex numLocalCols() const =0
Return the local number of columns on each processor.
virtual void open()=0
Open matrix for adding new entries.
virtual void apply(Vector const &src, Vector &dst) const override=0
Apply operator to a vector.
virtual void leftRightScale(Vector const &vecLeft, Vector const &vecRight)=0
Post-multiplies (right) with diagonal matrix consisting of the values in vecRight and pre-multiplies ...
virtual globalIndex numGlobalCols() const override=0
Returns the number of global columns.
virtual void insert(globalIndex const rowIndex, globalIndex const colIndex, real64 const value)=0
Insert one element.
virtual real64 norm1() const =0
Returns the one norm of the matrix.
virtual void create(CRSMatrixView< real64 const, globalIndex const > const &localMatrix, MPI_Comm const &comm)
Create parallel matrix from a local CRS matrix.
Definition: MatrixBase.hpp:225
virtual real64 normFrobenius() const =0
Returns the Frobenius norm of the matrix.
virtual globalIndex ilower() const =0
Returns the index of the first global row owned by that processor.
virtual void extractDiagonal(Vector &dst) const =0
Extract diagonal values into a vector.
bool m_closed
Flag indicating whether the matrix is currently open for adding new entries.
Definition: MatrixBase.hpp:944
virtual void multiply(EpetraMatrix const &src, EpetraMatrix &dst) const override
Matrix/Matrix multiplication.
virtual globalIndex numGlobalRows() const override
Returns the number of global rows.
virtual localIndex numLocalRows() const =0
Return the local number of columns on each processor.
virtual localIndex numLocalNonzeros() const =0
Returns the number of nonzeros in the local portion of the matrix.
virtual void createWithGlobalSize(globalIndex const globalSize, localIndex const maxEntriesPerRow, MPI_Comm const &comm)
Create a square matrix from global number of rows.
Definition: MatrixBase.hpp:182
virtual void addEntries(Matrix const &src, real64 const scale=1.0, bool const samePattern=true)=0
Add entries of another matrix to this.
~MatrixBase()=default
Destructor.
LAIOutputFormat
Definition: common.hpp:139
#define GEOSX_LAI_ASSERT_EQ(lhs, rhs)
Definition: common.hpp:47
virtual void multiplyRAP(Matrix const &R, Matrix const &P, Matrix &dst) const
Compute the triple product dst = R * this * P
Definition: MatrixBase.hpp:609
virtual globalIndex numGlobalCols() const override
Returns the number of global columns.
virtual void gemv(real64 const alpha, Vector const &x, real64 const beta, Vector &y, bool useTranspose=false) const =0
Compute gemv y = alpha*A*x + beta*y.
bool closed() const
Query matrix closed status.
Definition: MatrixBase.hpp:112
virtual void print(std::ostream &os=std::cout) const override
Print the matrix in Trilinos format to a stream.
virtual void scale(real64 const scalingFactor)=0
Multiply all elements by scalingFactor.
bool modifiable() const
Query matrix status.
Definition: MatrixBase.hpp:133
bool insertable() const
Query matrix status.
Definition: MatrixBase.hpp:140
std::ptrdiff_t localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:125
virtual void write(string const &filename, LAIOutputFormat const format=LAIOutputFormat::MATRIX_MARKET) const =0
Write the matrix to filename in a matlab-compatible format.
virtual void residual(Vector const &x, Vector const &b, Vector &r) const override
Compute residual r = Ax - b.
Definition: MatrixBase.hpp:546
virtual void zero()=0
Set all elements to zero.
#define GEOSX_LAI_ASSERT(expr)
Definition: common.hpp:33
virtual globalIndex jlower() const =0
Returns the index of the first global col owned by that processor.
virtual localIndex maxRowLength() const =0
Returns the number of nonzero entries in the longest row of the matrix.
virtual localIndex getLocalRowID(globalIndex const index) const =0
Map a global row index to local row index.
virtual void rightMultiplyTranspose(Matrix const &src, Matrix &dst) const =0
Matrix/Matrix transpose multiplication.
virtual real64 clearRow(globalIndex const row, bool const keepDiag=false, real64 const diagValue=0.0)=0
Clear a row, and optionally set diagonal element to diagValue.
virtual void reset()
Reset the matrix to default state.
Definition: MatrixBase.hpp:277
This class creates and provides basic support for the Epetra_FEVector vector object type used in Tril...
Abstract base class for linear operators.
constexpr std::enable_if_t< std::is_arithmetic< T >::value, T > max(T const a, T const b)
Definition: math.hpp:46
virtual void copy(EpetraVector const &x) override
Update vector y as y = x.
virtual globalIndex numGlobalRows() const override=0
Returns the number of global rows.
virtual MPI_Comm getComm() const =0
Get the MPI communicator the matrix was created with.
virtual void rightScale(Vector const &vec)=0
Post-multiplies (right) with diagonal matrix consisting of the values in vec.
Common base template for all matrix wrapper types.
Definition: MatrixBase.hpp:49
This class provides a view into a compressed row storage matrix.
virtual globalIndex getGlobalRowID(localIndex const index) const =0
Map a local row index to global row index.