Linear algebra changes

Author: klevzoff

.devcontainer/devcontainer.json

@@ -2,7 +2,7 @@
     "build": {
         "dockerfile": "Dockerfile",
         "args": {
-            "GEOS_TPL_TAG": "312-717"
+            "GEOS_TPL_TAG": "314-733"
         }
     },
     "runArgs": [

src/cmake/thirdparty/SetupGeosxThirdParty.cmake

@@ -581,27 +581,17 @@ endif()
 if(DEFINED SCOTCH_DIR)
     message(STATUS "SCOTCH_DIR = ${SCOTCH_DIR}")
 
-    find_and_import(NAME scotch
-                      INCLUDE_DIRECTORIES ${SCOTCH_DIR}/include
-                      LIBRARY_DIRECTORIES ${SCOTCH_DIR}/lib
-                      HEADER scotch.h
-                      LIBRARIES scotch scotcherr )
-
-    find_and_import(NAME ptscotch
-                      INCLUDE_DIRECTORIES ${SCOTCH_DIR}/include
-                      LIBRARY_DIRECTORIES ${SCOTCH_DIR}/lib
-                      DEPENDS scotch
-                      HEADER ptscotch.h
-                      LIBRARIES ptscotch ptscotcherr )
-
-    extract_version_from_header( NAME scotch
-                                 HEADER "${SCOTCH_DIR}/include/scotch.h"
-                                 MAJOR_VERSION_STRING "SCOTCH_VERSION"
-                                 MINOR_VERSION_STRING "SCOTCH_RELEASE"
-                                 PATCH_VERSION_STRING "SCOTCH_PATCHLEVEL")
-
-    set(ENABLE_SCOTCH ON CACHE BOOL "")
-    set(thirdPartyLibs ${thirdPartyLibs} scotch ptscotch)
+    find_package(SCOTCH REQUIRED
+                 PATHS ${SCOTCH_DIR}
+                 NO_DEFAULT_PATH)
+
+    message( " ----> scotch_VERSION = ${SCOTCH_VERSION}")
+
+    foreach(_scotch_target SCOTCH::scotch SCOTCH::ptscotch SCOTCH::scotcherr )
+        get_target_property(SCOTCH_INCLUDE_DIRS ${_scotch_target} INTERFACE_INCLUDE_DIRECTORIES)
+        set_target_properties(${_scotch_target} PROPERTIES INTERFACE_SYSTEM_INCLUDE_DIRECTORIES "${SCOTCH_INCLUDE_DIRS}")
+        set(thirdPartyLibs ${thirdPartyLibs} ${_scotch_target})
+    endforeach()
 else()
     if(ENABLE_SCOTCH)
         message(WARNING "ENABLE_SCOTCH is ON but SCOTCH_DIR isn't defined.")

src/coreComponents/codingUtilities/UnitTestUtilities.hpp

@@ -122,13 +122,11 @@ inline void checkRelativeError( real64 const v1, real64 const v2, real64 const r
   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 )
@@ -153,17 +151,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 >
@@ -178,45 +190,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

@@ -301,12 +301,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

@@ -98,6 +98,12 @@
 /// Enables use of PETSc library (CMake option ENABLE_PETSC)
 #cmakedefine GEOS_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 GEOS_USE_SCOTCH
 

src/coreComponents/common/MpiWrapper.hpp

@@ -382,6 +382,11 @@ struct MpiWrapper
                         array1d< T > & recvbuf,
                         MPI_Comm comm = MPI_COMM_GEOS );
 
+  template< typename T >
+  static int allGatherv( arrayView1d< T const > const & sendbuf,
+                         array1d< T > & recvbuf,
+                         MPI_Comm comm = MPI_COMM_GEOS );
+
   /**
    * @brief Convenience wrapper for the MPI_Allreduce function.
    * @tparam T type of data to reduce. Must correspond to a valid MPI_Datatype.
@@ -619,6 +624,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.
@@ -985,6 +1006,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 GEOS_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,
@@ -1239,6 +1292,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 GEOS_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

@@ -35,7 +35,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 );
   }
 }
@@ -96,10 +96,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

@@ -44,6 +44,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
@@ -70,6 +72,7 @@ set( linearAlgebra_sources
      solvers/CgSolver.cpp
      solvers/GmresSolver.cpp
      solvers/KrylovSolver.cpp
+     solvers/RichardsonSolver.cpp
      solvers/SeparateComponentPreconditioner.cpp
      utilities/ReverseCutHillMcKeeOrdering.cpp )