diff --git a/include/util.h b/include/util.h
index b0e67004ec8f1111825f9862724890ab9fa7dccb..4c2e7497c4531a252b6ebfd5a52fe4709b998ff0 100644
--- a/include/util.h
+++ b/include/util.h
@@ -32,42 +32,6 @@
 
 typedef enum { MPF_FALSE, MPF_TRUE } Logical;
 
-/**************************************************************************************/
-/***   Scalar Types Definitions   *****************************************************/
-/**************************************************************************************/
-/*! \brief Simple real type.
-
-  This scalar type corresponds to the SIMPLE_PRECISION value in the ScalarType enumeration.
-*/
-typedef float S_type ;
-/**************************************************************************************/
-/*! \brief Double real type.
-
-  This scalar type corresponds to the DOUBLE_PRECISION value in the ScalarType enumeration.
-*/
-typedef double D_type ;
-/**************************************************************************************/
-/*! \brief Simple complex type.
-
-  This scalar type corresponds to the SIMPLE_COMPLEX value in the ScalarType enumeration.
-*/
-typedef struct {
-  /*! \brief real part */
-  float r ;
-  /*! \brief imaginary part */
-  float i ;
-} C_type;
-/**************************************************************************************/
-/*! \brief Double complex type.
-
-  This scalar type corresponds to the DOUBLE_COMPLEX value in the ScalarType enumeration.
-*/
-typedef struct {
-  /*! \brief real part */
-  double r ;
-  /*! \brief imaginary part */
-  double i ;
-} Z_type;
 /**************************************************************************************/
 /*! \brief All the scalar types
 
@@ -75,224 +39,24 @@ typedef struct {
   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, S_type in MPF) */
