matrices.c File Reference

Matrix multiplication and related routines. More...

#include <assert.h>
#include <complex.h>
#include <matrices.h>
#include <string.h>

Include dependency graph for matrices.c:

Go to the source code of this file.

Functions
int	Dslash (Complex *phi, Complex *r, Complex *u11t, Complex *u12t, unsigned int *iu, unsigned int *id, Complex *gamval, int *gamin, double *dk4m, double *dk4p, Complex_f jqq, float akappa)
	Evaluates \(\Phi=M r\) in double precision.
int	Dslashd (Complex *phi, Complex *r, Complex *u11t, Complex *u12t, unsigned int *iu, unsigned int *id, Complex *gamval, int *gamin, double *dk4m, double *dk4p, Complex_f jqq, float akappa)
	Evaluates \(\Phi=M^\dagger r\) in double precision.
int	Hdslash (Complex *phi, Complex *r, Complex *u11t, Complex *u12t, unsigned int *iu, unsigned int *id, Complex *gamval, int *gamin, double *dk4m, double *dk4p, float akappa)
	Evaluates \(\Phi=M r\) in double precision.
int	Hdslashd (Complex *phi, Complex *r, Complex *u11t, Complex *u12t, unsigned int *iu, unsigned int *id, Complex *gamval, int *gamin, double *dk4m, double *dk4p, float akappa)
	Evaluates \(\Phi=M^\dagger r\) in double precision.
int	Dslash_f (Complex_f *phi, Complex_f *r, Complex_f *u11t_f, Complex_f *u12t_f, unsigned int *iu, unsigned int *id, Complex_f *gamval_f, int *gamin, float *dk4m_f, float *dk4p_f, Complex_f jqq, float akappa)
	Evaluates \(\Phi=M r\) in single precision.
int	Dslashd_f (Complex_f *phi, Complex_f *r, Complex_f *u11t_f, Complex_f *u12t_f, unsigned int *iu, unsigned int *id, Complex_f *gamval_f, int *gamin, float *dk4m_f, float *dk4p_f, Complex_f jqq, float akappa)
	Evaluates \(\Phi=M^\dagger r\) in single precision.
int	Hdslash_f (Complex_f *phi, Complex_f *r, Complex_f *u11t_f, Complex_f *u12t_f, unsigned int *iu, unsigned int *id, Complex_f *gamval_f, int *gamin, float *dk4m_f, float *dk4p_f, float akappa)
	Evaluates \(\Phi=M r\) in single precision.
int	Hdslashd_f (Complex_f *phi, Complex_f *r, Complex_f *u11t_f, Complex_f *u12t_f, unsigned int *iu, unsigned int *id, Complex_f *gamval_f, int *gamin, float *dk4m_f, float *dk4p_f, float akappa)
	Evaluates \(\Phi=M^\dagger r\) in single precision.
int	Reunitarise (Complex *u11t, Complex *u12t)
	Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.

Detailed Description

Matrix multiplication and related routines.

There are two four matrix mutiplication routines, and each had a double and single (_f) version The Hdslash? routines are called when acting on half of the fermions (up/down flavour partitioning) The Dslash routines act on everything

Any routine ending in a d is the daggered multiplication

Definition in file matrices.c.

Function Documentation

Dslash()

int Dslash	(	Complex *	phi,
		Complex *	r,
		Complex *	u11t,
		Complex *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex *	gamval,
		int *	gamin,
		double *	dk4m,
		double *	dk4p,
		Complex_f	jqq,
		float	akappa )

