From 1a6a3a54cc5a82ca7f4be14de74a3e75091d1e2d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?FEL=C5=A0=C3=96CI=20Marek?= <marek.felsoci@inria.fr>
Date: Sun, 5 Nov 2023 18:06:51 +0100
Subject: [PATCH] Keep only the linear system solving test
---
include/main.h | 11 --
src/help.c | 5 -
src/main.c | 51 +-------
src/testHMAT.c | 316 ++++++++++++++++++-------------------------------
4 files changed, 120 insertions(+), 263 deletions(-)
diff --git a/include/main.h b/include/main.h
index 3cc290f..38ba123 100644
--- a/include/main.h
+++ b/include/main.h
@@ -45,17 +45,6 @@ extern ScalarType stype;
/*! \brief Wavelength (for oscillatory kernels). */
extern double lambda;
-/*! \brief List of algorithms that we can test. */
-enum algo {
- _undefined,
- _gemvHMAT,
- _solveHMAT,
- _inverseHMAT
-} ;
-
-/*! \brief Algorithms that we want to test. */
-extern enum algo testedAlgo;
-
/*! \brief Flag to use simple assembly routine for HMAT (default is 0).
It can be changed with --use_simple_assembly (read in readFlagsTestHMAT() in hmat.c)
diff --git a/src/help.c b/src/help.c
index 6a71934..1fd7402 100644
--- a/src/help.c
+++ b/src/help.c
@@ -24,11 +24,6 @@ int printHelp() {
" -height: Height of the cylinder. Default is 4.\n"
" For FEM/BEM, the BEM is done on the outer surface of the cylinder (without the upper and lower disks).\n"
"\n"
- " Choice of the tested algorithm (no default):\n"
- " -gemvhmat: gemv HMAT.\n"
- " -solvehmat: solve HMAT.\n"
- " -inversehmat: inverse HMAT.\n"
- "\n"
" Options for HMAT tests\n"
" --use_simple_assembly: Using simple assembly routine for HMAT. Default: NO.\n"
" --hmat_divider: Value of the Hmat clustering divider (number of children per box). Default: 2.\n"
diff --git a/src/main.c b/src/main.c
index 37d3c6a..83f8bd4 100644
--- a/src/main.c
+++ b/src/main.c
@@ -19,7 +19,6 @@ int simplePrec = 0;
int complexALGO = 1;
ScalarType stype = DOUBLE_COMPLEX;
double lambda;
-enum algo testedAlgo = _undefined;
int use_simple_assembly = 0;
int divider = 2;
int hmat_residual = 0;
@@ -44,31 +43,6 @@ int main(int argc, char **argv) {
ierr=printHelp() ;
return ierr;
}
- /* --- Choix de l'algo a tester (_undefined par defaut) --- */
- if (MpfArgHasName(&argc, argv, 1, "-gemvhmat") > 0) {
- testedAlgo = _gemvHMAT ;
- printf("Testing : gemv HMAT.\n") ;
- }
- if (MpfArgHasName(&argc, argv, 1, "-solvehmat") > 0) {
- testedAlgo = _solveHMAT ;
- printf("Testing : solve HMAT.\n") ;
- }
- if (MpfArgHasName(&argc, argv, 1, "-inversehmat") > 0) {
- testedAlgo = _inverseHMAT ;
- printf("Testing : inverse HMAT.\n") ;
- }
- /* Check the test */
- switch(testedAlgo) {
- case _gemvHMAT:
- case _solveHMAT:
- case _inverseHMAT:
- // nothing
- break;
- case _undefined:
- default:
- SETERRQ(1, "No algorithm selected. Please choose one (use -h/--help for usage).");
- break;
- }
/* ------------------------------------------------------------------------- */
/* ---------------------- INITIATION DE MPF (PARALLELE) -------------------*/
@@ -76,18 +50,8 @@ int main(int argc, char **argv) {
ierr=SCAB_Init(&argc, &argv) ; CHKERRQ(ierr) ;
- switch(testedAlgo) {
- case _gemvHMAT:
- case _solveHMAT:
- case _inverseHMAT:
- ierr=readFlagsTestHMAT(&argc, &argv) ; CHKERRQ(ierr) ;
- ierr = init_hmat_interface(); CHKERRQ(ierr);
- break;
- case _undefined:
- default:
- SETERRQ(1, "No algorithm selected. Please choose one (use -h/--help for usage).");
- break;
- }
+ ierr=readFlagsTestHMAT(&argc, &argv) ; CHKERRQ(ierr) ;
+ ierr = init_hmat_interface(); CHKERRQ(ierr);
/* ------------------------------------------------------------------------- */
/* --------------------- DECODAGE DE LA LIGNE D'OPTIONS -------------------*/
@@ -153,16 +117,7 @@ int main(int argc, char **argv) {
ierr = prepareTEST() ; CHKERRQ(ierr) ;
double relative_error;
/* Run the test */
- switch(testedAlgo) {
- case _gemvHMAT:
- case _solveHMAT:
- case _inverseHMAT:
- ierr = testHMAT(&relative_error) ; CHKERRQ(ierr) ;
- break;
- default:
- SETERRQ(1, "Unknown algorithm=%d", testedAlgo);
- break;
- }
+ ierr = testHMAT(&relative_error) ; CHKERRQ(ierr) ;
int error_exit = 0;
if(relative_error > max_error) {
error_exit = 1;
diff --git a/src/testHMAT.c b/src/testHMAT.c
index 55e01f8..bd989a4 100644
--- a/src/testHMAT.c
+++ b/src/testHMAT.c
@@ -57,219 +57,137 @@ int testHMAT(double * relative_error) {
printf("Ratio: %g\n", (double)info.compressed_size/info.uncompressed_size);
}
- /* Switch on the various HMAT tests */
+ /* Perform a HMAT solve */
+ printf("\n") ;
+
+ hmat_matrix_t* hmatrix_orig = NULL;
+ void *solCLA_orig = NULL;
+ void *solCLA_sum = NULL;
+ double normRes, normOrig, temps_refine=-getTime();
+ int nbIter = 0;
+ if (hmat_residual) {
+ hmatrix_orig = hi->copy(hmatrix); // ideally, we do a copy with lower accuracy for factorisation
+ /* Vector to duplicate solCLA */
+ solCLA_orig = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_orig) ;
+ /* Norm of the originaal RHS */
+ ierr = computeVecNorm(solCLA, &normOrig); CHKERRQ(ierr);
+ }
+ if (hmat_refine) {
+ /* Vector to accumulate solCLA */
+ solCLA_sum = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_sum) ;
+ }
+
+ /* Factorize the H-matrix */
+
+ printf("\n**** Factorizing HMAT...\n") ;
+ temps_initial = getTime ();
+
+ hmat_factorization_context_t ctx;
+ hmat_factorization_context_init(&ctx);
+ ctx.progress = &mpf_hmat_settings.progress;
+ ctx.progress->user_data = NULL;
+ ctx.factorization = mpf_hmat_settings.factorization_type;
+ ierr = hi->factorize_generic(hmatrix, &ctx); CHKERRQ(ierr);
+
+ temps_final = getTime ();
+ temps_cpu = (temps_final - temps_initial) ;
+ printf("<PERFTESTS> TpsCpuFacto%s = %f \n", postfix_async, temps_cpu) ;
- switch(testedAlgo) {
- case _gemvHMAT:
+ temps_refine -= temps_cpu; // We do not count facto in the iterative refinement timer
+
+ /* display some informations (synchronous) */
+ if (hmat_get_sync_exec()) {
+ hmat_info_t info;
+ hi->get_info(hmatrix, &info);
+ printf("Compressed size: %ld\n", info.compressed_size);
+ printf("<PERFTESTS> FactorizedSizeMb = %f \n", (double)info.compressed_size*Mpf_stype_size[stype]/1024./1024.) ;
+ printf("Uncompressed size: %ld\n", info.uncompressed_size);
+ printf("Ratio: %g\n", (double)info.compressed_size/info.uncompressed_size);
+ }
- /* Appelle le produit mat.vec : A.rhs -> solHMAT */
+ /* Solve the system : A-1.solCLA -> solCLA */
- printf("\n**** Computing HMAT product...\n") ;
- void *solHMAT = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solHMAT) ;
- temps_initial = getTime ();
+ printf("\n**** Solve HMAT%s...\n", hmat_refine ? " with Iterative Refinement" : "") ;
- if (simplePrec)
- doubleToSimple(rhs, vectorSize);
+_refineLoop: // Iterative Refinement loop
+ nbIter++ ;
+ if (hmat_residual) {
+ /* Duplicate the RHS solCLA into solCLA_orig */
+ memcpy(solCLA_orig, solCLA, vectorSize*sizeof(D_type));
+ }
- /* Compute the mat.vec product A.rhs -> solHMAT */
- ierr = hi->gemv('N', Mpf_pone[stype], hmatrix, rhs, Mpf_zero[stype], solHMAT, nbRHS);
+ temps_initial = getTime ();
- if (simplePrec) {
- simpleToDouble(solHMAT, vectorSize);
- simpleToDouble(rhs, vectorSize);
- }
+ if (simplePrec)
+ doubleToSimple(solCLA, vectorSize);
- if (solHMAT) MpfFree(solHMAT);
- solHMAT=NULL;
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuGEMV%s = %f \n", postfix_async, temps_cpu) ;
+ ierr = hi->solve_systems(hmatrix, solCLA, nbRHS); CHKERRQ(ierr);
- /* Compare les 2 produits matrice-vecteur */
+ if (simplePrec)
+ simpleToDouble(solCLA, vectorSize);
- printf("\n**** Comparing results...\n") ;
- ierr=computeRelativeError(solHMAT, solCLA, relative_error); CHKERRQ(ierr);
+ temps_final = getTime ();
+ temps_cpu = (temps_final - temps_initial) ;
+ printf("<PERFTESTS> TpsCpuSolve%s = %f \n", postfix_async, temps_cpu) ;
- printf("<PERFTESTS> Error = %.4e \n", *relative_error);
- if (solHMAT) MpfFree(solHMAT);
- solHMAT=NULL ;
- break;
+ /* Compute the residual solCLA_orig-A_orig.solCLA */
- case _solveHMAT:
+ if (hmat_residual) {
+ temps_initial = getTime ();
+ if (simplePrec) {
+ doubleToSimple(solCLA, vectorSize);
+ doubleToSimple(solCLA_orig, vectorSize);
+ }
+ ierr = hi->gemv('N', Mpf_mone[stype], hmatrix_orig, solCLA, Mpf_pone[stype], solCLA_orig, nbRHS);
+ if (simplePrec) {
+ simpleToDouble(solCLA, vectorSize);
+ simpleToDouble(solCLA_orig, vectorSize);
+ }
+ temps_final = getTime ();
+ temps_cpu = (temps_final - temps_initial) ;
+ printf("<PERFTESTS> TpsCpuGEMV%s = %f \n", postfix_async, temps_cpu) ;
+ ierr = computeVecNorm(solCLA_orig, &normRes); CHKERRQ(ierr);
+ printf("<PERFTESTS> Residual_%d = %.4e \n", nbIter, normRes/normOrig);
+ // Test if convergence is achieved or if iteration number is too large
+ if (normRes/normOrig < mpf_hmat_settings.acaEpsilon/10 || nbIter > 20) {
+ // TODO stop the computation if the residual norm increases
printf("\n") ;
+ printf("<PERFTESTS> NbIterRefinement = %d \n", nbIter);
+ hmat_refine=0;
+ }
+ }
+ if (solCLA_sum) {
+ /* Compute the updated solution solCLA_sum += solCLA */
+ size_t k;
+ for (k=0 ; k<vectorSize ; k++)
+ ((double*)solCLA_sum)[k] += ((double*)solCLA)[k];
+ ierr=computeRelativeError(solCLA_sum, rhs, relative_error) ; CHKERRQ(ierr) ;
+ printf("<PERFTESTS> Error_%d = %.4e \n", nbIter, *relative_error);
+ }
+ if (hmat_refine) {
+ /* Copy the updated RHS of the system solCLA_orig into solCLA */
+ memcpy(solCLA, solCLA_orig, vectorSize*sizeof(D_type));
+ goto _refineLoop;
+ }
- hmat_matrix_t* hmatrix_orig = NULL;
- void *solCLA_orig = NULL;
- void *solCLA_sum = NULL;
- double normRes, normOrig, temps_refine=-getTime();
- int nbIter = 0;
- if (hmat_residual) {
- hmatrix_orig = hi->copy(hmatrix); // ideally, we do a copy with lower accuracy for factorisation
- /* Vector to duplicate solCLA */
- solCLA_orig = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_orig) ;
- /* Norm of the originaal RHS */
- ierr = computeVecNorm(solCLA, &normOrig); CHKERRQ(ierr);
- }
- if (hmat_refine) {
- /* Vector to accumulate solCLA */
- solCLA_sum = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA_sum) ;
- }
-
- /* Factorize the H-matrix */
-
- printf("\n**** Factorizing HMAT...\n") ;
- temps_initial = getTime ();
-
- hmat_factorization_context_t ctx;
- hmat_factorization_context_init(&ctx);
- ctx.progress = &mpf_hmat_settings.progress;
- ctx.progress->user_data = NULL;
- ctx.factorization = mpf_hmat_settings.factorization_type;
- ierr = hi->factorize_generic(hmatrix, &ctx); CHKERRQ(ierr);
-
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuFacto%s = %f \n", postfix_async, temps_cpu) ;
-
- temps_refine -= temps_cpu; // We do not count facto in the iterative refinement timer
-
- /* display some informations (synchronous) */
- if (hmat_get_sync_exec()) {
- hmat_info_t info;
- hi->get_info(hmatrix, &info);
- printf("Compressed size: %ld\n", info.compressed_size);
- printf("<PERFTESTS> FactorizedSizeMb = %f \n", (double)info.compressed_size*Mpf_stype_size[stype]/1024./1024.) ;
- printf("Uncompressed size: %ld\n", info.uncompressed_size);
- printf("Ratio: %g\n", (double)info.compressed_size/info.uncompressed_size);
- }
-
- /* Solve the system : A-1.solCLA -> solCLA */
-
- printf("\n**** Solve HMAT%s...\n", hmat_refine ? " with Iterative Refinement" : "") ;
+ if (hmatrix_orig) hi->destroy(hmatrix_orig);
+ hmatrix_orig=NULL;
+ if (solCLA_orig) MpfFree(solCLA_orig);
+ solCLA_orig=NULL ;
+ if (solCLA_sum) { // If solCLA_sum has been allocated, I use it as a replacement for solCLA
+ MpfFree(solCLA) ;
+ solCLA=solCLA_sum ;
+ temps_refine += getTime();
+ printf("<PERFTESTS> TpsCpuRefinement%s = %f \n", postfix_async, temps_refine) ;
-_refineLoop: // Iterative Refinement loop
- nbIter++ ;
- if (hmat_residual) {
- /* Duplicate the RHS solCLA into solCLA_orig */
- memcpy(solCLA_orig, solCLA, vectorSize*sizeof(D_type));
- }
-
- temps_initial = getTime ();
-
- if (simplePrec)
- doubleToSimple(solCLA, vectorSize);
-
- ierr = hi->solve_systems(hmatrix, solCLA, nbRHS); CHKERRQ(ierr);
-
- if (simplePrec)
- simpleToDouble(solCLA, vectorSize);
-
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuSolve%s = %f \n", postfix_async, temps_cpu) ;
-
- /* Compute the residual solCLA_orig-A_orig.solCLA */
-
- if (hmat_residual) {
- temps_initial = getTime ();
- if (simplePrec) {
- doubleToSimple(solCLA, vectorSize);
- doubleToSimple(solCLA_orig, vectorSize);
- }
- ierr = hi->gemv('N', Mpf_mone[stype], hmatrix_orig, solCLA, Mpf_pone[stype], solCLA_orig, nbRHS);
- if (simplePrec) {
- simpleToDouble(solCLA, vectorSize);
- simpleToDouble(solCLA_orig, vectorSize);
- }
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuGEMV%s = %f \n", postfix_async, temps_cpu) ;
- ierr = computeVecNorm(solCLA_orig, &normRes); CHKERRQ(ierr);
- printf("<PERFTESTS> Residual_%d = %.4e \n", nbIter, normRes/normOrig);
- // Test if convergence is achieved or if iteration number is too large
- if (normRes/normOrig < mpf_hmat_settings.acaEpsilon/10 || nbIter > 20) {
- // TODO stop the computation if the residual norm increases
- printf("\n") ;
- printf("<PERFTESTS> NbIterRefinement = %d \n", nbIter);
- hmat_refine=0;
- }
- }
- if (solCLA_sum) {
- /* Compute the updated solution solCLA_sum += solCLA */
- size_t k;
- for (k=0 ; k<vectorSize ; k++)
- ((double*)solCLA_sum)[k] += ((double*)solCLA)[k];
- ierr=computeRelativeError(solCLA_sum, rhs, relative_error) ; CHKERRQ(ierr) ;
- printf("<PERFTESTS> Error_%d = %.4e \n", nbIter, *relative_error);
- }
- if (hmat_refine) {
- /* Copy the updated RHS of the system solCLA_orig into solCLA */
- memcpy(solCLA, solCLA_orig, vectorSize*sizeof(D_type));
- goto _refineLoop;
- }
-
- if (hmatrix_orig) hi->destroy(hmatrix_orig);
- hmatrix_orig=NULL;
- if (solCLA_orig) MpfFree(solCLA_orig);
- solCLA_orig=NULL ;
- if (solCLA_sum) { // If solCLA_sum has been allocated, I use it as a replacement for solCLA
- MpfFree(solCLA) ;
- solCLA=solCLA_sum ;
- temps_refine += getTime();
- printf("<PERFTESTS> TpsCpuRefinement%s = %f \n", postfix_async, temps_refine) ;
-
- }
-
- /* Compare the two vectors solCLA and rhs */
-
- printf("\n**** Comparing results...\n") ;
- ierr=computeRelativeError(solCLA, rhs, relative_error); CHKERRQ(ierr);
-
- printf("<PERFTESTS> Error = %.4e \n", *relative_error);
- break;
- case _inverseHMAT:
-
- /* Inverse the H-matrix */
-
- printf("\n**** Inversing HMAT...\n") ;
- temps_initial = getTime ();
-
- hi->inverse(hmatrix);
-
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuInversion%s = %f \n", postfix_async, temps_cpu) ;
-
- /* Compute the gemv : A-1.solCLA -> solCLA */
-
- printf("\n**** GEMV HMAT...\n") ;
- temps_initial = getTime ();
-
- if (simplePrec)
- doubleToSimple(solCLA, vectorSize);
-
- void *solCLA2 = MpfCalloc(vectorSize, sizeof(D_type)) ; CHKPTRQ(solCLA) ;
- memcpy(solCLA2, solCLA, vectorSize*sizeof(D_type));
- ierr = hi->gemv('N', Mpf_pone[stype], hmatrix, solCLA2, Mpf_zero[stype], solCLA, nbRHS);
- MpfFree(solCLA2);
-
- if (simplePrec)
- simpleToDouble(solCLA, vectorSize);
-
- temps_final = getTime ();
- temps_cpu = (temps_final - temps_initial) ;
- printf("<PERFTESTS> TpsCpuGEMV%s = %f \n", postfix_async, temps_cpu) ;
-
- /* Compare the two vectors solCLA and rhs */
-
- printf("\n**** Comparing results...\n") ;
- ierr=computeRelativeError(solCLA, rhs, relative_error); CHKERRQ(ierr);
-
- printf("<PERFTESTS> Error = %.4e \n", *relative_error);
- break;
- default:
- break;
- } /* switch(testedAlgo) */
+ }
+
+ /* Compare the two vectors solCLA and rhs */
+
+ printf("\n**** Comparing results...\n") ;
+ ierr=computeRelativeError(solCLA, rhs, relative_error); CHKERRQ(ierr);
+
+ printf("<PERFTESTS> Error = %.4e \n", *relative_error);
/* Destroys the HMAT structure */
#ifdef TESTFEMBEM_DUMP_HMAT
--
GitLab