Linear algebra changes

Author: klevzoff

src/coreComponents/codingUtilities/UnitTestUtilities.hpp

@@ -121,13 +121,11 @@ void checkRelativeError( real64 const v1, real64 const v2, real64 const relTol,
   EXPECT_PRED_FORMAT4( checkRelativeErrorFormat, v1, v2, relTol, absTol );
 }
 
-template< typename ROW_INDEX, typename COL_INDEX, typename VALUE >
-void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE const absTol,
+template< typename COL_INDEX, typename VALUE >
+void compareMatrixRow( VALUE const relTol, VALUE const absTol,
                        localIndex const length1, COL_INDEX const * const indices1, VALUE const * const values1,
                        localIndex const length2, COL_INDEX const * const indices2, VALUE const * const values2 )
 {
-  SCOPED_TRACE( "Row " + std::to_string( rowNumber ));
-
   EXPECT_EQ( length1, length2 );
 
   for( localIndex j1 = 0, j2 = 0; j1 < length1 && j2 < length2; ++j1, ++j2 )
@@ -152,17 +150,31 @@ void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE cons
   }
 }
 
-template< typename ROW_INDEX, typename COL_INDEX, typename VALUE >
-void compareMatrixRow( ROW_INDEX const rowNumber, VALUE const relTol, VALUE const absTol,
-                       arraySlice1d< COL_INDEX const > indices1, arraySlice1d< VALUE const > values1,
-                       arraySlice1d< COL_INDEX const > indices2, arraySlice1d< VALUE const > values2 )
+template< typename T, typename COL_INDEX >
+void compareLocalMatrices( CRSMatrixView< T const, COL_INDEX const > const & matrix1,
+                           CRSMatrixView< T const, COL_INDEX const > const & matrix2,
+                           real64 const relTol = DEFAULT_REL_TOL,
+                           real64 const absTol = DEFAULT_ABS_TOL,
+                           globalIndex const rowOffset = 0 )
 {
-  ASSERT_EQ( indices1.size(), values1.size() );
-  ASSERT_EQ( indices2.size(), values2.size() );
+  ASSERT_EQ( matrix1.numRows(), matrix2.numRows() );
+  ASSERT_EQ( matrix1.numColumns(), matrix2.numColumns() );
+
+  matrix1.move( hostMemorySpace, false );
+  matrix2.move( hostMemorySpace, false );
 
-  compareMatrixRow( rowNumber, relTol, absTol,
-                    indices1.size(), indices1.dataIfContiguous(), values1.dataIfContiguous(),
-                    indices2.size(), indices2.dataIfContiguous(), values2.dataIfContiguous() );
+  // check the accuracy across local rows
+  for( localIndex i = 0; i < matrix1.numRows(); ++i )
+  {
+    SCOPED_TRACE( GEOS_FMT( "Row {}", i + rowOffset ) );
+    compareMatrixRow( relTol, absTol,
+                      matrix1.numNonZeros( i ),
+                      matrix1.getColumns( i ).dataIfContiguous(),
+                      matrix1.getEntries( i ).dataIfContiguous(),
+                      matrix2.numNonZeros( i ),
+                      matrix2.getColumns( i ).dataIfContiguous(),
+                      matrix2.getEntries( i ).dataIfContiguous() );
+  }
 }
 
 template< typename MATRIX >
@@ -177,45 +189,10 @@ void compareMatrices( MATRIX const & matrix1,
   ASSERT_EQ( matrix1.numLocalRows(), matrix2.numLocalRows() );
   ASSERT_EQ( matrix1.numLocalCols(), matrix2.numLocalCols() );
 
-  array1d< globalIndex > indices1, indices2;
-  array1d< real64 > values1, values2;
-
-  // check the accuracy across local rows
-  for( globalIndex i = matrix1.ilower(); i < matrix1.iupper(); ++i )
-  {
-    indices1.resize( matrix1.rowLength( i ) );
-    values1.resize( matrix1.rowLength( i ) );
-    matrix1.getRowCopy( i, indices1, values1 );
-
-    indices2.resize( matrix2.rowLength( i ) );
-    values2.resize( matrix2.rowLength( i ) );
-    matrix2.getRowCopy( i, indices2, values2 );
+  CRSMatrix< real64, globalIndex > const mat1 = matrix1.extract();
+  CRSMatrix< real64, globalIndex > const mat2 = matrix2.extract();
 
-    compareMatrixRow( i, relTol, absTol,
-                      indices1.size(), indices1.data(), values1.data(),
-                      indices2.size(), indices2.data(), values2.data() );
-  }
-}
-
-template< typename T, typename COL_INDEX >
-void compareLocalMatrices( CRSMatrixView< T const, COL_INDEX const > const & matrix1,
-                           CRSMatrixView< T const, COL_INDEX const > const & matrix2,
-                           real64 const relTol = DEFAULT_REL_TOL,
-                           real64 const absTol = DEFAULT_ABS_TOL )
-{
-  ASSERT_EQ( matrix1.numRows(), matrix2.numRows() );
-  ASSERT_EQ( matrix1.numColumns(), matrix2.numColumns() );
-
-  matrix1.move( hostMemorySpace, false );
-  matrix2.move( hostMemorySpace, false );
-
-  // check the accuracy across local rows
-  for( localIndex i = 0; i < matrix1.numRows(); ++i )
-  {
-    compareMatrixRow( i, relTol, absTol,
-                      matrix1.getColumns( i ), matrix1.getEntries( i ),
-                      matrix2.getColumns( i ), matrix2.getEntries( i ) );
-  }
+  compareLocalMatrices( mat1.toViewConst(), mat2.toViewConst(), relTol, absTol, matrix1.ilower() );
 }
 
 } // namespace testing

