Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 1a6a3a54 authored by FELŠÖCI Marek's avatar FELŠÖCI Marek
Browse files

Keep only the linear system solving test

parent b589b872
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
include/main.h
+ 0
− 11
View file @ 1a6a3a54
......@@ -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)
......
src/help.c
+ 0
− 5
View file @ 1a6a3a54
......@@ -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"
......
src/main.c
+ 3
− 48
View file @ 1a6a3a54
......@@ -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;
......
src/testHMAT.c
+ 117
− 199
View file @ 1a6a3a54
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment