From c4e7e77fc744940c4658999e5eb1f8243b92c3a8 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 19:03:05 +0100
Subject: [PATCH] Force double precision arithmetics
---
include/main.h | 15 +-------
include/util.h | 29 ----------------
src/classic.c | 17 +++------
src/compare.c | 20 ++---------
src/help.c | 3 --
src/hmat.c | 88 ++---------------------------------------------
src/kernel.c | 10 ++----
src/main.c | 32 -----------------
src/prepareTEST.c | 2 +-
src/rhs.c | 66 +++--------------------------------
src/testHMAT.c | 11 ++----
src/util.c | 34 ------------------
12 files changed, 22 insertions(+), 305 deletions(-)
diff --git a/include/main.h b/include/main.h
index e8b8c4c..3cd3e67 100644
--- a/include/main.h
+++ b/include/main.h
@@ -31,15 +31,6 @@ extern int nbPts;
/*! \brief RHS for my tests */
extern void *rhs;
-/*! \brief Simple precision accuracy flag */
-extern int simplePrec;
-
-/*! \brief Complex scalar type flag */
-extern int complexALGO;
-
-/*! \brief Scalar type used by the tests */
-extern ScalarType stype;
-
/*! \brief Wavelength (for oscillatory kernels). */
extern double lambda;
@@ -56,7 +47,7 @@ typedef struct {
int colDim;
} contextTestFEMBEM;
-int computeKernelBEM(double *coord1, double *coord2, int self, double _Complex *kernel) ;
+int computeKernelBEM(double *coord1, double *coord2, int self, double *kernel) ;
int produitClassiqueBEM(int ffar, void *sol) ;
int produitClassique(void **sol) ;
int computeRelativeError(void *sol, void *ref, double *eps) ;
@@ -66,16 +57,12 @@ int initCylinder(int *argc, char ***argv) ;
int getMeshStep(double *m) ;
int main(int argc, char **argv) ;
int computeRhs(void) ;
-void computeDenseBlockData (int *LigInf, int *LigSup, int *ColInf, int *ColSup, int *___Inf, int *___Sup, int *iloc, int *jloc, int *tailleS, int *tailleL, void *bloc, void *context);
double getTime() ;
int displayArray(char *s, void *f, int m, int n) ;
int testHMAT(double * relative_error) ;
int prepareTEST(void);
double* createCylinder(void) ;
-void interactionKernel(void* data, int i, int j, void* result) ;
int printHelp() ;
-void computeDenseBlockFEMBEM(int *LigInf, int *LigSup, int *ColInf, int *ColSup, int *___Inf, int *___Sup, int *iloc, int *jloc, int *tailleS, int *tailleL, void *bloc, void *context);
-void computeDenseBlockFEMBEM_Cprec(int *LigInf, int *LigSup, int *ColInf, int *ColSup, int *___Inf, int *___Sup, int *iloc, int *jloc, int *tailleS, int *tailleL, void *bloc, void *context);
void prepare_hmat(int, int, int, int, int*, int*, int*, int*, void*, hmat_block_info_t *);
void advanced_compute_hmat(struct hmat_block_compute_context_t*);
int init_hmat_interface() ;
diff --git a/include/util.h b/include/util.h
index 4c2e749..8b4742b 100644
--- a/include/util.h
+++ b/include/util.h
@@ -32,32 +32,6 @@
typedef enum { MPF_FALSE, MPF_TRUE } Logical;
-/**************************************************************************************/
-/*! \brief All the scalar types
-
- This enumeration defines the different scalar types available in MPF.
- They are transmitted to the vectors and matrices creation routine such as MatCreate().
-*/
-typedef enum {
- /*! \brief Simple real type (float in C, REAL in fortran, float in MPF) */
- SIMPLE_PRECISION=0,
- /*! \brief Double real type (double in C, REAL*8 in fortran, double in MPF) */
- DOUBLE_PRECISION=1,
- /*! \brief Simple complex type (doesn't exist in C, COMPLEX in fortran, float _Complex in MPF) */
- SIMPLE_COMPLEX=2,
- /*! \brief Double complex type (doesn't exist in C, DOUBLE COMPLEX in fortran, double _Complex in MPF) */
- DOUBLE_COMPLEX=3,
- /*! \brief Number of scalar types available. */
- nbScalarType=4
-} ScalarType;
-
-/*! \brief Name of the scalar types */
-_EXTERN_TEST_FEMBEM_ char *MPF_scalName[4]
-#ifdef ___TEST_FEMBEM___
-={"SIMPLE REAL", "DOUBLE REAL", "SIMPLE COMPLEX", "DOUBLE COMPLEX"}
- #endif
- ;
-#ifdef HAVE_HMAT
enum mpf_hmat_engine { mpf_hmat_seq, mpf_hmat_starpu, mpf_hmat_toyrt };
struct mpf_hmat_settings_t {
/** hmat interfaces for all scalar types */
@@ -80,7 +54,6 @@ struct mpf_hmat_create_compression_args_t {
};
int hmat_get_sync_exec(void) ;
void mpf_hmat_compression(struct mpf_hmat_create_compression_args_t *);
-#endif
#if !(defined (_WIN32)) /* It's a UNIX system, I know this ! */
/* Try to detect a debugger. This relies on the fact that GDB and IDB
@@ -180,7 +153,5 @@ Time time_interval(Time t1, Time t2);
Time add_times(Time t1, Time t2);
double time_in_s(Time t);
double time_interval_in_seconds(Time t1, Time t2);
-void simpleToDouble(void *data, size_t n) ;
-void doubleToSimple(void *data, size_t n) ;
int SCAB_Init(int* argc, char*** argv) ;
int SCAB_Exit(int* argc, char** argv) ;
diff --git a/src/classic.c b/src/classic.c
index dd13ec2..1938dce 100644
--- a/src/classic.c
+++ b/src/classic.c
@@ -23,8 +23,7 @@ int produitClassiqueBEM(int ffar, void *sol) {
for (int j=0 ; j<nbPts ; j+=sparseRHS) {
int i, k, ierr;
double coord_i[4], coord_j[3];
- double _Complex Aij ;
- double _Complex *z_sol=sol, *z_rhs=rhs ;
+ double Aij ;
double *d_sol=sol, *d_rhs=rhs ;
ierr=computeCoordCylinder(j, &(coord_j[0])) ; CHKERRA(ierr);
@@ -32,15 +31,9 @@ int produitClassiqueBEM(int ffar, void *sol) {
for (i=0 ; i<nbPts ; i++) {
ierr=computeCoordCylinder(i, &(coord_i[0])) ; CHKERRA(ierr);
{
- ierr=computeKernelBEM(&(coord_i[0]), &(coord_j[0]), i==j, (double _Complex *)&Aij) ;
- if (complexALGO) {
- for (k=0 ; k<nbRHS ; k++) {
- z_sol[(size_t)k*nbPts+i] += Aij * z_rhs[(size_t)k*nbPts+j];
- }
- } else {
- for (k=0 ; k<nbRHS ; k++) {
- d_sol[(size_t)k*nbPts+i] += creal(Aij) * d_rhs[(size_t)k*nbPts+j] ;
- }
+ ierr=computeKernelBEM(&(coord_i[0]), &(coord_j[0]), i==j, &Aij) ;
+ for (k=0 ; k<nbRHS ; k++) {
+ d_sol[(size_t)k*nbPts+i] += Aij * d_rhs[(size_t)k*nbPts+j] ;
}
} /* if testDofNeighbour(.. */
} /* for (i=0 */
@@ -61,7 +54,7 @@ int produitClassique(void **sol) {
int ierr;
double temps_initial, temps_final, temps_cpu;
- void *solCLA = MpfCalloc((size_t)nbPts*nbRHS, complexALGO ? sizeof(double _Complex) : sizeof(double)) ; CHKPTRQ(solCLA) ;
+ void *solCLA = MpfCalloc((size_t)nbPts*nbRHS, sizeof(double)) ; CHKPTRQ(solCLA) ;
temps_initial = getTime ();
ierr = produitClassiqueBEM(2, solCLA) ; CHKERRQ(ierr) ;
diff --git a/src/compare.c b/src/compare.c
index fa777b7..9eae7dd 100644
--- a/src/compare.c
+++ b/src/compare.c
@@ -19,7 +19,6 @@ frobenius_update( double *scale, double *sumsq, const double value )
int computeRelativeError(void *sol, void *ref, double *eps) {
double normRef[2], normDelta[2];
- double _Complex *z_sol=sol, *z_ref=ref ;
double *d_sol=sol, *d_ref=ref ;
*eps=0.;
@@ -28,16 +27,8 @@ int computeRelativeError(void *sol, void *ref, double *eps) {
normRef[0] = 1.;
normRef[1] = 0.;
for (size_t i=0 ; i<(size_t)nbPts*nbRHS ; i++) {
- if (complexALGO) {
- frobenius_update( normDelta, normDelta+1, (creal(z_sol[i]) - creal(z_ref[i])) );
- frobenius_update( normDelta, normDelta+1, (cimag(z_sol[i]) - cimag(z_ref[i])) );
-
- frobenius_update( normRef, normRef+1, creal(z_ref[i]) );
- frobenius_update( normRef, normRef+1, cimag(z_ref[i]) );
- } else {
- frobenius_update( normDelta, normDelta+1, (d_sol[i] - d_ref[i]) );
- frobenius_update( normRef, normRef+1, d_ref[i] );
- }
+ frobenius_update( normDelta, normDelta+1, (d_sol[i] - d_ref[i]) );
+ frobenius_update( normRef, normRef+1, d_ref[i] );
}
if ( normDelta[1] == 0. ) {/* To handle correctly the case where both vectors are null */
*eps = 0.;
@@ -52,17 +43,12 @@ int computeRelativeError(void *sol, void *ref, double *eps) {
int computeVecNorm(void *ref, double *norm) {
double normRef ;
- double _Complex *z_ref=ref ;
double *d_ref=ref ;
*norm=0.;
normRef=0. ;
for (size_t i=0 ; i<(size_t)nbPts*nbRHS ; i++) {
- if (complexALGO) {
- normRef += z_ref[i]*z_ref[i];
- } else {
- normRef += d_ref[i]*d_ref[i] ;
- }
+ normRef += d_ref[i]*d_ref[i] ;
}
*norm = sqrt(normRef) ;
diff --git a/src/help.c b/src/help.c
index 67991ff..9e3b356 100644
--- a/src/help.c
+++ b/src/help.c
@@ -8,11 +8,8 @@ int printHelp() {
printf("\ntest_FEMBEM usage:\n-----------------\n\n"
" -nbpts: Number of unknowns of the global linear system. Default: 16000 (can be in floating point notation, like 1.2e4).\n"
" -nbrhs: Number of right-hand sides of the global linear system. Default: 1.\n"
- " -s/-d/-c/-z: Scalar type of the computation (real computation uses 1/R kernel in BEM, complex computation uses oscillatory exp(ikr)/r kernel). Default: z.\n"
" -lambda: Wavelength (for oscillatory kernels). Mutually exclusive with -ptsperlambda. Default lambda is set with 10 points per wavelength.\n"
" -ptsperlambda: Points per wavelength on the main cylinder (for oscillatory kernels). Default: 10.\n"
- " --write-mesh: Write a VTK file 'mesh.vtu' of the mesh. Default is OFF.\n"
- " --write-mesh-unv: Write a UNV file 'mesh.unv' of the mesh. Default is OFF.\n"
"\n"
" BEM tests use a surfacic cylinder with an optionnal cylinder detail, defined with:\n"
" -radius: Radius of the cylinder. Default is 2.\n"
diff --git a/src/hmat.c b/src/hmat.c
index 4095dcd..fc63912 100644
--- a/src/hmat.c
+++ b/src/hmat.c
@@ -19,40 +19,6 @@ int init_hmat_interface() {
return 0;
}
-/* Simple assembly function, as needed by hmat */
-void interactionKernel(void* user_context, int i, int j, void* result) {
- int ierr;
- double coord_i[3], coord_j[3];
- double _Complex Aij = 0;
- contextTestFEMBEM *myCtx = (contextTestFEMBEM *)user_context;
-
- // Update indices (if we compute only a sub-part of the whole matrix)
- i += myCtx->rowOffset;
- j += myCtx->colOffset;
-
- // BEM part
- if (i<nbPts && j<nbPts) { // if testedModel is not _bem or _fembem, nbPtsBEM is 0
- ierr=computeCoordCylinder(i, coord_i); CHKERRA(ierr);
- ierr=computeCoordCylinder(j, coord_j); CHKERRA(ierr);
- ierr=computeKernelBEM(&(coord_i[0]), &(coord_j[0]), i==j, &Aij); CHKERRA(ierr);
- }
-
- switch(stype) {
- case SIMPLE_PRECISION:
- case DOUBLE_PRECISION:
- *((double*)result) = Aij;
- break;
- case SIMPLE_COMPLEX:
- case DOUBLE_COMPLEX:
- *((double _Complex*)result) = Aij;
- break;
- default:
- SETERRA(1, "Unknown stype");
- }
-
- return;
-}
-
/**
This structure is a glue between hmat library and application.
*/
@@ -182,7 +148,7 @@ advanced_compute_hmat(struct hmat_block_compute_context_t *b) {
int i, j, ierr;
double *dValues = (double *) b->block;
block_data_t* bdata = (block_data_t*) b->user_data;
- int step = complexALGO ? 2 : 1;
+ int step = 1;
double (*row_point)[3];
double (*col_point)[3];
@@ -198,26 +164,14 @@ advanced_compute_hmat(struct hmat_block_compute_context_t *b) {
for (i = 0; i < b->row_count; ++i, dValues+=step, row_point++) { // Idem with dof_row_used[]
int row = bdata->row_hmat2client ? bdata->row_hmat2client[i + b->row_start + bdata->row_start] : i + b->row_start + bdata->row_start;
int row_g = row + myCtx->rowOffset;
- double _Complex Aij = 0;
+ double Aij = 0;
// BEM part
if (row_g < nbPts && col_g < nbPts) {
ierr=computeKernelBEM((double*)row_point, (double*)col_point, row_g==col_g, &Aij); CHKERRA(ierr);
}
- switch(stype) {
- case SIMPLE_PRECISION:
- case DOUBLE_PRECISION:
- *((double*)dValues) = Aij;
- break;
- case SIMPLE_COMPLEX:
- case DOUBLE_COMPLEX:
- *(double _Complex*)dValues = Aij;
- break;
- default:
- SETERRA(1, "Unknown stype");
- }
-
+ *((double*)dValues) = Aij;
} /* for (i = 0; ... */
// Advance to the start of the next column
@@ -277,41 +231,6 @@ void computeDenseBlockFEMBEM(int *LigInf, int *LigSup, int *ColInf, int *ColSup,
free_block_data(block_info.user_data);
}
-/*! \brief Computes a block for the matrix of interactions.
-
- This is exactly like computeDenseBlockFEMBEM() but in simple complex or simple precision arithmetics.
-*/
-void computeDenseBlockFEMBEM_Cprec(int *LigInf, int *LigSup, int *ColInf, int *ColSup,
- int *___Inf, int *___Sup, int *iloc, int *jloc,
- int *tailleS, int *tailleL, void *bloc, void *context) {
- int id, ic, il, inxBloc ;
- ASSERTA(context);
-
- void * zbloc = MpfCalloc((size_t)(*tailleS)*(*___Sup - *___Inf), sizeof(double)*(complexALGO?2:1)); CHKPTRA(zbloc) ;
-
- /* Appelle la routine en double precision */
- computeDenseBlockFEMBEM( LigInf, LigSup, ColInf, ColSup,
- ___Inf, ___Sup, iloc, jloc,
- tailleS, tailleL, zbloc, context) ;
-
- for (id = 0 ; id < *___Sup - *___Inf ; id++) {
- for (ic = ColInf[0]-1; ic < ColSup[0] ; ic++) {
- for (il = LigInf[0]-1; il < LigSup[0] ; il++) {
- inxBloc = (*tailleL)*(ic-ColInf[0] +1) + il - LigInf[0] + 1 + (*iloc - 1) + (*jloc - 1) * (*tailleL) +
- (*tailleS)*id;
- if (stype == SIMPLE_COMPLEX) {
- ((float _Complex*)bloc)[inxBloc] = (float)(((double _Complex*)zbloc)[inxBloc]) ;
- } else {
- ((float*)bloc)[inxBloc] = (float)(((double*)zbloc)[inxBloc]) ;
- }
- }
- }
- }
-
- MpfFree(zbloc);
-}
-
-#ifdef HAVE_HMAT
HMAT_desc_t *HMAT_generate_matrix( hmat_interface_t *hi ) {
HMAT_desc_t *hdesc = MpfCalloc( 1, sizeof(HMAT_desc_t) );
@@ -393,4 +312,3 @@ void HMAT_destroy_matrix( hmat_interface_t *hi,
MpfFree(hdesc);
}
}
-#endif
diff --git a/src/kernel.c b/src/kernel.c
index 70ab6f2..24f57bb 100644
--- a/src/kernel.c
+++ b/src/kernel.c
@@ -1,7 +1,7 @@
#include "main.h"
extern double meshStep;
-int computeKernelBEM(double *coord1, double *coord2, int self, double complex *kernel) {
- double r, v[3], k=0. ;
+int computeKernelBEM(double *coord1, double *coord2, int self, double *kernel) {
+ double r, v[3];
v[0]=coord1[0]-coord2[0] ;
v[1]=coord1[1]-coord2[1] ;
@@ -9,11 +9,7 @@ int computeKernelBEM(double *coord1, double *coord2, int self, double complex *k
r=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]) ;
if (r<meshStep/2.) r=meshStep/2.;
- if (complexALGO) {
- k=2*M_PI/lambda;
- *kernel = cexp(I * k * r) / (4 * M_PI * r);
- } else
- *kernel = 1 / (4 * M_PI * r);
+ *kernel = 1 / (4 * M_PI * r);
return 0;
}
diff --git a/src/main.c b/src/main.c
index 3158032..836b3fd 100644
--- a/src/main.c
+++ b/src/main.c
@@ -15,9 +15,6 @@ int nbPtsPerLeaf = 50;
int nbAlgoRuns = 1;
char *sTYPE_ELEMENT_ARRAY_test_fembem[] = {NULL};
int coupled = 0;
-int simplePrec = 0;
-int complexALGO = 1;
-ScalarType stype = DOUBLE_COMPLEX;
double lambda;
/*! \brief Main routine
@@ -63,42 +60,13 @@ int main(int argc, char **argv) {
lambda = 10.*step ;
printf(" Setting lambda = %f (with 10 points per wavelength)\n", lambda) ;
- /* --- Choix de l'arithmetique de calcul (default is '-z') --- */
- if (MpfArgHasName(&argc, argv, 1, "-s") > 0) {
- simplePrec = 1 ;
- stype=SIMPLE_PRECISION ; /* For the fmm, it will be changed below depending on the FMM algo selected */
- complexALGO=0;
- printf("Simple Precision\n");
- }
- if (MpfArgHasName(&argc, argv, 1, "-d") > 0) {
- simplePrec = 0 ;
- stype=DOUBLE_PRECISION ; /* For the fmm, it will be changed below depending on the FMM algo selected */
- complexALGO=0;
- printf("Double Precision\n");
- }
- if (MpfArgHasName(&argc, argv, 1, "-c") > 0) {
- simplePrec = 1 ;
- stype=SIMPLE_COMPLEX ; /* For the fmm, it will be changed below depending on the FMM algo selected */
- complexALGO=1;
- printf("Simple Complex\n");
- }
- if (MpfArgHasName(&argc, argv, 1, "-z") > 0) {
- simplePrec = 0 ;
- stype=DOUBLE_COMPLEX ; /* For the fmm, it will be changed below depending on the FMM algo selected */
- complexALGO=1;
- printf("Double Complex\n");
- }
-
/* Setting remaining variables */
sparseRHS=(double)nbPts/80./log10((double)nbPts) ;
if (sparseRHS<1) sparseRHS = 1;
printf("Setting sparseRHS = %d\n", sparseRHS) ;
printf("<PERFTESTS> TEST_FEMBEM_Version = %s\n" , PACKAGE_VERSION);
-#ifdef HAVE_HMAT
printf("<PERFTESTS> HMAT_Version = %s\n" , hmat_get_version() );
-#endif
- printf("<PERFTESTS> ScalarType = %s\n" , MPF_scalName[stype]);
/* Prepare the mesh and RHS used for the FMM/HMAT tests */
ierr = prepareTEST() ; CHKERRQ(ierr) ;
diff --git a/src/prepareTEST.c b/src/prepareTEST.c
index 9a39ecd..175453f 100644
--- a/src/prepareTEST.c
+++ b/src/prepareTEST.c
@@ -17,7 +17,7 @@ int prepareTEST(void) {
/* Cree le second membre du produit mat-vec */
printf("Computing RHS...\n") ;
- rhs = MpfCalloc((size_t)nbPts*nbRHS, complexALGO ? sizeof(double _Complex) : sizeof(double)) ; CHKPTRQ(rhs) ;
+ rhs = MpfCalloc((size_t)nbPts*nbRHS, sizeof(double)) ; CHKPTRQ(rhs) ;
ierr = computeRhs() ; CHKERRQ(ierr) ;
leave_context();
diff --git a/src/rhs.c b/src/rhs.c
index b70d337..5e54a22 100644
--- a/src/rhs.c
+++ b/src/rhs.c
@@ -3,83 +3,25 @@
extern double height;
int computeRhs(void) {
- double coord[3], OP[3] ;
+ double coord[3];
int i, j, ierr ;
- double complex *z_rhs=rhs ;
double *d_rhs=rhs ;
- double small = simplePrec ? 1.e-20 : 1.e-100;
- double k=2.*M_PI/lambda;
+ double small = 1.e-100;
for (i=0 ; i<nbPts ; i++)
if (i % sparseRHS == 0) {
// Insert non-zero values on 1 value every 'sparseRHS'
ierr = computeCoordCylinder(i, &(coord[0])) ; CHKERRQ(ierr);
for (j=0 ; j<nbRHS ; j++) {
- if (complexALGO) {
- // OP = fake incident plane wave coming from direction rotating around the object
- ierr = computeCoordCylinder( (size_t)j*nbPts/nbRHS%nbPts , &(OP[0])) ; CHKERRQ(ierr);
- OP[2] -= height/2.; // to have plane waves coming from below too
- z_rhs[(size_t)j * nbPts + i] = cexp(
- I * k * (coord[0] * OP[0] + coord[1] * OP[1] + coord[2] * OP[2]));
- } else {
- d_rhs[(size_t)j*nbPts+i]=cos(coord[0]/(j+1)+coord[1]+sin(j)*coord[2])+sin(j) ;
- }
+ d_rhs[(size_t)j*nbPts+i]=cos(coord[0]/(j+1)+coord[1]+sin(j)*coord[2])+sin(j) ;
}
} else {
// Insert zero everywhere else (in fact we put a 'very small' value
// instead of zero, otherwise run-time optimisations produce unrealistic timers with FMM)
for (j=0 ; j<nbRHS ; j++) {
- if (complexALGO) {
- z_rhs[(size_t)j * nbPts + i] = small + I * small;
- } else {
- d_rhs[(size_t)j*nbPts+i]=small ;
- }
+ d_rhs[(size_t)j*nbPts+i]=small ;
}
}
return 0 ;
}
-
-void computeDenseBlockData (int *LigInf, int *LigSup, int *ColInf, int *ColSup, int *___Inf, int *___Sup, int *iloc, int *jloc, int *tailleS, int *tailleL, void *bloc, void *context) {
- int tailleL__ ;
- int iloc_____ ;
- int jloc_____ ;
- int il,ic,inxBloc;
- size_t inxGlob ;
- float *s_bloc=bloc ;
- double *d_bloc=bloc ;
- float _Complex *c_bloc=bloc ;
- double _Complex *z_bloc=bloc ;
- ASSERTA(context);
- contextTestFEMBEM *myCtx = (contextTestFEMBEM *)context;
- double _Complex *z_data=(double _Complex *)myCtx->dataVec ;
- double *d_data=(double *)myCtx->dataVec ;
-
- tailleL__ = *tailleL ;
- iloc_____ = *iloc - 1;
- jloc_____ = *jloc - 1;
-
- for (ic = ColInf[0]-1; ic < ColSup[0] ; ic++) {
- for (il = LigInf[0]-1; il < LigSup[0] ; il++) {
- inxBloc = tailleL__*(ic-ColInf[0]+1) + il-LigInf[0]+1 + iloc_____ + jloc_____ * tailleL__ ;
- inxGlob = (size_t)(ic + myCtx->colOffset) * nbPts + (il + myCtx->rowOffset);
-
- switch (stype) {
- case SIMPLE_PRECISION:
- s_bloc[inxBloc] = (float)d_data[inxGlob];
- break;
- case DOUBLE_PRECISION:
- d_bloc[inxBloc] = d_data[inxGlob];
- break;
- case SIMPLE_COMPLEX:
- c_bloc[inxBloc] = (float _Complex)z_data[inxGlob];
- break;
- case DOUBLE_COMPLEX:
- z_bloc[inxBloc] = z_data[inxGlob];
- break;
- default:
- break;
- } /* switch (stype) ... */
- }
- }
-}
diff --git a/src/testHMAT.c b/src/testHMAT.c
index 2ce85c2..2b535da 100644
--- a/src/testHMAT.c
+++ b/src/testHMAT.c
@@ -7,7 +7,6 @@ int testHMAT(double * relative_error) {
enter_context("testHMAT()");
int ierr;
double temps_initial, temps_final, temps_cpu;
- size_t vectorSize = (size_t)nbPts*nbRHS*(complexALGO?2:1);
printf("<PERFTESTS> HAcaAccuracy = %.4e \n", mpf_hmat_settings.acaEpsilon);
printf("<PERFTESTS> HRecompressionAccuracy = %.4e \n", mpf_hmat_settings.epsilon);
@@ -27,13 +26,13 @@ int testHMAT(double * relative_error) {
printf("\n**** Creating HMAT...\n") ;
/* We use the hmat interface initialized in MPF (default is starpu, --hmat-seq or --hmat-toyrt to change it) */
- hmat_interface_t *hi = mpf_hmat_settings.interfaces[stype];
+ hmat_interface_t *hi = mpf_hmat_settings.interfaces[HMAT_DOUBLE_PRECISION];
// hmat_factorization_none mean the user did not chose a factorization so we
// must choose one for him
if(mpf_hmat_settings.factorization_type == hmat_factorization_none) {
// use LDLT for real matrices
- mpf_hmat_settings.factorization_type = complexALGO ? hmat_factorization_llt : hmat_factorization_ldlt;
+ mpf_hmat_settings.factorization_type = hmat_factorization_ldlt;
}
HMAT_desc_t *hdesc;
@@ -88,14 +87,8 @@ int testHMAT(double * relative_error) {
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) ;
diff --git a/src/util.c b/src/util.c
index f5af67d..c8825a4 100644
--- a/src/util.c
+++ b/src/util.c
@@ -463,39 +463,6 @@ double time_interval_in_seconds(Time t1, Time t2) {
return time_in_s(interval);
}
/* ================================================================================== */
-/* Convertit de simple vers double precision. Attention : Le tableau 'data' doit etre suffisamment grand */
-void simpleToDouble(void *data, size_t n) {
- ssize_t i ;
- double *dest ;
- float tmp, *src ;
-
- dest =(double*)(data) ;
- src =(float*)(data) ;
- /* Commence la conversion par la fin pour ne pas ecraser les donnees */
- for (i=n-1 ; i>=0 ; i--) {
- tmp=src[i] ;
- dest[i]=0. ;
- dest[i]=(double)tmp ;
- }
-
-}
-/* ================================================================================== */
-/* Convertit de double vers simple precision. */
-void doubleToSimple(void *data, size_t n) {
- size_t i ;
- double tmp, *src =(double*)data ;
- float *dest =(float*)data ;
-
- for (i=0 ; i<n ; i++) {
- tmp=src[i] ;
- dest[i]=0. ;
- dest[i]=(float)tmp ;
- }
-
-}
-
-#ifdef HAVE_HMAT
-/* ================================================================================== */
int hmat_get_sync_exec(void) {
return 1;
}
@@ -525,7 +492,6 @@ void mpf_hmat_compression(struct mpf_hmat_create_compression_args_t *args) {
default: SETWARN(1, "unknown compression method");
}
}
-#endif
/* ================================================================================== */
int SCAB_Init(int* argc, char*** argv) {
int ierr;
--
GitLab