Evaluates \(\Phi=M r\) in double precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
jqq	Diquark source
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 18 of file matrices.c.

                                                                                           {
   /*
    * @brief Evaluates @f(\Phi=M r@f) in double precision.
    *
    * @param   phi:        The product
    * @param   r:          The array being acted on by M
    * @param   u11t:    First colour trial field
    * @param   u12t:    Second colour trial field
    * @param   iu:         Upper halo indices
    * @param   id:         Lower halo indices
    * @param   gamval:  Gamma matrices
    * @param   gamin:      Indices for dirac terms
    * @param   dk4m:    
    * @param   dk4p:    
    * @param   jqq:        Diquark source
    * @param   akappa:     Hopping parameter
    *
    * @see ZHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Dslash";
   //Get the halos in order
#if(nproc>1)
   ZHalo_swap_all(r, 16);
#endif
 
   //Mass term
   //Diquark Term (antihermitian)
#ifdef __NVCC__
   cuDslash(phi,r,u11t,u12t,iu,id,gamval,gamin,dk4m,dk4p,jqq,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm*sizeof(Complex));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#pragma omp simd aligned(phi,r,gamval:AVX)
      for(int idirac = 0; idirac<ndirac; idirac++){
         int igork = idirac+4;
         Complex a_1, a_2;
         a_1=conj(jqq)*gamval[4*ndirac+idirac];
         //We subtract a_2, hence the minus
         a_2=-jqq*gamval[4*ndirac+idirac];
         phi[(i*ngorkov+idirac)*nc]+=a_1*r[(i*ngorkov+igork)*nc+0];
         phi[(i*ngorkov+idirac)*nc+1]+=a_1*r[(i*ngorkov+igork)*nc+1];
         phi[(i*ngorkov+igork)*nc+0]+=a_2*r[(i*ngorkov+idirac)*nc];
         phi[(i*ngorkov+igork)*nc+1]+=a_2*r[(i*ngorkov+idirac)*nc+1];
      }
 
      //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t,u12t,gamval:AVX)
         for(int igorkov=0; igorkov<ngorkov; igorkov++){
            //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing in the dirac term.
            int idirac=igorkov%4;      
            int igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            phi[(i*ngorkov+igorkov)*nc]+=-akappa*(u11t[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc]+\
                  u12t[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc+1]+\
                  conj(u11t[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]-\
                  u12t[did*ndim+mu]*r[(did*ngorkov+igorkov)*nc+1])+\
                                         //Dirac term. Reminder! gamval was rescaled by kappa when we defined it
                                         gamval[mu*ndirac+idirac]*(u11t[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc]+\
                                               u12t[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc+1]-\
                                               conj(u11t[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]+\
                                               u12t[did*ndim+mu]*r[(did*ngorkov+igork1)*nc+1]);
 
            phi[(i*ngorkov+igorkov)*nc+1]+=-akappa*(-conj(u12t[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc]+\
                  conj(u11t[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc+1]+\
                  conj(u12t[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]+\
                  u11t[did*ndim+mu]*r[(did*ngorkov+igorkov)*nc+1])+\
                                           //Dirac term
                                           gamval[mu*ndirac+idirac]*(-conj(u12t[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc]+\
                                                 conj(u11t[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc+1]-\
                                                 conj(u12t[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]-\
                                                 u11t[did*ndim+mu]*r[(did*ngorkov+igork1)*nc+1]);
         }
      }
      //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
      //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7  separately
      //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
#endif
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t,u12t,dk4m,dk4p:AVX)
      for(int igorkov=0; igorkov<4; igorkov++){
         int igorkovPP=igorkov+4;   //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
                                    //the FORTRAN code did it.
         int igork1 = gamin[3*ndirac+igorkov];  int igork1PP = igork1+4;
 
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         phi[(i*ngorkov+igorkov)*nc]+=
            -dk4p[i]*(u11t[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc]-r[(uid*ngorkov+igork1)*nc])
                  +u12t[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc+1]-r[(uid*ngorkov+igork1)*nc+1]))
            -dk4m[did]*(conj(u11t[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]+r[(did*ngorkov+igork1)*nc])
                  -u12t[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]+r[(did*ngorkov+igork1)*nc+1]));
         phi[(i*ngorkov+igorkov)*nc+1]+=
            -dk4p[i]*(-conj(u12t[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc]-r[(uid*ngorkov+igork1)*nc])
                  +conj(u11t[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc+1]-r[(uid*ngorkov+igork1)*nc+1]))
            -dk4m[did]*(conj(u12t[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]+r[(did*ngorkov+igork1)*nc])
                  +u11t[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]+r[(did*ngorkov+igork1)*nc+1]));
 
         //And the +4 terms. Note that dk4p and dk4m swap positions compared to the above          
         phi[(i*ngorkov+igorkovPP)*nc]+=-dk4m[i]*(u11t[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc]-r[(uid*ngorkov+igork1PP)*nc])+\
               u12t[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc+1]-r[(uid*ngorkov+igork1PP)*nc+1]))-\
                                        dk4p[did]*(conj(u11t[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]+r[(did*ngorkov+igork1PP)*nc])-\
                                              u12t[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]+r[(did*ngorkov+igork1PP)*nc+1]));
 
         phi[(i*ngorkov+igorkovPP)*nc+1]+=-dk4m[i]*(conj(-u12t[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc]-r[(uid*ngorkov+igork1PP)*nc])+\
               conj(u11t[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc+1]-r[(uid*ngorkov+igork1PP)*nc+1]))-\
                                          dk4p[did]*(conj(u12t[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]+r[(did*ngorkov+igork1PP)*nc])+\
                                                u11t[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]+r[(did*ngorkov+igork1PP)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References Complex, kferm, kvol, nc, ndim, ndirac, ngorkov, and ZHalo_swap_all().

Here is the call graph for this function:

Dslash_f()

int Dslash_f	(	Complex_f *	phi,
		Complex_f *	r,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	gamval,
		int *	gamin,
		float *	dk4m,
		float *	dk4p,
		Complex_f	jqq,
		float	akappa )

Evaluates \(\Phi=M r\) in single precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
jqq	Diquark source
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 462 of file matrices.c.

                                                                                                 {
   /*
    * @brief Evaluates @f(\Phi=M r@f) in single precision.
    *
    * @param   phi:        The product
    * @param   r:          The array being acted on by M
    * @param   u11t_f:     First colour trial field
    * @param   u12t_f:     Second colour trial field
    * @param   iu:         Upper halo indices
    * @param   id:         Lower halo indices
    * @param   gamval_f:   Gamma matrices
    * @param   gamin:      Indices for dirac terms
    * @param   dk4m_f:     
    * @param   dk4p_f:     
    * @param   jqq:        Diquark source
    * @param   akappa:     Hopping parameter
    *
    * @see CHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Dslash_f";
   //Get the halos in order
#if(nproc>1)
   CHalo_swap_all(r, 16);
#endif
 
   //Mass term
   //Diquark Term (antihermitian)
#ifdef __NVCC__
   cuDslash_f(phi,r,u11t_f,u12t_f,iu,id,gamval_f,gamin,dk4m_f,dk4p_f,jqq,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm*sizeof(Complex_f));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#pragma omp simd aligned(phi,r,gamval_f:AVX)
      for(int idirac = 0; idirac<ndirac; idirac++){
         int igork = idirac+4;
         Complex_f a_1, a_2;
         a_1=conj(jqq)*gamval_f[4*ndirac+idirac];
         //We subtract a_2, hence the minus
         a_2=-jqq*gamval_f[4*ndirac+idirac];
         phi[(i*ngorkov+idirac)*nc]+=a_1*r[(i*ngorkov+igork)*nc+0];
         phi[(i*ngorkov+idirac)*nc+1]+=a_1*r[(i*ngorkov+igork)*nc+1];
         phi[(i*ngorkov+igork)*nc+0]+=a_2*r[(i*ngorkov+idirac)*nc];
         phi[(i*ngorkov+igork)*nc+1]+=a_2*r[(i*ngorkov+idirac)*nc+1];
      }
 
      //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,gamval_f,gamin:AVX)
         for(int igorkov=0; igorkov<ngorkov; igorkov++){
            //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing in the dirac term.
            int idirac=igorkov%4;      
            int igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            phi[(i*ngorkov+igorkov)*nc]+=-akappa*(u11t_f[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc]+\
                  u12t_f[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc+1]+\
                  conj(u11t_f[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]-\
                  u12t_f[did*ndim+mu]*r[(did*ngorkov+igorkov)*nc+1])+\
                                         //Dirac term. Reminder! gamval was rescaled by kappa when we defined it
                                         gamval_f[mu*ndirac+idirac]*(u11t_f[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc]+\
                                               u12t_f[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc+1]-\
                                               conj(u11t_f[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]+\
                                               u12t_f[did*ndim+mu]*r[(did*ngorkov+igork1)*nc+1]);
 
            phi[(i*ngorkov+igorkov)*nc+1]+=-akappa*(-conj(u12t_f[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc]+\
                  conj(u11t_f[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc+1]+\
                  conj(u12t_f[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]+\
                  u11t_f[did*ndim+mu]*r[(did*ngorkov+igorkov)*nc+1])+\
                                           //Dirac term
                                           gamval_f[mu*ndirac+idirac]*(-conj(u12t_f[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc]+\
                                                 conj(u11t_f[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc+1]-\
                                                 conj(u12t_f[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]-\
                                                 u11t_f[did*ndim+mu]*r[(did*ngorkov+igork1)*nc+1]);
         }
      }
      //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
      //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7  separately
      //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
#endif
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,dk4m_f,dk4p_f,gamin:AVX)
      for(int igorkov=0; igorkov<4; igorkov++){
         int igorkovPP=igorkov+4;   //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
                                    //the FORTRAN code did it.
         int igork1 = gamin[3*ndirac+igorkov];  int igork1PP = igork1+4;
 
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         phi[(i*ngorkov+igorkov)*nc]+=
            -dk4p_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc]-r[(uid*ngorkov+igork1)*nc])
                  +u12t_f[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc+1]-r[(uid*ngorkov+igork1)*nc+1]))
            -dk4m_f[did]*(conj(u11t_f[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]+r[(did*ngorkov+igork1)*nc])
                  -u12t_f[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]+r[(did*ngorkov+igork1)*nc+1]));
         phi[(i*ngorkov+igorkov)*nc+1]+=
            -dk4p_f[i]*(-conj(u12t_f[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc]-r[(uid*ngorkov+igork1)*nc])
                  +conj(u11t_f[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc+1]-r[(uid*ngorkov+igork1)*nc+1]))
            -dk4m_f[did]*(conj(u12t_f[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]+r[(did*ngorkov+igork1)*nc])
                  +u11t_f[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]+r[(did*ngorkov+igork1)*nc+1]));
 
         //And the +4 terms. Note that dk4p_f and dk4m_f swap positions compared to the above            
         phi[(i*ngorkov+igorkovPP)*nc]+=-dk4m_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc]-r[(uid*ngorkov+igork1PP)*nc])
               +u12t_f[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc+1]-r[(uid*ngorkov+igork1PP)*nc+1]))
            -dk4p_f[did]*(conj(u11t_f[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]+r[(did*ngorkov+igork1PP)*nc])
                  -u12t_f[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]+r[(did*ngorkov+igork1PP)*nc+1]));
 
         phi[(i*ngorkov+igorkovPP)*nc+1]+=-dk4m_f[i]*(conj(-u12t_f[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc]-r[(uid*ngorkov+igork1PP)*nc])
               +conj(u11t_f[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc+1]-r[(uid*ngorkov+igork1PP)*nc+1]))
            -dk4p_f[did]*(conj(u12t_f[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]+r[(did*ngorkov+igork1PP)*nc])
                  +u11t_f[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]+r[(did*ngorkov+igork1PP)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References CHalo_swap_all(), Complex_f, kferm, kvol, nc, ndim, ndirac, and ngorkov.

Here is the call graph for this function:

Here is the caller graph for this function:

Dslashd()

int Dslashd	(	Complex *	phi,
		Complex *	r,
		Complex *	u11t,
		Complex *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex *	gamval,
		int *	gamin,
		double *	dk4m,
		double *	dk4p,
		Complex_f	jqq,
		float	akappa )

Evaluates \(\Phi=M^\dagger r\) in double precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
jqq	Diquark source
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 140 of file matrices.c.

                                                                                           {
   /*
    * @brief Evaluates @f(\Phi=M^\dagger r@f) in double precision.
    *
    * @param   phi:        The product
    * @param   r:          The array being acted on by M
    * @param   u11t:    First colour trial field
    * @param   u12t:    Second colour trial field
    * @param   iu:         Upper halo indices
    * @param   id:         Lower halo indices
    * @param   gamval:  Gamma matrices
    * @param   gamin:      Indices for dirac terms
    * @param   dk4m:    
    * @param   dk4p:    
    * @param   jqq:        Diquark source
    * @param   akappa:     Hopping parameter
    *
    * @see ZHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Dslashd";
   //Get the halos in order
#if(nproc>1)
   ZHalo_swap_all(r, 16);
#endif
 
   //Mass term
#ifdef __NVCC__
   cuDslashd(phi,r,u11t,u12t,iu,id,gamval,gamin,dk4m,dk4p,jqq,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm*sizeof(Complex));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#pragma omp simd aligned(phi,r,gamval:AVX)
      //Diquark Term (antihermitian) The signs of a_1 and a_2 below flip under dagger
      for(int idirac = 0; idirac<ndirac; idirac++){
         int igork = idirac+4;
         Complex a_1, a_2;
         //We subtract a_1, hence the minus
         a_1=-conj(jqq)*gamval[4*ndirac+idirac];
         a_2=jqq*gamval[4*ndirac+idirac];
         phi[(i*ngorkov+idirac)*nc]+=a_1*r[(i*ngorkov+igork)*nc];
         phi[(i*ngorkov+idirac)*nc+1]+=a_1*r[(i*ngorkov+igork)*nc+1];
         phi[(i*ngorkov+igork)*nc]+=a_2*r[(i*ngorkov+idirac)*nc];
         phi[(i*ngorkov+igork)*nc+1]+=a_2*r[(i*ngorkov+idirac)*nc+1];
      }
 
      //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t,u12t,gamval:AVX)
         for(int igorkov=0; igorkov<ngorkov; igorkov++){
            //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing.
            int idirac=igorkov%4;      
            int igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            //Reminder! gamval was rescaled by kappa when we defined it
            phi[(i*ngorkov+igorkov)*nc]+=
               -akappa*(      u11t[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc]
                     +u12t[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc+1]
                     +conj(u11t[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]
                     -u12t[did*ndim+mu] *r[(did*ngorkov+igorkov)*nc+1])
               -gamval[mu*ndirac+idirac]*
               (          u11t[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc]
                          +u12t[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc+1]
                          -conj(u11t[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]
                          +u12t[did*ndim+mu] *r[(did*ngorkov+igork1)*nc+1]);
 
            phi[(i*ngorkov+igorkov)*nc+1]+=
               -akappa*(-conj(u12t[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc]
                     +conj(u11t[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc+1]
                     +conj(u12t[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]
                     +u11t[did*ndim+mu] *r[(did*ngorkov+igorkov)*nc+1])
               -gamval[mu*ndirac+idirac]*
               (-conj(u12t[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc]
                +conj(u11t[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc+1]
                -conj(u12t[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]
                -u11t[did*ndim+mu] *r[(did*ngorkov+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
      //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7  separately
      //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
      //Under dagger, dk4p and dk4m get swapped and the dirac component flips sign.
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t,u12t,dk4m,dk4p:AVX)
      for(int igorkov=0; igorkov<4; igorkov++){
         //the FORTRAN code did it.
         int igork1 = gamin[3*ndirac+igorkov];  
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         phi[(i*ngorkov+igorkov)*nc]+=
            -dk4m[i]*(u11t[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc]+r[(uid*ngorkov+igork1)*nc])
                  +u12t[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc+1]+r[(uid*ngorkov+igork1)*nc+1]))
            -dk4p[did]*(conj(u11t[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]-r[(did*ngorkov+igork1)*nc])
                  -u12t[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]-r[(did*ngorkov+igork1)*nc+1]));
         phi[(i*ngorkov+igorkov)*nc+1]+=
            -dk4m[i]*(-conj(u12t[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc]+r[(uid*ngorkov+igork1)*nc])
                  +conj(u11t[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc+1]+r[(uid*ngorkov+igork1)*nc+1]))
            -dk4p[did]*(conj(u12t[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]-r[(did*ngorkov+igork1)*nc])
                  +u11t[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]-r[(did*ngorkov+igork1)*nc+1]));
 
 
         int igorkovPP=igorkov+4;   //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
         int igork1PP = igork1+4;
         //And the +4 terms. Note that dk4p and dk4m swap positions compared to the above          
         phi[(i*ngorkov+igorkovPP)*nc]+=-dk4p[i]*(u11t[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc]+r[(uid*ngorkov+igork1PP)*nc])+\
               u12t[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc+1]+r[(uid*ngorkov+igork1PP)*nc+1]))-\
                                        dk4m[did]*(conj(u11t[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]-r[(did*ngorkov+igork1PP)*nc])-\
                                              u12t[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]-r[(did*ngorkov+igork1PP)*nc+1]));
 
         phi[(i*ngorkov+igorkovPP)*nc+1]+=dk4p[i]*(conj(u12t[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc]+r[(uid*ngorkov+igork1PP)*nc])-\
               conj(u11t[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc+1]+r[(uid*ngorkov+igork1PP)*nc+1]))-\
                                          dk4m[did]*(conj(u12t[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]-r[(did*ngorkov+igork1PP)*nc])+
                                                u11t[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]-r[(did*ngorkov+igork1PP)*nc+1]));
 
      }
#endif
   }
#endif
   return 0;
}

References Complex, kferm, kvol, nc, ndim, ndirac, ngorkov, and ZHalo_swap_all().

Here is the call graph for this function:

Dslashd_f()

int Dslashd_f	(	Complex_f *	phi,
		Complex_f *	r,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	gamval,
		int *	gamin,
		float *	dk4m,
		float *	dk4p,
		Complex_f	jqq,
		float	akappa )

Evaluates \(\Phi=M^\dagger r\) in single precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
jqq	Diquark source
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 584 of file matrices.c.

                                                                                                 {
   /*
    * @brief Evaluates @f(\Phi=M^\dagger r@f) in single precision.
    *
    * @param   phi:        The product
    * @param   r:          The array being acted on by M
    * @param   u11t_f:     First colour trial field
    * @param   u12t_f:     Second colour trial field
    * @param   iu:         Upper halo indices
    * @param   id:         Lower halo indices
    * @param   gamval_f:   Gamma matrices
    * @param   gamin:      Indices for dirac terms
    * @param   dk4m_f:     
    * @param   dk4p_f:     
    * @param   jqq:        Diquark source
    * @param   akappa:     Hopping parameter
    *
    * @see CHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Dslashd_f";
   //Get the halos in order
#if(nproc>1)
   CHalo_swap_all(r, 16);
#endif
 
   //Mass term
#ifdef __NVCC__
   cuDslashd_f(phi,r,u11t_f,u12t_f,iu,id,gamval_f,gamin,dk4m_f,dk4p_f,jqq,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm*sizeof(Complex_f));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#pragma omp simd aligned(phi,r,gamval_f:AVX)
      //Diquark Term (antihermitian) The signs of a_1 and a_2 below flip under dagger
      for(int idirac = 0; idirac<ndirac; idirac++){
         int igork = idirac+4;
         Complex_f a_1, a_2;
         //We subtract a_1, hence the minus
         a_1=-conj(jqq)*gamval_f[4*ndirac+idirac];
         a_2=jqq*gamval_f[4*ndirac+idirac];
         phi[(i*ngorkov+idirac)*nc]+=a_1*r[(i*ngorkov+igork)*nc];
         phi[(i*ngorkov+idirac)*nc+1]+=a_1*r[(i*ngorkov+igork)*nc+1];
         phi[(i*ngorkov+igork)*nc]+=a_2*r[(i*ngorkov+idirac)*nc];
         phi[(i*ngorkov+igork)*nc+1]+=a_2*r[(i*ngorkov+idirac)*nc+1];
      }
 
      //Spacelike terms. Here's hoping I haven't put time as the zeroth component somewhere!
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,gamval_f:AVX)
         for(int igorkov=0; igorkov<ngorkov; igorkov++){
            //FORTRAN had mod((igorkov-1),4)+1 to prevent issues with non-zero indexing.
            int idirac=igorkov%4;      
            int igork1 = (igorkov<4) ? gamin[mu*ndirac+idirac] : gamin[mu*ndirac+idirac]+4;
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            //Reminder! gamval was rescaled by kappa when we defined it
            phi[(i*ngorkov+igorkov)*nc]+=
               -akappa*(      u11t_f[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc]
                     +u12t_f[i*ndim+mu]*r[(uid*ngorkov+igorkov)*nc+1]
                     +conj(u11t_f[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]
                     -u12t_f[did*ndim+mu] *r[(did*ngorkov+igorkov)*nc+1])
               -gamval_f[mu*ndirac+idirac]*
               (          u11t_f[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc]
                          +u12t_f[i*ndim+mu]*r[(uid*ngorkov+igork1)*nc+1]
                          -conj(u11t_f[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]
                          +u12t_f[did*ndim+mu] *r[(did*ngorkov+igork1)*nc+1]);
 
            phi[(i*ngorkov+igorkov)*nc+1]+=
               -akappa*(-conj(u12t_f[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc]
                     +conj(u11t_f[i*ndim+mu])*r[(uid*ngorkov+igorkov)*nc+1]
                     +conj(u12t_f[did*ndim+mu])*r[(did*ngorkov+igorkov)*nc]
                     +u11t_f[did*ndim+mu] *r[(did*ngorkov+igorkov)*nc+1])
               -gamval_f[mu*ndirac+idirac]*
               (-conj(u12t_f[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc]
                +conj(u11t_f[i*ndim+mu])*r[(uid*ngorkov+igork1)*nc+1]
                -conj(u12t_f[did*ndim+mu])*r[(did*ngorkov+igork1)*nc]
                -u11t_f[did*ndim+mu] *r[(did*ngorkov+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms next. These run from igorkov=0..3 and 4..7 with slightly different rules for each
      //We can fit it into a single loop by declaring igorkovPP=igorkov+4 instead of looping igorkov=4..7  separately
      //Note that for the igorkov 4..7 loop idirac=igorkov-4, so we don't need to declare idiracPP separately
      //Under dagger, dk4p_f and dk4m_f get swapped and the dirac component flips sign.
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,dk4m_f,dk4p_f:AVX)
      for(int igorkov=0; igorkov<4; igorkov++){
         //the FORTRAN code did it.
         int igork1 = gamin[3*ndirac+igorkov];  
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         phi[(i*ngorkov+igorkov)*nc]+=
            -dk4m_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc]+r[(uid*ngorkov+igork1)*nc])
                  +u12t_f[i*ndim+3]*(r[(uid*ngorkov+igorkov)*nc+1]+r[(uid*ngorkov+igork1)*nc+1]))
            -dk4p_f[did]*(conj(u11t_f[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]-r[(did*ngorkov+igork1)*nc])
                  -u12t_f[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]-r[(did*ngorkov+igork1)*nc+1]));
         phi[(i*ngorkov+igorkov)*nc+1]+=
            -dk4m_f[i]*(-conj(u12t_f[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc]+r[(uid*ngorkov+igork1)*nc])
                  +conj(u11t_f[i*ndim+3])*(r[(uid*ngorkov+igorkov)*nc+1]+r[(uid*ngorkov+igork1)*nc+1]))
            -dk4p_f[did]*(conj(u12t_f[did*ndim+3])*(r[(did*ngorkov+igorkov)*nc]-r[(did*ngorkov+igork1)*nc])
                  +u11t_f[did*ndim+3] *(r[(did*ngorkov+igorkov)*nc+1]-r[(did*ngorkov+igork1)*nc+1]));
 
 
         int igorkovPP=igorkov+4;   //idirac = igorkov; It is a bit redundant but I'll mention it as that's how
         int igork1PP = igork1+4;
         //And the +4 terms. Note that dk4p_f and dk4m_f swap positions compared to the above            
         phi[(i*ngorkov+igorkovPP)*nc]+=-dk4p_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc]+r[(uid*ngorkov+igork1PP)*nc])
               +u12t_f[i*ndim+3]*(r[(uid*ngorkov+igorkovPP)*nc+1]+r[(uid*ngorkov+igork1PP)*nc+1]))
            -dk4m_f[did]*(conj(u11t_f[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]-r[(did*ngorkov+igork1PP)*nc])
                  -u12t_f[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]-r[(did*ngorkov+igork1PP)*nc+1]));
 
         phi[(i*ngorkov+igorkovPP)*nc+1]+=dk4p_f[i]*(conj(u12t_f[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc]+r[(uid*ngorkov+igork1PP)*nc])
               -conj(u11t_f[i*ndim+3])*(r[(uid*ngorkov+igorkovPP)*nc+1]+r[(uid*ngorkov+igork1PP)*nc+1]))
            -dk4m_f[did]*(conj(u12t_f[did*ndim+3])*(r[(did*ngorkov+igorkovPP)*nc]-r[(did*ngorkov+igork1PP)*nc])
                  +u11t_f[did*ndim+3]*(r[(did*ngorkov+igorkovPP)*nc+1]-r[(did*ngorkov+igork1PP)*nc+1]));
 
      }
#endif
   }
#endif
   return 0;
}

References CHalo_swap_all(), Complex_f, kferm, kvol, nc, ndim, ndirac, and ngorkov.

Here is the call graph for this function:

Here is the caller graph for this function:

Hdslash()

int Hdslash	(	Complex *	phi,
		Complex *	r,
		Complex *	u11t,
		Complex *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex *	gamval,
		int *	gamin,
		double *	dk4m,
		double *	dk4p,
		float	akappa )

Evaluates \(\Phi=M r\) in double precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 267 of file matrices.c.

                                                                            {
   /*
    * @brief Evaluates @f(\Phi=M r@f) in double precision
    *
    * @param   phi:     The product
    * @param   r:       The array being acted on by M
    * @param   u11t: First colour trial field
    * @param   u12t: Second colour trial field
    * @param   iu:      Upper halo indices
    * @param   id:      Lower halo indices
    * @param   gamval:  Gamma matrices
    * @param   gamin:   Indices for dirac terms
    * @param   dk4m: 
    * @param   dk4p: 
    * @param   akappa:  Hopping parameter
    *
    * @see ZHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Hdslash";
   //Get the halos in order
#if(nproc>1)
   ZHalo_swap_all(r, 8);
#endif
 
   //Mass term
   //Spacelike term
#ifdef __NVCC__
   cuHdslash(phi,r,u11t,u12t,iu,id,gamval,gamin,dk4m,dk4p,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm2*sizeof(Complex));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t,u12t,gamval:AVX)
         for(int idirac=0; idirac<ndirac; idirac++){
            //FORTRAN had mod((idirac-1),4)+1 to prevent issues with non-zero indexing.
            int igork1 = gamin[mu*ndirac+idirac];
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            phi[(i*ndirac+idirac)*nc]+=-akappa*(u11t[i*ndim+mu]*r[(uid*ndirac+idirac)*nc]+\
                  u12t[i*ndim+mu]*r[(uid*ndirac+idirac)*nc+1]+\
                  conj(u11t[did*ndim+mu])*r[(did*ndirac+idirac)*nc]-\
                  u12t[did*ndim+mu]*r[(did*ndirac+idirac)*nc+1])+\
                                       //Dirac term
                                       gamval[mu*ndirac+idirac]*(u11t[i*ndim+mu]*r[(uid*ndirac+igork1)*nc]+\
                                             u12t[i*ndim+mu]*r[(uid*ndirac+igork1)*nc+1]-\
                                             conj(u11t[did*ndim+mu])*r[(did*ndirac+igork1)*nc]+\
                                             u12t[did*ndim+mu]*r[(did*ndirac+igork1)*nc+1]);
 
            phi[(i*ndirac+idirac)*nc+1]+=-akappa*(-conj(u12t[i*ndim+mu])*r[(uid*ndirac+idirac)*nc]+\
                  conj(u11t[i*ndim+mu])*r[(uid*ndirac+idirac)*nc+1]+\
                  conj(u12t[did*ndim+mu])*r[(did*ndirac+idirac)*nc]+\
                  u11t[did*ndim+mu]*r[(did*ndirac+idirac)*nc+1])+\
                                         //Dirac term
                                         gamval[mu*ndirac+idirac]*(-conj(u12t[i*ndim+mu])*r[(uid*ndirac+igork1)*nc]+\
                                               conj(u11t[i*ndim+mu])*r[(uid*ndirac+igork1)*nc+1]-\
                                               conj(u12t[did*ndim+mu])*r[(did*ndirac+igork1)*nc]-\
                                               u11t[did*ndim+mu]*r[(did*ndirac+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t,u12t,dk4m,dk4p:AVX)
      for(int idirac=0; idirac<ndirac; idirac++){
         int igork1 = gamin[3*ndirac+idirac];
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         //Reminder! gamval was rescaled by kappa when we defined it
         phi[(i*ndirac+idirac)*nc]+=
            -dk4p[i]*(u11t[i*ndim+3]*(r[(uid*ndirac+idirac)*nc]-r[(uid*ndirac+igork1)*nc])
                  +u12t[i*ndim+3]*(r[(uid*ndirac+idirac)*nc+1]-r[(uid*ndirac+igork1)*nc+1]))
            -dk4m[did]*(conj(u11t[did*ndim+3])*(r[(did*ndirac+idirac)*nc]+r[(did*ndirac+igork1)*nc])
                  -u12t[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]+r[(did*ndirac+igork1)*nc+1]));
         phi[(i*ndirac+idirac)*nc+1]+=
            -dk4p[i]*(-conj(u12t[i*ndim+3])*(r[(uid*ndirac+idirac)*nc]-r[(uid*ndirac+igork1)*nc])
                  +conj(u11t[i*ndim+3])*(r[(uid*ndirac+idirac)*nc+1]-r[(uid*ndirac+igork1)*nc+1]))
            -dk4m[did]*(conj(u12t[did*ndim+3])*(r[(did*ndirac+idirac)*nc]+r[(did*ndirac+igork1)*nc])
                  +u11t[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]+r[(did*ndirac+igork1)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References Complex, kferm2, kvol, nc, ndim, ndirac, and ZHalo_swap_all().

Here is the call graph for this function:

Here is the caller graph for this function:

Hdslash_f()

int Hdslash_f	(	Complex_f *	phi,
		Complex_f *	r,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	gamval,
		int *	gamin,
		float *	dk4m,
		float *	dk4p,
		float	akappa )

Evaluates \(\Phi=M r\) in single precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 711 of file matrices.c.

                                                                                  {
   /*
    * @brief Evaluates @f(\Phi=M r@f) in single precision.
    *
    * @param   phi:     The product
    * @param   r:       The array being acted on by M
    * @param   u11t_f:  First colour trial field
    * @param   u12t_f:  Second colour trial field
    * @param   iu:      Upper halo indices
    * @param   id:      Lower halo indices
    * @param   gamval_f:   Gamma matrices
    * @param   gamin:   Indices for dirac terms
    * @param   dk4m_f:  
    * @param   dk4p_f:  
    * @param   akappa:  Hopping parameter
    *
    * @see CHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Hdslash_f";
   //Get the halos in order
#if(nproc>1)
   CHalo_swap_all(r, 8);
#endif
#ifdef __NVCC__
   cuHdslash_f(phi,r,u11t_f,u12t_f,iu,id,gamval_f,gamin,dk4m_f,dk4p_f,akappa,dimGrid,dimBlock);
#else
   //Mass term
   memcpy(phi, r, kferm2*sizeof(Complex_f));
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
      //Spacelike term
#ifndef NO_SPACE
      for(int mu = 0; mu <3; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,gamval_f:AVX)
         for(int idirac=0; idirac<ndirac; idirac++){
            //FORTRAN had mod((idirac-1),4)+1 to prevent issues with non-zero indexing.
            int igork1 = gamin[mu*ndirac+idirac];
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
            phi[(i*ndirac+idirac)*nc]+=-akappa*(u11t_f[i*ndim+mu]*r[(uid*ndirac+idirac)*nc]+\
                  u12t_f[i*ndim+mu]*r[(uid*ndirac+idirac)*nc+1]+\
                  conjf(u11t_f[did*ndim+mu])*r[(did*ndirac+idirac)*nc]-\
                  u12t_f[did*ndim+mu]*r[(did*ndirac+idirac)*nc+1]);
            //Dirac term
            //Reminder! gamval was rescaled by kappa when we defined it
            phi[(i*ndirac+idirac)*nc]+=gamval_f[mu*ndirac+idirac]*(u11t_f[i*ndim+mu]*r[(uid*ndirac+igork1)*nc]+\
                  u12t_f[i*ndim+mu]*r[(uid*ndirac+igork1)*nc+1]-\
                  conjf(u11t_f[did*ndim+mu])*r[(did*ndirac+igork1)*nc]+\
                  u12t_f[did*ndim+mu]*r[(did*ndirac+igork1)*nc+1]);
 
            phi[(i*ndirac+idirac)*nc+1]+=-akappa*(-conjf(u12t_f[i*ndim+mu])*r[(uid*ndirac+idirac)*nc]+\
                  conjf(u11t_f[i*ndim+mu])*r[(uid*ndirac+idirac)*nc+1]+\
                  conjf(u12t_f[did*ndim+mu])*r[(did*ndirac+idirac)*nc]+\
                  u11t_f[did*ndim+mu]*r[(did*ndirac+idirac)*nc+1]);
            //Dirac term
            phi[(i*ndirac+idirac)*nc+1]+=gamval_f[mu*ndirac+idirac]*(-conjf(u12t_f[i*ndim+mu])*r[(uid*ndirac+igork1)*nc]+\
                  conjf(u11t_f[i*ndim+mu])*r[(uid*ndirac+igork1)*nc+1]-\
                  conjf(u12t_f[did*ndim+mu])*r[(did*ndirac+igork1)*nc]-\
                  u11t_f[did*ndim+mu]*r[(did*ndirac+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,dk4m_f,dk4p_f:AVX)
      for(int idirac=0; idirac<ndirac; idirac++){
         int igork1 = gamin[3*ndirac+idirac];
         //Factorising for performance, we get dk4?*(float)u1?*(+/-r_wilson -/+ r_dirac)
         phi[(i*ndirac+idirac)*nc]+=
            -dk4p_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ndirac+idirac)*nc]-r[(uid*ndirac+igork1)*nc])
                  +u12t_f[i*ndim+3]*(r[(uid*ndirac+idirac)*nc+1]-r[(uid*ndirac+igork1)*nc+1]))
            -dk4m_f[did]*(conjf(u11t_f[did*ndim+3])*(r[(did*ndirac+idirac)*nc]+r[(did*ndirac+igork1)*nc])
                  -u12t_f[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]+r[(did*ndirac+igork1)*nc+1]));
 
         phi[(i*ndirac+idirac)*nc+1]+=
            -dk4p_f[i]*(-conjf(u12t_f[i*ndim+3])*(r[(uid*ndirac+idirac)*nc]-r[(uid*ndirac+igork1)*nc])
                  +conjf(u11t_f[i*ndim+3])*(r[(uid*ndirac+idirac)*nc+1]-r[(uid*ndirac+igork1)*nc+1]))
            -dk4m_f[did]*(conjf(u12t_f[did*ndim+3])*(r[(did*ndirac+idirac)*nc]+r[(did*ndirac+igork1)*nc])
                  +u11t_f[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]+r[(did*ndirac+igork1)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References CHalo_swap_all(), Complex_f, kferm2, kvol, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

Hdslashd()

int Hdslashd	(	Complex *	phi,
		Complex *	r,
		Complex *	u11t,
		Complex *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex *	gamval,
		int *	gamin,
		double *	dk4m,
		double *	dk4p,
		float	akappa )

Evaluates \(\Phi=M^\dagger r\) in double precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 358 of file matrices.c.

                                                                            {
   /*
    * @brief Evaluates @f(\Phi=M^\dagger r@f) in double precision
    *
    * @param   phi:     The product
    * @param   r:       The array being acted on by M
    * @param   u11t: First colour trial field
    * @param   u12t: Second colour trial field
    * @param   iu:      Upper halo indices
    * @param   id:      Lower halo indices
    * @param   gamval:  Gamma matrices
    * @param   gamin:   Indices for dirac terms
    * @param   dk4m: 
    * @param   dk4p: 
    * @param   akappa:  Hopping parameter
    *
    * @see ZHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Hdslashd";
   //Get the halos in order. Because C is row major, we need to extract the correct
   //terms for each halo first. Changing the indices was considered but that caused
   //issues with the BLAS routines.
#if(nproc>1)
   ZHalo_swap_all(r, 8);
#endif
 
   //Looks like flipping the array ordering for C has meant a lot
   //of for loops. Sense we're jumping around quite a bit the cache is probably getting refreshed
   //anyways so memory access patterns mightn't be as big of an limiting factor here anyway
 
   //Mass term
#ifdef __NVCC__
   cuHdslashd(phi,r,u11t,u12t,iu,id,gamval,gamin,dk4m,dk4p,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm2*sizeof(Complex));
   //Spacelike term
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#ifndef NO_SPACE
      for(int mu = 0; mu <ndim-1; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t,u12t,gamval:AVX)
         for(int idirac=0; idirac<ndirac; idirac++){
            //FORTRAN had mod((idirac-1),4)+1 to prevent issues with non-zero indexing.
            int igork1 = gamin[mu*ndirac+idirac];
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
 
            //Reminder! gamval was rescaled by kappa when we defined it
            phi[(i*ndirac+idirac)*nc]+=
               -akappa*(u11t[i*ndim+mu]*r[(uid*ndirac+idirac)*nc]
                     +u12t[i*ndim+mu]*r[(uid*ndirac+idirac)*nc+1]
                     +conj(u11t[did*ndim+mu])*r[(did*ndirac+idirac)*nc]
                     -u12t[did*ndim+mu] *r[(did*ndirac+idirac)*nc+1])
               -gamval[mu*ndirac+idirac]*
               (          u11t[i*ndim+mu]*r[(uid*ndirac+igork1)*nc]
                          +u12t[i*ndim+mu]*r[(uid*ndirac+igork1)*nc+1]
                          -conj(u11t[did*ndim+mu])*r[(did*ndirac+igork1)*nc]
                          +u12t[did*ndim+mu] *r[(did*ndirac+igork1)*nc+1]);
 
            phi[(i*ndirac+idirac)*nc+1]+=
               -akappa*(-conj(u12t[i*ndim+mu])*r[(uid*ndirac+idirac)*nc]
                     +conj(u11t[i*ndim+mu])*r[(uid*ndirac+idirac)*nc+1]
                     +conj(u12t[did*ndim+mu])*r[(did*ndirac+idirac)*nc]
                     +u11t[did*ndim+mu] *r[(did*ndirac+idirac)*nc+1])
               -gamval[mu*ndirac+idirac]*
               (-conj(u12t[i*ndim+mu])*r[(uid*ndirac+igork1)*nc]
                +conj(u11t[i*ndim+mu])*r[(uid*ndirac+igork1)*nc+1]
                -conj(u12t[did*ndim+mu])*r[(did*ndirac+igork1)*nc]
                -u11t[did*ndim+mu] *r[(did*ndirac+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t,u12t,dk4m,dk4p:AVX)
      for(int idirac=0; idirac<ndirac; idirac++){
         int igork1 = gamin[3*ndirac+idirac];
         //Factorising for performance, we get dk4?*u1?*(+/-r_wilson -/+ r_dirac)
         //dk4m and dk4p swap under dagger
         phi[(i*ndirac+idirac)*nc]+=
            -dk4m[i]*(u11t[i*ndim+3]*(r[(uid*ndirac+idirac)*nc]+r[(uid*ndirac+igork1)*nc])
                  +u12t[i*ndim+3]*(r[(uid*ndirac+idirac)*nc+1]+r[(uid*ndirac+igork1)*nc+1]))
            -dk4p[did]*(conj(u11t[did*ndim+3])*(r[(did*ndirac+idirac)*nc]-r[(did*ndirac+igork1)*nc])
                  -u12t[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]-r[(did*ndirac+igork1)*nc+1]));
 
         phi[(i*ndirac+idirac)*nc+1]+=
            -dk4m[i]*(-conj(u12t[i*ndim+3])*(r[(uid*ndirac+idirac)*nc]+r[(uid*ndirac+igork1)*nc])
                  +conj(u11t[i*ndim+3])*(r[(uid*ndirac+idirac)*nc+1]+r[(uid*ndirac+igork1)*nc+1]))
            -dk4p[did]*(conj(u12t[did*ndim+3])*(r[(did*ndirac+idirac)*nc]-r[(did*ndirac+igork1)*nc])
                  +u11t[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]-r[(did*ndirac+igork1)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References Complex, kferm2, kvol, nc, ndim, ndirac, and ZHalo_swap_all().

Here is the call graph for this function:

Hdslashd_f()

int Hdslashd_f	(	Complex_f *	phi,
		Complex_f *	r,
		Complex_f *	u11t,
		Complex_f *	u12t,
		unsigned int *	iu,
		unsigned int *	id,
		Complex_f *	gamval,
		int *	gamin,
		float *	dk4m,
		float *	dk4p,
		float	akappa )

Evaluates \(\Phi=M^\dagger r\) in single precision.

Parameters

phi	The product
r	The array being acted on by M
u11t	First colour trial field
u12t	Second colour trial field
iu	Upper halo indices
id	Lower halo indices
gamval	Gamma matrices rescaled by kappa
gamin	Indices for dirac terms
dk4m	\(\left(1+\gamma_0\right)e^{-\mu}\)
dk4p	\(\left(1-\gamma_0\right)e^\mu\)
akappa	Hopping parameter

Returns

Zero on success, integer error code otherwise

Definition at line 802 of file matrices.c.

                                                                                  {
   /*
    * @brief Evaluates @f(\Phi=M^\dagger r@f) in single precision
    *
    * @param   phi:     The product
    * @param   r:       The array being acted on by M
    * @param   u11t_f:  First colour trial field
    * @param   u12t_f:  Second colour trial field
    * @param   iu:      Upper halo indices
    * @param   id:      Lower halo indices
    * @param   gamval_f:   Gamma matrices
    * @param   gamin:   Indices for dirac terms
    * @param   dk4m_f:  
    * @param   dk4p_f:  
    * @param   akappa:  Hopping parameter
    *
    * @see CHalo_swap_all (MPI only)
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Hdslashd_f";
   //Get the halos in order. Because C is row major, we need to extract the correct
   //terms for each halo first. Changing the indices was considered but that caused
   //issues with the BLAS routines.
#if(nproc>1)
   CHalo_swap_all(r, 8);
#endif
 
   //Looks like flipping the array ordering for C has meant a lot
   //of for loops. Sense we're jumping around quite a bit the cache is probably getting refreshed
   //anyways so memory access patterns mightn't be as big of an limiting factor here anyway
 
   //Mass term
#ifdef __NVCC__
   cuHdslashd_f(phi,r,u11t_f,u12t_f,iu,id,gamval_f,gamin,dk4m_f,dk4p_f,akappa,dimGrid,dimBlock);
#else
   memcpy(phi, r, kferm2*sizeof(Complex_f));
   //Spacelike term
#pragma omp parallel for
   for(int i=0;i<kvol;i++){
#ifndef NO_SPACE
      for(int mu = 0; mu <ndim-1; mu++){
         int did=id[mu+ndim*i]; int uid = iu[mu+ndim*i];
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,gamval_f:AVX)
         for(int idirac=0; idirac<ndirac; idirac++){
            //FORTRAN had mod((idirac-1),4)+1 to prevent issues with non-zero indexing.
            int igork1 = gamin[mu*ndirac+idirac];
            //Can manually vectorise with a pragma?
            //Wilson + Dirac term in that order. Definitely easier
            //to read when split into different loops, but should be faster this way
 
            //Reminder! gamval was rescaled by kappa when we defined it
            phi[(i*ndirac+idirac)*nc]+=
               -akappa*(u11t_f[i*ndim+mu]*r[(uid*ndirac+idirac)*nc]
                     +u12t_f[i*ndim+mu]*r[(uid*ndirac+idirac)*nc+1]
                     +conjf(u11t_f[did*ndim+mu])*r[(did*ndirac+idirac)*nc]
                     -u12t_f[did*ndim+mu] *r[(did*ndirac+idirac)*nc+1])
               -gamval_f[mu*ndirac+idirac]*
               (          u11t_f[i*ndim+mu]*r[(uid*ndirac+igork1)*nc]
                          +u12t_f[i*ndim+mu]*r[(uid*ndirac+igork1)*nc+1]
                          -conjf(u11t_f[did*ndim+mu])*r[(did*ndirac+igork1)*nc]
                          +u12t_f[did*ndim+mu] *r[(did*ndirac+igork1)*nc+1]);
 
            phi[(i*ndirac+idirac)*nc+1]+=
               -akappa*(-conjf(u12t_f[i*ndim+mu])*r[(uid*ndirac+idirac)*nc]
                     +conjf(u11t_f[i*ndim+mu])*r[(uid*ndirac+idirac)*nc+1]
                     +conjf(u12t_f[did*ndim+mu])*r[(did*ndirac+idirac)*nc]
                     +u11t_f[did*ndim+mu] *r[(did*ndirac+idirac)*nc+1])
               -gamval_f[mu*ndirac+idirac]*
               (-conjf(u12t_f[i*ndim+mu])*r[(uid*ndirac+igork1)*nc]
                +conjf(u11t_f[i*ndim+mu])*r[(uid*ndirac+igork1)*nc+1]
                -conjf(u12t_f[did*ndim+mu])*r[(did*ndirac+igork1)*nc]
                -u11t_f[did*ndim+mu] *r[(did*ndirac+igork1)*nc+1]);
         }
      }
#endif
      //Timelike terms
      int did=id[3+ndim*i]; int uid = iu[3+ndim*i];
#ifndef NO_TIME
#pragma omp simd aligned(phi,r,u11t_f,u12t_f,gamval_f:AVX)
      for(int idirac=0; idirac<ndirac; idirac++){
         int igork1 = gamin[3*ndirac+idirac];
         //Factorising for performance, we get (float)dk4?*(float)u1?*(+/-r_wilson -/+ r_dirac)
         //(float)dk4m and dk4p_f swap under dagger
         phi[(i*ndirac+idirac)*nc]+=
            -dk4m_f[i]*(u11t_f[i*ndim+3]*(r[(uid*ndirac+idirac)*nc]+r[(uid*ndirac+igork1)*nc])
                  +u12t_f[i*ndim+3]*(r[(uid*ndirac+idirac)*nc+1]+r[(uid*ndirac+igork1)*nc+1]))
            -dk4p_f[did]*(conjf(u11t_f[did*ndim+3])*(r[(did*ndirac+idirac)*nc]-r[(did*ndirac+igork1)*nc])
                  -u12t_f[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]-r[(did*ndirac+igork1)*nc+1]));
 
         phi[(i*ndirac+idirac)*nc+1]+=
            -dk4m_f[i]*(-conjf(u12t_f[i*ndim+3])*(r[(uid*ndirac+idirac)*nc]+r[(uid*ndirac+igork1)*nc])
                  +conjf(u11t_f[i*ndim+3])*(r[(uid*ndirac+idirac)*nc+1]+r[(uid*ndirac+igork1)*nc+1]))
            -dk4p_f[did]*(conjf(u12t_f[did*ndim+3])*(r[(did*ndirac+idirac)*nc]-r[(did*ndirac+igork1)*nc])
                  +u11t_f[did*ndim+3] *(r[(did*ndirac+idirac)*nc+1]-r[(did*ndirac+igork1)*nc+1]));
      }
#endif
   }
#endif
   return 0;
}

References CHalo_swap_all(), Complex_f, kferm2, kvol, nc, ndim, and ndirac.

Here is the call graph for this function:

Here is the caller graph for this function:

Reunitarise()

int Reunitarise ( Complex * u11t, Complex * u12t )

inline

Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1.

If you're looking at the FORTRAN code be careful. There are two header files for the /trial/ header. One with u11 u12 (which was included here originally) and the other with u11t and u12t.

Parameters

u11t,u12t

Trial fields to be reunitarised

Returns

Zero on success, integer error code otherwise

Definition at line 904 of file matrices.c.

                                                    {
   /*
    * @brief Reunitarises u11t and u12t as in conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]=1
    *
    * If you're looking at the FORTRAN code be careful. There are two header files
    * for the /trial/ header. One with u11 u12 (which was included here originally)
    * and the other with u11t and u12t.
    *
    * @see cuReunitarise (CUDA Wrapper)
    *
    * @param u11t, u12t Trial fields to be reunitarised
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Reunitarise";
#ifdef __NVCC__
   cuReunitarise(u11t,u12t,dimGrid,dimBlock);
#else
#pragma omp parallel for simd aligned(u11t,u12t:AVX)
   for(int i=0; i<kvol*ndim; i++){
      //Declaring anorm inside the loop will hopefully let the compiler know it
      //is safe to vectorise aggressively
      double anorm=sqrt(conj(u11t[i])*u11t[i]+conj(u12t[i])*u12t[i]);
      //    Exception handling code. May be faster to leave out as the exit prevents vectorisation.
      //    if(anorm==0){
      //       fprintf(stderr, "Error %i in %s on rank %i: anorm = 0 for μ=%i and i=%i.\nExiting...\n\n",
      //             DIVZERO, funcname, rank, mu, i);
      //       MPI_Finalise();
      //       exit(DIVZERO);
      //    }
      u11t[i]/=anorm;
      u12t[i]/=anorm;
   }
#endif
   return 0;
}

References kvol, and ndim.

Here is the caller graph for this function: