DON'T MERGE YET - Code cleaning and extension to general SIMD vectorization

.gitignore

@@ -1,4 +1,3 @@
-*
 !.gitignore
 
 

CREDITS

@@ -1,6 +1,6 @@
 This software is an outcome of the PhD thesis of Matthias Rottmann, University of Wuppertal.
 
-Code Designers: Matthias Rottmann, Artur Strebel, Simon Heybrock, Simone Bacchio, Bjoern Leder.
+Code Designers: Matthias Rottmann, Simone Bacchio, Artur Strebel, Simon Heybrock, Bjoern Leder.
 
 Contributions by: Andreas Frommer, Karsten Kahl, Stefan Krieg, Kalman Szabo, Wolfgang Soeldner,
                   Holger Arndt, Peter Georg.

Makefile

@@ -2,8 +2,8 @@
 CC = mpiicc 
 
 # --- CFLAGS -----------------------------------------
-CFLAGS_gnu = -std=gnu99 -Wall -pedantic -fopenmp -O3 -ffast-math -msse4.2
-CFLAGS_intel = -std=gnu99 -Wall -pedantic -qopenmp -O3  -xHOST
+CFLAGS_gnu = -std=gnu99 -Wall -pedantic -O3 -ffast-math -msse4.2 -fopenmp 
+CFLAGS_intel = -std=gnu99 -Wall -pedantic -O3  -xHOST -qopenmp 
 CFLAGS = $(CFLAGS_intel)
 
 # --- DO NOT CHANGE -----------------------------------
@@ -45,8 +45,8 @@ LIMELIB= -L$(LIMEDIR)/lib -llime
 # -DPARAMOUTPUT -DTRACK_RES -DFGMRES_RESTEST -DPROFILING
 # -DSINGLE_ALLREDUCE_ARNOLDI
 # -DCOARSE_RES -DSCHWARZ_RES -DTESTVECTOR_ANALYSIS -DDEBUG
-OPT_VERSION_FLAGS = $(CFLAGS) $(LIMEFLAGS) $(H5FLAGS) -DOPENMP -DSSE -DPARAMOUTPUT -DTRACK_RES
-DEVEL_VERSION_FLAGS = $(CFLAGS) $(LIMEFLAGS) -DOPENMP -DSSE -DDEBUG -DPARAMOUTPUT -DTRACK_RES -DFGMRES_RESTEST -DPROFILING -DCOARSE_RES -DSCHWARZ_RES -DTESTVECTOR_ANALYSIS
+OPT_VERSION_FLAGS = $(CFLAGS) $(LIMEFLAGS) $(H5FLAGS) -DPARAMOUTPUT -DTRACK_RES -DSSE -DOPENMP
+DEVEL_VERSION_FLAGS = $(CFLAGS) $(LIMEFLAGS) -DDEBUG -DPARAMOUTPUT -DTRACK_RES -DFGMRES_RESTEST -DPROFILING -DCOARSE_RES -DSCHWARZ_RES -DTESTVECTOR_ANALYSIS -DSSE -DOPENMP
 
 
 all: execs library exec-tests
