Mybranch

include/tmLQCD.h

@@ -41,8 +41,20 @@ extern "C"
   typedef struct {
     unsigned int nproc, nproc_t, nproc_x, nproc_y, nproc_z, cart_id, proc_id, time_rank, omp_num_threads;
     unsigned int proc_coords[4];
+#ifdef TM_USE_MPI
+    MPI_Comm cart_grid;
+#endif
   } tmLQCD_mpi_params;
 
+  typedef struct {
+    char type_name[100];
+    int  eoprec;
+    void *evecs;
+    void *evals;
+    int  prec;
+    int  nev;
+  } tmLQCD_deflator_params;
+
   int tmLQCD_invert_init(int argc, char *argv[], const int verbose, const int external_id);
   int tmLQCD_read_gauge(const int nconfig);
   int tmLQCD_invert(double * const propagator, double * const source,
@@ -58,6 +70,12 @@ extern "C"
                     const int op_id, const int gauge_persist);
 #endif
 
+  int tmLQCD_invert_eo(double * const propagator, double * const source, const int op_id);
+
+  int tmLQCD_get_deflator_params(tmLQCD_deflator_params*params, const int op_id);
+  int tmLQCD_init_deflator(const int op_id);
+  int tmLQCD_set_deflator_fields(const int op_id1, const int op_id2);
+
 #ifdef __cplusplus
 }
 #endif

invert_eo.h

@@ -28,6 +28,7 @@
 #define _INVERT_EO_H
 #include "global.h"
 #include "solver/solver_params.h"
+#include "solver/deflator.h"
 
 int invert_eo(spinor * const Even_new, spinor * const Odd_new, 
               spinor * const Even, spinor * const Odd,

linalg/arpack.h

@@ -0,0 +1,74 @@
+/*************************************************************************
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008 Carsten Urbach
+ *
+ * This file is part of tmLQCD.
+ *
+ * tmLQCD is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ * 
+ * tmLQCD is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with tmLQCD.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *
+ * Header file to interface with ARPACK package for solving large 
+ * eigenvalue problems
+ *
+ * Author: Abdou M. Abdel-Rehim (amabdelrehim@gmail.com) 2014
+ *************************************************************************/
+
+
+#ifndef _ARPACK_H
+#define _ARPACK_H
+
+#include <complex.h>
+#include "linalg/fortran.h"
+#ifdef TM_USE_MPI
+# include <mpi.h>
+#endif
+
+
+//ARPACK initlog and finilog routines for printing the ARPACK log (same for serial and parallel version)
+extern int _AFT(initlog) (int*, char*, int);
+extern int _AFT(finilog) (int*);
+
+
+//ARPACK driver routines for computing eigenvectors (serial version) 
+extern int _AFT(znaupd) (int *ido, char *bmat, int *n, char *which, int *nev, double *tol,
+                         _Complex double *resid, int *ncv, _Complex double *v, int *ldv, 
+                         int *iparam, int *ipntr, _Complex double *workd, _Complex double *workl, 
+                         int *lworkl, double *rwork, int *info, int bmat_size, int which_size );
+
+extern int _AFT(zneupd) (int *comp_evecs, char *howmany, int *select, _Complex double *evals, 
+			 _Complex double *v, int *ldv, _Complex double *sigma, _Complex double *workev, 
+			 char *bmat, int *n, char *which, int *nev, double *tol, _Complex double *resid, 
+                         int *ncv, _Complex double *v1, int *ldv1, int *iparam, int *ipntr, 
+                         _Complex double *workd, _Complex double *workl, int *lworkl, double *rwork, int *info,
+                         int howmany_size, int bmat_size, int which_size);
+
+extern int _AFT(mcinitdebug)(int*,int*,int*,int*,int*,int*,int*,int*);
+
+//PARPACK routines (parallel version)
+#ifdef TM_USE_MPI
+extern int _AFT(pznaupd) (int *comm, int *ido, char *bmat, int *n, char *which, int *nev, double *tol,
+                         _Complex double *resid, int *ncv, _Complex double *v, int *ldv, 
+                         int *iparam, int *ipntr, _Complex double *workd, _Complex double *workl, 
+                         int *lworkl, double *rwork, int *info, int bmat_size, int which_size );
+
+extern int _AFT(pzneupd) (int *comm, int *comp_evecs, char *howmany, int *select, _Complex double *evals, 
+                         _Complex double *v, int *ldv, _Complex double *sigma, _Complex double *workev, 
+                         char *bmat, int *n, char *which, int *nev, double *tol, _Complex double *resid, 
+                         int *ncv, _Complex double *v1, int *ldv1, int *iparam, int *ipntr, 
+                         _Complex double *workd, _Complex double *workl, int *lworkl, double *rwork, int *info,
+                         int howmany_size, int bmat_size, int which_size);
+
+extern int _AFT(pmcinitdebug)(int*,int*,int*,int*,int*,int*,int*,int*);
+#endif
+
+#endif

linalg/assign.c

@@ -34,6 +34,7 @@
 #include <math.h>
 #include <string.h>
 #include "su3.h"
+#include "global.h"
 #include "assign.h"
 
 
@@ -67,3 +68,164 @@ void assign_su3vect(su3_vector * const R, su3_vector * const S, const int N)
   }
 }
 #endif
