/*
* mri_exp_mult_template.c
* Template for mri_exp_mult_* to support both double and float.
* Sangwoo Lee, The University of Michigan, Jun, 2004
* 2004-6-21 modified by JF
*/

/*
in:
	A	[N,L]	complex matrix
	u	[N,1]	vector
	v	[M,1]	vector
			one (and only one) of u and v must be complex!
out:
	D	[L,M]	complex vector, D = A' * exp(-u * v.')
			D_lm = sum_n A_nl^* exp(-u_n v_m)
*/

/*
* static bool mri_exp_mult_*
*/
(
dtype *D_r,		/* [L,M] */
dtype *D_i,
double *work_r,		/* [N,1] work space */
double *work_i,
const dtype *A_r,	/* [N,L] */
const dtype *A_i,
const dtype *un_r,	/* [N,1] */
const dtype *un_i,
const dtype *vm_r,	/* [M,1] */
const dtype *vm_i,
const int LL,
const int MM,
const int NN)
{
	int ll, mm, nn;

	for (mm=0; mm < MM; ++mm) {
		const dtype *Ar = A_r;
		const dtype *Ai = A_i;

		/* real and imag part of exp(-u(:) v(i)) */
		if (!vm_i) { /* vm is real */
			const double nvmr = -vm_r[mm];
			for (nn=0; nn < NN; ++nn) {
				const double ee = exp(nvmr * un_r[nn]);
				work_r[nn] = ee * cos(nvmr * un_i[nn]);
				work_i[nn] = ee * sin(nvmr * un_i[nn]);
			}
		}

		else if (!un_i) { /* un is real */
			const double nvmr = -vm_r[mm];
			const double nvmi = -vm_i[mm];
			for (nn=0; nn < NN; ++nn) {
				const double ee = exp(nvmr * un_r[nn]);
				work_r[nn] = ee * cos(nvmi * un_r[nn]);
				work_i[nn] = ee * sin(nvmi * un_r[nn]);
			}
		}

		else
			Fail("bug")

		for (ll=0; ll < LL; ++ll) {
			const double *wr = work_r;
			const double *wi = work_i;
			double sumr = 0., sumi = 0.;

			for (nn=0; nn < NN; ++nn, ++Ar, ++Ai, ++wr, ++wi) {
				sumr += *Ar * *wr + *Ai * *wi;
				sumi += *Ar * *wi - *Ai * *wr;
			}

			*D_r++ = sumr;
			*D_i++ = sumi;
		}
	}

	return 1;
}