src/coreComponents/common/DataTypes.hpp

@@ -342,12 +342,12 @@ template< typename COL_INDEX, typename INDEX_TYPE=localIndex >
 using SparsityPatternView = LvArray::SparsityPatternView< COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer >;
 
 /// Alias for CRS Matrix class.
-template< typename T, typename COL_INDEX=globalIndex >
-using CRSMatrix = LvArray::CRSMatrix< T, COL_INDEX, localIndex, LvArray::ChaiBuffer >;
+template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex >
+using CRSMatrix = LvArray::CRSMatrix< T, COL_INDEX, INDEX_TYPE, LvArray::ChaiBuffer >;
 
 /// Alias for CRS Matrix View.
-template< typename T, typename COL_INDEX=globalIndex >
-using CRSMatrixView = LvArray::CRSMatrixView< T, COL_INDEX, localIndex const, LvArray::ChaiBuffer >;
+template< typename T, typename COL_INDEX=globalIndex, typename INDEX_TYPE=localIndex >
+using CRSMatrixView = LvArray::CRSMatrixView< T, COL_INDEX, INDEX_TYPE const, LvArray::ChaiBuffer >;
 
 ///@}
 

src/coreComponents/common/GeosxConfig.hpp.in

@@ -83,6 +83,12 @@
 /// Enables use of PETSc library (CMake option ENABLE_PETSC)
 #cmakedefine GEOSX_USE_PETSC
 
+///  Enables use of METIS library (CMake option ENABLE_METIS)
+#cmakedefine GEOSX_USE_METIS
+
+///  Enables use of ParMETIS library (CMake option ENABLE_PARMETIS)
+#cmakedefine GEOSX_USE_PARMETIS
+
 /// Enables use of Scotch library (CMake option ENABLE_SCOTCH)
 #cmakedefine GEOSX_USE_SCOTCH
 

src/coreComponents/common/MpiWrapper.hpp

@@ -22,6 +22,8 @@
 #include "common/Span.hpp"
 #include "mesh/ElementType.hpp"
 
+#include <numeric>
+
 #if defined(GEOSX_USE_MPI)
   #include <mpi.h>
 #define MPI_PARAM( x ) x
@@ -343,6 +345,11 @@ struct MpiWrapper
                         array1d< T > & recvbuf,
                         MPI_Comm comm = MPI_COMM_GEOSX );
 