+
+
+/* copy a complex double S of size N into a spinor R of size N/12  */
+/* S input, R output        */
+/* S and R must not overlap */
+void assign_complex_to_spinor(spinor * const R, _Complex double * const S, const int N)
+{
+  //check the dimension
+  if( (N%12) != 0){
+    if(g_proc_id==g_stdio_proc)
+      fprintf(stderr,"Error in assign_double_to_spinor. N is not a multiple of 12");
+  }
+
+  int k; //spinor index
+  spinor *r;
+  _Complex double *s;
+
+  k=0;
+  for(int ix=0; ix<N ; ix +=12)
+  {
+     s=S+ix;
+     r=R+k;
+
+     (r->s0).c0 = *s;
+     (r->s0).c1 = *(s+1);
+     (r->s0).c2 = *(s+2);
+
+     (r->s1).c0 = *(s+3);
+     (r->s1).c1 = *(s+4);
+     (r->s1).c2 = *(s+5);
+
+     (r->s2).c0 = *(s+6);
+     (r->s2).c1 = *(s+7);
+     (r->s2).c2 = *(s+8);
+
+
+     (r->s3).c0 = *(s+9);
+     (r->s3).c1 = *(s+10);
+     (r->s3).c2 = *(s+11);
+
+     k++;
+  } 
+
+  return;
+}
+
+
+void assign_complex_to_spinor_32(spinor32 * const R, _Complex float* const S, const int N)
+{
+  //check the dimension
+  if( (N%12) != 0){
+    if(g_proc_id==g_stdio_proc)
+      fprintf(stderr,"Error in assign_double_to_spinor. N is not a multiple of 12");
+  }
+
+  int k; //spinor index
+  spinor32 *r;
+  _Complex float *s;
+
+  k=0;
+  for(int ix=0; ix<N ; ix +=12)
+  {
+     s=S+ix;
+     r=R+k;
+
+     (r->s0).c0 = *s;
+     (r->s0).c1 = *(s+1);
+     (r->s0).c2 = *(s+2);
+
+     (r->s1).c0 = *(s+3);
+     (r->s1).c1 = *(s+4);
+     (r->s1).c2 = *(s+5);
+
+     (r->s2).c0 = *(s+6);
+     (r->s2).c1 = *(s+7);
+     (r->s2).c2 = *(s+8);
+
+
+     (r->s3).c0 = *(s+9);
+     (r->s3).c1 = *(s+10);
+     (r->s3).c2 = *(s+11);
+
+     k++;
+  } 
+
+  return;
+}
+
+/* copy a spinor S of size N into a complex double R of size 12*N  */
+/* S input, R output        */
+/* S and R must not overlap */
+void assign_spinor_to_complex(_Complex double * const R, spinor * const S, const int N)
+{
+
+  int k; //complex double index
+  _Complex double *r;
+  spinor *s;
+
+  k=0;
+  for(int ix=0; ix<N ; ix++)
+  {
+     s=S+ix;
+     r=R+k;
+
+     *r     =  (s->s0).c0;
+     *(r+1) =  (s->s0).c1;
+     *(r+2) =  (s->s0).c2;
+
+     *(r+3) =  (s->s1).c0;
+     *(r+4) =  (s->s1).c1;
+     *(r+5) =  (s->s1).c2;
+
+     *(r+6) =  (s->s2).c0;
+     *(r+7) =  (s->s2).c1;
+     *(r+8) =  (s->s2).c2;
+
+     *(r+9)  =  (s->s3).c0;
+     *(r+10) =  (s->s3).c1;
+     *(r+11) =  (s->s3).c2;
+
+     k +=12;
+  } 
+
+  return;
+}
+
+
+void assign_spinor_to_complex_32(_Complex float * const R, spinor32 * const S, const int N)
+{
+
+  int k; //complex double index
+  _Complex float *r;
+  spinor32 *s;
+
+  k=0;
+  for(int ix=0; ix<N ; ix++)
+  {
+     s=S+ix;
+     r=R+k;
+
+     *r     =  (s->s0).c0;
+     *(r+1) =  (s->s0).c1;
+     *(r+2) =  (s->s0).c2;
+
+     *(r+3) =  (s->s1).c0;
+     *(r+4) =  (s->s1).c1;
+     *(r+5) =  (s->s1).c2;
+
+     *(r+6) =  (s->s2).c0;
+     *(r+7) =  (s->s2).c1;
+     *(r+8) =  (s->s2).c2;
+
+     *(r+9)  =  (s->s3).c0;
+     *(r+10) =  (s->s3).c1;
+     *(r+11) =  (s->s3).c2;
+
+     k +=12;
+  } 
+
+  return;
+}

linalg/assign.h

@@ -27,4 +27,10 @@ void assign(spinor * const R, spinor * const S, const int N);
 void assign_32(spinor32 * const R, spinor32 * const S, const int N);
 void assign_su3vect(su3_vector * const R, su3_vector * const S, const int N);
 
+void assign_complex_to_spinor(spinor * const R, _Complex double * const S, const int N); /* N is the size of S */
+void assign_spinor_to_complex(_Complex double * const R, spinor * const S, const int N); /* N is the size of S */
+
+void assign_complex_to_spinor_32(spinor32 * const R, _Complex float* const S, const int N);
+void assign_spinor_to_complex_32(_Complex float* const R, spinor32 * const S, const int N); /* N is the size of S */
+
 #endif

linalg/fortran.h

@@ -25,4 +25,10 @@
 #define _FT(s) s ## _
 #endif
 
+#if (defined NOARPACKUNDERSCORE)
+#define _AFT(s) s
+#else
+#define _AFT(s) s ## _
+#endif
+
 #endif