+  /*! \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, D_type in MPF) */
+  /*! \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, C_type in MPF) */
+  /*! \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, Z_type in MPF) */
+  /*! \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;
 
-/**************************************************************************************/
-/***   Sparse Scalar Types Definitions   **********************************************/
-/**************************************************************************************/
-/*! \brief Simple real sparse type.
-
-  This sparse scalar type corresponds to the position (row, column) of a non-zero element
-  followed by its value (S_type).
-*/
-typedef struct {
-  /*! \brief Row position */
-  int i ;
-  /*! \brief Column position */
-  int j ;
-  /*! \brief Value */
-  S_type v;
-} SnonZero;
-/**************************************************************************************/
-/*! \brief Double real sparse type.
-
-  This sparse scalar type corresponds to the position (row, column) of a non-zero element
-  followed by its value (D_type).
-*/
-typedef struct {
-  /*! \brief Row position */
-  int i ;
-  /*! \brief Column position */
-  int j ;
-  /*! \brief Value */
-  D_type v;
-} DnonZero;
-/**************************************************************************************/
-/*! \brief Simple complex sparse type.
-
-  This sparse scalar type corresponds to the position (row, column) of a non-zero element
-  followed by its value (C_type).
-*/
-typedef struct {
-  /*! \brief Row position */
-  int i ;
-  /*! \brief Column position */
-  int j ;
-  union {
-    float _Complex cv;
-    /*! \brief Value */
-    C_type v;
-  };
-} CnonZero;
-/**************************************************************************************/
-/*! \brief Double complex sparse type.
-
-  This sparse scalar type corresponds to the position (row, column) of a non-zero element
-  followed by its value (Z_type).
-*/
-typedef struct {
-  /*! \brief Row position */
-  int i ;
-  /*! \brief Column position */
-  int j ;
-  union {
-    double _Complex cv;
-    /*! \brief Value */
-    Z_type v;
-  };
-} ZnonZero;
-/**************************************************************************************/
-
 /*! \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
  ;
-
-/* ================================================================================== */
-/*! \brief  Multi-arithmetics constant.
-
-This structure is used to store any floating point constants in all of the four possible
-arithmetics. This is an array of 4 void* pointers, pointing on the constant in S, D, C and Z precision
-(in that order). For instance, if X is an MpfConst, then X[SIMPLE_COMPLEX] is a pointer on a C_type.
-*/
-typedef void* MpfConst [nbScalarType] ;
-/* ================================================================================== */
-/* +1 in Floating point */
-/*! \brief +1 in simple real precision. */
-_EXTERN_TEST_FEMBEM_ S_type Mpf_s_pone
-#ifdef   ___TEST_FEMBEM___
-=1.
-    #endif
-    ;
-/*! \brief +1 in double real precision. */
-_EXTERN_TEST_FEMBEM_ D_type Mpf_d_pone
-#ifdef   ___TEST_FEMBEM___
-=1.
-    #endif
-    ;
-/*! \brief +1 in simple complex precision. */
-_EXTERN_TEST_FEMBEM_ C_type Mpf_c_pone
-#ifdef   ___TEST_FEMBEM___
-={1., 0.}
- #endif
- ;
-/*! \brief +1 in double complex precision. */
-_EXTERN_TEST_FEMBEM_ Z_type Mpf_z_pone
-#ifdef   ___TEST_FEMBEM___
-={1., 0.}
- #endif
- ;
-/*! \brief +1 in all precisions. */
-_EXTERN_TEST_FEMBEM_ MpfConst Mpf_pone
-#ifdef   ___TEST_FEMBEM___
-={&Mpf_s_pone, &Mpf_d_pone, &Mpf_c_pone, &Mpf_z_pone}
- #endif
- ;
-/* ================================================================================== */
-/* 0 in Floating point */
-/*! \brief Zero in simple real precision. */
-_EXTERN_TEST_FEMBEM_ S_type Mpf_s_zero
-#ifdef   ___TEST_FEMBEM___
-=0.
-    #endif
-    ;
-/*! \brief Zero in double real precision. */
-_EXTERN_TEST_FEMBEM_ D_type Mpf_d_zero
-#ifdef   ___TEST_FEMBEM___
-=0.
-    #endif
-    ;
-/*! \brief Zero in simple complex precision. */
-_EXTERN_TEST_FEMBEM_ C_type Mpf_c_zero
-#ifdef   ___TEST_FEMBEM___
-={0., 0.}
- #endif
- ;
-/*! \brief Zero in double complex precision. */
-_EXTERN_TEST_FEMBEM_ Z_type Mpf_z_zero
-#ifdef   ___TEST_FEMBEM___
-={0., 0.}
- #endif
- ;
-/*! \brief Zero in all precisions. */
-_EXTERN_TEST_FEMBEM_ MpfConst Mpf_zero
-#ifdef   ___TEST_FEMBEM___
-={&Mpf_s_zero, &Mpf_d_zero, &Mpf_c_zero, &Mpf_z_zero}
- #endif
- ;
-/* ================================================================================== */
-/* -1 in Floating point */
-/*! \brief -1 in simple real precision. */
-_EXTERN_TEST_FEMBEM_ S_type Mpf_s_mone
-#ifdef   ___TEST_FEMBEM___
-=-1.
-#endif
-;
-/*! \brief -1 in double real precision. */
-_EXTERN_TEST_FEMBEM_ D_type Mpf_d_mone
-#ifdef   ___TEST_FEMBEM___
-=-1.
-#endif
-;
-/*! \brief -1 in simple complex precision. */
-_EXTERN_TEST_FEMBEM_ C_type Mpf_c_mone
-#ifdef   ___TEST_FEMBEM___
-={-1., 0.}
-#endif
-;
-/*! \brief -1 in double complex precision. */
-_EXTERN_TEST_FEMBEM_ Z_type Mpf_z_mone
-#ifdef   ___TEST_FEMBEM___
-={-1., 0.}
-#endif
-;
-/*! \brief -1 in all precisions. */
-_EXTERN_TEST_FEMBEM_ MpfConst Mpf_mone
-#ifdef   ___TEST_FEMBEM___
-={&Mpf_s_mone, &Mpf_d_mone, &Mpf_c_mone, &Mpf_z_mone}
-#endif
-;
-
-/* ================================================================================== */
-/*! \brief 0 in integer type. */
-_EXTERN_TEST_FEMBEM_ int Mpf_i_zero
-#ifdef   ___TEST_FEMBEM___
-=0
-#endif
-;
-/*! \brief 1 in integer type. */
-_EXTERN_TEST_FEMBEM_ int Mpf_i_one
-#ifdef   ___TEST_FEMBEM___
-=1
-#endif
-;
-
-/* ================================================================================== */
-/*! \brief Size (in bytes) of the different scalar types. */
-_EXTERN_TEST_FEMBEM_ size_t Mpf_stype_size[nbScalarType]
-#ifdef   ___TEST_FEMBEM___
-={sizeof(S_type), sizeof(D_type), sizeof(C_type), sizeof(Z_type)}
-#endif
-;
-/*! \brief Size (in bytes) of the different non-zero scalar types. */
-_EXTERN_TEST_FEMBEM_ size_t Mpf_nonzerostype_size[nbScalarType]
-#ifdef   ___TEST_FEMBEM___
-={sizeof(SnonZero), sizeof(DnonZero), sizeof(CnonZero), sizeof(ZnonZero)}
- #endif
- ;
-
 #ifdef HAVE_HMAT
 enum mpf_hmat_engine { mpf_hmat_seq, mpf_hmat_starpu, mpf_hmat_toyrt };
 struct mpf_hmat_settings_t {
@@ -411,15 +175,6 @@ int Mpf_progressBar(int currentValue, int maximumValue) ;
 int MpfArgGetDouble( int *Argc, char **argv, int rflag, const char *name, double *val ) ;
 int MpfArgHasName( int *Argc, char **argv, int rflag, const char *name );
 int MpfArgGetInt( int  *Argc, char **argv, int rflag, const char *name, int  *val );
-int S_sortAndCompact(SnonZero *mat, int *size);
-int sortAndCompact(DnonZero *mat, int *size);
-int C_sortAndCompact(CnonZero *mat, int *size);
-int Z_sortAndCompact(ZnonZero *mat, int *size);
-int Z_sortAndCompactL(ZnonZero *mat, size_t *size);
-int compare_SnonZero_Line(const void *pt1, const void *pt2);
-int compare_DnonZero_Line(const void *pt1, const void *pt2);
-int compare_CnonZero_Line(const void *pt1, const void *pt2);
-int compare_ZnonZero_Line(const void *pt1, const void *pt2);
 Time get_time(void);
 Time time_interval(Time t1, Time t2);
 Time add_times(Time t1, Time t2);
diff --git a/src/classic.c b/src/classic.c
index 5670b1508da5cbf259b0b9038c365c5a7c20a9b2..dd13ec214e8e864289af536742e69ff3b79f3815 100644
--- a/src/classic.c
+++ b/src/classic.c
@@ -23,9 +23,9 @@ 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];
-    Z_type Aij ;
-    Z_type *z_sol=sol, *z_rhs=rhs ;
-    D_type *d_sol=sol, *d_rhs=rhs ;
+    double _Complex Aij ;
+    double _Complex *z_sol=sol, *z_rhs=rhs ;
+    double *d_sol=sol, *d_rhs=rhs ;
 
     ierr=computeCoordCylinder(j, &(coord_j[0])) ; CHKERRA(ierr);
     /* Boucle sur les lignes */
