/****************  m_matvec_s.c  (in su3.a) *****************************
*									*
* void mult_su3_mat_vec_sum( su3_matrix *a, su3_vector *b,*c )		*
* su3_matrix times su3_vector multiply and add to another su3_vector 	*
* C  <-  C + A*B 							*
*/
#include "complex.h"
#include "su3.h"

#ifndef FAST
/* su3_matrix times su3_vector multiply and add to another su3_vector */
/* c  <-  A*b+c */
void mult_su3_mat_vec_sum( a,b,c )su3_matrix *a; su3_vector *b,*c;{
register int i,j;
register complex x,y;
    for(i=0;i<3;i++){
	x.real=x.imag=0.0;
	for(j=0;j<3;j++){
	    CMUL( a->e[i][j] , b->c[j] , y )
	    CSUM( x , y );
	}
	c->c[i].real += x.real;
	c->c[i].imag += x.imag;
    }
}

#else
void mult_su3_mat_vec_sum(a,b,c) su3_matrix *a; su3_vector *b,*c; {
int i,j,k;
register float t,ar,ai,br,bi,cr,ci;
    for(i=0;i<3;i++){

	ar=a->e[i][0].real; ai=a->e[i][0].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[i][1].real; ai=a->e[i][1].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[i][2].real; ai=a->e[i][2].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 FAST" */