Add FEM with compose

src/testCOMPOSE.cpp

@@ -7,6 +7,8 @@
 char * C_PACKAGE_NAME = const_cast<char *>(PACKAGE_NAME) ;
 char * C_FILE = const_cast<char *>(__FILE__);
 
+// Redefine the macros to be C++ compatible
+
 /*! \brief Displays an error message and that's all. */
 #define C_SETWARN(n,s, ...)   {MpfWarning(__LINE__,C_PACKAGE_NAME,C_FILE,__func__,n,s, ## __VA_ARGS__);}
 /*! \brief Displays an error message and quits the current routine. */
@@ -26,20 +28,176 @@ char * C_FILE = const_cast<char *>(__FILE__);
 /*! \brief Evaluates the expression x, if it is false then it displays an error message and aborts the code. */
 #define C_ASSERTA(x) { if (!(x)) { C_SETERRA(1, "Assert failure %s", #x); }}
 
-
 #include <maphys.hpp>
 #include <memory>
 #include <cstring>
 #include <string>
+#include <span>
 
 using namespace maphys;
 
 template<typename Scalar>
-int testCOMPOSE_Scalar(double * relative_error){
+struct nonzero_type { using type = std::false_type; };
+template<> struct nonzero_type<float> { using type = SnonZero; };
+template<> struct nonzero_type<double> { using type = DnonZero; };
+template<> struct nonzero_type<std::complex<float>> { using type = CnonZero; };
+template<> struct nonzero_type<std::complex<double>> { using type = ZnonZero; };
+
+template<typename Scalar>
+int testCOMPOSE_FEM(double * relative_error){
+
+  using NonZero = typename nonzero_type<Scalar>::type;
 
   int ierr = 0;
   double temps_initial, temps_final, temps_cpu;
-  size_t vectorSize = (size_t)nbPts*nbRHS*(complexALGO?2:1);
+  size_t vectorSize = (size_t) nbPts*nbRHS*(complexALGO?2:1);
+
+  /* Realise le produit matrice vecteur classique : A.rhs -> solCLA */
+
+  Mpf_printf(MPI_COMM_WORLD, "\n**** Computing Classical product...\n") ;
+  void *solCLA=NULL;
+  ierr = produitClassique(&solCLA, MPI_COMM_WORLD) ; C_CHKERRQ(ierr) ;
+
+  const int N = nbPts;
+
+  auto myCtx = std::make_unique<contextTestFEMBEM>();
+  std::memset(myCtx.get(), 0, sizeof(contextTestFEMBEM));
+  myCtx->colDim = N;
+
+  temps_initial = getTime ();
+  Mpf_printf( MPI_COMM_WORLD, "\n**** Creating COMPOSE ...\n" );
+
+  int row_min = 1;
+  int row_max = N;
+  int col_min = 1;
+  int col_max = N;
+  int nnz;
+  NonZero * matComp = (NonZero *) malloc(sizeof(NonZero)); // For some reason, realloc is used on this
+
+  freeFemMatrix();
+  createFemMatrix();
+  computeSparseBlockFEMBEM(&row_min, &row_max, &col_min, &col_max,
+                           &Mpf_i_zero, &Mpf_i_one, &Mpf_i_one, &Mpf_i_one,
+                           &nnz, (void **) &matComp, (void *) myCtx.get());
+  freeFemMatrix();
+
+  std::span<NonZero> matspan(matComp, nnz);
+
+  IndexArray<int> Ai(nnz), Aj(nnz);
+  IndexArray<Scalar> Av(nnz);
+  int k = 0;
+  for(auto elt : matspan){
+    Ai[k] = elt.i; Aj[k] = elt.j;
+    if constexpr(is_complex<Scalar>::value) {
+      Av[k] =  Scalar{elt.v.r, elt.v.i};
+    }
+    else{
+      Av[k] = elt.v;
+    }
+    k++;
+  }
+
+  SparseMatrixCOO<Scalar> A(N, N, static_cast<size_t>(nnz), std::move(Ai), std::move(Aj), std::move(Av));
+
+  temps_final = getTime();
+  temps_cpu = temps_final - temps_initial ;
+  Mpf_printf(MPI_COMM_WORLD,"<PERFTESTS> TpsCpuCOMPOSECreation = %f \n", temps_cpu) ;
+
+  DenseMatrix<Scalar> B;
+  DenseMatrix<Scalar> X;
+  void *solCOMPOSE = nullptr;
+
+  switch(testedAlgo){
+  case _gemvCOMPOSE:
+
+    Mpf_printf(MPI_COMM_WORLD, "\n**** Computing COMPOSE product...\n") ;
+    temps_initial = getTime();
+
+    if (simplePrec)
+      doubleToSimple(rhs, vectorSize);
+
+    B = DenseMatrix<Scalar>::view(N, nbRHS, (Scalar*) rhs, N);
+
+    solCOMPOSE = MpfCalloc(vectorSize, sizeof(D_type)) ; C_CHKPTRQ(solCOMPOSE) ;
+    X = DenseMatrix<Scalar>::view(N, nbRHS, (Scalar*) solCOMPOSE, N);
+
+    X = A * B;
+
+    temps_final = getTime();
+    temps_cpu = temps_final - temps_initial ;
+
+    if (simplePrec) {
+      simpleToDouble(solCOMPOSE, vectorSize);
+      simpleToDouble(rhs, vectorSize);
+    }
+
+    Mpf_printf(MPI_COMM_WORLD,"<PERFTESTS> TpsCpuGEMVCOMPOSE= %f \n", temps_cpu) ;
+
+    // Compare les 2 produits matrice-vecteur
+
+    Mpf_printf(MPI_COMM_WORLD, "\n**** Comparing results...\n") ;
+    ierr=computeRelativeError(solCOMPOSE, solCLA, relative_error); C_CHKERRQ(ierr);
+
+    Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> Error = %.4e \n", *relative_error);
+
+    if(solCOMPOSE) MpfFree(solCOMPOSE);
+
+    break;
+
+  case _solveCOMPOSE:
+    {
+      Mpf_printf(MPI_COMM_WORLD, "\n\n**** Factorizing COMPOSE Mat...\n") ;
+      temps_initial = getTime();
+
+      Pastix<SparseMatrixCOO<Scalar>, DenseMatrix<Scalar>> solver(A);
+      //Mumps<SparseMatrixCOO<Scalar>, DenseMatrix<Scalar>> solver(A);
+      solver.factorize();
+
+      temps_final = getTime();
+      temps_cpu = temps_final - temps_initial;
+
+      Mpf_printf(MPI_COMM_WORLD,"<PERFTESTS> TpsCpuCOMPOSEFacto = %f \n", temps_cpu) ;
+      temps_initial = getTime();
+
+      if (simplePrec)
+	doubleToSimple(solCLA, vectorSize);
+
+      B = DenseMatrix<Scalar>::view(N, nbRHS, (Scalar*) solCLA, N);
+      //B.display("B");
+
+      solCOMPOSE = MpfCalloc(vectorSize, sizeof(D_type)) ; C_CHKPTRQ(solCOMPOSE) ;
+      X = DenseMatrix<Scalar>::view(N, nbRHS, (Scalar*) solCOMPOSE, N);
+
+      X = solver * B;
+
+      temps_final = getTime();
+      temps_cpu = temps_final - temps_initial;
+
+      if (simplePrec) {
+	simpleToDouble(solCOMPOSE, vectorSize);
+      }
+
+      Mpf_printf(MPI_COMM_WORLD,"<PERFTESTS> TpsCpuSolveCOMPOSE = %f \n", temps_cpu);
+      // Compare the two vectors solCOMPOSE and rhs
+      Mpf_printf(MPI_COMM_WORLD, "\n**** Comparing results...\n") ;
+      ierr=computeRelativeError(solCOMPOSE, rhs, relative_error) ; C_CHKERRQ(ierr) ;
+      Mpf_printf(MPI_COMM_WORLD, "<PERFTESTS> Error = %.4e \n", *relative_error);
+    }
+
+    break;
+  default:
+    break;
+  }
+
+  return 0;
+}
+
+template<typename Scalar>
+int testCOMPOSE_BEM(double * relative_error){
+
+  int ierr = 0;
+  double temps_initial, temps_final, temps_cpu;
+  size_t vectorSize = (size_t) nbPts*nbRHS*(complexALGO?2:1);
 
   /* Realise le produit matrice vecteur classique : A.rhs -> solCLA */
 
@@ -67,8 +225,8 @@ int testCOMPOSE_Scalar(double * relative_error){
   (simplePrec ?
    computeDenseBlockFEMBEM_Cprec :
    computeDenseBlockFEMBEM)(&row_min, &row_max, &col_min, &col_max,
-				  &Mpf_i_zero, &Mpf_i_one, &Mpf_i_one, &Mpf_i_one,
-				  &size_of_buffer, &ld, (void *) get_ptr(A), (void *) myCtx.get());
+			    &Mpf_i_zero, &Mpf_i_one, &Mpf_i_one, &Mpf_i_one,
+			    &size_of_buffer, &ld, (void *) get_ptr(A), (void *) myCtx.get());
 
   if(symMatSolver){
     if(symMatContent == 1){
@@ -188,22 +346,44 @@ int testCOMPOSE_Scalar(double * relative_error){
 int testCOMPOSE(double * relative_error){
 
   int ierr = 0;
-  switch(stype) {
-  case (SIMPLE_PRECISION) :
-    ierr = testCOMPOSE_Scalar<float>(relative_error);
-    break ;
-  case (DOUBLE_PRECISION) :
-    ierr = testCOMPOSE_Scalar<double>(relative_error);
-    break ;
-  case (SIMPLE_COMPLEX) :
-    ierr = testCOMPOSE_Scalar<std::complex<float>>(relative_error);
-    break ;
-  case (DOUBLE_COMPLEX) :
-    ierr = testCOMPOSE_Scalar<std::complex<double>>(relative_error);
-    break ;
-  default :
-    C_SETERRQ(1, "testCOMPOSE : unknown scalar type\n") ;
-    break ;
+
+  if(testedModel == _bem){
+    switch(stype) {
+    case (SIMPLE_PRECISION) :
+      ierr = testCOMPOSE_BEM<float>(relative_error);
+      break ;
+    case (DOUBLE_PRECISION) :
+      ierr = testCOMPOSE_BEM<double>(relative_error);
+      break ;
+    case (SIMPLE_COMPLEX) :
+      ierr = testCOMPOSE_BEM<std::complex<float>>(relative_error);
+      break ;
+    case (DOUBLE_COMPLEX) :
+      ierr = testCOMPOSE_BEM<std::complex<double>>(relative_error);
+      break ;
+    default :
+      C_SETERRQ(1, "testCOMPOSE : unknown scalar type\n") ;
+      break ;
+    }
+  }
+  else{ // if(testedModel == _fem)
+    switch(stype) {
+    case (SIMPLE_PRECISION) :
+      ierr = testCOMPOSE_FEM<float>(relative_error);
+      break ;
+    case (DOUBLE_PRECISION) :
+      ierr = testCOMPOSE_FEM<double>(relative_error);
+      break ;
+    case (SIMPLE_COMPLEX) :
+      ierr = testCOMPOSE_FEM<std::complex<float>>(relative_error);
+      break ;
+    case (DOUBLE_COMPLEX) :
+      ierr = testCOMPOSE_FEM<std::complex<double>>(relative_error);
+      break ;
+    default :
+      C_SETERRQ(1, "testCOMPOSE : unknown scalar type\n") ;
+      break ;
+    }
   }
 
   return 0;