/**************  m_amat_hwvec.c  (in su3.a) **********************
*									*
*  void mult_adj_su3_mat_hwvec( su3_matrix *mat,			*
*	half_wilson_vector *src,*dest )					*
*  multiply a Wilson half-vector by the adjoint of a matrix		*
*/
#include "../include/config.h"
#include "../include/complex.h"
#include "../include/su3.h"
#include "../include/macros.h"

extern __targetConst__ int t_sites_on_node;

__targetHost__ __target__ void mult_adj_su3_mat_hwvec( su3_matrix * mat, half_wilson_vector * src, half_wilson_vector * dest )
{

#ifdef NATIVEDOUBLE
    register double a0r, a0i, a1r, a1i, a2r, a2i;
    register double b0r, b0i, b1r, b1i, b2r, b2i;
#else
    register double a0r, a0i, a1r, a1i, a2r, a2i;
    register double b0r, b0i, b1r, b1i, b2r, b2i;
#endif

    /*    mult_adj_su3_mat_vec_KE(mat, &(src->h[0]), &(dest->h[0]) ); */

    a0r = mat->ROWCOL( 0, 0 ).real;
    a0i = mat->ROWCOL( 0, 0 ).imag;
    b0r = src->h[0].c[0].real;
    b0i = src->h[0].c[0].imag;
    a1r = mat->ROWCOL( 1, 0 ).real;
    a1i = mat->ROWCOL( 1, 0 ).imag;
    b1r = src->h[0].c[1].real;
    b1i = src->h[0].c[1].imag;
    a2r = mat->ROWCOL( 2, 0 ).real;
    a2i = mat->ROWCOL( 2, 0 ).imag;
    b2r = src->h[0].c[2].real;
    b2i = src->h[0].c[2].imag;

    dest->h[0].c[0].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;



    dest->h[0].c[0].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;

    a0r = mat->ROWCOL( 0, 1 ).real;
    a0i = mat->ROWCOL( 0, 1 ).imag;
    b0r = src->h[0].c[0].real;
    b0i = src->h[0].c[0].imag;
    a1r = mat->ROWCOL( 1, 1 ).real;
    a1i = mat->ROWCOL( 1, 1 ).imag;
    b1r = src->h[0].c[1].real;
    b1i = src->h[0].c[1].imag;
    a2r = mat->ROWCOL( 2, 1 ).real;
    a2i = mat->ROWCOL( 2, 1 ).imag;
    b2r = src->h[0].c[2].real;
    b2i = src->h[0].c[2].imag;

    dest->h[0].c[1].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;
    dest->h[0].c[1].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;

    a0r = mat->ROWCOL( 0, 2 ).real;
    a0i = mat->ROWCOL( 0, 2 ).imag;
    b0r = src->h[0].c[0].real;
    b0i = src->h[0].c[0].imag;
    a1r = mat->ROWCOL( 1, 2 ).real;
    a1i = mat->ROWCOL( 1, 2 ).imag;
    b1r = src->h[0].c[1].real;
    b1i = src->h[0].c[1].imag;
    a2r = mat->ROWCOL( 2, 2 ).real;
    a2i = mat->ROWCOL( 2, 2 ).imag;
    b2r = src->h[0].c[2].real;
    b2i = src->h[0].c[2].imag;

    dest->h[0].c[2].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;
    dest->h[0].c[2].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;


    /*    mult_adj_su3_mat_vec_KE(mat, &(src->h[1]), &(dest->h[1]) ); */

    a0r = mat->ROWCOL( 0, 0 ).real;
    a0i = mat->ROWCOL( 0, 0 ).imag;
    b0r = src->h[1].c[0].real;
    b0i = src->h[1].c[0].imag;
    a1r = mat->ROWCOL( 1, 0 ).real;
    a1i = mat->ROWCOL( 1, 0 ).imag;
    b1r = src->h[1].c[1].real;
    b1i = src->h[1].c[1].imag;
    a2r = mat->ROWCOL( 2, 0 ).real;
    a2i = mat->ROWCOL( 2, 0 ).imag;
    b2r = src->h[1].c[2].real;
    b2i = src->h[1].c[2].imag;

    dest->h[1].c[0].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;
    dest->h[1].c[0].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;

    a0r = mat->ROWCOL( 0, 1 ).real;
    a0i = mat->ROWCOL( 0, 1 ).imag;
    b0r = src->h[1].c[0].real;
    b0i = src->h[1].c[0].imag;
    a1r = mat->ROWCOL( 1, 1 ).real;
    a1i = mat->ROWCOL( 1, 1 ).imag;
    b1r = src->h[1].c[1].real;
    b1i = src->h[1].c[1].imag;
    a2r = mat->ROWCOL( 2, 1 ).real;
    a2i = mat->ROWCOL( 2, 1 ).imag;
    b2r = src->h[1].c[2].real;
    b2i = src->h[1].c[2].imag;

    dest->h[1].c[1].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;
    dest->h[1].c[1].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;

    a0r = mat->ROWCOL( 0, 2 ).real;
    a0i = mat->ROWCOL( 0, 2 ).imag;
    b0r = src->h[1].c[0].real;
    b0i = src->h[1].c[0].imag;
    a1r = mat->ROWCOL( 1, 2 ).real;
    a1i = mat->ROWCOL( 1, 2 ).imag;
    b1r = src->h[1].c[1].real;
    b1i = src->h[1].c[1].imag;
    a2r = mat->ROWCOL( 2, 2 ).real;
    a2i = mat->ROWCOL( 2, 2 ).imag;
    b2r = src->h[1].c[2].real;
    b2i = src->h[1].c[2].imag;

    dest->h[1].c[2].real = a0r * b0r + a0i * b0i + a1r * b1r + a1i * b1i + a2r * b2r + a2i * b2i;
    dest->h[1].c[2].imag = a0r * b0i - a0i * b0r + a1r * b1i - a1i * b1r + a2r * b2i - a2i * b2r;

}


