1 /* scalevec.c */
2
3 #include "../SubMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 static void diagonal_scale3vec ( SubMtx *mtxA, double y0[], double y1[],
7 double y2[], double x0[], double x1[], double x2[] ) ;
8 static void diagonal_scale2vec ( SubMtx *mtxA, double y0[], double y1[],
9 double x0[], double x1[] ) ;
10 static void diagonal_scale1vec ( SubMtx *mtxA,
11 double y0[], double x0[] ) ;
12 static void block_diagonal_sym_scale3vec ( SubMtx *mtxA, double y0[],
13 double y1[], double y2[], double x0[], double x1[], double x2[] ) ;
14 static void block_diagonal_sym_scale2vec ( SubMtx *mtxA, double y0[],
15 double y1[], double x0[], double x1[] ) ;
16 static void block_diagonal_sym_scale1vec ( SubMtx *mtxA, double y0[],
17 double x0[] ) ;
18 static void block_diagonal_herm_scale3vec ( SubMtx *mtxA, double y0[],
19 double y1[], double y2[], double x0[], double x1[], double x2[] ) ;
20 static void block_diagonal_herm_scale2vec ( SubMtx *mtxA, double y0[],
21 double y1[], double x0[], double x1[] ) ;
22 static void block_diagonal_herm_scale1vec ( SubMtx *mtxA, double y0[],
23 double x0[] ) ;
24 /*--------------------------------------------------------------------*/
25 /*
26 -----------------------------------
27 purpose -- compute [y0] := A * [x0]
28
29 created -- 98apr17, cca
30 -----------------------------------
31 */
32 void
SubMtx_scale1vec(SubMtx * mtxA,double y0[],double x0[])33 SubMtx_scale1vec (
34 SubMtx *mtxA,
35 double y0[],
36 double x0[]
37 ) {
38 /*
39 ---------------
40 check the input
41 ---------------
42 */
43 if ( mtxA == NULL || y0 == NULL || x0 == NULL ) {
44 fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
45 "\n bad input\n", mtxA, y0, x0) ;
46 exit(-1) ;
47 }
48 if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
49 fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
50 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
51 mtxA, y0, x0, mtxA->type) ;
52 exit(-1) ;
53 }
54 switch ( mtxA->mode ) {
55 case SUBMTX_DIAGONAL :
56 diagonal_scale1vec(mtxA, y0, x0) ;
57 break ;
58 case SUBMTX_BLOCK_DIAGONAL_SYM :
59 block_diagonal_sym_scale1vec(mtxA, y0, x0) ;
60 break ;
61 case SUBMTX_BLOCK_DIAGONAL_HERM :
62 if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
63 fprintf(stderr, "\n fatal error in SubMtx_scale1vec(%p,%p,%p)"
64 "\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
65 mtxA, y0, x0, mtxA->type) ;
66 exit(-1) ;
67 }
68 block_diagonal_herm_scale1vec(mtxA, y0, x0) ;
69 break ;
70 default :
71 fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
72 "\n matrix mode not supported"
73 "\n must be SUBMTX_DIAGONAL,"
74 "\n or SUBMTX_BLOCK_DIAGONAL_SYM"
75 "\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
76 exit(-1) ;
77 }
78 return ; }
79
80 /*--------------------------------------------------------------------*/
81 /*
82 -------------------------------------
83 purpose -- compute [y0 y1] := [x0 x1]
84
85 created -- 98apr17, cca
86 -------------------------------------
87 */
88 void
SubMtx_scale2vec(SubMtx * mtxA,double y0[],double y1[],double x0[],double x1[])89 SubMtx_scale2vec (
90 SubMtx *mtxA,
91 double y0[],
92 double y1[],
93 double x0[],
94 double x1[]
95 ) {
96 /*
97 ---------------
98 check the input
99 ---------------
100 */
101 if ( mtxA == NULL || y0 == NULL || y1 == NULL
102 || x0 == NULL || x1 == NULL ) {
103 fprintf(stderr, "\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
104 "\n bad input\n", mtxA, y0, y1, x0, x1) ;
105 exit(-1) ;
106 }
107 if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
108 fprintf(stderr, "\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
109 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
110 mtxA, y0, y1, x0, x1, mtxA->type) ;
111 exit(-1) ;
112 }
113 switch ( mtxA->mode ) {
114 case SUBMTX_DIAGONAL :
115 diagonal_scale2vec(mtxA, y0, y1, x0, x1) ;
116 break ;
117 case SUBMTX_BLOCK_DIAGONAL_SYM :
118 block_diagonal_sym_scale2vec(mtxA, y0, y1, x0, x1) ;
119 break ;
120 case SUBMTX_BLOCK_DIAGONAL_HERM :
121 if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
122 fprintf(stderr,
123 "\n fatal error in SubMtx_scale2vec(%p,%p,%p,%p,%p)"
124 "\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
125 mtxA, y0, y1, x0, x1, mtxA->type) ;
126 exit(-1) ;
127 }
128 block_diagonal_herm_scale2vec(mtxA, y0, y1, x0, x1) ;
129 break ;
130 default :
131 fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
132 "\n matrix type not supported"
133 "\n must be SUBMTX_DIAGONAL,"
134 "\n or SUBMTX_BLOCK_DIAGONAL_SYM"
135 "\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
136 exit(-1) ;
137 }
138 return ; }
139
140 /*--------------------------------------------------------------------*/
141 /*
142 -----------------------------------------------
143 purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
144
145 created -- 98apr17, cca
146 -----------------------------------------------
147 */
148 void
SubMtx_scale3vec(SubMtx * mtxA,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[])149 SubMtx_scale3vec (
150 SubMtx *mtxA,
151 double y0[],
152 double y1[],
153 double y2[],
154 double x0[],
155 double x1[],
156 double x2[]
157 ) {
158 /*
159 ---------------
160 check the input
161 ---------------
162 */
163 if ( mtxA == NULL || y0 == NULL || y1 == NULL || y2 == NULL
164 || x0 == NULL || x1 == NULL || x2 == NULL ) {
165 fprintf(stderr,
166 "\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
167 "\n bad input\n", mtxA, y0, y1, y2, x0, x1, x2) ;
168 exit(-1) ;
169 }
170 if ( ! (SUBMTX_IS_REAL(mtxA) || SUBMTX_IS_COMPLEX(mtxA)) ) {
171 fprintf(stderr,
172 "\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
173 "\n bad type %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
174 mtxA, y0, y1, y2, x0, x1, x2, mtxA->type) ;
175 exit(-1) ;
176 }
177 switch ( mtxA->mode ) {
178 case SUBMTX_DIAGONAL :
179 diagonal_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
180 break ;
181 case SUBMTX_BLOCK_DIAGONAL_SYM :
182 block_diagonal_sym_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
183 break ;
184 case SUBMTX_BLOCK_DIAGONAL_HERM :
185 if ( ! SUBMTX_IS_COMPLEX(mtxA) ) {
186 fprintf(stderr,
187 "\n fatal error in SubMtx_scale3vec(%p,%p,%p,%p,%p,%p,%p)"
188 "\n hermitian matrix, type %d is not SPOOLES_COMPLEX\n",
189 mtxA, y0, y1, y2, x0, x1, x2, mtxA->type) ;
190 exit(-1) ;
191 }
192 block_diagonal_herm_scale3vec(mtxA, y0, y1, y2, x0, x1, x2) ;
193 break ;
194 default :
195 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
196 "\n matrix type not supported"
197 "\n must be SUBMTX_DIAGONAL,"
198 "\n or SUBMTX_BLOCK_DIAGONAL_SYM"
199 "\n or SUBMTX_BLOCK_DIAGONAL_HERM\n") ;
200 exit(-1) ;
201 }
202 return ; }
203
204 /*--------------------------------------------------------------------*/
205 /*
206 -----------------------------------------------
207 purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
208 where A is diagonal
209
210 created -- 98apr17, cca
211 -----------------------------------------------
212 */
213 static void
diagonal_scale3vec(SubMtx * mtxA,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[])214 diagonal_scale3vec (
215 SubMtx *mtxA,
216 double y0[],
217 double y1[],
218 double y2[],
219 double x0[],
220 double x1[],
221 double x2[]
222 ) {
223 double *entriesA ;
224 int nrowA ;
225
226 SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
227 if ( SUBMTX_IS_REAL(mtxA) ) {
228 double a ;
229 int irowA ;
230 for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
231 a = entriesA[irowA] ;
232 y0[irowA] = a*x0[irowA] ;
233 y1[irowA] = a*x1[irowA] ;
234 y2[irowA] = a*x2[irowA] ;
235 }
236 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
237 double ai, ar, xi0, xi1, xi2, xr0, xr1, xr2 ;
238 int iloc, irowA, rloc ;
239 for ( irowA = rloc = 0, iloc = 1 ;
240 irowA < nrowA ;
241 irowA++, rloc += 2, iloc += 2 ) {
242 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
243 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
244 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
245 ar = entriesA[rloc] ; ai = entriesA[iloc] ;
246 y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
247 y1[rloc] = ar*xr1 - ai*xi1 ; y1[iloc] = ar*xi1 + ai*xr1 ;
248 y2[rloc] = ar*xr2 - ai*xi2 ; y2[iloc] = ar*xi2 + ai*xr2 ;
249 }
250 }
251 return ; }
252
253 /*--------------------------------------------------------------------*/
254 /*
255 -----------------------------------------
256 purpose -- compute [y0 y1] := A * [x0 x1]
257 where A is diagonal
258
259 created -- 98apr17, cca
260 -----------------------------------------
261 */
262 static void
diagonal_scale2vec(SubMtx * mtxA,double y0[],double y1[],double x0[],double x1[])263 diagonal_scale2vec (
264 SubMtx *mtxA,
265 double y0[],
266 double y1[],
267 double x0[],
268 double x1[]
269 ) {
270 double *entriesA ;
271 int nrowA ;
272
273 SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
274 if ( SUBMTX_IS_REAL(mtxA) ) {
275 double a ;
276 int irowA ;
277 for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
278 a = entriesA[irowA] ;
279 y0[irowA] = a*x0[irowA] ;
280 y1[irowA] = a*x1[irowA] ;
281 }
282 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
283 double ai, ar, xi0, xi1, xr0, xr1 ;
284 int iloc, irowA, rloc ;
285 for ( irowA = rloc = 0, iloc = 1 ;
286 irowA < nrowA ;
287 irowA++, rloc += 2, iloc += 2 ) {
288 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
289 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
290 ar = entriesA[rloc] ; ai = entriesA[iloc] ;
291 y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
292 y1[rloc] = ar*xr1 - ai*xi1 ; y1[iloc] = ar*xi1 + ai*xr1 ;
293 }
294 }
295 return ; }
296
297 /*--------------------------------------------------------------------*/
298 /*
299 -----------------------------------
300 purpose -- compute [y0] := A * [x0]
301 where A is diagonal
302
303 created -- 98apr17, cca
304 -----------------------------------
305 */
306 static void
diagonal_scale1vec(SubMtx * mtxA,double y0[],double x0[])307 diagonal_scale1vec (
308 SubMtx *mtxA,
309 double y0[],
310 double x0[]
311 ) {
312 double *entriesA ;
313 int nrowA ;
314
315 SubMtx_diagonalInfo(mtxA, &nrowA, &entriesA) ;
316 if ( SUBMTX_IS_REAL(mtxA) ) {
317 double a ;
318 int irowA ;
319 for ( irowA = 0 ; irowA < nrowA ; irowA++ ) {
320 a = entriesA[irowA] ;
321 y0[irowA] = a*x0[irowA] ;
322 }
323 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
324 double ai, ar, xi0, xr0 ;
325 int iloc, irowA, rloc ;
326 for ( irowA = rloc = 0, iloc = 1 ;
327 irowA < nrowA ;
328 irowA++, rloc += 2, iloc += 2 ) {
329 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
330 ar = entriesA[rloc] ; ai = entriesA[iloc] ;
331 y0[rloc] = ar*xr0 - ai*xi0 ; y0[iloc] = ar*xi0 + ai*xr0 ;
332 }
333 }
334 return ; }
335
336 /*--------------------------------------------------------------------*/
337 /*
338 --------------------------------------------------
339 purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
340 where A is block diagonal and symmetric
341
342 created -- 98apr17, cca
343 --------------------------------------------------
344 */
345 static void
block_diagonal_sym_scale3vec(SubMtx * mtxA,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[])346 block_diagonal_sym_scale3vec (
347 SubMtx *mtxA,
348 double y0[],
349 double y1[],
350 double y2[],
351 double x0[],
352 double x1[],
353 double x2[]
354 ) {
355 double *entries ;
356 int nentA, nrowA ;
357 int *pivotsizes ;
358
359 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
360 if ( SUBMTX_IS_REAL(mtxA) ) {
361 int ipivot, irowA, kk ;
362
363 for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
364 if ( pivotsizes[ipivot] == 1 ) {
365 double a00, x00, x01, x02 ;
366
367 a00 = entries[kk] ;
368 x00 = x0[irowA] ;
369 x01 = x1[irowA] ;
370 x02 = x2[irowA] ;
371 y0[irowA] = a00*x00 ;
372 y1[irowA] = a00*x01 ;
373 y2[irowA] = a00*x02 ;
374 kk++, irowA++ ;
375 } else if ( pivotsizes[ipivot] == 2 ) {
376 double a00, a01, a11,
377 x00, x01, x02, x10, x11, x12 ;
378
379 a00 = entries[kk] ;
380 a01 = entries[kk+1] ;
381 a11 = entries[kk+2] ;
382 x00 = x0[irowA] ;
383 x01 = x1[irowA] ;
384 x02 = x2[irowA] ;
385 x10 = x0[irowA+1] ;
386 x11 = x1[irowA+1] ;
387 x12 = x2[irowA+1] ;
388 y0[irowA] = a00*x00 + a01*x10 ;
389 y1[irowA] = a00*x01 + a01*x11 ;
390 y2[irowA] = a00*x02 + a01*x12 ;
391 y0[irowA+1] = a01*x00 + a11*x10 ;
392 y1[irowA+1] = a01*x01 + a11*x11 ;
393 y2[irowA+1] = a01*x02 + a11*x12 ;
394 kk += 3, irowA += 2 ;
395 } else {
396 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
397 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
398 exit(-1) ;
399 }
400 }
401 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
402 int iloc, ipivot, irowA, kk, rloc ;
403
404 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
405 irowA < nrowA ;
406 ipivot++ ) {
407 if ( pivotsizes[ipivot] == 1 ) {
408 double ai00, ar00, xi0, xi1, xi2, xr0, xr1, xr2 ;
409
410 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
411 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
412 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
413 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
414 y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
415 y1[rloc] = ar00*xr1 - ai00*xi1; y1[iloc] = ar00*xi1 + ai00*xr1;
416 y2[rloc] = ar00*xr2 - ai00*xi2; y2[iloc] = ar00*xi2 + ai00*xr2;
417 kk += 2, irowA++, rloc += 2, iloc += 2 ;
418 } else if ( pivotsizes[ipivot] == 2 ) {
419 double ai00, ai01, ai11, ar00, ar01, ar11,
420 xi00, xi01, xi02, xi10, xi11, xi12,
421 xr00, xr01, xr02, xr10, xr11, xr12 ;
422 int iloc0, iloc1, rloc0, rloc1 ;
423
424 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
425 ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
426 ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
427 rloc0 = rloc ; iloc0 = iloc ;
428 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
429 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
430 xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
431 xr02 = x2[rloc0] ; xi02 = x2[iloc0] ;
432 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
433 xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
434 xr12 = x2[rloc1] ; xi12 = x2[iloc1] ;
435 y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
436 y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
437 y1[rloc0] = ar00*xr01 - ai00*xi01 + ar01*xr11 - ai01*xi11 ;
438 y1[iloc0] = ar00*xi01 + ai00*xr01 + ar01*xi11 + ai01*xr11 ;
439 y2[rloc0] = ar00*xr02 - ai00*xi02 + ar01*xr12 - ai01*xi12 ;
440 y2[iloc0] = ar00*xi02 + ai00*xr02 + ar01*xi12 + ai01*xr12 ;
441 y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
442 y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
443 y1[rloc1] = ar01*xr01 - ai01*xi01 + ar11*xr11 - ai11*xi11 ;
444 y1[iloc1] = ar01*xi01 + ai01*xr01 + ar11*xi11 + ai11*xr11 ;
445 y2[rloc1] = ar01*xr02 - ai01*xi02 + ar11*xr12 - ai11*xi12 ;
446 y2[iloc1] = ar01*xi02 + ai01*xr02 + ar11*xi12 + ai11*xr12 ;
447 kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
448 } else {
449 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
450 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
451 exit(-1) ;
452 }
453 }
454 }
455 return ; }
456
457 /*--------------------------------------------------------------------*/
458 /*
459 --------------------------------------------------
460 purpose -- compute [y0 y1] := A * [x0 x1]
461 where A is block diagonal and symmetric
462
463 created -- 98apr17, cca
464 --------------------------------------------------
465 */
466 static void
block_diagonal_sym_scale2vec(SubMtx * mtxA,double y0[],double y1[],double x0[],double x1[])467 block_diagonal_sym_scale2vec (
468 SubMtx *mtxA,
469 double y0[],
470 double y1[],
471 double x0[],
472 double x1[]
473 ) {
474 double *entries ;
475 int nentA, nrowA ;
476 int *pivotsizes ;
477
478 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
479 if ( SUBMTX_IS_REAL(mtxA) ) {
480 int ipivot, irowA, kk ;
481
482 for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
483 if ( pivotsizes[ipivot] == 1 ) {
484 double a00, x00, x01 ;
485
486 a00 = entries[kk] ;
487 x00 = x0[irowA] ;
488 x01 = x1[irowA] ;
489 y0[irowA] = a00*x00 ;
490 y1[irowA] = a00*x01 ;
491 kk++, irowA++ ;
492 } else if ( pivotsizes[ipivot] == 2 ) {
493 double a00, a01, a11, x00, x01, x10, x11 ;
494
495 a00 = entries[kk] ;
496 a01 = entries[kk+1] ;
497 a11 = entries[kk+2] ;
498 x00 = x0[irowA] ;
499 x01 = x1[irowA] ;
500 x10 = x0[irowA+1] ;
501 x11 = x1[irowA+1] ;
502 y0[irowA] = a00*x00 + a01*x10 ;
503 y1[irowA] = a00*x01 + a01*x11 ;
504 y0[irowA+1] = a01*x00 + a11*x10 ;
505 y1[irowA+1] = a01*x01 + a11*x11 ;
506 kk += 3, irowA += 2 ;
507 } else {
508 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
509 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
510 exit(-1) ;
511 }
512 }
513 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
514 int iloc, ipivot, irowA, kk, rloc ;
515
516 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
517 irowA < nrowA ;
518 ipivot++ ) {
519 if ( pivotsizes[ipivot] == 1 ) {
520 double ai00, ar00, xi0, xi1, xr0, xr1 ;
521
522 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
523 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
524 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
525 y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
526 y1[rloc] = ar00*xr1 - ai00*xi1; y1[iloc] = ar00*xi1 + ai00*xr1;
527 kk += 2, irowA++, rloc += 2, iloc += 2 ;
528 } else if ( pivotsizes[ipivot] == 2 ) {
529 double ai00, ai01, ai11, ar00, ar01, ar11,
530 xi00, xi01, xi10, xi11, xr00, xr01, xr10, xr11 ;
531 int iloc0, iloc1, rloc0, rloc1 ;
532
533 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
534 ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
535 ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
536 rloc0 = rloc ; iloc0 = iloc ;
537 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
538 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
539 xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
540 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
541 xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
542 y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
543 y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
544 y1[rloc0] = ar00*xr01 - ai00*xi01 + ar01*xr11 - ai01*xi11 ;
545 y1[iloc0] = ar00*xi01 + ai00*xr01 + ar01*xi11 + ai01*xr11 ;
546 y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
547 y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
548 y1[rloc1] = ar01*xr01 - ai01*xi01 + ar11*xr11 - ai11*xi11 ;
549 y1[iloc1] = ar01*xi01 + ai01*xr01 + ar11*xi11 + ai11*xr11 ;
550 kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
551 } else {
552 fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
553 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
554 exit(-1) ;
555 }
556 }
557 }
558 return ; }
559
560 /*--------------------------------------------------------------------*/
561 /*
562 --------------------------------------------------
563 purpose -- compute [y0] := A * [x0]
564 where A is block diagonal and symmetric
565
566 created -- 98apr17, cca
567 --------------------------------------------------
568 */
569 static void
block_diagonal_sym_scale1vec(SubMtx * mtxA,double y0[],double x0[])570 block_diagonal_sym_scale1vec (
571 SubMtx *mtxA,
572 double y0[],
573 double x0[]
574 ) {
575 double *entries ;
576 int nentA, nrowA ;
577 int *pivotsizes ;
578
579 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
580 if ( SUBMTX_IS_REAL(mtxA) ) {
581 int ipivot, irowA, kk ;
582
583 for ( irowA = ipivot = kk = 0 ; irowA < nrowA ; ipivot++ ) {
584 if ( pivotsizes[ipivot] == 1 ) {
585 double a00, x00 ;
586
587 a00 = entries[kk] ;
588 x00 = x0[irowA] ;
589 y0[irowA] = a00*x00 ;
590 kk++, irowA++ ;
591 } else if ( pivotsizes[ipivot] == 2 ) {
592 double a00, a01, a11, x00, x10 ;
593
594 a00 = entries[kk] ;
595 a01 = entries[kk+1] ;
596 a11 = entries[kk+2] ;
597 x00 = x0[irowA] ;
598 x10 = x0[irowA+1] ;
599 y0[irowA] = a00*x00 + a01*x10 ;
600 y0[irowA+1] = a01*x00 + a11*x10 ;
601 kk += 3, irowA += 2 ;
602 } else {
603 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
604 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
605 exit(-1) ;
606 }
607 }
608 } else if ( SUBMTX_IS_COMPLEX(mtxA) ) {
609 int iloc, ipivot, irowA, kk, rloc ;
610
611 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
612 irowA < nrowA ;
613 ipivot++ ) {
614 if ( pivotsizes[ipivot] == 1 ) {
615 double ai00, ar00, xi0, xr0 ;
616
617 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
618 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
619 y0[rloc] = ar00*xr0 - ai00*xi0; y0[iloc] = ar00*xi0 + ai00*xr0;
620 kk += 2, irowA++, rloc += 2, iloc += 2 ;
621 } else if ( pivotsizes[ipivot] == 2 ) {
622 double ai00, ai01, ai11, ar00, ar01, ar11,
623 xi00, xi10, xr00, xr10 ;
624 int iloc0, iloc1, rloc0, rloc1 ;
625
626 ar00 = entries[kk] ; ai00 = entries[kk+1] ;
627 ar01 = entries[kk+2] ; ai01 = entries[kk+3] ;
628 ar11 = entries[kk+4] ; ai11 = entries[kk+5] ;
629 rloc0 = rloc ; iloc0 = iloc ;
630 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
631 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
632 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
633 y0[rloc0] = ar00*xr00 - ai00*xi00 + ar01*xr10 - ai01*xi10 ;
634 y0[iloc0] = ar00*xi00 + ai00*xr00 + ar01*xi10 + ai01*xr10 ;
635 y0[rloc1] = ar01*xr00 - ai01*xi00 + ar11*xr10 - ai11*xi10 ;
636 y0[iloc1] = ar01*xi00 + ai01*xr00 + ar11*xi10 + ai11*xr10 ;
637 kk += 6, irowA += 2 ; rloc += 4 ; iloc += 4 ;
638 } else {
639 fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
640 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
641 exit(-1) ;
642 }
643 }
644 }
645 return ; }
646
647 /*--------------------------------------------------------------------*/
648 /*
649 --------------------------------------------------
650 purpose -- compute [y0 y1 y2] := A * [x0 x1 x2]
651 where A is block diagonal and hermitian
652
653 created -- 98apr17, cca
654 --------------------------------------------------
655 */
656 static void
block_diagonal_herm_scale3vec(SubMtx * mtxA,double y0[],double y1[],double y2[],double x0[],double x1[],double x2[])657 block_diagonal_herm_scale3vec (
658 SubMtx *mtxA,
659 double y0[],
660 double y1[],
661 double y2[],
662 double x0[],
663 double x1[],
664 double x2[]
665 ) {
666 double *entries ;
667 int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
668 int *pivotsizes ;
669
670 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
671 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
672 irowA < nrowA ;
673 ipivot++ ) {
674 if ( pivotsizes[ipivot] == 1 ) {
675 double ar00, ai00, xi0, xi1, xi2, xr0, xr1, xr2 ;
676
677 /*
678 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
679 */
680 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
681 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
682 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
683 xr2 = x2[rloc] ; xi2 = x2[iloc] ;
684 y0[rloc] = ar00*xr0 ; y0[iloc] = ar00*xi0 ;
685 y1[rloc] = ar00*xr1 ; y1[iloc] = ar00*xi1 ;
686 y2[rloc] = ar00*xr2 ; y2[iloc] = ar00*xi2 ;
687 irowA++, rloc += 2, iloc += 2 ;
688 } else if ( pivotsizes[ipivot] == 2 ) {
689 double ai00, ai01, ai11, ar00, ar01, ar11,
690 xi00, xi01, xi02, xi10, xi11, xi12,
691 xr00, xr01, xr02, xr10, xr11, xr12 ;
692 int iloc0, iloc1, rloc0, rloc1 ;
693
694 rloc0 = rloc ; iloc0 = iloc ;
695 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
696 /*
697 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
698 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
699 ar11 = entries[kk++] ; ai11 = entries[kk++] ;
700 */
701 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
702 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
703 ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
704 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
705 xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
706 xr02 = x2[rloc0] ; xi02 = x2[iloc0] ;
707 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
708 xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
709 xr12 = x2[rloc1] ; xi12 = x2[iloc1] ;
710 y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
711 y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
712 y1[rloc0] = ar00*xr01 + ar01*xr11 - ai01*xi11 ;
713 y1[iloc0] = ar00*xi01 + ar01*xi11 + ai01*xr11 ;
714 y2[rloc0] = ar00*xr02 + ar01*xr12 - ai01*xi12 ;
715 y2[iloc0] = ar00*xi02 + ar01*xi12 + ai01*xr12 ;
716 y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
717 y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
718 y1[rloc1] = ar01*xr01 + ai01*xi01 + ar11*xr11 ;
719 y1[iloc1] = ar01*xi01 - ai01*xr01 + ar11*xi11 ;
720 y2[rloc1] = ar01*xr02 + ai01*xi02 + ar11*xr12 ;
721 y2[iloc1] = ar01*xi02 - ai01*xr02 + ar11*xi12 ;
722 irowA += 2 ; rloc += 4 ; iloc += 4 ;
723 } else {
724 fprintf(stderr, "\n fatal error in SubMtx_scale3vec()"
725 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
726 exit(-1) ;
727 }
728 }
729 return ; }
730
731 /*--------------------------------------------------------------------*/
732 /*
733 --------------------------------------------------
734 purpose -- compute [y0 y1] := A * [x0 x1]
735 where A is block diagonal and hermitian
736
737 created -- 98apr17, cca
738 --------------------------------------------------
739 */
740 static void
block_diagonal_herm_scale2vec(SubMtx * mtxA,double y0[],double y1[],double x0[],double x1[])741 block_diagonal_herm_scale2vec (
742 SubMtx *mtxA,
743 double y0[],
744 double y1[],
745 double x0[],
746 double x1[]
747 ) {
748 double *entries ;
749 int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
750 int *pivotsizes ;
751
752 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
753 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
754 irowA < nrowA ;
755 ipivot++ ) {
756 if ( pivotsizes[ipivot] == 1 ) {
757 double ai00, ar00, xi0, xi1, xr0, xr1 ;
758
759 /*
760 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
761 */
762 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
763 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
764 xr1 = x1[rloc] ; xi1 = x1[iloc] ;
765 y0[rloc] = ar00*xr0 - ai00*xi0 ; y0[iloc] = ar00*xi0 + ai00*xr0 ;
766 y1[rloc] = ar00*xr1 - ai00*xi1 ; y1[iloc] = ar00*xi1 + ai00*xr1 ;
767 irowA++, rloc += 2, iloc += 2 ;
768 } else if ( pivotsizes[ipivot] == 2 ) {
769 double ai00, ai01, ai11, ar00, ar01, ar11,
770 xi00, xi01, xi10, xi11, xr00, xr01, xr10, xr11 ;
771 int iloc0, iloc1, rloc0, rloc1 ;
772
773 rloc0 = rloc ; iloc0 = iloc ;
774 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
775 /*
776 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
777 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
778 ar11 = entries[kk++] ; ai11 = entries[kk++] ;
779 */
780 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
781 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
782 ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
783 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
784 xr01 = x1[rloc0] ; xi01 = x1[iloc0] ;
785 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
786 xr11 = x1[rloc1] ; xi11 = x1[iloc1] ;
787 y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
788 y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
789 y1[rloc0] = ar00*xr01 + ar01*xr11 - ai01*xi11 ;
790 y1[iloc0] = ar00*xi01 + ar01*xi11 + ai01*xr11 ;
791 y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
792 y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
793 y1[rloc1] = ar01*xr01 + ai01*xi01 + ar11*xr11 ;
794 y1[iloc1] = ar01*xi01 - ai01*xr01 + ar11*xi11 ;
795 irowA += 2 ; rloc += 4 ; iloc += 4 ;
796 } else {
797 fprintf(stderr, "\n fatal error in SubMtx_scale2vec()"
798 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
799 exit(-1) ;
800 }
801 }
802 return ; }
803
804 /*--------------------------------------------------------------------*/
805 /*
806 --------------------------------------------------
807 purpose -- compute [y0] := A * [x0]
808 where A is block diagonal and hermitian
809
810 created -- 98apr17, cca
811 --------------------------------------------------
812 */
813 static void
block_diagonal_herm_scale1vec(SubMtx * mtxA,double y0[],double x0[])814 block_diagonal_herm_scale1vec (
815 SubMtx *mtxA,
816 double y0[],
817 double x0[]
818 ) {
819 double *entries ;
820 int iloc, ipivot, irowA, kk, nentA, nrowA, rloc ;
821 int *pivotsizes ;
822
823 SubMtx_blockDiagonalInfo(mtxA, &nrowA, &nentA, &pivotsizes, &entries) ;
824 for ( irowA = ipivot = kk = rloc = 0, iloc = 1 ;
825 irowA < nrowA ;
826 ipivot++ ) {
827 if ( pivotsizes[ipivot] == 1 ) {
828 double ai00, ar00, xi0, xr0 ;
829
830 /*
831 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
832 */
833 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
834 xr0 = x0[rloc] ; xi0 = x0[iloc] ;
835 y0[rloc] = ar00*xr0 - ai00*xi0 ; y0[iloc] = ar00*xi0 + ai00*xr0 ;
836 irowA++, rloc += 2, iloc += 2 ;
837 } else if ( pivotsizes[ipivot] == 2 ) {
838 double ai00, ai01, ai11, ar00, ar01, ar11,
839 xi00, xi10, xr00, xr10 ;
840 int iloc0, iloc1, rloc0, rloc1 ;
841
842 rloc0 = rloc ; iloc0 = iloc ;
843 rloc1 = rloc + 2 ; iloc1 = iloc + 2 ;
844 /*
845 ar00 = entries[kk++] ; ai00 = entries[kk++] ;
846 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
847 ar11 = entries[kk++] ; ai11 = entries[kk++] ;
848 */
849 ar00 = entries[kk++] ; ai00 = 0.0 ; kk++ ;
850 ar01 = entries[kk++] ; ai01 = entries[kk++] ;
851 ar11 = entries[kk++] ; ai11 = 0.0 ; kk++ ;
852 xr00 = x0[rloc0] ; xi00 = x0[iloc0] ;
853 xr10 = x0[rloc1] ; xi10 = x0[iloc1] ;
854 y0[rloc0] = ar00*xr00 + ar01*xr10 - ai01*xi10 ;
855 y0[iloc0] = ar00*xi00 + ar01*xi10 + ai01*xr10 ;
856 y0[rloc1] = ar01*xr00 + ai01*xi00 + ar11*xr10 ;
857 y0[iloc1] = ar01*xi00 - ai01*xr00 + ar11*xi10 ;
858 irowA += 2 ; rloc += 4 ; iloc += 4 ;
859 } else {
860 fprintf(stderr, "\n fatal error in SubMtx_scale1vec()"
861 "\n pivotsizes[%d] = %d", ipivot, pivotsizes[ipivot]) ;
862 exit(-1) ;
863 }
864 }
865 return ; }
866
867 /*--------------------------------------------------------------------*/
868