@@ -92,10 +92,10 @@ $(INCDIR)/%: $(SRCDIR)/%
 	cp $(SRCDIR)/`basename $@` $@
 
 $(BUILDDIR)/%.o: $(GSRCDIR)/%.c $(SRCDIR)/*.h
-	$(CC) $(CFLAGS) $(OPT_VERSION_FLAGS) -c $< -o $@
+	$(CC) $(OPT_VERSION_FLAGS) -c $< -o $@
 
 $(BUILDDIR)/%_devel.o: $(GSRCDIR)/%.c $(SRCDIR)/*.h
-	$(CC) -g $(CFLAGS) $(DEVEL_VERSION_FLAGS) -c $< -o $@
+	$(CC) -g $(DEVEL_VERSION_FLAGS) -c $< -o $@
 
 $(GSRCDIR)/%.h: $(SRCDIR)/%.h $(firstword $(MAKEFILE_LIST))
 	cp $< $@

NEWS

@@ -1,9 +1,25 @@
-Version v1606-sbacchio/master:
-  -Twisted mass fermions operator for Nf=2
-  -Support of different shift on even and odd sites
+
+
+Version v1610 - TM Nf=2+1+1:
+
+  - Twisted mass fermions operator for Nf=1+1
+  - Aupport of different shift on even and odd sites
   -Personalized version of the library
+  - Following parameters have been added to input file:
+    -- //TODO
+
+
+Version v1606 - TM Nf=2:
+
+  - Twisted mass fermions operator for Nf=2
+  - A different TM shift can be applied on even and odd sites
+  - Personalized version of the library
+  - Following parameters have been added to input file:
+    -- //TODO
+
 
 Version v1606 - first release:
+
   This is the first release of the DDalphaAMG solver library.
   For an overview of its features included, please consult the
   user documentation in /doc.

README

@@ -22,6 +22,6 @@ INSTALL:
 
 HOWTO:
 
-  After having compiled the the user documentation via
+  After having compiled the user documentation via
   "make documentation" please consult the compiled PDF in /doc for
   further information.

sample.ini

@@ -33,48 +33,52 @@ right hand side: 2
 | 2 - twisted boundary cond.                   |
 |         (M_PI,M_PI,M_PI,M_PI)*t.b.c.         |
 |----------------------------------------------|
-boundary conditions: 2
+boundary conditions: 1
 twisted boundary conditions: 1 0 0 0
 number of levels: 3
-number of openmp threads: 2
+number of openmp threads: 1
 
 |--- depth 0 ----------------------------------|
 d0 global lattice: 8 8 8 8  // global lattice size
-d0 local lattice: 4 8 8 8   // lattice size on each process
-
+d0 local lattice: 8 8 4 4   // lattice size on each process
                              // nproc = prod(global lattice)/prod(local lattice)
-d0 block lattice: 2 2 2 2    // Schwarz blocks
+d0 block lattice: 4 4 4 4    // Schwarz blocks
+  d0 block lattice: 2 2 2 2    // Schwarz blocks
 d0 post smooth iter: 2       // number of Schwarz cycles per MG iteration
 d0 block iter: 4
 d0 test vectors: 24          // number of test vectors used to construct the interpolation
-d0 setup iter: 4             // number of bootstrap setup iteration (excluding the initial step)
+d0 setup iter: 3             // number of bootstrap setup iteration (excluding the initial step)
 
-d2 mu factor: 4.0
-d3 mu factor: 8.0
+d1 mu factor: 5.0
+  d2 mu factor: 5.0
 
 #wilson_param                // parameters for the inverter
-#the following OR kappa: ...
-m0: -0.5
-csw: 1.0
-#the following OR 2KappaMu: ...
-mu: 0.5
-mu odd shift: 0.0
-mu even shift: 0.0
+m0: -0.4
+csw: 1.6
+mu: 0.005
+setup mu: 0.001
+  mu odd shift: 0.0
+  mu even shift: 0.0
+epsbar: 0.11
+  epsbar odd shift: 0.0
+  epsbar even shift: 0.0
 addDownPropagator: 1
 
 tolerance for relative residual: 1E-10
-iterations between restarts: 50         // should be increased for ill-conditioned cnfgs    
-maximum of restarts: 20                 // should be increased for ill-conditioned cnfgs
-coarse grid tolerance: 5E-2
-coarse grid iterations: 100             // should be increased for ill-conditioned cnfgs
-coarse grid restarts: 5                 // should be increased for ill-conditioned cnfgs
+iterations between restarts: 20
+maximum of restarts: 50        
+coarse grid tolerance: 1E-2
+coarse grid iterations: 25     
+coarse grid restarts: 20       
 
 
 #general_param
 print mode: 1
 method: 2
-mixed precision: 1
+interpolation: 2
+mixed precision: 2
 randomize test vectors: 0               // initialize random number generator with time(0) ? 0=no/1=yes
+odd even preconditioning: 1
 
 // for further information, please read the user documentation in doc/
 // developers version of an input file in sample_devel.ini

src/DDalphaAMG.h

@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2016, Matthias Rottmann, Artur Strebel, Simon Heybrock, Simone Bacchio, Bjoern Leder.
+ * Copyright (C) 2016, Simone Bacchio.
  * 
  * This file is part of the DDalphaAMG solver library.
  * 
@@ -145,43 +145,82 @@
    **  mg_status.success = 0: not converged, 1: converged
    **  mg_status.info = final residual
    **/