//version of the above ported to targetDP

__target__ void mult_adj_su3_mat_hwvec_tdp(  double* matstart, int idirmat,   double* srcstart,   double* deststart, int isite )
{

  int i,j;
  int iv=0;

    /*    mult_adj_su3_mat_vec_KE(mat, &(src->h[0]), &(dest->h[0]) ); */


  // load matrix into temporary structure
  double matloc[3][3][2][VVL];
  for(i=0;i<3;i++){
    for(j=0;j<3;j++){
      __targetILP__(iv) matloc[i][j][REPART][iv]=matstart[SU3MI(isite+iv,i,j,idirmat,REPART)];
      __targetILP__(iv) matloc[i][j][IMPART][iv]=matstart[SU3MI(isite+iv,i,j,idirmat,IMPART)];
    }
  }


  // perform computation
  __targetILP__(iv) deststart[HWVI(isite+iv,0,0,REPART)] = 
    matloc[0][0][REPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[0][0][IMPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    +  matloc[1][0][REPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +   matloc[1][0][IMPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    +   matloc[2][0][REPART][iv] * srcstart[HWVLI(iv,2,0,REPART)]
    +  matloc[2][0][IMPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)];
  
					
  
  __targetILP__(iv) deststart[HWVI(isite+iv,0,0,IMPART)] = 
    matloc[0][0][REPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    -  matloc[0][0][IMPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[1][0][REPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    -   matloc[1][0][IMPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +  matloc[2][0][REPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)]
    -  matloc[2][0][IMPART][iv] * srcstart[HWVLI(iv,2,0,REPART)];
  
  
  
  __targetILP__(iv) deststart[HWVI(isite+iv,1,0,REPART)] = 
    matloc[0][1][REPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[0][1][IMPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    +  matloc[1][1][REPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +   matloc[1][1][IMPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    +   matloc[2][1][REPART][iv] * srcstart[HWVLI(iv,2,0,REPART)]
    +  matloc[2][1][IMPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)];
  
  __targetILP__(iv) deststart[HWVI(isite+iv,1,0,IMPART)] = 
    matloc[0][1][REPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    -  matloc[0][1][IMPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[1][1][REPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    -   matloc[1][1][IMPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +  matloc[2][1][REPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)]
    -  matloc[2][1][IMPART][iv] * srcstart[HWVLI(iv,2,0,REPART)];
  
  
  __targetILP__(iv) deststart[HWVI(isite+iv,2,0,REPART)] = 
    matloc[0][2][REPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[0][2][IMPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    +  matloc[1][2][REPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +   matloc[1][2][IMPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    +   matloc[2][2][REPART][iv] * srcstart[HWVLI(iv,2,0,REPART)]
    +  matloc[2][2][IMPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)];
  
  __targetILP__(iv) deststart[HWVI(isite+iv,2,0,IMPART)] = 
    matloc[0][2][REPART][iv] * srcstart[HWVLI(iv,0,0,IMPART)]
    -  matloc[0][2][IMPART][iv] * srcstart[HWVLI(iv,0,0,REPART)]
    +  matloc[1][2][REPART][iv] *  srcstart[HWVLI(iv,1,0,IMPART)]
    -   matloc[1][2][IMPART][iv] * srcstart[HWVLI(iv,1,0,REPART)]
    +  matloc[2][2][REPART][iv] * srcstart[HWVLI(iv,2,0,IMPART)]
    -  matloc[2][2][IMPART][iv] * srcstart[HWVLI(iv,2,0,REPART)];
  
  
  
  __targetILP__(iv) deststart[HWVI(isite+iv,0,1,REPART)] = 
    matloc[0][0][REPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[0][0][IMPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    +  matloc[1][0][REPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +   matloc[1][0][IMPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    +   matloc[2][0][REPART][iv]* srcstart[HWVLI(iv,2,1,REPART)]
    +  matloc[2][0][IMPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)];
  
  __targetILP__(iv) deststart[HWVI(isite+iv,0,1,IMPART)] = 
    matloc[0][0][REPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    -  matloc[0][0][IMPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[1][0][REPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    -   matloc[1][0][IMPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +  matloc[2][0][REPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)]
    -  matloc[2][0][IMPART][iv] * srcstart[HWVLI(iv,2,1,REPART)];
  
  
  __targetILP__(iv) deststart[HWVI(isite+iv,1,1,REPART)] = 
    matloc[0][1][REPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[0][1][IMPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    +  matloc[1][1][REPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +   matloc[1][1][IMPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    +   matloc[2][1][REPART][iv] * srcstart[HWVLI(iv,2,1,REPART)]
    +  matloc[2][1][IMPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)];
  
  __targetILP__(iv) deststart[HWVI(isite+iv,1,1,IMPART)] = 
    matloc[0][1][REPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    -  matloc[0][1][IMPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[1][1][REPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    -   matloc[1][1][IMPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +  matloc[2][1][REPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)]
    -  matloc[2][1][IMPART][iv] * srcstart[HWVLI(iv,2,1,REPART)];
  
  
  __targetILP__(iv) deststart[HWVI(isite+iv,2,1,REPART)] = 
    matloc[0][2][REPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[0][2][IMPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    +  matloc[1][2][REPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +   matloc[1][2][IMPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    +   matloc[2][2][REPART][iv] * srcstart[HWVLI(iv,2,1,REPART)]
    +  matloc[2][2][IMPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)];
  
  __targetILP__(iv) deststart[HWVI(isite+iv,2,1,IMPART)] = 
    matloc[0][2][REPART][iv] * srcstart[HWVLI(iv,0,1,IMPART)]
    -  matloc[0][2][IMPART][iv] * srcstart[HWVLI(iv,0,1,REPART)]
    +  matloc[1][2][REPART][iv] *  srcstart[HWVLI(iv,1,1,IMPART)]
    -   matloc[1][2][IMPART][iv] * srcstart[HWVLI(iv,1,1,REPART)]
    +  matloc[2][2][REPART][iv] * srcstart[HWVLI(iv,2,1,IMPART)]
    -  matloc[2][2][IMPART][iv] * srcstart[HWVLI(iv,2,1,REPART)];
  
}


