/*****************  m_amv_4dir.c  (in su3.a) *****************************
*									*
*  void mult_adj_su3_mat_vec_4dir( su3_matrix *mat,			*
*  su3_vector *src, su3_vector *dest )					*
*  Multiply an su3_vector by an array of four adjoint su3_matrices,	*
*  result in an array of four su3_vectors.				*
*  dest[i]  <-  A_adjoint[i] * src					*
*/
#include "complex.h"
#include "su3.h"

#ifndef FAST
void mult_adj_su3_mat_vec_4dir( mat, src, dest )
su3_matrix *mat;
su3_vector *src, *dest;
{
    mult_adj_su3_mat_vec( mat+0, src, dest+0 );
    mult_adj_su3_mat_vec( mat+1, src, dest+1 );
    mult_adj_su3_mat_vec( mat+2, src, dest+2 );
    mult_adj_su3_mat_vec( mat+3, src, dest+3 );
}

#else
/* Fast code, with subroutines inlined */
#ifdef _IBMR2 /* IBM RS6000 version */
void mult_adj_su3_mat_vec_4dir( mat, src, dest )
su3_matrix *mat;
su3_vector *src, *dest;
{
  register int n;
  register double a0r,a0i,a1r,a1i,a2r,a2i;
  register double b0r,b0i,b1r,b1i,b2r,b2i;
  register su3_matrix *a;
  register su3_vector *b,*c;

  a = mat; c = dest ; b = src;
  for(n=0;n<4;n++,a++,c++){

  a0r=a->e[0][0].real; a0i=a->e[0][0].imag;
  b0r=b->c[0].real;    b0i=b->c[0].imag;
  a1r=a->e[1][0].real; a1i=a->e[1][0].imag;
  b1r=b->c[1].real;    b1i=b->c[1].imag;
  a2r=a->e[2][0].real; a2i=a->e[2][0].imag;
  b2r=b->c[2].real;    b2i=b->c[2].imag;

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

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

  a0r=a->e[0][2].real; a0i=a->e[0][2].imag;
  b0r=b->c[0].real;    b0i=b->c[0].imag;
  a1r=a->e[1][2].real; a1i=a->e[1][2].imag;
  b1r=b->c[1].real;    b1i=b->c[1].imag;
  a2r=a->e[2][2].real; a2i=a->e[2][2].imag;
  b2r=b->c[2].real;    b2i=b->c[2].imag;

  c->c[2].real = a0r*b0r + a0i*b0i + a1r*b1r + a1i*b1i + a2r*b2r + a2i*b2i;
  c->c[2].imag = a0r*b0i - a0i*b0r + a1r*b1i - a1i*b1r + a2r*b2i - a2i*b2r;
  }
}
#else
void mult_adj_su3_mat_vec_4dir( mat, src, dest )
su3_matrix *mat;
su3_vector *src, *dest;
{
int i,n;
register float t,ar,ai,br,bi,cr,ci;
register su3_matrix *a;
register su3_vector *b,*c;
  a = mat; c = dest ; b = src;
  for(n=0;n<4;n++,a++,c++){
    for(i=0;i<3;i++){
	ar=a->e[0][i].real; ai=a->e[0][i].imag;
	br=b->c[0].real; bi=b->c[0].imag;
	cr=ar*br; t=ai*bi; cr += t;
	ci=ar*bi; t=ai*br; ci -= t;

	ar=a->e[1][i].real; ai=a->e[1][i].imag;
	br=b->c[1].real; bi=b->c[1].imag;
	t=ar*br; cr += t; t=ai*bi; cr += t;
	t=ar*bi; ci += t; t=ai*br; ci -= t;

	ar=a->e[2][i].real; ai=a->e[2][i].imag;
	br=b->c[2].real; bi=b->c[2].imag;
	t=ar*br; cr += t; t=ai*bi; cr += t;
	t=ar*bi; ci += t; t=ai*br; ci -= t;

	c->c[i].real=cr;
	c->c[i].imag=ci;
    }
  }
}
#endif	/* End of "#ifdef _IBMR2" */
#endif	/* End of "#ifndef FAST" */