From 28811b14cac36b4bf6c4203525f517c65aab9ce6 Mon Sep 17 00:00:00 2001
From: MARAIT Gilles <gilles.marait@inria.fr>
Date: Fri, 2 Feb 2024 11:40:50 +0100
Subject: [PATCH] Add FEM with compose
---
src/testCOMPOSE.cpp | 222 +++++++++++++++++++++++++++++++++++++++-----
1 file changed, 201 insertions(+), 21 deletions(-)
diff --git a/src/testCOMPOSE.cpp b/src/testCOMPOSE.cpp
index 0c3c9f8..69736c5 100644
--- a/src/testCOMPOSE.cpp
+++ b/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;
--
GitLab