Matrix class extension

src/linearAlgebra/AbstractMatrix.h

@@ -40,6 +40,93 @@ namespace dftefe
     template <typename ValueType, dftefe::utils::MemorySpace memorySpace>
     class AbstractMatrix
     {
+    public:
+      /**
+       * @brief The value setter interface for the matrix class. This setter
+       *        assumes the given pointer to ValueType to be a serially repeated
+       *        matrix and copy the data to the corresponding local owned
+       *        location trivially.
+       * @warning There is no boundary check in this function. The user should
+       *          be responsible that the size of data array should have d_m*d_n
+       *          size.
+       * @param data The pointer to the data to be copied from.
+       */
+      virtual void
+      setValues(const ValueType *data) = 0;
+
+      /**
+       * @brief The value setter interface for the matrix class. This setter
+       *        assumes the given pointer to ValueType to be a serially repeated
+       *        matrix and copy the data to the corresponding local owned
+       *        sub-matrix (i1:i2-1, j1:j2-1) trivially.
+       * @warning There is no boundary check in this function. The user should
+       *          be responsible that the size of data array should have
+       *          (j2-j1)*(i2-i1) size.
+       * @param i1 the global index of the first row of the submatrix.
+       * @param i2 one passes the global index of the last row of the submatrix.
+       * @param j1 the global index of the first column of the submatrix.
+       * @param j2 one passes the global index of the last column of the
+       *           submatrix.
+       * @param data The pointer to the data to be copied from. The user should
+       *             be responsible that the size of data array should have
+       *             (j2-j1)*(i2-i1) size.
+       */
+      virtual void
+      setValues(size_t           i1,
+                size_t           i2,
+                size_t           j1,
+                size_t           j2,
+                const ValueType *data) = 0;
+
+      //      /**
+      //       * @brief The value setter interface for the matrix class. This setter
+      //       *        allows data to be different on each processor. The
+      //       locally owned
+      //       *        data contains two parts: (1) values belongs to the
+      //       locally owned
+      //       *        portion of the matrix and (2) values belong to the
+      //       off-processor
+      //       *        portion of the matrix. This setter will distribute the
+      //       *        off-processor values correspondingly.
+      //       * @warning This routine currently is not optimized and fully
+      //       *          tested. Using this routing to assign values can
+      //       *          deteriorate the performance dramatically. This should
+      //       *          only be used for debugging purpose or not performance
+      //       *          relevant part of the code.
+      //       * @param i1 the global index of the first row of the submatrix.
+      //       * @param i2 one passes the global index of the last row of the submatrix.
+      //       * @param j1 the global index of the first column of the submatrix.
+      //       * @param j2 one passes the global index of the last column of the
+      //       *           submatrix.
+      //       * @param data The pointer to the data to be copied from. The user should
+      //       *             be responsible that the size of data array should
+      //       have
+      //       *             (j2-j1)*(i2-i1) size.
+      //       */
+      //      virtual void
+      //      setDistributedValues(size_t           i1,
+      //                           size_t           i2,
+      //                           size_t           j1,
+      //                           size_t           j2,
+      //                           const ValueType *data) = 0;
+      //
+      //      /**
+      //       * @brief The value setter interface for inserting a single value.
+      //       * @warning This routine is sub-optimal and can seriously deteriorate the
+      //       *          the performance. It should be only used when
+      //       necessary. Only
+      //       *          the processor which owns (i, j) element will be
+      //       inserting
+      //       *          value.
+      //       * @param i the row index of the value.
+      //       * @param j the column index of the value.
+      //       * @param d the value.
+      //       */
+      //      virtual void
+      //      setValue(size_t i, size_t j, ValueType d) = 0;
+
+      AbstractMatrix() = delete;
+
     protected:
       AbstractMatrix(size_t   m,
                      size_t   n,
@@ -51,7 +138,24 @@ namespace dftefe
       size_t                        d_m, d_n;
       size_t                        d_mb, d_nb, d_p, d_q;
       MPI_Comm                      d_comm;
+      int                           d_num_ranks;
       slate::BaseMatrix<ValueType> *d_baseMatrix = nullptr;
+
+      void
+      setValueSlateMatrix(slate::BaseMatrix<ValueType> *matrix,
+                          const ValueType *             data);
+
+      void
+      setValueSlateMatrix(slate::BaseMatrix<ValueType> *matrix,
+                          const ValueType *             data,
+                          int                           i1,
+                          int                           i2,
+                          int                           j1,
+                          int                           j2);
+
+
+      void
+      initSlateMatrix(slate::BaseMatrix<ValueType> *matrix, ValueType val = 0);
     };
 
   } // namespace linearAlgebra

src/linearAlgebra/AbstractMatrix.t.cpp

@@ -23,6 +23,8 @@
  * @author Ian C. Lin.
  */
 
