/**************  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 "complex.h"
#include "su3.h"

#ifndef FAST

void mult_adj_su3_mat_hwvec( su3_matrix *mat,
       half_wilson_vector *src, half_wilson_vector *dest ){
    mult_adj_su3_mat_vec(mat, &(src->h[0]), &(dest->h[0]) );
    mult_adj_su3_mat_vec(mat, &(src->h[1]), &(dest->h[1]) );
}

#else  /* Fast version */

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 float a0r,a0i,a1r,a1i,a2r,a2i;
  register float b0r,b0i,b1r,b1i,b2r,b2i;
#endif

/*    mult_adj_su3_mat_vec(mat, &(src->h[0]), &(dest->h[0]) ); */
  
  a0r=mat->e[0][0].real;   a0i=mat->e[0][0].imag;
  b0r=src->h[0].c[0].real; b0i=src->h[0].c[0].imag;
  a1r=mat->e[1][0].real;   a1i=mat->e[1][0].imag;
  b1r=src->h[0].c[1].real; b1i=src->h[0].c[1].imag;
  a2r=mat->e[2][0].real;   a2i=mat->e[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->e[0][1].real;   a0i=mat->e[0][1].imag;
  b0r=src->h[0].c[0].real; b0i=src->h[0].c[0].imag;  
  a1r=mat->e[1][1].real;   a1i=mat->e[1][1].imag;
  b1r=src->h[0].c[1].real; b1i=src->h[0].c[1].imag;
  a2r=mat->e[2][1].real;   a2i=mat->e[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->e[0][2].real;   a0i=mat->e[0][2].imag;
  b0r=src->h[0].c[0].real; b0i=src->h[0].c[0].imag;
  a1r=mat->e[1][2].real;   a1i=mat->e[1][2].imag;
  b1r=src->h[0].c[1].real; b1i=src->h[0].c[1].imag;
  a2r=mat->e[2][2].real;   a2i=mat->e[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(mat, &(src->h[1]), &(dest->h[1]) ); */

  a0r=mat->e[0][0].real;   a0i=mat->e[0][0].imag;
  b0r=src->h[1].c[0].real; b0i=src->h[1].c[0].imag;
  a1r=mat->e[1][0].real;   a1i=mat->e[1][0].imag;
  b1r=src->h[1].c[1].real; b1i=src->h[1].c[1].imag;
  a2r=mat->e[2][0].real;   a2i=mat->e[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->e[0][1].real;   a0i=mat->e[0][1].imag;
  b0r=src->h[1].c[0].real; b0i=src->h[1].c[0].imag;  
  a1r=mat->e[1][1].real;   a1i=mat->e[1][1].imag;
  b1r=src->h[1].c[1].real; b1i=src->h[1].c[1].imag;
  a2r=mat->e[2][1].real;   a2i=mat->e[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->e[0][2].real;   a0i=mat->e[0][2].imag;
  b0r=src->h[1].c[0].real; b0i=src->h[1].c[0].imag;
  a1r=mat->e[1][2].real;   a1i=mat->e[1][2].imag;
  b1r=src->h[1].c[1].real; b1i=src->h[1].c[1].imag;
  a2r=mat->e[2][2].real;   a2i=mat->e[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;

}
#endif /* "ifndef FAST" */