1 /* ZV.c */
2
3 #include "../Utilities.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 --------------------------------------------------------
8 purpose -- to return the absolute value of (areal,aimag)
9
10 created -- 98jan24, cca
11 --------------------------------------------------------
12 */
13 double
Zabs(double real,double imag)14 Zabs (
15 double real,
16 double imag
17 ) {
18 double abs, val ;
19
20 if ( real == 0.0 ) {
21 abs = fabs(imag) ;
22 } else if ( imag == 0.0 ) {
23 abs = fabs(real) ;
24 } else if ( real >= imag ) {
25 val = imag/real ;
26 abs = fabs(real)*sqrt(1 + val*val) ;
27 } else {
28 val = real/imag ;
29 abs = fabs(imag)*sqrt(1 + val*val) ;
30 }
31 return(abs) ; }
32
33 /*--------------------------------------------------------------------*/
34 /*
35 ---------------------------------------------------
36 purpose -- given (areal,aimag),
37 compute (breal,bimag) = 1/(areal,aimag),
38
39 put breal into *pbreal
40 put bimag into *pbimag
41
42 created -- 98jan23, cca
43 ---------------------------------------------------
44 */
45 int
Zrecip(double areal,double aimag,double * pbreal,double * pbimag)46 Zrecip (
47 double areal,
48 double aimag,
49 double *pbreal,
50 double *pbimag
51 ) {
52 double bimag, breal, fac ;
53 if ( areal == 0.0 && aimag == 0.0 ) {
54 return(0) ;
55 }
56 if ( fabs(areal) >= fabs(aimag) ) {
57 fac = aimag/areal ;
58 breal = 1./(areal + fac*aimag) ;
59 bimag = -fac * breal ;
60 } else {
61 fac = areal/aimag ;
62 bimag = -1./(aimag + fac*areal) ;
63 breal = -fac * bimag ;
64 }
65 *pbreal = breal ;
66 *pbimag = bimag ;
67
68 return(1) ; }
69
70 /*--------------------------------------------------------------------*/
71 /*
72 ---------------------------------------------------------------------
73 given [ (areal,aimag) (breal,bimag) ]
74 [ (creal,cimag) (dreal,dimag) ]
75 compute [ (ereal,eimag) (freal,fimag) ]
76 [ (greal,gimag) (hreal,himag) ]
77 where
78 I = [ (areal,aimag) (breal,bimag) ] * [ (ereal,eimag) (freal,fimag) ]
79 [ (creal,cimag) (dreal,dimag) ] [ (greal,gimag) (hreal,himag) ]
80
81 note, any of {pereal, peimag, pfreal, pfimag, pgreal, pgimag,
82 phreal, phimag} can be NULL. if not NULL, their value is filled.
83 this is useful when the input matrix is symmetric or hermitian.
84
85 created -- 98jan23, cca
86 ---------------------------------------------------------------------
87 */
88 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)89 Zrecip2 (
90 double areal,
91 double aimag,
92 double breal,
93 double bimag,
94 double creal,
95 double cimag,
96 double dreal,
97 double dimag,
98 double *pereal,
99 double *peimag,
100 double *pfreal,
101 double *pfimag,
102 double *pgreal,
103 double *pgimag,
104 double *phreal,
105 double *phimag
106 ) {
107 double xreal, ximag, yreal, yimag ;
108 /*
109 -------------------------------
110 step one, compute x = a*d - b*c
111 -------------------------------
112 */
113 xreal = areal*dreal - aimag*dimag - breal*creal + bimag*cimag ;
114 ximag = areal*dimag + aimag*dreal - breal*cimag - bimag*creal ;
115 /*
116 -------------------------
117 step two, compute y = 1/x
118 -------------------------
119 */
120 Zrecip(xreal, ximag, &yreal, &yimag) ;
121 /*
122 -----------------------
123 step three,
124 [ e f ] = y * [ d -b ]
125 [ g h ] [ -c a ]
126 -----------------------
127 */
128 if ( pereal != NULL ) {
129 *pereal = dreal*yreal - dimag*yimag ;
130 }
131 if ( peimag != NULL ) {
132 *peimag = dreal*yimag + dimag*yreal ;
133 }
134 if ( pfreal != NULL ) {
135 *pfreal = -breal*yreal + bimag*yimag ;
136 }
137 if ( pfimag != NULL ) {
138 *pfimag = -breal*yimag - bimag*yreal ;
139 }
140 if ( pgreal != NULL ) {
141 *pgreal = -creal*yreal + cimag*yimag ;
142 }
143 if ( pgimag != NULL ) {
144 *pgimag = -creal*yimag - cimag*yreal ;
145 }
146 if ( phreal != NULL ) {
147 *phreal = areal*yreal - aimag*yimag ;
148 }
149 if ( phimag != NULL ) {
150 *phimag = areal*yimag + aimag*yreal ;
151 }
152 return(1) ; }
153
154 /*--------------------------------------------------------------------*/
155 /*
156 ------------------------------------------------
157 purpose -- to allocate, initialize and
158 return a complex vector x[]
159
160 n -- length of the complex vector (2*n double's)
161
162 x[ii] = (real, imag) for 0 <= ii < n
163
164 created -- 98jan23, cca
165 ------------------------------------------------
166 */
167 double *
ZVinit(int n,double real,double imag)168 ZVinit (
169 int n,
170 double real,
171 double imag
172 ) {
173 double *x ;
174 int ii, jj ;
175
176 if ( n <= 0 ) {
177 fprintf(stderr, "\n fatal error in ZVinit(%d,%f,%f)"
178 "\n bad input\n", n, real, imag) ;
179 exit(-1) ;
180 }
181 ALLOCATE(x, double, 2*n) ;
182 for ( ii = jj = 0 ; ii < n ; ii++, jj += 2 ) {
183 x[jj] = real ;
184 x[jj+1] = imag ;
185 }
186 return(x) ; }
187
188 /*--------------------------------------------------------------------*/
189 /*
190 -------------------------------------------
191 purpose -- to perform a complex dot product
192
193 (*prdot,*pidot) = y^T x
194
195 where y and x are complex
196
197 created -- 98apr15, cca
198 -------------------------------------------
199 */
200 void
ZVdotU(int size,double y[],double x[],double * prdot,double * pidot)201 ZVdotU (
202 int size,
203 double y[],
204 double x[],
205 double *prdot,
206 double *pidot
207 ) {
208 double isum, rsum, ximag, xreal, yimag, yreal ;
209 int ii, jj ;
210
211 if ( size < 0 || y == NULL || x == NULL
212 || prdot == NULL || pidot == NULL ) {
213 fprintf(stderr, "\n fatal error in ZVdotU(%d,%p,%p,%p,%p)"
214 "\n bad input\n", size, y, x, prdot, pidot) ;
215 exit(-1) ;
216 }
217 isum = rsum = 0.0 ;
218 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
219 xreal = x[jj] ;
220 ximag = x[jj+1] ;
221 yreal = y[jj] ;
222 yimag = y[jj+1] ;
223 rsum += xreal*yreal - ximag*yimag ;
224 isum += xreal*yimag + ximag*yreal ;
225 }
226 *prdot = rsum ;
227 *pidot = isum ;
228
229 return ; }
230
231 /*--------------------------------------------------------------------*/
232 /*
233 ------------------------------------------------------
234 purpose -- to perform a conjugated complex dot product
235
236 (*prdot,*pidot) = (conjugate(y))^T x = y^H x
237
238 where y and x are complex
239
240 created -- 98apr15, cca
241 ------------------------------------------------------
242 */
243 void
ZVdotC(int size,double y[],double x[],double * prdot,double * pidot)244 ZVdotC (
245 int size,
246 double y[],
247 double x[],
248 double *prdot,
249 double *pidot
250 ) {
251 double isum, rsum, ximag, xreal, yimag, yreal ;
252 int ii, jj ;
253
254 if ( size < 0 || y == NULL || x == NULL
255 || prdot == NULL || pidot == NULL ) {
256 fprintf(stderr, "\n fatal error in ZVdotC(%d,%p,%p,%p,%p)"
257 "\n bad input\n", size, y, x, prdot, pidot) ;
258 exit(-1) ;
259 }
260 isum = rsum = 0.0 ;
261 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
262 xreal = x[jj] ;
263 ximag = x[jj+1] ;
264 yreal = y[jj] ;
265 yimag = y[jj+1] ;
266 rsum += xreal*yreal + ximag*yimag ;
267 isum += - xreal*yimag + ximag*yreal ;
268 }
269 *prdot = rsum ;
270 *pidot = isum ;
271
272 return ; }
273
274 /*--------------------------------------------------------------------*/
275 /*
276 ---------------------------------------------------
277 purpose -- to perform a indexed complex dot product
278
279 (*prdot,*pidot) = sum_k y[index[k]]*x[k]
280
281 where y and x are complex
282
283 created -- 98apr15, cca
284 ---------------------------------------------------
285 */
286 void
ZVdotiU(int size,double y[],int index[],double x[],double * prdot,double * pidot)287 ZVdotiU (
288 int size,
289 double y[],
290 int index[],
291 double x[],
292 double *prdot,
293 double *pidot
294 ) {
295 double isum, rsum, ximag, xreal, yimag, yreal ;
296 int ii, jj ;
297
298 if ( size < 0 || y == NULL || index == NULL || x == NULL
299 || prdot == NULL || pidot == NULL ) {
300 fprintf(stderr, "\n fatal error in ZVdotiU(%d,%p,%p,%p,%p,%p)"
301 "\n bad input\n", size, y, index, x, prdot, pidot) ;
302 exit(-1) ;
303 }
304 isum = rsum = 0.0 ;
305 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
306 /*
307 fprintf(stdout,
308 "\n %% ii = %d, jj = %d, kk = %d", ii, jj, index[ii]) ;
309 fflush(stdout) ;
310 */
311 xreal = x[jj] ;
312 ximag = x[jj+1] ;
313 yreal = y[2*index[ii]] ;
314 yimag = y[2*index[ii]+1] ;
315 rsum += xreal*yreal - ximag*yimag ;
316 isum += xreal*yimag + ximag*yreal ;
317 }
318 *prdot = rsum ;
319 *pidot = isum ;
320
321 return ; }
322
323 /*--------------------------------------------------------------------*/
324 /*
325 -------------------------------------------------------------
326 purpose -- to perform a indexed conjugate complex dot product
327
328 (*prdot,*pidot) = sum_k conjugate(y[index[k]])*x[k]
329
330 where y and x are complex
331
332 created -- 98apr15, cca
333 -------------------------------------------------------------
334 */
335 void
ZVdotiC(int size,double y[],int index[],double x[],double * prdot,double * pidot)336 ZVdotiC (
337 int size,
338 double y[],
339 int index[],
340 double x[],
341 double *prdot,
342 double *pidot
343 ) {
344 double isum, rsum, ximag, xreal, yimag, yreal ;
345 int ii, jj ;
346
347 if ( size < 0 || y == NULL || index == NULL || x == NULL
348 || prdot == NULL || pidot == NULL ) {
349 fprintf(stderr, "\n fatal error in ZVdotiU(%d,%p,%p,%p,%p,%p)"
350 "\n bad input\n", size, y, index, x, prdot, pidot) ;
351 exit(-1) ;
352 }
353 isum = rsum = 0.0 ;
354 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
355 xreal = x[jj] ;
356 ximag = x[jj+1] ;
357 yreal = y[2*index[ii]] ;
358 yimag = y[2*index[ii]+1] ;
359 rsum += xreal*yreal + ximag*yimag ;
360 isum += - xreal*yimag + ximag*yreal ;
361 }
362 *prdot = rsum ;
363 *pidot = isum ;
364
365 return ; }
366
367 /*--------------------------------------------------------------------*/
368 /*
369 ------------------------------------
370 purpose -- to perform a complex axpy
371
372 y := y + (areal, aimag) * x
373
374 where y and x are complex
375
376 created -- 98jan22, cca
377 ------------------------------------
378 */
379 void
ZVaxpy(int size,double y[],double areal,double aimag,double x[])380 ZVaxpy (
381 int size,
382 double y[],
383 double areal,
384 double aimag,
385 double x[]
386 ) {
387 double ximag, xreal ;
388 int ii, jj ;
389
390 if ( size < 0 || y == NULL || x == NULL ) {
391 fprintf(stderr, "\n fatal error in ZVaxpy(%d,%p,%f,%f,%p)"
392 "\n bad input\n", size, y, areal, aimag, x) ;
393 exit(-1) ;
394 }
395 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
396 xreal = x[jj] ;
397 ximag = x[jj+1] ;
398 /*
399 fprintf(stdout, "\n ii = %d, xreal = %20.12e, ximag = %20.12e",
400 ii, xreal, ximag) ;
401 fprintf(stdout, "\n before yreal = %20.12e, yimag = %20.12e",
402 y[jj], y[jj+1]) ;
403 */
404 y[jj] += areal*xreal - aimag*ximag ;
405 y[jj+1] += areal*ximag + aimag*xreal ;
406 /*
407 fprintf(stdout, "\n after yreal = %20.12e, yimag = %20.12e",
408 y[jj], y[jj+1]) ;
409 */
410 }
411 return ; }
412
413 /*--------------------------------------------------------------------*/
414 /*
415 ------------------------------------------------------------
416 purpose -- to scale a double complex vector by (xreal,ximag)
417
418 y := y * (areal, aimag)
419
420 created -- 98jan22, cca
421 ------------------------------------------------------------
422 */
423 void
ZVscale(int size,double y[],double areal,double aimag)424 ZVscale (
425 int size,
426 double y[],
427 double areal,
428 double aimag
429 ) {
430 double yimag, yreal ;
431 int ii, jj ;
432
433 if ( size < 0 || y == NULL ) {
434 fprintf(stderr, "\n fatal error in ZVscale(%d,%p,%f,%f)"
435 "\n bad input\n", size, y, areal, aimag) ;
436 exit(-1) ;
437 }
438 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
439 yreal = y[jj] ;
440 yimag = y[jj+1] ;
441 y[jj] = areal*yreal - aimag*yimag ;
442 y[jj+1] = areal*yimag + aimag*yreal ;
443 }
444 return ; }
445
446 /*--------------------------------------------------------------------*/
447 /*
448 -----------------------------------------------
449 purpose --- to print a complex vector to a file
450
451 created -- 98jan23, cca
452 -----------------------------------------------
453 */
454 void
ZVfprintf(FILE * fp,int size,double y[])455 ZVfprintf (
456 FILE *fp,
457 int size,
458 double y[]
459 ) {
460 int ii, jj ;
461
462 if ( size < 0 || y == NULL ) {
463 fprintf(stderr, "\n fatal error in ZVfprintf(%p,%d,%p)"
464 "\n bad input\n", fp, size, y) ;
465 exit(-1) ;
466 }
467 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
468 /*
469 if ( ii % 2 == 0 ) {
470 fprintf(fp, "\n") ;
471 }
472 */
473 fprintf(fp, "\n < %12.4e, %12.4e >", y[jj], y[jj+1]) ;
474 }
475
476 return ; }
477
478 /*--------------------------------------------------------------------*/
479 /*
480 ----------------------------------
481 return the minimum absolute value
482 of the entries in a complex vector
483
484 created -- 98jan23, cca
485 ----------------------------------
486 */
487 double
ZVminabs(int size,double x[])488 ZVminabs (
489 int size,
490 double x[]
491 ) {
492 double abs, imag, minabs, real, val ;
493 int ii, jj ;
494
495 if ( size < 0 || x == NULL ) {
496 fprintf(stderr, "\n fatal error in ZVminabs(%d,%p)"
497 "\n bad input\n", size, x) ;
498 exit(-1) ;
499 }
500 minabs = 0.0 ;
501 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
502 real = fabs(x[jj]) ;
503 imag = fabs(x[jj+1]) ;
504 /*
505 fprintf(stdout, "\n i = %d, <%12.4e, %12.4e>", ii, real, imag) ;
506 */
507 if ( real == 0.0 ) {
508 abs = imag ;
509 } else if ( imag == 0.0 ) {
510 abs = real ;
511 } else if ( real >= imag ) {
512 val = imag/real ;
513 abs = real*sqrt(1 + val*val) ;
514 } else {
515 val = real/imag ;
516 abs = imag*sqrt(1 + val*val) ;
517 }
518 /*
519 fprintf(stdout, " abs = %12.4e", abs) ;
520 */
521 if ( ii == 0 || minabs > abs ) {
522 minabs = abs ;
523 }
524 }
525 /*
526 fprintf(stdout, "\n minabs = %12.4e", minabs) ;
527 */
528
529 return(minabs) ; }
530
531 /*--------------------------------------------------------------------*/
532 /*
533 ----------------------------------
534 return the maximum absolute value
535 of the entries in a complex vector
536
537 created -- 98jan23, cca
538 ----------------------------------
539 */
540 double
ZVmaxabs(int size,double x[])541 ZVmaxabs (
542 int size,
543 double x[]
544 ) {
545 double abs, imag, maxabs, real, val ;
546 int ii, jj ;
547
548 if ( size < 0 || x == NULL ) {
549 fprintf(stderr, "\n fatal error in ZVmaxabs(%d,%p)"
550 "\n bad input\n", size, x) ;
551 exit(-1) ;
552 }
553 maxabs = 0.0 ;
554 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
555 real = fabs(x[jj]) ;
556 imag = fabs(x[jj+1]) ;
557 /*
558 fprintf(stdout, "\n i = %d, <%12.4e, %12.4e>", ii, real, imag) ;
559 */
560 if ( real == 0.0 ) {
561 abs = imag ;
562 } else if ( imag == 0.0 ) {
563 abs = real ;
564 } else if ( real >= imag ) {
565 val = imag/real ;
566 abs = real*sqrt(1 + val*val) ;
567 } else {
568 val = real/imag ;
569 abs = imag*sqrt(1 + val*val) ;
570 }
571 /*
572 fprintf(stdout, " abs = %12.4e", abs) ;
573 */
574 if ( ii == 0 || maxabs < abs ) {
575 maxabs = abs ;
576 }
577 }
578 /*
579 fprintf(stdout, "\n maxabs = %12.4e", maxabs) ;
580 */
581
582 return(maxabs) ; }
583
584 /*--------------------------------------------------------------------*/
585 /*
586 ----------------------------------
587 copy a complex vector into another
588 y[] := x[]
589
590 created -- 98jan23, cca
591 ----------------------------------
592 */
593 void
ZVcopy(int size,double y[],double x[])594 ZVcopy (
595 int size,
596 double y[],
597 double x[]
598 ) {
599 int ii, jj ;
600
601 if ( size < 0 || y == NULL || x == NULL ) {
602 fprintf(stderr, "\n fatal error in ZVcopy(%d,%p,%p)"
603 "\n bad input\n", size, y, x) ;
604 exit(-1) ;
605 }
606 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
607 y[jj] = x[jj] ;
608 y[jj+1] = x[jj+1] ;
609 }
610 return ; }
611
612 /*--------------------------------------------------------------------*/
613 /*
614 --------------------------------------
615 subtract a complex vector from another
616 y[] := y[] - x[]
617
618 created -- 98may25, cca
619 --------------------------------------
620 */
621 void
ZVsub(int size,double y[],double x[])622 ZVsub (
623 int size,
624 double y[],
625 double x[]
626 ) {
627 int ii, jj ;
628
629 if ( size < 0 || y == NULL || x == NULL ) {
630 fprintf(stderr, "\n fatal error in ZVsub(%d,%p,%p)"
631 "\n bad input\n", size, y, x) ;
632 exit(-1) ;
633 }
634 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
635 y[jj] -= x[jj] ;
636 y[jj+1] -= x[jj+1] ;
637 }
638 return ; }
639
640 /*--------------------------------------------------------------------*/
641 /*
642 -----------------------------------------------------
643 purpose -- to perform a complex axpy with two vectors
644
645 z := z + (areal, aimag)*x + (breal, bimag)*y
646
647 where y and x are complex
648
649 created -- 98jan23, cca
650 ----------------------------------------------------
651 */
652 void
ZVaxpy2(int size,double z[],double areal,double aimag,double x[],double breal,double bimag,double y[])653 ZVaxpy2 (
654 int size,
655 double z[],
656 double areal,
657 double aimag,
658 double x[],
659 double breal,
660 double bimag,
661 double y[]
662 ) {
663 double ximag, xreal, yimag, yreal ;
664 int ii, jj ;
665
666 if ( size < 0 || y == NULL || x == NULL ) {
667 fprintf(stderr, "\n fatal error in ZVaxpy(%d,%p,%f,%f,%p)"
668 "\n bad input\n", size, y, areal, aimag, x) ;
669 exit(-1) ;
670 }
671 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
672 xreal = x[jj] ;
673 ximag = x[jj+1] ;
674 yreal = y[jj] ;
675 yimag = y[jj+1] ;
676 /*
677 fprintf(stdout, "\n ii = %d, xreal = %20.12e, ximag = %20.12e",
678 ii, xreal, ximag) ;
679 fprintf(stdout, "\n yreal = %20.12e, yimag = %20.12e",
680 y[jj], y[jj+1]) ;
681 fprintf(stdout, "\n before zreal = %20.12e, zimag = %20.12e",
682 z[jj], z[jj+1]) ;
683 */
684 z[jj] += areal*xreal - aimag*ximag + breal*yreal - bimag*yimag ;
685 z[jj+1] += areal*ximag + aimag*xreal + breal*yimag + bimag*yreal ;
686 /*
687 fprintf(stdout, "\n after zreal = %20.12e, zimag = %20.12e",
688 z[jj], z[jj+1]) ;
689 */
690 }
691 return ; }
692
693 /*--------------------------------------------------------------------*/
694 /*
695 ------------------------------------------------------------
696 purpose -- to scale a double complex vector by a 2x2 matrix
697
698 [ x ] := [ a b ] [ x ]
699 [ y ] [ c d ] [ y ]
700
701 created -- 98jan23, cca
702 ------------------------------------------------------------
703 */
704 void
ZVscale2(int size,double x[],double y[],double areal,double aimag,double breal,double bimag,double creal,double cimag,double dreal,double dimag)705 ZVscale2 (
706 int size,
707 double x[],
708 double y[],
709 double areal,
710 double aimag,
711 double breal,
712 double bimag,
713 double creal,
714 double cimag,
715 double dreal,
716 double dimag
717 ) {
718 double ximag, xreal, yimag, yreal ;
719 int ii, jj ;
720
721 if ( size < 0 || x == NULL || y == NULL ) {
722 fprintf(stderr, "\n fatal error in ZVscale2(%d,%p,%p,...)"
723 "\n bad input\n", size, x, y) ;
724 exit(-1) ;
725 }
726 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
727 xreal = x[jj] ;
728 ximag = x[jj+1] ;
729 yreal = y[jj] ;
730 yimag = y[jj+1] ;
731 x[jj] = areal*xreal - aimag*ximag + breal*yreal - bimag*yimag ;
732 x[jj+1] = areal*ximag + aimag*xreal + breal*yimag + bimag*yreal ;
733 y[jj] = creal*xreal - cimag*ximag + dreal*yreal - dimag*yimag ;
734 y[jj+1] = creal*ximag + cimag*xreal + dreal*yimag + dimag*yreal ;
735 }
736 return ; }
737
738 /*--------------------------------------------------------------------*/
739 /*
740 ---------------------------------------
741 purpose -- to gather y[*] = x[index[*]]
742
743 created -- 98apr15, cca
744 ---------------------------------------
745 */
746 void
ZVgather(int size,double y[],double x[],int index[])747 ZVgather (
748 int size,
749 double y[],
750 double x[],
751 int index[]
752 ) {
753 if ( size > 0 ) {
754 if ( y == NULL || x == NULL || index == NULL ) {
755 fprintf(stderr, "\n fatal error in ZVgather, invalid input"
756 "\n size = %d, y = %p, x = %p, index = %p\n",
757 size, y, x, index) ;
758 exit(-1) ;
759 } else {
760 int i, j, k ;
761 for ( i = j = 0 ; i < size ; i++, j += 2 ) {
762 k = 2*index[i] ;
763 y[j] = x[k] ;
764 y[j+1] = x[k+1] ;
765 }
766 }
767 }
768 return ; }
769
770 /*--------------------------------------------------------------------*/
771 /*
772 ----------------------------------------
773 purpose -- to scatter y[index[*]] = x[*]
774
775 created -- 98apr15, cca
776 ----------------------------------------
777 */
778 void
ZVscatter(int size,double y[],int index[],double x[])779 ZVscatter (
780 int size,
781 double y[],
782 int index[],
783 double x[]
784 ) {
785 if ( size > 0 ) {
786 if ( y == NULL || x == NULL || index == NULL ) {
787 fprintf(stderr, "\n fatal error in ZVscatter, invalid data"
788 "\n size = %d, y = %p, index = %p, x = %p\n",
789 size, y, index, x) ;
790 exit(-1) ;
791 } else {
792 int i, j, k ;
793 for ( i = j = 0 ; i < size ; i++, j += 2 ) {
794 k = 2*index[i] ;
795 y[k] = x[j] ;
796 y[k+1] = x[j+1] ;
797 }
798 }
799 }
800 return ; }
801
802 /*--------------------------------------------------------------------*/
803 /*
804 ----------------------------------------------
805 purpose -- to compute the multiple dot product
806
807 [ y0 y1 y2 ]^T [ x0 x1 x2]
808
809 created -- 98apr17, cca
810 ----------------------------------------------
811 */
812 void
ZVdotU33(int n,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[],double sums[])813 ZVdotU33 (
814 int n,
815 double y0[],
816 double y1[],
817 double y2[],
818 double x0[],
819 double x1[],
820 double x2[],
821 double sums[]
822 ) {
823 double i00, i01, i02, i10, i11, i12, i20, i21, i22,
824 r00, r01, r02, r10, r11, r12, r20, r21, r22,
825 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yi2, yr0, yr1, yr2 ;
826 int ii, iloc, rloc ;
827
828 i00 = i01 = i02 = i10 = i11 = i12 = i20 = i21 = i22
829 = r00 = r01 = r02 = r10 = r11 = r12 = r20 = r21 = r22 = 0.0 ;
830 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
831 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
832 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
833 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
834 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
835 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
836 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
837 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
838 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
839 r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
840 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
841 r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
842 r12 += yr1*xr2 - yi1*xi2 ; i12 += yr1*xi2 + yi1*xr2 ;
843 r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
844 r21 += yr2*xr1 - yi2*xi1 ; i21 += yr2*xi1 + yi2*xr1 ;
845 r22 += yr2*xr2 - yi2*xi2 ; i22 += yr2*xi2 + yi2*xr2 ;
846 }
847 sums[ 0] = r00 ; sums[ 1] = i00 ;
848 sums[ 2] = r01 ; sums[ 3] = i01 ;
849 sums[ 4] = r02 ; sums[ 5] = i02 ;
850 sums[ 6] = r10 ; sums[ 7] = i10 ;
851 sums[ 8] = r11 ; sums[ 9] = i11 ;
852 sums[10] = r12 ; sums[11] = i12 ;
853 sums[12] = r20 ; sums[13] = i20 ;
854 sums[14] = r21 ; sums[15] = i21 ;
855 sums[16] = r22 ; sums[17] = i22 ;
856
857 return ; }
858
859 /*--------------------------------------------------------------------*/
860 /*
861 ----------------------------------------------
862 purpose -- to compute the multiple dot product
863
864 [ y0 y1 y2 ]^T [ x0 x1 ]
865
866 created -- 98apr17, cca
867 ----------------------------------------------
868 */
869 void
ZVdotU32(int n,double y0[],double y1[],double y2[],double x0[],double x1[],double sums[])870 ZVdotU32 (
871 int n,
872 double y0[],
873 double y1[],
874 double y2[],
875 double x0[],
876 double x1[],
877 double sums[]
878 ) {
879 double i00, i01, i10, i11, i20, i21,
880 r00, r01, r10, r11, r20, r21,
881 xi0, xi1, xr0, xr1, yi0, yi1, yi2, yr0, yr1, yr2 ;
882 int ii, iloc, rloc ;
883
884 i00 = i01 = i10 = i11 = i20 = i21
885 = r00 = r01 = r10 = r11 = r20 = r21 = 0.0 ;
886 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
887 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
888 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
889 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
890 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
891 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
892 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
893 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
894 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
895 r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
896 r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
897 r21 += yr2*xr1 - yi2*xi1 ; i21 += yr2*xi1 + yi2*xr1 ;
898 }
899 sums[ 0] = r00 ; sums[ 1] = i00 ;
900 sums[ 2] = r01 ; sums[ 3] = i01 ;
901 sums[ 4] = r10 ; sums[ 5] = i10 ;
902 sums[ 6] = r11 ; sums[ 7] = i11 ;
903 sums[ 8] = r20 ; sums[ 9] = i20 ;
904 sums[10] = r21 ; sums[11] = i21 ;
905
906 return ; }
907
908 /*--------------------------------------------------------------------*/
909 /*
910 ----------------------------------------------
911 purpose -- to compute the multiple dot product
912
913 [ y0 y1 y2 ]^T [ x0 ]
914
915 created -- 98apr17, cca
916 ----------------------------------------------
917 */
918 void
ZVdotU31(int n,double y0[],double y1[],double y2[],double x0[],double sums[])919 ZVdotU31 (
920 int n,
921 double y0[],
922 double y1[],
923 double y2[],
924 double x0[],
925 double sums[]
926 ) {
927 double i00, i10, i20, r00, r10, r20,
928 xi0, xr0, yi0, yi1, yi2, yr0, yr1, yr2 ;
929 int ii, iloc, rloc ;
930
931 i00 = i10 = i20 = r00 = r10 = r20 = 0.0 ;
932 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
933 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
934 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
935 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
936 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
937 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
938 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
939 r20 += yr2*xr0 - yi2*xi0 ; i20 += yr2*xi0 + yi2*xr0 ;
940 }
941 sums[ 0] = r00 ; sums[ 1] = i00 ;
942 sums[ 2] = r10 ; sums[ 3] = i10 ;
943 sums[ 4] = r20 ; sums[ 5] = i20 ;
944
945 return ; }
946
947 /*--------------------------------------------------------------------*/
948 /*
949 ----------------------------------------------
950 purpose -- to compute the multiple dot product
951
952 [ y0 y1 ]^T [ x0 x1 x2]
953
954 created -- 98apr17, cca
955 ----------------------------------------------
956 */
957 void
ZVdotU23(int n,double y0[],double y1[],double x0[],double x1[],double x2[],double sums[])958 ZVdotU23 (
959 int n,
960 double y0[],
961 double y1[],
962 double x0[],
963 double x1[],
964 double x2[],
965 double sums[]
966 ) {
967 double i00, i01, i02, i10, i11, i12,
968 r00, r01, r02, r10, r11, r12,
969 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yr0, yr1 ;
970 int ii, iloc, rloc ;
971
972 i00 = i01 = i02 = i10 = i11 = i12
973 = r00 = r01 = r02 = r10 = r11 = r12 = 0.0 ;
974 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
975 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
976 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
977 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
978 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
979 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
980 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
981 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
982 r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
983 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
984 r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
985 r12 += yr1*xr2 - yi1*xi2 ; i12 += yr1*xi2 + yi1*xr2 ;
986 }
987 sums[ 0] = r00 ; sums[ 1] = i00 ;
988 sums[ 2] = r01 ; sums[ 3] = i01 ;
989 sums[ 4] = r02 ; sums[ 5] = i02 ;
990 sums[ 6] = r10 ; sums[ 7] = i10 ;
991 sums[ 8] = r11 ; sums[ 9] = i11 ;
992 sums[10] = r12 ; sums[11] = i12 ;
993
994 return ; }
995
996 /*--------------------------------------------------------------------*/
997 /*
998 ----------------------------------------------
999 purpose -- to compute the multiple dot product
1000
1001 [ y0 y1 ]^T [ x0 x1 ]
1002
1003 created -- 98apr17, cca
1004 ----------------------------------------------
1005 */
1006 void
ZVdotU22(int n,double y0[],double y1[],double x0[],double x1[],double sums[])1007 ZVdotU22 (
1008 int n,
1009 double y0[],
1010 double y1[],
1011 double x0[],
1012 double x1[],
1013 double sums[]
1014 ) {
1015 double i00, i01, i10, i11, r00, r01, r10, r11,
1016 xi0, xi1, xr0, xr1, yi0, yi1, yr0, yr1 ;
1017 int ii, iloc, rloc ;
1018
1019 i00 = i01 = i10 = i11 = r00 = r01 = r10 = r11 = 0.0 ;
1020 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1021 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1022 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1023 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1024 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1025 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
1026 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
1027 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
1028 r11 += yr1*xr1 - yi1*xi1 ; i11 += yr1*xi1 + yi1*xr1 ;
1029 }
1030 sums[ 0] = r00 ; sums[ 1] = i00 ;
1031 sums[ 2] = r01 ; sums[ 3] = i01 ;
1032 sums[ 4] = r10 ; sums[ 5] = i10 ;
1033 sums[ 6] = r11 ; sums[ 7] = i11 ;
1034
1035 return ; }
1036
1037 /*--------------------------------------------------------------------*/
1038 /*
1039 ----------------------------------------------
1040 purpose -- to compute the multiple dot product
1041
1042 [ y0 y1 ]^T [ x0 ]
1043
1044 created -- 98apr17, cca
1045 ----------------------------------------------
1046 */
1047 void
ZVdotU21(int n,double y0[],double y1[],double x0[],double sums[])1048 ZVdotU21 (
1049 int n,
1050 double y0[],
1051 double y1[],
1052 double x0[],
1053 double sums[]
1054 ) {
1055 double i00, i10, r00, r10, xi0, xr0, yi0, yi1, yr0, yr1 ;
1056 int ii, iloc, rloc ;
1057
1058 i00 = i10 = r00 = r10 = 0.0 ;
1059 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1060 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1061 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1062 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1063 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
1064 r10 += yr1*xr0 - yi1*xi0 ; i10 += yr1*xi0 + yi1*xr0 ;
1065 }
1066 sums[ 0] = r00 ; sums[ 1] = i00 ;
1067 sums[ 2] = r10 ; sums[ 3] = i10 ;
1068
1069 return ; }
1070
1071 /*--------------------------------------------------------------------*/
1072 /*
1073 ----------------------------------------------
1074 purpose -- to compute the multiple dot product
1075
1076 [ y0 ]^T [ x0 x1 x2]
1077
1078 created -- 98apr17, cca
1079 ----------------------------------------------
1080 */
1081 void
ZVdotU13(int n,double y0[],double x0[],double x1[],double x2[],double sums[])1082 ZVdotU13 (
1083 int n,
1084 double y0[],
1085 double x0[],
1086 double x1[],
1087 double x2[],
1088 double sums[]
1089 ) {
1090 double i00, i01, i02, r00, r01, r02,
1091 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yr0 ;
1092 int ii, iloc, rloc ;
1093
1094 i00 = i01 = i02 = r00 = r01 = r02 = 0.0 ;
1095 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1096 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1097 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1098 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1099 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
1100 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
1101 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
1102 r02 += yr0*xr2 - yi0*xi2 ; i02 += yr0*xi2 + yi0*xr2 ;
1103 }
1104 sums[ 0] = r00 ; sums[ 1] = i00 ;
1105 sums[ 2] = r01 ; sums[ 3] = i01 ;
1106 sums[ 4] = r02 ; sums[ 5] = i02 ;
1107
1108 return ; }
1109
1110 /*--------------------------------------------------------------------*/
1111 /*
1112 ----------------------------------------------
1113 purpose -- to compute the multiple dot product
1114
1115 [ y0 ]^T [ x0 x1 ]
1116
1117 created -- 98apr17, cca
1118 ----------------------------------------------
1119 */
1120 void
ZVdotU12(int n,double y0[],double x0[],double x1[],double sums[])1121 ZVdotU12 (
1122 int n,
1123 double y0[],
1124 double x0[],
1125 double x1[],
1126 double sums[]
1127 ) {
1128 double i00, i01, r00, r01, xi0, xi1, xr0, xr1, yi0, yr0 ;
1129 int ii, iloc, rloc ;
1130
1131 i00 = i01 = r00 = r01 = 0.0 ;
1132 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1133 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1134 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1135 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1136 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
1137 r01 += yr0*xr1 - yi0*xi1 ; i01 += yr0*xi1 + yi0*xr1 ;
1138 }
1139 sums[ 0] = r00 ; sums[ 1] = i00 ;
1140 sums[ 2] = r01 ; sums[ 3] = i01 ;
1141
1142 return ; }
1143
1144 /*--------------------------------------------------------------------*/
1145 /*
1146 ----------------------------------------------
1147 purpose -- to compute the multiple dot product
1148
1149 [ y0 ]^T [ x0 ]
1150
1151 created -- 98apr17, cca
1152 ----------------------------------------------
1153 */
1154 void
ZVdotU11(int n,double y0[],double x0[],double sums[])1155 ZVdotU11 (
1156 int n,
1157 double y0[],
1158 double x0[],
1159 double sums[]
1160 ) {
1161 double i00, r00, xi0, xr0, yi0, yr0 ;
1162 int ii, iloc, rloc ;
1163
1164 i00 = r00 = 0.0 ;
1165 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1166 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1167 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1168 r00 += yr0*xr0 - yi0*xi0 ; i00 += yr0*xi0 + yi0*xr0 ;
1169 }
1170 sums[ 0] = r00 ; sums[ 1] = i00 ;
1171
1172 return ; }
1173
1174 /*--------------------------------------------------------------------*/
1175 /*
1176 ----------------------------------------------
1177 purpose -- to compute the multiple dot product
1178
1179 [ y0 y1 y2 ]^H [ x0 x1 x2]
1180
1181 created -- 98apr17, cca
1182 ----------------------------------------------
1183 */
1184 void
ZVdotC33(int n,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[],double sums[])1185 ZVdotC33 (
1186 int n,
1187 double y0[],
1188 double y1[],
1189 double y2[],
1190 double x0[],
1191 double x1[],
1192 double x2[],
1193 double sums[]
1194 ) {
1195 double i00, i01, i02, i10, i11, i12, i20, i21, i22,
1196 r00, r01, r02, r10, r11, r12, r20, r21, r22,
1197 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yi2, yr0, yr1, yr2 ;
1198 int ii, iloc, rloc ;
1199
1200 i00 = i01 = i02 = i10 = i11 = i12 = i20 = i21 = i22
1201 = r00 = r01 = r02 = r10 = r11 = r12 = r20 = r21 = r22 = 0.0 ;
1202 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1203 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1204 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1205 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
1206 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1207 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1208 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
1209 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1210 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1211 r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
1212 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1213 r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
1214 r12 += yr1*xr2 + yi1*xi2 ; i12 += yr1*xi2 - yi1*xr2 ;
1215 r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
1216 r21 += yr2*xr1 + yi2*xi1 ; i21 += yr2*xi1 - yi2*xr1 ;
1217 r22 += yr2*xr2 + yi2*xi2 ; i22 += yr2*xi2 - yi2*xr2 ;
1218 }
1219 sums[ 0] = r00 ; sums[ 1] = i00 ;
1220 sums[ 2] = r01 ; sums[ 3] = i01 ;
1221 sums[ 4] = r02 ; sums[ 5] = i02 ;
1222 sums[ 6] = r10 ; sums[ 7] = i10 ;
1223 sums[ 8] = r11 ; sums[ 9] = i11 ;
1224 sums[10] = r12 ; sums[11] = i12 ;
1225 sums[12] = r20 ; sums[13] = i20 ;
1226 sums[14] = r21 ; sums[15] = i21 ;
1227 sums[16] = r22 ; sums[17] = i22 ;
1228
1229 return ; }
1230
1231 /*--------------------------------------------------------------------*/
1232 /*
1233 ----------------------------------------------
1234 purpose -- to compute the multiple dot product
1235
1236 [ y0 y1 y2 ]^H [ x0 x1 ]
1237
1238 created -- 98apr17, cca
1239 ----------------------------------------------
1240 */
1241 void
ZVdotC32(int n,double y0[],double y1[],double y2[],double x0[],double x1[],double sums[])1242 ZVdotC32 (
1243 int n,
1244 double y0[],
1245 double y1[],
1246 double y2[],
1247 double x0[],
1248 double x1[],
1249 double sums[]
1250 ) {
1251 double i00, i01, i10, i11, i20, i21,
1252 r00, r01, r10, r11, r20, r21,
1253 xi0, xi1, xr0, xr1, yi0, yi1, yi2, yr0, yr1, yr2 ;
1254 int ii, iloc, rloc ;
1255
1256 i00 = i01 = i10 = i11 = i20 = i21
1257 = r00 = r01 = r10 = r11 = r20 = r21 = 0.0 ;
1258 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1259 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1260 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1261 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
1262 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1263 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1264 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1265 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1266 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1267 r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
1268 r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
1269 r21 += yr2*xr1 + yi2*xi1 ; i21 += yr2*xi1 - yi2*xr1 ;
1270 }
1271 sums[ 0] = r00 ; sums[ 1] = i00 ;
1272 sums[ 2] = r01 ; sums[ 3] = i01 ;
1273 sums[ 4] = r10 ; sums[ 5] = i10 ;
1274 sums[ 6] = r11 ; sums[ 7] = i11 ;
1275 sums[ 8] = r20 ; sums[ 9] = i20 ;
1276 sums[10] = r21 ; sums[11] = i21 ;
1277
1278 return ; }
1279
1280 /*--------------------------------------------------------------------*/
1281 /*
1282 ----------------------------------------------
1283 purpose -- to compute the multiple dot product
1284
1285 [ y0 y1 y2 ]^H [ x0 ]
1286
1287 created -- 98apr17, cca
1288 ----------------------------------------------
1289 */
1290 void
ZVdotC31(int n,double y0[],double y1[],double y2[],double x0[],double sums[])1291 ZVdotC31 (
1292 int n,
1293 double y0[],
1294 double y1[],
1295 double y2[],
1296 double x0[],
1297 double sums[]
1298 ) {
1299 double i00, i10, i20, r00, r10, r20,
1300 xi0, xr0, yi0, yi1, yi2, yr0, yr1, yr2 ;
1301 int ii, iloc, rloc ;
1302
1303 i00 = i10 = i20 = r00 = r10 = r20 = 0.0 ;
1304 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1305 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1306 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1307 yr2 = y2[rloc] ; yi2 = y2[iloc] ;
1308 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1309 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1310 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1311 r20 += yr2*xr0 + yi2*xi0 ; i20 += yr2*xi0 - yi2*xr0 ;
1312 }
1313 sums[ 0] = r00 ; sums[ 1] = i00 ;
1314 sums[ 2] = r10 ; sums[ 3] = i10 ;
1315 sums[ 4] = r20 ; sums[ 5] = i20 ;
1316
1317 return ; }
1318
1319 /*--------------------------------------------------------------------*/
1320 /*
1321 ----------------------------------------------
1322 purpose -- to compute the multiple dot product
1323
1324 [ y0 y1 ]^H [ x0 x1 x2]
1325
1326 created -- 98apr17, cca
1327 ----------------------------------------------
1328 */
1329 void
ZVdotC23(int n,double y0[],double y1[],double x0[],double x1[],double x2[],double sums[])1330 ZVdotC23 (
1331 int n,
1332 double y0[],
1333 double y1[],
1334 double x0[],
1335 double x1[],
1336 double x2[],
1337 double sums[]
1338 ) {
1339 double i00, i01, i02, i10, i11, i12,
1340 r00, r01, r02, r10, r11, r12,
1341 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yi1, yr0, yr1 ;
1342 int ii, iloc, rloc ;
1343
1344 i00 = i01 = i02 = i10 = i11 = i12
1345 = r00 = r01 = r02 = r10 = r11 = r12 = 0.0 ;
1346 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1347 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1348 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1349 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1350 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1351 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
1352 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1353 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1354 r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
1355 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1356 r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
1357 r12 += yr1*xr2 + yi1*xi2 ; i12 += yr1*xi2 - yi1*xr2 ;
1358 }
1359 sums[ 0] = r00 ; sums[ 1] = i00 ;
1360 sums[ 2] = r01 ; sums[ 3] = i01 ;
1361 sums[ 4] = r02 ; sums[ 5] = i02 ;
1362 sums[ 6] = r10 ; sums[ 7] = i10 ;
1363 sums[ 8] = r11 ; sums[ 9] = i11 ;
1364 sums[10] = r12 ; sums[11] = i12 ;
1365
1366 return ; }
1367
1368 /*--------------------------------------------------------------------*/
1369 /*
1370 ----------------------------------------------
1371 purpose -- to compute the multiple dot product
1372
1373 [ y0 y1 ]^H [ x0 x1 ]
1374
1375 created -- 98apr17, cca
1376 ----------------------------------------------
1377 */
1378 void
ZVdotC22(int n,double y0[],double y1[],double x0[],double x1[],double sums[])1379 ZVdotC22 (
1380 int n,
1381 double y0[],
1382 double y1[],
1383 double x0[],
1384 double x1[],
1385 double sums[]
1386 ) {
1387 double i00, i01, i10, i11, r00, r01, r10, r11,
1388 xi0, xi1, xr0, xr1, yi0, yi1, yr0, yr1 ;
1389 int ii, iloc, rloc ;
1390
1391 i00 = i01 = i10 = i11 = r00 = r01 = r10 = r11 = 0.0 ;
1392 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1393 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1394 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1395 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1396 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1397 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1398 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1399 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1400 r11 += yr1*xr1 + yi1*xi1 ; i11 += yr1*xi1 - yi1*xr1 ;
1401 }
1402 sums[ 0] = r00 ; sums[ 1] = i00 ;
1403 sums[ 2] = r01 ; sums[ 3] = i01 ;
1404 sums[ 4] = r10 ; sums[ 5] = i10 ;
1405 sums[ 6] = r11 ; sums[ 7] = i11 ;
1406
1407 return ; }
1408
1409 /*--------------------------------------------------------------------*/
1410 /*
1411 ----------------------------------------------
1412 purpose -- to compute the multiple dot product
1413
1414 [ y0 y1 ]^H [ x0 ]
1415
1416 created -- 98apr17, cca
1417 ----------------------------------------------
1418 */
1419 void
ZVdotC21(int n,double y0[],double y1[],double x0[],double sums[])1420 ZVdotC21 (
1421 int n,
1422 double y0[],
1423 double y1[],
1424 double x0[],
1425 double sums[]
1426 ) {
1427 double i00, i10, r00, r10, xi0, xr0, yi0, yi1, yr0, yr1 ;
1428 int ii, iloc, rloc ;
1429
1430 i00 = i10 = r00 = r10 = 0.0 ;
1431 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1432 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1433 yr1 = y1[rloc] ; yi1 = y1[iloc] ;
1434 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1435 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1436 r10 += yr1*xr0 + yi1*xi0 ; i10 += yr1*xi0 - yi1*xr0 ;
1437 }
1438 sums[ 0] = r00 ; sums[ 1] = i00 ;
1439 sums[ 2] = r10 ; sums[ 3] = i10 ;
1440
1441 return ; }
1442
1443 /*--------------------------------------------------------------------*/
1444 /*
1445 ----------------------------------------------
1446 purpose -- to compute the multiple dot product
1447
1448 [ y0 ]^H [ x0 x1 x2]
1449
1450 created -- 98apr17, cca
1451 ----------------------------------------------
1452 */
1453 void
ZVdotC13(int n,double y0[],double x0[],double x1[],double x2[],double sums[])1454 ZVdotC13 (
1455 int n,
1456 double y0[],
1457 double x0[],
1458 double x1[],
1459 double x2[],
1460 double sums[]
1461 ) {
1462 double i00, i01, i02, r00, r01, r02,
1463 xi0, xi1, xi2, xr0, xr1, xr2, yi0, yr0 ;
1464 int ii, iloc, rloc ;
1465
1466 i00 = i01 = i02 = r00 = r01 = r02 = 0.0 ;
1467 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1468 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1469 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1470 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1471 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
1472 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1473 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1474 r02 += yr0*xr2 + yi0*xi2 ; i02 += yr0*xi2 - yi0*xr2 ;
1475 }
1476 sums[ 0] = r00 ; sums[ 1] = i00 ;
1477 sums[ 2] = r01 ; sums[ 3] = i01 ;
1478 sums[ 4] = r02 ; sums[ 5] = i02 ;
1479
1480 return ; }
1481
1482 /*--------------------------------------------------------------------*/
1483 /*
1484 ----------------------------------------------
1485 purpose -- to compute the multiple dot product
1486
1487 [ y0 ]^H [ x0 x1 ]
1488
1489 created -- 98apr17, cca
1490 ----------------------------------------------
1491 */
1492 void
ZVdotC12(int n,double y0[],double x0[],double x1[],double sums[])1493 ZVdotC12 (
1494 int n,
1495 double y0[],
1496 double x0[],
1497 double x1[],
1498 double sums[]
1499 ) {
1500 double i00, i01, r00, r01, xi0, xi1, xr0, xr1, yi0, yr0 ;
1501 int ii, iloc, rloc ;
1502
1503 i00 = i01 = r00 = r01 = 0.0 ;
1504 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1505 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1506 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1507 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
1508 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1509 r01 += yr0*xr1 + yi0*xi1 ; i01 += yr0*xi1 - yi0*xr1 ;
1510 }
1511 sums[ 0] = r00 ; sums[ 1] = i00 ;
1512 sums[ 2] = r01 ; sums[ 3] = i01 ;
1513
1514 return ; }
1515
1516 /*--------------------------------------------------------------------*/
1517 /*
1518 ----------------------------------------------
1519 purpose -- to compute the multiple dot product
1520
1521 [ y0 ]^H [ x0 ]
1522
1523 created -- 98apr17, cca
1524 ----------------------------------------------
1525 */
1526 void
ZVdotC11(int n,double y0[],double x0[],double sums[])1527 ZVdotC11 (
1528 int n,
1529 double y0[],
1530 double x0[],
1531 double sums[]
1532 ) {
1533 double i00, r00, xi0, xr0, yi0, yr0 ;
1534 int ii, iloc, rloc ;
1535
1536 i00 = r00 = 0.0 ;
1537 for ( ii = rloc = 0, iloc = 1 ; ii < n ; ii++, rloc += 2, iloc += 2 ) {
1538 yr0 = y0[rloc] ; yi0 = y0[iloc] ;
1539 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
1540 r00 += yr0*xr0 + yi0*xi0 ; i00 += yr0*xi0 - yi0*xr0 ;
1541 }
1542 sums[ 0] = r00 ; sums[ 1] = i00 ;
1543
1544 return ; }
1545
1546 /*--------------------------------------------------------------------*/
1547 /*
1548 -----------------------------
1549 purpose -- to zero the vector
1550
1551 y := 0
1552
1553 where y is complex
1554
1555 created -- 98apr25, cca
1556 -----------------------------
1557 */
1558 void
ZVzero(int size,double y[])1559 ZVzero (
1560 int size,
1561 double y[]
1562 ) {
1563 int ii, jj ;
1564
1565 if ( size < 0 || y == NULL ) {
1566 fprintf(stderr, "\n fatal error in ZVzero(%d,%p)"
1567 "\n bad input\n", size, y) ;
1568 exit(-1) ;
1569 }
1570 for ( ii = jj = 0 ; ii < size ; ii++, jj += 2 ) {
1571 y[jj] = y[jj+1] = 0.0 ;
1572 }
1573 return ; }
1574
1575 /*--------------------------------------------------------------------*/
1576