+#include "AbstractMatrix.h"
+
 namespace dftefe
 {
   namespace linearAlgebra
@@ -44,6 +46,95 @@ namespace dftefe
       , d_mb(mb)
     {}
 
+    template <typename ValueType, dftefe::utils::MemorySpace memorySpace>
+    void
+    AbstractMatrix<ValueType, memorySpace>::setValues(const ValueType *data)
+    {
+      setValueSlateMatrix(d_baseMatrix, data);
+    }
+
+    template <typename ValueType, dftefe::utils::MemorySpace memorySpace>
+    void
+    AbstractMatrix<ValueType, memorySpace>::initSlateMatrix(
+      slate::BaseMatrix<ValueType> *matrix,
+      ValueType                     val)
+    {
+      for (int64_t j = 0, j_offset = 0; j < matrix->nt();
+           j_offset += matrix->tileNb(j++))
+        {
+          for (int64_t i = 0, i_offset = 0; i < matrix->mt();
+               i_offset += matrix->tileMb(i++))
+            {
+              if (matrix->tileIsLocal(i, j))
+                {
+                  matrix->at(i, j).set(val);
+                }
+            }
+        }
+    }
 
+    template <typename ValueType, dftefe::utils::MemorySpace memorySpace>
+    void
+    AbstractMatrix<ValueType, memorySpace>::setValueSlateMatrix(
+      slate::BaseMatrix<ValueType> *matrix,
+      const ValueType *             data)
+    {
+      for (int64_t j = 0, j_offset = 0; j < matrix->nt();
+           j_offset += matrix->tileNb(j++))
+        {
+          for (int64_t i = 0, i_offset = 0; i < matrix->mt();
+               i_offset += matrix->tileMb(i++))
+            {
+              if (matrix->tileIsLocal(i, j))
+                {
+                  slate::Tile<double> T = (*matrix)(i, j);
+                  // todo: check for transpose case (d_m and d_n)
+                  int64_t mb = T.mb(), nb = T.nb(),
+                          offset = i_offset + j_offset * d_m;
+                  lapack::lacpy(lapack::MatrixType::General,
+                                mb,
+                                nb,
+                                &data[offset],
+                                d_m,
+                                T.data(),
+                                mb);
+                }
+            }
+        }
+    }
+
+    template <typename ValueType, dftefe::utils::MemorySpace memorySpace>
+    void
+    AbstractMatrix<ValueType, memorySpace>::setValueSlateMatrix(
+      slate::BaseMatrix<ValueType> *matrix,
+      const ValueType *             data,
+      int                           i1,
+      int                           i2,
+      int                           j1,
+      int                           j2)
+    {
+      for (int64_t j = 0, j_offset = 0; j < matrix->nt();
+           j_offset += matrix->tileNb(j++))
+        {
+          for (int64_t i = 0, i_offset = 0; i < matrix->mt();
+               i_offset += matrix->tileMb(i++))
+            {
+              if (matrix->tileIsLocal(i, j))
+                {
+                  slate::Tile<double> T = (*matrix)(i, j);
+                  // todo: check for transpose case (d_m and d_n)
+                  int64_t mb = T.mb(), nb = T.nb(),
+                          offset = i_offset + j_offset * d_m;
+                  int64_t ii1    = lapack::lacpy(lapack::MatrixType::General,
+                                              mb,
+                                              nb,
+                                              &data[offset],
+                                              d_m,
+                                              T.data(),
+                                              mb);
+                }
+            }
+        }
+    }
   } // namespace linearAlgebra
 } // namespace dftefe

src/linearAlgebra/CMakeLists.txt

@@ -13,20 +13,29 @@ if (NOT TARGET dft-efe-linalg)
             QueueManager.cpp
             LinAlgOpContext.cpp
             VectorAttributes.cpp
-	    Vector.cpp
-	    MultiVector.cpp
+	          Vector.cpp
+	          MultiVector.cpp
             AbstractMatrix.cpp
             GeneralMatrix.cpp
             HermitianMatrix.cpp
             TriangularMatrix.cpp
             MatrixKernels.cu
             MatrixKernels.cpp
-            MatrixOperations.cpp)
+            MatrixOperations.cpp
+            MatVecOperations.cpp)
 
     add_library(dft-efe-linalg SHARED ${DFT-EFE-LINALG-SOURCES})
 
     include_directories(../)
 
+    find_package(deal.II 9.3.0 REQUIRED HINTS ${DEALII_PATH})
+    target_include_directories(dft-efe-linalg PUBLIC ${DEAL_II_INCLUDE_DIRS})
+    IF("${CMAKE_BUILD_TYPE}" STREQUAL "Release")
+        target_link_libraries (dft-efe-linalg PUBLIC ${DEAL_II_LIBRARIES_RELEASE})
+    ELSE()
+        target_link_libraries (dft-efe-linalg PUBLIC ${DEAL_II_LIBRARIES_DEBUG})
+    ENDIF()
+
     if (ENABLE_MPI)
         add_compile_definitions(DFTEFE_WITH_MPI)
         if (NOT MPI_FOUND)
@@ -55,5 +64,5 @@ if (NOT TARGET dft-efe-linalg)
     target_include_directories(dft-efe-linalg PUBLIC ${SLATE_DIR}/include)
     target_link_directories(dft-efe-linalg PUBLIC ${SLATE_DIR}/lib64)
 
-    target_link_libraries(dft-efe-linalg PUBLIC blaspp slate dft-efe-utils ${DFTEFE_LINEARALGEBRA_MPI_LIBRARIES})
+    target_link_libraries(dft-efe-linalg PUBLIC blaspp lapackpp slate gomp dft-efe-utils ${DFTEFE_LINEARALGEBRA_MPI_LIBRARIES})
 endif ()

src/linearAlgebra/DO_NOT_USE_ScalapckWrapper/README

@@ -0,0 +1,4 @@
+This directory is adapted from Scalapack wrapper in DFT-FE and Deal.ii code. It contains the scalapack wrapper when we
+really need a fall back option if SLATE does not work.
+
+This project currently does NOT support Scalapck.