+  template< typename T >
+  static int allGatherv( arrayView1d< T const > const & sendbuf,
+                         array1d< T > & recvbuf,
+                         MPI_Comm comm = MPI_COMM_GEOSX );
+
   /**
    * @brief Strongly typed wrapper around MPI_Allreduce.
    * @param[in] sendbuf The pointer to the sending buffer.
@@ -522,6 +529,22 @@ struct MpiWrapper
                     MPI_Comm comm,
                     MPI_Request * request );
 
+  /**
+   * @brief Strongly typed wrapper around MPI_Send()
+   * @param[in] buf The pointer to the buffer that contains the data to be sent.
+   * @param[in] count The number of elements in \p buf.
+   * @param[in] dest The rank of the destination process within \p comm.
+   * @param[in] tag The message tag that is be used to distinguish different types of messages.
+   * @param[in] comm The handle to the MPI_Comm.
+   * @return
+   */
+  template< typename T >
+  static int send( T const * const buf,
+                   int count,
+                   int dest,
+                   int tag,
+                   MPI_Comm comm );
+
   /**
    * @brief Strongly typed wrapper around MPI_Isend()
    * @param[in] buf The pointer to the buffer that contains the data to be sent.
@@ -754,6 +777,38 @@ int MpiWrapper::allGather( arrayView1d< T const > const & sendValues,
 #endif
 }
 
+template< typename T >
+int MpiWrapper::allGatherv( arrayView1d< T const > const & sendValues,
+                            array1d< T > & allValues,
+                            MPI_Comm MPI_PARAM( comm ) )
+{
+  int const sendSize = LvArray::integerConversion< int >( sendValues.size() );
+#ifdef GEOSX_USE_MPI
+  int const mpiSize = commSize( comm );
+  array1d< int > counts;
+  allGather( sendSize, counts, comm );
+  array1d< int > displs( mpiSize + 1 );
+  std::partial_sum( counts.begin(), counts.end(), displs.begin() + 1 );
+  allValues.resize( displs.back() );
+  return MPI_Allgatherv( sendValues.data(),
+                         sendSize,
+                         internal::getMpiType< T >(),
+                         allValues.data(),
+                         counts.data(),
+                         displs.data(),
+                         internal::getMpiType< T >(),
+                         comm );
+
+#else
+  allValues.resize( sendSize );
+  for( localIndex a=0; a<sendSize; ++a )
+  {
+    allValues[a] = sendValues[a];
+  }
+  return 0;
+#endif
+}
+
 template< typename T >
 int MpiWrapper::allReduce( T const * const sendbuf,
                            T * const recvbuf,
@@ -989,6 +1044,20 @@ int MpiWrapper::iSend( arrayView1d< T > const & buf,
 #endif
 }
 
+template< typename T >
+int MpiWrapper::send( T const * const buf,
+                      int count,
+                      int dest,
+                      int tag,
+                      MPI_Comm comm )
+{
+#ifdef GEOSX_USE_MPI
+  return MPI_Send( buf, count, internal::getMpiType< T >(), dest, tag, comm );
+#else
+  GEOS_ERROR( "Not implemented without MPI" );
+#endif
+}
+
 template< typename T >
 int MpiWrapper::iSend( T const * const buf,
                        int count,

src/coreComponents/common/TimingMacros.hpp

@@ -34,7 +34,7 @@ namespace timingHelpers
   {
     std::string input(prettyFunction);
     std::string::size_type const end = input.find_first_of( '(' );
-    std::string::size_type const beg = input.find_last_of( ' ', end)+1;
+    std::string::size_type const beg = input.find_last_of( ' ', end ) + 1;
     return input.substr( beg, end-beg );
   }
 }
@@ -95,10 +95,10 @@ namespace timingHelpers
 #include <iostream>
 
 /// Mark a function or scope for timing with a given name
-#define GEOS_CALIPER_MARK_SCOPE(name) cali::Function __cali_ann##__LINE__(STRINGIZE_NX(name))
+#define GEOS_CALIPER_MARK_SCOPE(name) cali::Function GEOS_CONCAT(_cali_ann, __LINE__)(STRINGIZE_NX(name))
 
 /// Mark a function for timing using a compiler-provided name
-#define GEOS_CALIPER_MARK_FUNCTION cali::Function __cali_ann##__func__(timingHelpers::stripPF(__PRETTY_FUNCTION__).c_str())
+#define GEOS_CALIPER_MARK_FUNCTION cali::Function _cali_ann_func(timingHelpers::stripPF(__PRETTY_FUNCTION__).c_str())
 
 /// Mark the beginning of timed statement group
 #define GEOS_CALIPER_MARK_BEGIN(name) CALI_MARK_BEGIN(STRINGIZE(name))

src/coreComponents/linearAlgebra/CMakeLists.txt

@@ -19,6 +19,8 @@ set( linearAlgebra_headers
      solvers/PreconditionerBlockJacobi.hpp
      solvers/PreconditionerIdentity.hpp
      solvers/PreconditionerJacobi.hpp
+     solvers/PreconditionerNull.hpp
+     solvers/RichardsonSolver.hpp
      solvers/SeparateComponentPreconditioner.hpp
      utilities/Arnoldi.hpp
      utilities/BlockOperator.hpp
@@ -44,6 +46,7 @@ set( linearAlgebra_sources
      solvers/CgSolver.cpp
      solvers/GmresSolver.cpp
      solvers/KrylovSolver.cpp
+     solvers/RichardsonSolver.cpp
      solvers/SeparateComponentPreconditioner.cpp
      utilities/ReverseCutHillMcKeeOrdering.cpp )