-  void DDalphaAMG_solve( double *vector_out, double *vector_in, double tol,
-		     DDalphaAMG_status *mg_status );
+  void DDalphaAMG_solve( double *vector_out, double *vector_in,
+                         double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_doublet( double *vector1_out, double *vector1_in, 
+                                 double *vector2_out, double *vector2_in,
+                                 double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_ms_doublet( double **vector1_out, double *vector1_in, 
+                                    double **vector2_out, double *vector2_in,
+                                    double  *even_shifts, double *odd_shifts, int n_shifts,
+                                    double  *tol, DDalphaAMG_status *mg_status );
 
   /**
    ** Optional - Solve squared operator performing two inversions: 
    **   vector_out = D_{d/u}^{-1} * \Gamma_5 * D_{u/d}^{-1} * \Gamma_5 * vector_in.
    **  mg_status.success = 0: not converged, 1: converged
    **  mg_status.info = final residual
    **/
-  void DDalphaAMG_solve_squared( double *vector_out, double *vector_in, double tol,
-			     DDalphaAMG_status *mg_status );
+  void DDalphaAMG_solve_squared( double *vector_out, double *vector_in,
+                                 double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_doublet_squared( double *vector1_out, double *vector1_in, 
+                                         double *vector2_out, double *vector2_in,
+                                         double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_ms_doublet_squared( double **vector1_out, double *vector1_in, 
+                                            double **vector2_out, double *vector2_in,
+                                            double  *even_shifts, double *odd_shifts, int n_shifts,
+                                            double  *tol, DDalphaAMG_status *mg_status );
 
   /**
    ** Optional - Solve squared operator against the odd compoments performing two inversions: 
    **   vector_out = D_{d/u}^{-1} * \Gamma_5 * P_{odd} * D_{u/d}^{-1} * \Gamma_5 * vector_in.
    **  mg_status.success = 0: not converged, 1: converged
    **  mg_status.info = final residual
    **/
- void DDalphaAMG_solve_squared_odd( double *vector_out, double *vector_in, double tol,
-				DDalphaAMG_status *mg_status );
+  void DDalphaAMG_solve_squared_odd( double *vector_out, double *vector_in,
+                                     double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_doublet_squared_odd( double *vector1_out, double *vector1_in, 
+                                             double *vector2_out, double *vector2_in,
+                                             double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_ms_doublet_squared_odd( double **vector1_out, double *vector1_in, 
+                                                double **vector2_out, double *vector2_in,
+                                                double  *even_shifts, double *odd_shifts, int n_shifts,
+                                                double  *tol, DDalphaAMG_status *mg_status );
 
   /**
    ** Optional - Solve squared operator against the even compoments performing two inversions:
    **   vector_out = D_{d/u}^{-1} * \Gamma_5 * P_{even} * D_{u/d}^{-1} * \Gamma_5 * vector_in.
    **  mg_status.success = 0: not converged, 1: converged
    **  mg_status.info = final residual
    **/
- void DDalphaAMG_solve_squared_even( double *vector_out, double *vector_in, double tol,
-				 DDalphaAMG_status *mg_status );
+  void DDalphaAMG_solve_squared_even( double *vector_out, double *vector_in,
+                                      double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_doublet_squared_even( double *vector1_out, double *vector1_in, 
+                                              double *vector2_out, double *vector2_in,
+                                              double tol, DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_solve_ms_doublet_squared_even( double **vector1_out, double *vector1_in, 
+                                                 double **vector2_out, double *vector2_in,
+                                                 double  *even_shifts, double *odd_shifts, int n_shifts,
+                                                 double  *tol, DDalphaAMG_status *mg_status );
 
   /**
    ** Optional - Apply the operator:
    **   vector_out = D * vector_in.
    **  mg_status.success = 1
    **/
   void DDalphaAMG_apply_operator( double *vector_out, double *vector_in,
-			      DDalphaAMG_status *mg_status );
+                                  DDalphaAMG_status *mg_status );
+
+  void DDalphaAMG_apply_operator_doublet( double *vector1_out, double *vector1_in,
+                                          double *vector2_out, double *vector2_in, DDalphaAMG_status *mg_status );
 
   /**
    ** Optional - Apply a preconditioner step:
@@ -190,7 +229,9 @@
    **  mg_status.info = residual after preconditioning
    **/
   void DDalphaAMG_preconditioner( double *vector_out, double *vector_in,
-			      DDalphaAMG_status *mg_status );
+                                  DDalphaAMG_status *mg_status );
+  void DDalphaAMG_preconditioner_doublet( double *vector1_out, double *vector1_in,
+                                          double *vector2_out, double *vector2_in, DDalphaAMG_status *mg_status );
 
   /*
    *  Concluding the following functions have to be call for freeing the memory and finalizing
@@ -228,11 +269,11 @@
    **    -> mg_params.conf_index_fct = NULL, mg_params.vector_index_fct = NULL;
    **/
   void DDalphaAMG_read_configuration( double *gauge_field, char *filename, int format,
-				      DDalphaAMG_status *mg_status );
+                                      DDalphaAMG_status *mg_status );
   void DDalphaAMG_read_vector( double *vector_in, char *filename, int format,
-			       DDalphaAMG_status *mg_status );
+                               DDalphaAMG_status *mg_status );
   void DDalphaAMG_write_vector( double *vector_out, char *filename, int format,
-				DDalphaAMG_status *mg_status );
+                                DDalphaAMG_status *mg_status );
 
   /**
    ** Extra - Define vector with constant or random components
@@ -414,23 +455,39 @@
      ** Hopping parameter
      **/
     double kappa;
-    
+
     /**
      ** Twisted mass parameter and shifts on even/odd sites
      **/
     double mu;
     double mu_odd_shift;
     double mu_even_shift;
-    
+
     /**
-     ** Twisted mass factor for the preconditioner on each level.
+     ** Twisted mass factor for the preconditioner on each level, l.
      **  Default 6 on the coarsest level 
      ** 
-     **   -> mu_o[l] = (mu + mu_odd_shift)  * mu_factor
-     **   -> mu_e[l] = (mu + mu_even_shift) * mu_factor
+     **   -> mu_o[l] = (mu + mu_odd_shift)  * mu_factor[l]
+     **   -> mu_e[l] = (mu + mu_even_shift) * mu_factor[l]
      **/
     double mu_factor[MAX_MG_LEVELS];  
 
+    /**
+     ** Twisted mass doublet parameter and shifts on even/odd sites
+     **/
+    double epsbar;
+    double epsbar_ig5_odd_shift;
+    double epsbar_ig5_even_shift;
+
+    /**
+     ** Twisted mass doublet factor for the preconditioner on each level, l.
+     **  Default 6 on the coarsest level 
+     ** 
+     **   -> epsbar_o[l] = ( epsbar + i * gamma_5 * epsbar_ig5_odd_shift ) * epsbar_factor[l]
+     **   -> epsbar_e[l] = ( epsbar + i * gamma_5 * epsbar_ig5_even_shift ) * epsbar_factor[l]
+     **/
+    double epsbar_factor[MAX_MG_LEVELS];  
+
     /**
      ** Function returning the index of a element at the corresponding
      **    position (t,z,y,x are local position w.r.t the process ).