/* ZV.c */
#include "../Utilities.h"
/*--------------------------------------------------------------------*/
/*
--------------------------------------------------------
purpose -- to return the absolute value of (areal,aimag)
created -- 98jan24, cca
--------------------------------------------------------
*/
double
Zabs (
double real,
double imag
) {
double abs, val ;
if ( real == 0.0 ) {
abs = fabs(imag) ;
} else if ( imag == 0.0 ) {
abs = fabs(real) ;
} else if ( real >= imag ) {
val = imag/real ;
abs = fabs(real)*sqrt(1 + val*val) ;
} else {
val = real/imag ;
abs = fabs(imag)*sqrt(1 + val*val) ;
}
return(abs) ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------
purpose -- given (areal,aimag),
compute (breal,bimag) = 1/(areal,aimag),
put breal into *pbreal
put bimag into *pbimag
created -- 98jan23, cca
---------------------------------------------------
*/
int
Zrecip (
double areal,
double aimag,
double *pbreal,
double *pbimag
) {
double bimag, breal, fac ;
if ( areal == 0.0 && aimag == 0.0 ) {
return(0) ;
}
if ( fabs(areal) >= fabs(aimag) ) {
fac = aimag/areal ;
breal = 1./(areal + fac*aimag) ;
bimag = -fac * breal ;
} else {
fac = areal/aimag ;
bimag = -1./(aimag + fac*areal) ;
breal = -fac * bimag ;
}
*pbreal = breal ;
*pbimag = bimag ;
return(1) ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------------------------
given [ (areal,aimag) (breal,bimag) ]
[ (creal,cimag) (dreal,dimag) ]
compute [ (ereal,eimag) (freal,fimag) ]
[ (greal,gimag) (hreal,himag) ]
where
I = [ (areal,aimag) (breal,bimag) ] * [ (ereal,eimag) (freal,fimag) ]
[ (creal,cimag) (dreal,dimag) ] [ (greal,gimag) (hreal,himag) ]
note, any of {pereal, peimag, pfreal, pfimag, pgreal, pgimag,
phreal, phimag} can be NULL. if not NULL, their value is filled.
this is useful when the input matrix is symmetric or hermitian.
created -- 98jan23, cca
---------------------------------------------------------------------
*/
int
Zrecip2 (
double areal,
double aimag,
double breal,
double bimag,
double creal,
double cimag,
double dreal,
double dimag,
double *pereal,
double *peimag,
double *pfreal,
double *pfimag,
double *pgreal,
double *pgimag,
double *phreal,
double *phimag
) {
double xreal, ximag, yreal, yimag ;
/*
-------------------------------
step one, compute x = a*d - b*c
-------------------------------
*/
xreal = areal*dreal - aimag*dimag - breal*creal + bimag*cimag ;
ximag = areal*dimag + aimag*dreal - breal*cimag - bimag*creal ;
/*
-------------------------
step two, compute y = 1/x
-------------------------
*/
Zrecip(xreal, ximag, &yreal, &yimag) ;
/*
-----------------------
step three,
[ e f ] = y * [ d -b ]
[ g h ] [ -c a ]
-----------------------
*/
if ( pereal != NULL ) {
*pereal = dreal*yreal - dimag*yimag ;
}
if ( peimag != NULL ) {
*peimag = dreal*yimag + dimag*yreal ;
}
if ( pfreal != NULL ) {
*pfreal = -breal*yreal + bimag*yimag ;
}
if ( pfimag != NULL ) {
*pfimag = -breal*yimag - bimag*yreal ;
}
if ( pgreal != NULL ) {
*pgreal = -creal*yreal + cimag*yimag ;
}
if ( pgimag != NULL ) {
*pgimag = -creal*yimag - cimag*yreal ;
}
if ( phreal != NULL ) {
*phreal = areal*yreal - aimag*yimag ;
}
if ( phimag != NULL ) {
*phimag = areal*yimag + aimag*yreal ;
}
return(1) ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------
purpose -- to allocate, initialize and
return a complex vector x[]
n -- length of the complex vector (2*n double's)
x[ii] = (real, imag) for 0 <= ii < n
created -- 98jan23, cca
------------------------------------------------
*/
double *
ZVinit (
int n,
double real,
double imag
) {
double *x ;
int ii, jj ;
if ( n <= 0 ) {
fprintf(stderr, "\n fatal error in ZVinit(%d,%f,%f)"
"\n bad input\n", n, real, imag) ;
exit(-1) ;
}
ALLOCATE(x, double, 2*n) ;
for ( ii = jj = 0 ; ii < n ; ii++, jj += 2 ) {
x[jj] = real ;
x[jj+1] = imag ;
}
return(x) ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------------
purpose -- to perform a complex dot product
(*prdot,*pidot) = y^T x
where y and x are complex
created -- 98apr15, cca
-------------------------------------------
*/
void
ZVdotU (
int size,
double y[],
double x[],
double *prdot,
double *pidot
) {
double isum, rsum, ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL
|| prdot == NULL || pidot == NULL ) {
fprintf(stderr, "\n fatal error in ZVdotU(%d,%p,%p,%p,%p)"
"\n bad input\n", size, y, x, prdot, pidot) ;
exit(-1) ;
}
isum = rsum = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[jj] ;
yimag = y[jj+1] ;
rsum += xreal*yreal - ximag*yimag ;
isum += xreal*yimag + ximag*yreal ;
}
*prdot = rsum ;
*pidot = isum ;
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------
purpose -- to perform a conjugated complex dot product
(*prdot,*pidot) = (conjugate(y))^T x = y^H x
where y and x are complex
created -- 98apr15, cca
------------------------------------------------------
*/
void
ZVdotC (
int size,
double y[],
double x[],
double *prdot,
double *pidot
) {
double isum, rsum, ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL
|| prdot == NULL || pidot == NULL ) {
fprintf(stderr, "\n fatal error in ZVdotC(%d,%p,%p,%p,%p)"
"\n bad input\n", size, y, x, prdot, pidot) ;
exit(-1) ;
}
isum = rsum = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[jj] ;
yimag = y[jj+1] ;
rsum += xreal*yreal + ximag*yimag ;
isum += - xreal*yimag + ximag*yreal ;
}
*prdot = rsum ;
*pidot = isum ;
return ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------------------
purpose -- to perform a indexed complex dot product
(*prdot,*pidot) = sum_k y[index[k]]*x[k]
where y and x are complex
created -- 98apr15, cca
---------------------------------------------------
*/
void
ZVdotiU (
int size,
double y[],
int index[],
double x[],
double *prdot,
double *pidot
) {
double isum, rsum, ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL || index == NULL || x == NULL
|| prdot == NULL || pidot == NULL ) {
fprintf(stderr, "\n fatal error in ZVdotiU(%d,%p,%p,%p,%p,%p)"
"\n bad input\n", size, y, index, x, prdot, pidot) ;
exit(-1) ;
}
isum = rsum = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
/*
fprintf(stdout,
"\n %% ii = %d, jj = %d, kk = %d", ii, jj, index[ii]) ;
fflush(stdout) ;
*/
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[2*index[ii]] ;
yimag = y[2*index[ii]+1] ;
rsum += xreal*yreal - ximag*yimag ;
isum += xreal*yimag + ximag*yreal ;
}
*prdot = rsum ;
*pidot = isum ;
return ; }
/*--------------------------------------------------------------------*/
/*
-------------------------------------------------------------
purpose -- to perform a indexed conjugate complex dot product
(*prdot,*pidot) = sum_k conjugate(y[index[k]])*x[k]
where y and x are complex
created -- 98apr15, cca
-------------------------------------------------------------
*/
void
ZVdotiC (
int size,
double y[],
int index[],
double x[],
double *prdot,
double *pidot
) {
double isum, rsum, ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL || index == NULL || x == NULL
|| prdot == NULL || pidot == NULL ) {
fprintf(stderr, "\n fatal error in ZVdotiU(%d,%p,%p,%p,%p,%p)"
"\n bad input\n", size, y, index, x, prdot, pidot) ;
exit(-1) ;
}
isum = rsum = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[2*index[ii]] ;
yimag = y[2*index[ii]+1] ;
rsum += xreal*yreal + ximag*yimag ;
isum += - xreal*yimag + ximag*yreal ;
}
*prdot = rsum ;
*pidot = isum ;
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------
purpose -- to perform a complex axpy
y := y + (areal, aimag) * x
where y and x are complex
created -- 98jan22, cca
------------------------------------
*/
void
ZVaxpy (
int size,
double y[],
double areal,
double aimag,
double x[]
) {
double ximag, xreal ;
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVaxpy(%d,%p,%f,%f,%p)"
"\n bad input\n", size, y, areal, aimag, x) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
/*
fprintf(stdout, "\n ii = %d, xreal = %20.12e, ximag = %20.12e",
ii, xreal, ximag) ;
fprintf(stdout, "\n before yreal = %20.12e, yimag = %20.12e",
y[jj], y[jj+1]) ;
*/
y[jj] += areal*xreal - aimag*ximag ;
y[jj+1] += areal*ximag + aimag*xreal ;
/*
fprintf(stdout, "\n after yreal = %20.12e, yimag = %20.12e",
y[jj], y[jj+1]) ;
*/
}
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------
purpose -- to scale a double complex vector by (xreal,ximag)
y := y * (areal, aimag)
created -- 98jan22, cca
------------------------------------------------------------
*/
void
ZVscale (
int size,
double y[],
double areal,
double aimag
) {
double yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL ) {
fprintf(stderr, "\n fatal error in ZVscale(%d,%p,%f,%f)"
"\n bad input\n", size, y, areal, aimag) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
yreal = y[jj] ;
yimag = y[jj+1] ;
y[jj] = areal*yreal - aimag*yimag ;
y[jj+1] = areal*yimag + aimag*yreal ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------
purpose --- to print a complex vector to a file
created -- 98jan23, cca
-----------------------------------------------
*/
void
ZVfprintf (
FILE *fp,
int size,
double y[]
) {
int ii, jj ;
if ( size < 0 || y == NULL ) {
fprintf(stderr, "\n fatal error in ZVfprintf(%p,%d,%p)"
"\n bad input\n", fp, size, y) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
/*
if ( ii % 2 == 0 ) {
fprintf(fp, "\n") ;
}
*/
fprintf(fp, "\n < %12.4e, %12.4e >", y[jj], y[jj+1]) ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------
return the minimum absolute value
of the entries in a complex vector
created -- 98jan23, cca
----------------------------------
*/
double
ZVminabs (
int size,
double x[]
) {
double abs, imag, minabs, real, val ;
int ii, jj ;
if ( size < 0 || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVminabs(%d,%p)"
"\n bad input\n", size, x) ;
exit(-1) ;
}
minabs = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
real = fabs(x[jj]) ;
imag = fabs(x[jj+1]) ;
/*
fprintf(stdout, "\n i = %d, <%12.4e, %12.4e>", ii, real, imag) ;
*/
if ( real == 0.0 ) {
abs = imag ;
} else if ( imag == 0.0 ) {
abs = real ;
} else if ( real >= imag ) {
val = imag/real ;
abs = real*sqrt(1 + val*val) ;
} else {
val = real/imag ;
abs = imag*sqrt(1 + val*val) ;
}
/*
fprintf(stdout, " abs = %12.4e", abs) ;
*/
if ( ii == 0 || minabs > abs ) {
minabs = abs ;
}
}
/*
fprintf(stdout, "\n minabs = %12.4e", minabs) ;
*/
return(minabs) ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------
return the maximum absolute value
of the entries in a complex vector
created -- 98jan23, cca
----------------------------------
*/
double
ZVmaxabs (
int size,
double x[]
) {
double abs, imag, maxabs, real, val ;
int ii, jj ;
if ( size < 0 || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVmaxabs(%d,%p)"
"\n bad input\n", size, x) ;
exit(-1) ;
}
maxabs = 0.0 ;
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
real = fabs(x[jj]) ;
imag = fabs(x[jj+1]) ;
/*
fprintf(stdout, "\n i = %d, <%12.4e, %12.4e>", ii, real, imag) ;
*/
if ( real == 0.0 ) {
abs = imag ;
} else if ( imag == 0.0 ) {
abs = real ;
} else if ( real >= imag ) {
val = imag/real ;
abs = real*sqrt(1 + val*val) ;
} else {
val = real/imag ;
abs = imag*sqrt(1 + val*val) ;
}
/*
fprintf(stdout, " abs = %12.4e", abs) ;
*/
if ( ii == 0 || maxabs < abs ) {
maxabs = abs ;
}
}
/*
fprintf(stdout, "\n maxabs = %12.4e", maxabs) ;
*/
return(maxabs) ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------
copy a complex vector into another
y[] := x[]
created -- 98jan23, cca
----------------------------------
*/
void
ZVcopy (
int size,
double y[],
double x[]
) {
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVcopy(%d,%p,%p)"
"\n bad input\n", size, y, x) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
y[jj] = x[jj] ;
y[jj+1] = x[jj+1] ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
--------------------------------------
subtract a complex vector from another
y[] := y[] - x[]
created -- 98may25, cca
--------------------------------------
*/
void
ZVsub (
int size,
double y[],
double x[]
) {
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVsub(%d,%p,%p)"
"\n bad input\n", size, y, x) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
y[jj] -= x[jj] ;
y[jj+1] -= x[jj+1] ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------------------------------
purpose -- to perform a complex axpy with two vectors
z := z + (areal, aimag)*x + (breal, bimag)*y
where y and x are complex
created -- 98jan23, cca
----------------------------------------------------
*/
void
ZVaxpy2 (
int size,
double z[],
double areal,
double aimag,
double x[],
double breal,
double bimag,
double y[]
) {
double ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || y == NULL || x == NULL ) {
fprintf(stderr, "\n fatal error in ZVaxpy(%d,%p,%f,%f,%p)"
"\n bad input\n", size, y, areal, aimag, x) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[jj] ;
yimag = y[jj+1] ;
/*
fprintf(stdout, "\n ii = %d, xreal = %20.12e, ximag = %20.12e",
ii, xreal, ximag) ;
fprintf(stdout, "\n yreal = %20.12e, yimag = %20.12e",
y[jj], y[jj+1]) ;
fprintf(stdout, "\n before zreal = %20.12e, zimag = %20.12e",
z[jj], z[jj+1]) ;
*/
z[jj] += areal*xreal - aimag*ximag + breal*yreal - bimag*yimag ;
z[jj+1] += areal*ximag + aimag*xreal + breal*yimag + bimag*yreal ;
/*
fprintf(stdout, "\n after zreal = %20.12e, zimag = %20.12e",
z[jj], z[jj+1]) ;
*/
}
return ; }
/*--------------------------------------------------------------------*/
/*
------------------------------------------------------------
purpose -- to scale a double complex vector by a 2x2 matrix
[ x ] := [ a b ] [ x ]
[ y ] [ c d ] [ y ]
created -- 98jan23, cca
------------------------------------------------------------
*/
void
ZVscale2 (
int size,
double x[],
double y[],
double areal,
double aimag,
double breal,
double bimag,
double creal,
double cimag,
double dreal,
double dimag
) {
double ximag, xreal, yimag, yreal ;
int ii, jj ;
if ( size < 0 || x == NULL || y == NULL ) {
fprintf(stderr, "\n fatal error in ZVscale2(%d,%p,%p,...)"
"\n bad input\n", size, x, y) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
xreal = x[jj] ;
ximag = x[jj+1] ;
yreal = y[jj] ;
yimag = y[jj+1] ;
x[jj] = areal*xreal - aimag*ximag + breal*yreal - bimag*yimag ;
x[jj+1] = areal*ximag + aimag*xreal + breal*yimag + bimag*yreal ;
y[jj] = creal*xreal - cimag*ximag + dreal*yreal - dimag*yimag ;
y[jj+1] = creal*ximag + cimag*xreal + dreal*yimag + dimag*yreal ;
}
return ; }
/*--------------------------------------------------------------------*/
/*
---------------------------------------
purpose -- to gather y[*] = x[index[*]]
created -- 98apr15, cca
---------------------------------------
*/
void
ZVgather (
int size,
double y[],
double x[],
int index[]
) {
if ( size > 0 ) {
if ( y == NULL || x == NULL || index == NULL ) {
fprintf(stderr, "\n fatal error in ZVgather, invalid input"
"\n size = %d, y = %p, x = %p, index = %p\n",
size, y, x, index) ;
exit(-1) ;
} else {
int i, j, k ;
for ( i = j = 0 ; i < size ; i++, j += 2 ) {
k = 2*index[i] ;
y[j] = x[k] ;
y[j+1] = x[k+1] ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------
purpose -- to scatter y[index[*]] = x[*]
created -- 98apr15, cca
----------------------------------------
*/
void
ZVscatter (
int size,
double y[],
int index[],
double x[]
) {
if ( size > 0 ) {
if ( y == NULL || x == NULL || index == NULL ) {
fprintf(stderr, "\n fatal error in ZVscatter, invalid data"
"\n size = %d, y = %p, index = %p, x = %p\n",
size, y, index, x) ;
exit(-1) ;
} else {
int i, j, k ;
for ( i = j = 0 ; i < size ; i++, j += 2 ) {
k = 2*index[i] ;
y[k] = x[j] ;
y[k+1] = x[j+1] ;
}
}
}
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^T [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU33 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, i10, i11, i12, i20, i21, i22,
r00, r01, r02, r10, r11, r12, r20, r21, r22,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = i10 = i11 = i12 = i20 = i21 = i22
= r00 = r01 = r02 = r10 = r11 = r12 = r20 = r21 = r22 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
r12 += yr1*xr2 - yi1*xi2 ; i12 += yr1*xi2 + yi1*xr2 ;
r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
r21 += yr2*xr1 - yi2*xi1 ; i21 += yr2*xi1 + yi2*xr1 ;
r22 += yr2*xr2 - yi2*xi2 ; i22 += yr2*xi2 + yi2*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
sums[ 6] = r10 ; sums[ 7] = i10 ;
sums[ 8] = r11 ; sums[ 9] = i11 ;
sums[10] = r12 ; sums[11] = i12 ;
sums[12] = r20 ; sums[13] = i20 ;
sums[14] = r21 ; sums[15] = i21 ;
sums[16] = r22 ; sums[17] = i22 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^T [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU32 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, i10, i11, i20, i21,
r00, r01, r10, r11, r20, r21,
xi0, xi1, xr0, xr1, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i01 = i10 = i11 = i20 = i21
= r00 = r01 = r10 = r11 = r20 = r21 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
r21 += yr2*xr1 - yi2*xi1 ; i21 += yr2*xi1 + yi2*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r10 ; sums[ 5] = i10 ;
sums[ 6] = r11 ; sums[ 7] = i11 ;
sums[ 8] = r20 ; sums[ 9] = i20 ;
sums[10] = r21 ; sums[11] = i21 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^T [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU31 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double sums[]
) {
double i00, i10, i20, r00, r10, r20,
xi0, xr0, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i10 = i20 = r00 = r10 = r20 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r10 ; sums[ 3] = i10 ;
sums[ 4] = r20 ; sums[ 5] = i20 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^T [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU23 (
int n,
double y0[],
double y1[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, i10, i11, i12,
r00, r01, r02, r10, r11, r12,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = i10 = i11 = i12
= r00 = r01 = r02 = r10 = r11 = r12 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
r12 += yr1*xr2 - yi1*xi2 ; i12 += yr1*xi2 + yi1*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
sums[ 6] = r10 ; sums[ 7] = i10 ;
sums[ 8] = r11 ; sums[ 9] = i11 ;
sums[10] = r12 ; sums[11] = i12 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^T [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU22 (
int n,
double y0[],
double y1[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, i10, i11, r00, r01, r10, r11,
xi0, xi1, xr0, xr1, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i01 = i10 = i11 = r00 = r01 = r10 = r11 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r10 ; sums[ 5] = i10 ;
sums[ 6] = r11 ; sums[ 7] = i11 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^T [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU21 (
int n,
double y0[],
double y1[],
double x0[],
double sums[]
) {
double i00, i10, r00, r10, xi0, xr0, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i10 = r00 = r10 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r10 ; sums[ 3] = i10 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^T [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU13 (
int n,
double y0[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, r00, r01, r02,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = r00 = r01 = r02 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^T [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU12 (
int n,
double y0[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, r00, r01, xi0, xi1, xr0, xr1, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = i01 = r00 = r01 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^T [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotU11 (
int n,
double y0[],
double x0[],
double sums[]
) {
double i00, r00, xi0, xr0, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = r00 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^H [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC33 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, i10, i11, i12, i20, i21, i22,
r00, r01, r02, r10, r11, r12, r20, r21, r22,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = i10 = i11 = i12 = i20 = i21 = i22
= r00 = r01 = r02 = r10 = r11 = r12 = r20 = r21 = r22 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
r12 += yr1*xr2 + yi1*xi2 ; i12 += yr1*xi2 - yi1*xr2 ;
r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
r21 += yr2*xr1 + yi2*xi1 ; i21 += yr2*xi1 - yi2*xr1 ;
r22 += yr2*xr2 + yi2*xi2 ; i22 += yr2*xi2 - yi2*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
sums[ 6] = r10 ; sums[ 7] = i10 ;
sums[ 8] = r11 ; sums[ 9] = i11 ;
sums[10] = r12 ; sums[11] = i12 ;
sums[12] = r20 ; sums[13] = i20 ;
sums[14] = r21 ; sums[15] = i21 ;
sums[16] = r22 ; sums[17] = i22 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^H [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC32 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, i10, i11, i20, i21,
r00, r01, r10, r11, r20, r21,
xi0, xi1, xr0, xr1, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i01 = i10 = i11 = i20 = i21
= r00 = r01 = r10 = r11 = r20 = r21 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
r21 += yr2*xr1 + yi2*xi1 ; i21 += yr2*xi1 - yi2*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r10 ; sums[ 5] = i10 ;
sums[ 6] = r11 ; sums[ 7] = i11 ;
sums[ 8] = r20 ; sums[ 9] = i20 ;
sums[10] = r21 ; sums[11] = i21 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 y2 ]^H [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC31 (
int n,
double y0[],
double y1[],
double y2[],
double x0[],
double sums[]
) {
double i00, i10, i20, r00, r10, r20,
xi0, xr0, yi0, yi1, yi2, yr0, yr1, yr2 ;
int ii, iloc, rloc ;
i00 = i10 = i20 = r00 = r10 = r20 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
yr2 = y2[rloc] ; yi2 = y2[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r10 ; sums[ 3] = i10 ;
sums[ 4] = r20 ; sums[ 5] = i20 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^H [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC23 (
int n,
double y0[],
double y1[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, i10, i11, i12,
r00, r01, r02, r10, r11, r12,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = i10 = i11 = i12
= r00 = r01 = r02 = r10 = r11 = r12 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
r12 += yr1*xr2 + yi1*xi2 ; i12 += yr1*xi2 - yi1*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
sums[ 6] = r10 ; sums[ 7] = i10 ;
sums[ 8] = r11 ; sums[ 9] = i11 ;
sums[10] = r12 ; sums[11] = i12 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^H [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC22 (
int n,
double y0[],
double y1[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, i10, i11, r00, r01, r10, r11,
xi0, xi1, xr0, xr1, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i01 = i10 = i11 = r00 = r01 = r10 = r11 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r10 ; sums[ 5] = i10 ;
sums[ 6] = r11 ; sums[ 7] = i11 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 y1 ]^H [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC21 (
int n,
double y0[],
double y1[],
double x0[],
double sums[]
) {
double i00, i10, r00, r10, xi0, xr0, yi0, yi1, yr0, yr1 ;
int ii, iloc, rloc ;
i00 = i10 = r00 = r10 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
yr1 = y1[rloc] ; yi1 = y1[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r10 ; sums[ 3] = i10 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^H [ x0 x1 x2]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC13 (
int n,
double y0[],
double x0[],
double x1[],
double x2[],
double sums[]
) {
double i00, i01, i02, r00, r01, r02,
xi0, xi1, xi2, xr0, xr1, xr2, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = i01 = i02 = r00 = r01 = r02 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
xr2 = x2[rloc] ; xi2 = x2[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
sums[ 4] = r02 ; sums[ 5] = i02 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^H [ x0 x1 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC12 (
int n,
double y0[],
double x0[],
double x1[],
double sums[]
) {
double i00, i01, r00, r01, xi0, xi1, xr0, xr1, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = i01 = r00 = r01 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
xr1 = x1[rloc] ; xi1 = x1[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
sums[ 2] = r01 ; sums[ 3] = i01 ;
return ; }
/*--------------------------------------------------------------------*/
/*
----------------------------------------------
purpose -- to compute the multiple dot product
[ y0 ]^H [ x0 ]
created -- 98apr17, cca
----------------------------------------------
*/
void
ZVdotC11 (
int n,
double y0[],
double x0[],
double sums[]
) {
double i00, r00, xi0, xr0, yi0, yr0 ;
int ii, iloc, rloc ;
i00 = r00 = 0.0 ;
for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
yr0 = y0[rloc] ; yi0 = y0[iloc] ;
xr0 = x0[rloc] ; xi0 = x0[iloc] ;
r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
}
sums[ 0] = r00 ; sums[ 1] = i00 ;
return ; }
/*--------------------------------------------------------------------*/
/*
-----------------------------
purpose -- to zero the vector
y := 0
where y is complex
created -- 98apr25, cca
-----------------------------
*/
void
ZVzero (
int size,
double y[]
) {
int ii, jj ;
if ( size < 0 || y == NULL ) {
fprintf(stderr, "\n fatal error in ZVzero(%d,%p)"
"\n bad input\n", size, y) ;
exit(-1) ;
}
for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
y[jj] = y[jj+1] = 0.0 ;
}
return ; }
/*--------------------------------------------------------------------*/
syntax highlighted by Code2HTML, v. 0.9.1