@@ -35,12 +35,11 @@ int produitClassiqueBEM(int ffar, void *sol) {
         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].r += Aij.r * z_rhs[(size_t)k*nbPts+j].r - Aij.i * z_rhs[(size_t)k*nbPts+j].i ;
-            z_sol[(size_t)k*nbPts+i].i += Aij.r * z_rhs[(size_t)k*nbPts+j].i + Aij.i * z_rhs[(size_t)k*nbPts+j].r ;
+            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] += Aij.r * d_rhs[(size_t)k*nbPts+j] ;
+            d_sol[(size_t)k*nbPts+i] += creal(Aij) * d_rhs[(size_t)k*nbPts+j] ;
           }
         }
       } /* if testDofNeighbour(.. */
@@ -62,7 +61,7 @@ int produitClassique(void **sol) {
   int ierr;
   double temps_initial, temps_final, temps_cpu;
 
-  void *solCLA = MpfCalloc((size_t)nbPts*nbRHS, complexALGO ? sizeof(Z_type) : sizeof(D_type)) ; CHKPTRQ(solCLA) ;
+  void *solCLA = MpfCalloc((size_t)nbPts*nbRHS, complexALGO ? sizeof(double _Complex) : sizeof(double)) ; CHKPTRQ(solCLA) ;
 
   temps_initial = getTime ();
   ierr = produitClassiqueBEM(2, solCLA) ; CHKERRQ(ierr) ;
diff --git a/src/compare.c b/src/compare.c
index 39c0f23016b97e149c36b67bf48d731856df6bbf..fa777b78459920275eaf0960a990e98358da2f67 100644
--- a/src/compare.c
+++ b/src/compare.c
@@ -19,8 +19,8 @@ frobenius_update( double *scale, double *sumsq, const double value )
 
 int computeRelativeError(void *sol, void *ref, double *eps) {
   double normRef[2], normDelta[2];
-  Z_type *z_sol=sol, *z_ref=ref ;
-  D_type *d_sol=sol, *d_ref=ref ;
+  double _Complex *z_sol=sol, *z_ref=ref ;
+  double *d_sol=sol, *d_ref=ref ;
 
   *eps=0.;
   normDelta[0] = 1.;
@@ -29,11 +29,11 @@ int computeRelativeError(void *sol, void *ref, double *eps) {
   normRef[1]   = 0.;
   for (size_t i=0 ; i<(size_t)nbPts*nbRHS ; i++) {
     if (complexALGO) {
-      frobenius_update( normDelta, normDelta+1, (z_sol[i].r - z_ref[i].r) );
-      frobenius_update( normDelta, normDelta+1, (z_sol[i].i - z_ref[i].i) );
+      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, z_ref[i].r );
-      frobenius_update( normRef, normRef+1, z_ref[i].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] );
@@ -52,14 +52,14 @@ int computeRelativeError(void *sol, void *ref, double *eps) {
 
 int computeVecNorm(void *ref, double *norm) {
   double normRef ;
-  Z_type *z_ref=ref ;
-  D_type *d_ref=ref ;
+  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].r*z_ref[i].r+z_ref[i].i*z_ref[i].i ;
+      normRef += z_ref[i]*z_ref[i];
     } else {
       normRef += d_ref[i]*d_ref[i] ;
     }
diff --git a/src/hmat.c b/src/hmat.c
index 941365e2984d878d5969e0638a44360f275feeeb..4095dcdc0aa4179f0ba0187eee29e170bf18f1d2 100644
--- a/src/hmat.c
+++ b/src/hmat.c
@@ -300,8 +300,7 @@ void computeDenseBlockFEMBEM_Cprec(int *LigInf, int *LigSup, int *ColInf, int *C
         inxBloc = (*tailleL)*(ic-ColInf[0] +1)  + il - LigInf[0] + 1 + (*iloc   - 1) + (*jloc   - 1) * (*tailleL) +
             (*tailleS)*id;
         if (stype == SIMPLE_COMPLEX) {
-          ((C_type*)bloc)[inxBloc].r = (float)(((Z_type*)zbloc)[inxBloc].r) ;
-          ((C_type*)bloc)[inxBloc].i = (float)(((Z_type*)zbloc)[inxBloc].i) ;
+          ((float _Complex*)bloc)[inxBloc] = (float)(((double _Complex*)zbloc)[inxBloc]) ;
         } else {
           ((float*)bloc)[inxBloc]   = (float)(((double*)zbloc)[inxBloc]) ;
         }
diff --git a/src/prepareTEST.c b/src/prepareTEST.c
index c67d5f96d9a629eed9b66d71075aa8a84eca1119..9a39ecdaaaa2933a459984812cd6da3fef5f8532 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(Z_type) : sizeof(D_type)) ;  CHKPTRQ(rhs) ;
+  rhs = MpfCalloc((size_t)nbPts*nbRHS, complexALGO ? sizeof(double _Complex) : sizeof(double)) ;  CHKPTRQ(rhs) ;
   ierr = computeRhs() ; CHKERRQ(ierr) ;
 
   leave_context();
@@ -30,14 +30,14 @@ double getTime() {
 
 int displayArray(char *s, void *f, int m, int n) {
   int i, j ;
-  Z_type *z_f=f ;
-  D_type *d_f=f ;
+  //double _Complex *z_f=f ;
+  double *d_f=f ;
 
   for (j=0 ; j<n ; j++)
     for (i=0 ; i<m ; i++)
-      if (complexALGO)
-        printf("%s[%d, %d]=(%e, %e)\n", s, i, j, z_f[(size_t)j*m+i].r, z_f[(size_t)j*m+i].i) ;
-      else
+      //if (complexALGO)
+        //printf("%s[%d, %d]=(%e, %e)\n", s, i, j, z_f[(size_t)j*m+i].r, z_f[(size_t)j*m+i].i) ;
+      //else
         printf("%s[%d, %d]=%e\n", s, i, j, d_f[(size_t)j*m+i]) ;
   return 0 ;
 }
diff --git a/src/rhs.c b/src/rhs.c
index 27b0d1e503a79ce85620712eec1fe0d03768b881..b70d33709ef6697bc4884bc856e0511e7b931b71 100644
--- a/src/rhs.c
+++ b/src/rhs.c
@@ -6,7 +6,7 @@ int computeRhs(void) {
   double coord[3], OP[3] ;
   int i, j, ierr ;
   double complex *z_rhs=rhs ;
-  D_type *d_rhs=rhs ;
+  double *d_rhs=rhs ;
   double small = simplePrec ? 1.e-20 : 1.e-100;
   double k=2.*M_PI/lambda;
 
@@ -46,14 +46,14 @@ void computeDenseBlockData (int *LigInf, int *LigSup, int *ColInf, int *ColSup,
   int jloc_____  ;
   int il,ic,inxBloc;
   size_t inxGlob ;
-  S_type *s_bloc=bloc ;
-  D_type *d_bloc=bloc ;
-  C_type *c_bloc=bloc ;
-  Z_type *z_bloc=bloc ;
+  float *s_bloc=bloc ;
+  double *d_bloc=bloc ;
+  float _Complex *c_bloc=bloc ;
+  double _Complex *z_bloc=bloc ;
   ASSERTA(context);
   contextTestFEMBEM *myCtx = (contextTestFEMBEM *)context;
-  Z_type *z_data=(Z_type *)myCtx->dataVec ;
-  D_type *d_data=(D_type *)myCtx->dataVec ;
+  double _Complex *z_data=(double _Complex *)myCtx->dataVec ;
+  double *d_data=(double *)myCtx->dataVec ;
 
   tailleL__ = *tailleL   ;
   iloc_____ = *iloc   - 1;
@@ -66,14 +66,13 @@ void computeDenseBlockData (int *LigInf, int *LigSup, int *ColInf, int *ColSup,
 
       switch (stype) {
         case SIMPLE_PRECISION:
-          s_bloc[inxBloc] = (S_type)d_data[inxGlob];
+          s_bloc[inxBloc] = (float)d_data[inxGlob];
           break;
         case DOUBLE_PRECISION:
           d_bloc[inxBloc] = d_data[inxGlob];
           break;
         case SIMPLE_COMPLEX:
-          c_bloc[inxBloc].r = (S_type)z_data[inxGlob].r;
-          c_bloc[inxBloc].i = (S_type)z_data[inxGlob].i;
+          c_bloc[inxBloc] = (float _Complex)z_data[inxGlob];
           break;
         case DOUBLE_COMPLEX:
           z_bloc[inxBloc] = z_data[inxGlob];
diff --git a/src/testHMAT.c b/src/testHMAT.c
index b61febed96f66897373bd621495cb89694b594fe..2ce85c24dfcba3957688ca06700d9190664f7f8b 100644
--- a/src/testHMAT.c
+++ b/src/testHMAT.c
@@ -48,7 +48,7 @@ int testHMAT(double * relative_error) {
     hmat_info_t info;
     hi->get_info(hmatrix, &info);
     printf("Compressed size: %ld\n", info.compressed_size);
-    printf("<PERFTESTS> AssembledSizeMb = %f \n", (double)info.compressed_size*Mpf_stype_size[stype]/1024./1024.) ;
+    printf("<PERFTESTS> AssembledSizeMb = %f \n", (double)info.compressed_size*sizeof(double)/1024./1024.) ;
     printf("Uncompressed size: %ld\n", info.uncompressed_size);
     printf("Ratio: %g\n", (double)info.compressed_size/info.uncompressed_size);
   }
@@ -77,7 +77,7 @@ int testHMAT(double * relative_error) {
     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("<PERFTESTS> FactorizedSizeMb = %f \n", (double)info.compressed_size*sizeof(double)/1024./1024.) ;
     printf("Uncompressed size: %ld\n", info.uncompressed_size);
     printf("Ratio: %g\n", (double)info.compressed_size/info.uncompressed_size);
   }
diff --git a/src/util.c b/src/util.c
index ee00d11b21d7832588d022ce86589f95c9cec607..f5af67d141139647d2b6921e1c93a27cf4d14a67 100644
--- a/src/util.c
+++ b/src/util.c
@@ -391,241 +391,6 @@ int MpfArgGetString( int *Argc, char **argv, int rflag, const char *name, char *
   return 1;
 }
 /* ================================================================================== */
-/*! \brief Sort the array mat, remove the zeros and sum the elements placed at the same position
-
-  For simple precision only (S)
-\param mat the sparse matrix to sort
-\param size the size (modified on output)
-\return 0 for success
-*/
-int S_sortAndCompact(SnonZero *mat, int *size) {
-  int ie, il ;
-
-  /* Trie le tableau */
-  qsort(mat, *size, sizeof(SnonZero), compare_SnonZero_Line) ;
-
-  /* somme les valeurs situees a la meme position et enleve les zeros */
-  ie=-1 ; /* ie = position d'ecriture */
-  for (il=0 ; il<*size ; il++) /* il = position de lecture */
-    if (mat[il].v) {
-      if (ie>=0 && compare_SnonZero_Line( &(mat[ie]), &(mat[il]) ) == 0) /* s'ils sont a la meme position, on somme */
-        mat[ie].v += mat[il].v ;
-      else { /* sinon, on garde l'element 'il' en le recopiant dans une nouvelle case */
-        ie ++ ;
-        if (il>ie) mat[ie] = mat[il] ;
-      }
-    }
-  *size=ie+1 ;
-  return 0 ;
-}
-/* ================================================================================== */
-/*! \brief Sort the array mat, remove the zeros and sum the elements placed at the same position
-
-  For double precision only (D)
-\param mat the sparse matrix to sort
-\param size the size (modified on output)
-\return 0 for success
-*/
-int sortAndCompact(DnonZero *mat, int *size) {
-  int ie, il ;
-
-  /* Trie le tableau */
-  qsort(mat, *size, sizeof(DnonZero), compare_DnonZero_Line) ;
-
-  /* somme les valeurs situees a la meme position et enleve les zeros */
-  ie=-1 ; /* ie = position d'ecriture */
-  for (il=0 ; il<*size ; il++) /* il = position de lecture */
-    if (mat[il].v) {
-      if (ie>=0 && compare_DnonZero_Line( &(mat[ie]), &(mat[il]) ) == 0) /* s'ils sont a la meme position, on somme */
-        mat[ie].v += mat[il].v ;
-      else { /* sinon, on garde l'element 'il' en le recopiant dans une nouvelle case */
-        ie ++ ;
-        if (il>ie) mat[ie] = mat[il] ;
-      }
-    }
-  *size=ie+1 ;
-  return 0 ;
-}
-/* ================================================================================== */
-/*! \brief Sort the array mat, remove the zeros and sum the elements placed at the same position
-
-  For simple complex only (C)
-\param mat the sparse matrix to sort
-\param size the size (modified on output)
-\return 0 for success
-*/
-int C_sortAndCompact(CnonZero *mat, int *size) {
-  int ie, il ;
-
-  /* Trie le tableau */
-  qsort(mat, *size, sizeof(CnonZero), compare_CnonZero_Line) ;
-
-  /* somme les valeurs situees a la meme position et enleve les zeros */
-  ie=-1 ; /* ie = position d'ecriture */
-  for (il=0 ; il<*size ; il++) /* il = position de lecture */
-    if (mat[il].v.r || mat[il].v.i) {
-      if (ie>=0 && compare_CnonZero_Line( &(mat[ie]), &(mat[il]) ) == 0) { /* s'ils sont a la meme position, on somme */
-        mat[ie].v.r += mat[il].v.r ;
-        mat[ie].v.i += mat[il].v.i ;
-      } else { /* sinon, on garde l'element 'il' en le recopiant dans une nouvelle case */
-        ie ++ ;
-        if (il>ie) mat[ie] = mat[il] ;
-      }
-    }
-  *size=ie+1 ;
-  return 0 ;
-}
-/* ================================================================================== */
-/*! \brief Sort the array mat, remove the zeros and sum the elements placed at the same position
-
-  For double complex only (Z)
-\param mat the sparse matrix to sort
-\param size the size (modified on output)
-\return 0 for success
-*/
-int Z_sortAndCompact(ZnonZero *mat, int *size) {
-  int ie, il ;
-
-  /* Trie le tableau */
-  qsort(mat, *size, sizeof(ZnonZero), compare_ZnonZero_Line) ;
-
-  /* somme les valeurs situees a la meme position et enleve les zeros */
-  ie=-1 ; /* ie = position d'ecriture */
-  for (il=0 ; il<*size ; il++) /* il = position de lecture */
-    if (mat[il].v.r || mat[il].v.i) {
-      if (ie>=0 && compare_ZnonZero_Line( &(mat[ie]), &(mat[il]) ) == 0) { /* s'ils sont a la meme position, on somme */
-        mat[ie].v.r += mat[il].v.r ;
-        mat[ie].v.i += mat[il].v.i ;
-      } else { /* sinon, on garde l'element 'il' en le recopiant dans une nouvelle case */
-        ie ++ ;
-        if (il>ie) mat[ie] = mat[il] ;
-      }
-    }
-  *size=ie+1 ;
-  return 0 ;
-}
-/* ================================================================================== */
-/*! \brief Sort the array mat, remove the zeros and sum the elements placed at the same position
-
-  For double complex only (Z) with size_t indices
-\param mat the sparse matrix to sort
-\param size the size (modified on output)
-\return 0 for success
-*/
-int Z_sortAndCompactL(ZnonZero *mat, size_t *size) {
-  ssize_t ie, il ;
-
-  /* Trie le tableau */
-  qsort(mat, *size, sizeof(ZnonZero), compare_ZnonZero_Line) ;
-
-  /* somme les valeurs situees a la meme position et enleve les zeros */
-  ie=-1 ; /* ie = position d'ecriture */
-  for (il=0 ; il<*size ; il++) /* il = position de lecture */
-    if (mat[il].v.r || mat[il].v.i) {
-      if (ie>=0 && compare_ZnonZero_Line( &(mat[ie]), &(mat[il]) ) == 0) { /* s'ils sont a la meme position, on somme */
-        mat[ie].v.r += mat[il].v.r ;
-        mat[ie].v.i += mat[il].v.i ;
-      } else { /* sinon, on garde l'element 'il' en le recopiant dans une nouvelle case */
-        ie ++ ;
-        if (il>ie) mat[ie] = mat[il] ;
-      }
-    }
-  *size=ie+1 ;
-  return 0 ;
-}
-/* ================================================================================== */
-/*! \brief compares two SnonZeros by lines
-
-\param pt1 pointer on the first SnonZero
-\param pt2 pointer on the second SnonZero
-\return -1, 0 or +1 if \a pt1 is before, equal or after \a pt2.
-*/
-int compare_SnonZero_Line(const void *pt1, const void *pt2) {
-  SnonZero *s1=(SnonZero*)pt1 ;
-  SnonZero *s2=(SnonZero*)pt2 ;
-
-  if (s1->i > s2->i)
-    return(+1) ;
-  if (s1->i < s2->i)
-    return(-1) ;
-
-  if (s1->j > s2->j)
-    return(+1) ;
-  if (s1->j < s2->j)
-    return(-1) ;
-
-  return(0) ;
-}
-/* ================================================================================== */
-/*! \brief compares two DnonZeros by lines
-
-\param pt1 pointer on the first DnonZero
-\param pt2 pointer on the second DnonZero
-\return -1, 0 or +1 if \a pt1 is before, equal or after \a pt2.
-*/
-int compare_DnonZero_Line(const void *pt1, const void *pt2) {
-  DnonZero *d1=(DnonZero*)pt1 ;
-  DnonZero *d2=(DnonZero*)pt2 ;
-
-  if (d1->i > d2->i)
-    return(+1) ;
-  if (d1->i < d2->i)
-    return(-1) ;
-
-  if (d1->j > d2->j)
-    return(+1) ;
-  if (d1->j < d2->j)
-    return(-1) ;
-
-  return(0) ;
-}
-/* ================================================================================== */
-/*! \brief compares two CnonZeros by lines
-
-\param pt1 pointer on the first CnonZero
-\param pt2 pointer on the second CnonZero
-\return -1, 0 or +1 if \a pt1 is before, equal or after \a pt2.
-*/
-int compare_CnonZero_Line(const void *pt1, const void *pt2) {
-  CnonZero *c1=(CnonZero*)pt1 ;
-  CnonZero *c2=(CnonZero*)pt2 ;
-
-  if (c1->i > c2->i)
-    return(+1) ;
-  if (c1->i < c2->i)
-    return(-1) ;
-
-  if (c1->j > c2->j)
-    return(+1) ;
-  if (c1->j < c2->j)
-    return(-1) ;
-
-  return(0) ;
-}
-/* ================================================================================== */
-/*! \brief compares two ZnonZeros by lines
-
-\param pt1 pointer on the first ZnonZero
-\param pt2 pointer on the second ZnonZero
-\return -1, 0 or +1 if \a pt1 is before, equal or after \a pt2.
-*/
-int compare_ZnonZero_Line(const void *pt1, const void *pt2) {
-  ZnonZero *z1=(ZnonZero*)pt1 ;
-  ZnonZero *z2=(ZnonZero*)pt2 ;
-
-  if (z1->i > z2->i)
-    return(+1) ;
-  if (z1->i < z2->i)
-    return(-1) ;
-
-  if (z1->j > z2->j)
-    return(+1) ;
-  if (z1->j < z2->j)
-    return(-1) ;
-
-  return(0) ;
-}
-/* ================================================================================== */
 Time get_time() {
   Time result ;
 #ifdef _WIN32