1 /* assemble.c */
2
3 #include "../Chv.h"
4
5 #define MYDEBUG 0
6
7 /*--------------------------------------------------------------------*/
8 /*
9 -----------------------------------------------------------------
10 add a scaled multiple of a simple chevron to a Chv object.
11 the indices are offsets.
12 note: for this purpose, (assembling original entries into the
13 matrix), the row and column indices of the chevron are identical.
14 also, the indices of both the Chv object and the chvind[]
15 vector are assumed to be in ascending order.
16
17 created -- 98apr30, cca
18 -----------------------------------------------------------------
19 */
20 void
Chv_addChevron(Chv * chv,double alpha[],int ichv,int chvsize,int chvind[],double chvent[])21 Chv_addChevron (
22 Chv *chv,
23 double alpha[],
24 int ichv,
25 int chvsize,
26 int chvind[],
27 double chvent[]
28 ) {
29 int ii, iloc, jcol, jj, jjfirst, jjlast,
30 ncol, nD, nL, nU, offset ;
31 int *colind ;
32 double *diag ;
33 /*
34 ---------------
35 check the input
36 ---------------
37 */
38 if ( chv == NULL || ichv < 0 || chvsize < 0
39 || chvind == NULL || chvent == NULL ) {
40 fprintf(stderr,
41 "\n fatal error in Chv_addChevron(%p,%p,%d,%d,%p,%p)"
42 "\n bad input\n",
43 chv, alpha, ichv, chvsize, chvind, chvent) ;
44 exit(-1) ;
45 }
46 switch ( chv->type ) {
47 case SPOOLES_REAL :
48 switch ( chv->symflag ) {
49 case SPOOLES_SYMMETRIC :
50 case SPOOLES_NONSYMMETRIC :
51 break ;
52 default :
53 fprintf(stderr, "\n fatal error in Chv_addChevron()"
54 "\n type is SPOOLES_REAL, symflag = %d"
55 "\n must be SPOOLES_SYMMETRIC or SPOOLES_NONSYMMETRIC\n",
56 chv->symflag) ;
57 exit(-1) ;
58 break ;
59 }
60 case SPOOLES_COMPLEX :
61 switch ( chv->symflag ) {
62 case SPOOLES_SYMMETRIC :
63 case SPOOLES_HERMITIAN :
64 case SPOOLES_NONSYMMETRIC :
65 break ;
66 default :
67 fprintf(stderr, "\n fatal error in Chv_addChevron()"
68 "\n type is SPOOLES_REAL, symflag = %d"
69 "\n must be SPOOLES_SYMMETRIC, SPOOLES_HERMITIAN"
70 "\n or SPOOLES_NONSYMMETRIC\n",
71 chv->symflag) ;
72 exit(-1) ;
73 break ;
74 }
75 break ;
76 default :
77 fprintf(stderr, "\n fatal error in Chv_addChevron()"
78 "\n type is %d, must be SPOOLES_REAL or SPOOLES_COMPLEX\n",
79 chv->type) ;
80 exit(-1) ;
81 break ;
82 }
83 #if MYDEBUG > 0
84 fprintf(stdout, "\n alpha = %f, ichv = %d, chvsize = %d",
85 alpha, ichv, chvsize) ;
86 #endif
87 if ( chvsize == 0
88 || (CHV_IS_REAL(chv) && alpha[0] == 0.0)
89 || (CHV_IS_COMPLEX(chv) && (alpha[0] == 0.0 && alpha[1] == 0.0)) ) {
90 /*
91 ----------------------------
92 quick return, nothing to add
93 ----------------------------
94 */
95 return ;
96 }
97 #if MYDEBUG > 0
98 fprintf(stdout, "\n\n Chv_addChevron(%d): ", ichv) ;
99 IVfprintf(stdout, chvsize, chvind) ;
100 DVfprintf(stdout, chvsize, chvent) ;
101 fflush(stdout) ;
102 #endif
103 Chv_dimensions(chv, &nD, &nL, &nU) ;
104 Chv_columnIndices(chv, &ncol, &colind) ;
105 /*
106 -------------------------------------
107 locate the chevron in the Chv object
108 that will accumulate these entries
109 -------------------------------------
110 */
111 for ( iloc = 0 ; iloc < nD ; iloc++ ) {
112 if ( colind[iloc] == ichv ) {
113 break ;
114 }
115 }
116 if ( iloc == nD ) {
117 /*
118 --------------------------------------------
119 unable to assemble these entries, error exit
120 --------------------------------------------
121 */
122 fprintf(stderr, "\n fatal error in Chv_addChevron(%p,%d,%d,%p,%p)"
123 "\n chevron id %d not found in colind[]",
124 chv, ichv, chvsize, chvind, chvent, ichv) ;
125 exit(-1) ;
126 }
127 if ( CHV_IS_SYMMETRIC(chv) || CHV_IS_HERMITIAN(chv) ) {
128 /*
129 ---------------------
130 symmetric chevron
131 get the local indices
132 ---------------------
133 */
134 jjlast = nD + nU - 1 ;
135 for ( ii = 0, jj = iloc ; ii < chvsize ; ii++ ) {
136 if ( (offset = chvind[ii]) < 0 ) {
137 fprintf(stderr,
138 "\n fatal error in Chv_addChevron(%p,%d,%d,%p,%p)"
139 "\n ii %d, negative offset %d\n",
140 chv, ichv, chvsize, chvind, chvent, ii, chvind[ii]) ;
141 IVfprintf(stderr, chvsize, chvind) ;
142 exit(-1) ;
143 }
144 jcol = ichv + offset ;
145 #if MYDEBUG > 0
146 fprintf(stdout, "\n ii = %d, offset = %d, jcol = %d",
147 ii, offset, jcol) ;
148 fflush(stdout) ;
149 #endif
150 while ( jj <= jjlast && jcol != colind[jj] ) {
151 jj++ ;
152 }
153 #if MYDEBUG > 0
154 fprintf(stdout, ", jj = %d", jj) ;
155 fflush(stdout) ;
156 #endif
157 if ( jj > jjlast ) {
158 fprintf(stderr,
159 "\n fatal error in Chv_addChevron(%p,%d,%d,%p,%p)"
160 "\n jcol %d not found in colind[]\n",
161 chv, ichv, chvsize, chvind, chvent, jcol) ;
162 fprintf(stderr, "\n colind") ;
163 IVfprintf(stderr, ncol, colind) ;
164 fprintf(stderr, "\n chvind") ;
165 IVfprintf(stderr, chvsize, chvind) ;
166 exit(-1) ;
167 }
168 chvind[ii] = jj ;
169 }
170 #if MYDEBUG > 0
171 fprintf(stdout, "\n local indices") ;
172 IVfprintf(stdout, chvsize, chvind) ;
173 fflush(stdout) ;
174 #endif
175 /*
176 --------------------
177 assemble the chevron
178 --------------------
179 */
180 if ( CHV_IS_REAL(chv) ) {
181 diag = Chv_diagLocation(chv, iloc) - iloc ;
182 #if MYDEBUG > 0
183 fprintf(stdout, "\n ichv = %d, iloc = %d, diag = %p"
184 "\n chv->entries = %p, diag - chv->entries = %d",
185 ichv, iloc, diag, chv->entries,
186 diag - chv->entries) ;
187 fflush(stdout) ;
188 #endif
189 if ( alpha[0] == 1.0 ) {
190 for ( ii = 0 ; ii < chvsize ; ii++ ) {
191 #if MYDEBUG > 0
192 fprintf(stdout, "\n location %d",
193 &diag[chvind[ii]] - chv->entries) ;
194 fflush(stdout) ;
195 #endif
196 diag[chvind[ii]] += chvent[ii] ;
197 }
198 } else {
199 for ( ii = 0 ; ii < chvsize ; ii++ ) {
200 diag[chvind[ii]] += alpha[0]*chvent[ii] ;
201 }
202 }
203 } else if ( CHV_IS_COMPLEX(chv) ) {
204 diag = Chv_diagLocation(chv, iloc) - 2*iloc ;
205 #if MYDEBUG > 0
206 fprintf(stdout, "\n ichv = %d, iloc = %d, diag = %p"
207 "\n chv->entries = %p, diag - chv->entries = %d",
208 ichv, iloc, diag, chv->entries,
209 diag - chv->entries) ;
210 fflush(stdout) ;
211 #endif
212 if ( alpha[0] == 1.0 && alpha[1] == 0.0 ) {
213 for ( ii = 0 ; ii < chvsize ; ii++ ) {
214 #if MYDEBUG > 0
215 fprintf(stdout, "\n location %d",
216 &diag[chvind[ii]] - chv->entries) ;
217 fflush(stdout) ;
218 #endif
219 diag[2*chvind[ii]] += chvent[2*ii] ;
220 diag[2*chvind[ii]+1] += chvent[2*ii+1] ;
221 }
222 } else if ( alpha[0] != 0.0 && alpha[1] == 0.0 ) {
223 for ( ii = 0 ; ii < chvsize ; ii++ ) {
224 diag[2*chvind[ii]] += alpha[0]*chvent[2*ii] ;
225 diag[2*chvind[ii]+1] += alpha[0]*chvent[2*ii+1] ;
226 }
227 } else if ( CHV_IS_SYMMETRIC(chv) ) {
228 double alphareal, alphaimag, xreal, ximag ;
229 for ( ii = 0 ; ii < chvsize ; ii++ ) {
230 alphareal = alpha[0] ; alphaimag = alpha[1] ;
231 xreal = chvent[2*ii] ; ximag = chvent[2*ii+1] ;
232 diag[2*chvind[ii]] += alphareal*xreal - alphaimag*ximag ;
233 diag[2*chvind[ii]+1] += alphareal*ximag + alphaimag*xreal ;
234 }
235 } else {
236 fprintf(stderr, "\n fatal error in Chv_addChevron()"
237 "\n chevron is hermitian, but the scalar has nonzero imaginary part"
238 "\n sum is no longer hermitian\n") ;
239 exit(-1) ;
240 }
241 }
242 /*
243 ------------------------------------
244 restore the indices to their offsets
245 ------------------------------------
246 */
247 for ( ii = 0 ; ii < chvsize ; ii++ ) {
248 chvind[ii] = colind[chvind[ii]] - ichv ;
249 }
250 #if MYDEBUG > 0
251 fprintf(stdout, "\n restored indices") ;
252 IVfprintf(stdout, chvsize, chvind) ;
253 fflush(stdout) ;
254 #endif
255 } else if ( CHV_IS_NONSYMMETRIC(chv) ) {
256 /*
257 -----------------------------------------
258 nonsymmetric chevron, symmetric structure
259 overwrite chvind[] with local indices
260 -----------------------------------------
261 */
262 jjfirst = iloc ;
263 jjlast = nD + nU - 1 ;
264 for ( ii = 0, jj = jjlast ; ii < chvsize ; ii++ ) {
265 if ( (offset = chvind[ii]) >= 0 ) {
266 break ;
267 }
268 jcol = ichv - offset ;
269 while ( jj >= jjfirst && jcol != colind[jj] ) {
270 jj-- ;
271 }
272 if ( jj < jjfirst ) {
273 fprintf(stderr,
274 "\n fatal error in Chv_addChevron(%p,%d,%d,%p,%p)"
275 "\n jcol %d not found in colind[]\n",
276 chv, ichv, chvsize, chvind, chvent, jcol) ;
277 exit(-1) ;
278 }
279 chvind[ii] = -jj + iloc ;
280 }
281 for ( jj = jjfirst ; ii < chvsize ; ii++ ) {
282 jcol = ichv + chvind[ii] ;
283 while ( jj <= jjlast && jcol != colind[jj] ) {
284 jj++ ;
285 }
286 if ( jj > jjlast ) {
287 fprintf(stderr,
288 "\n fatal error in Chv_addChevron(%p,%d,%d,%p,%p)"
289 "\n jcol %d not found in colind[]\n",
290 chv, ichv, chvsize, chvind, chvent, jcol) ;
291 exit(-1) ;
292 }
293 chvind[ii] = jj - iloc ;
294 }
295 #if MYDEBUG > 0
296 fprintf(stdout, "\n local indices") ;
297 IVfprintf(stdout, chvsize, chvind) ;
298 fflush(stdout) ;
299 #endif
300 /*
301 --------------------
302 assemble the chevron
303 --------------------
304 */
305 diag = Chv_diagLocation(chv, iloc) ;
306 #if MYDEBUG > 0
307 fprintf(stdout, "\n ichv = %d, iloc = %d, diag = %p"
308 "\n chv->entries = %p, diag - chv->entries = %d",
309 ichv, iloc, diag, chv->entries,
310 diag - chv->entries) ;
311 fflush(stdout) ;
312 #endif
313 if ( CHV_IS_REAL(chv) ) {
314 if ( alpha[0] == 1.0 ) {
315 for ( ii = 0 ; ii < chvsize ; ii++ ) {
316 #if MYDEBUG > 0
317 fprintf(stdout, "\n diag[%d] += %12.4e, ii = %d",
318 chvind[ii], chvent[ii], ii) ;
319 fflush(stdout) ;
320 #endif
321 diag[chvind[ii]] += chvent[ii] ;
322 }
323 } else {
324 for ( ii = 0 ; ii < chvsize ; ii++ ) {
325 #if MYDEBUG > 0
326 fprintf(stdout, "\n diag[%d] += %12.4e, ii = %d",
327 chvind[ii], chvent[ii], ii) ;
328 fflush(stdout) ;
329 #endif
330 diag[chvind[ii]] += alpha[0] * chvent[ii] ;
331 }
332 }
333 } else if ( CHV_IS_COMPLEX(chv) ) {
334 if ( alpha[0] == 1.0 && alpha[1] == 0.0 ) {
335 for ( ii = 0 ; ii < chvsize ; ii++ ) {
336 #if MYDEBUG > 0
337 fprintf(stdout, "\n diag[%d] += %12.4e, ii = %d",
338 chvind[ii], chvent[ii], ii) ;
339 fflush(stdout) ;
340 #endif
341 diag[2*chvind[ii]] += chvent[2*ii] ;
342 diag[2*chvind[ii]+1] += chvent[2*ii+1] ;
343 }
344 } else if ( alpha[0] != 1.0 && alpha[1] == 0.0 ) {
345 for ( ii = 0 ; ii < chvsize ; ii++ ) {
346 diag[2*chvind[ii]] += alpha[0] * chvent[2*ii] ;
347 diag[2*chvind[ii]+1] += alpha[0] * chvent[2*ii+1] ;
348 }
349 } else {
350 double alphareal, alphaimag, xreal, ximag ;
351 for ( ii = 0 ; ii < chvsize ; ii++ ) {
352 alphareal = alpha[0] ; alphaimag = alpha[1] ;
353 xreal = chvent[2*ii] ; ximag = chvent[2*ii+1] ;
354 diag[2*chvind[ii]] += alphareal*xreal - alphaimag*ximag ;
355 diag[2*chvind[ii]+1] += alphareal*ximag + alphaimag*xreal ;
356 }
357 }
358 }
359 /*
360 ------------------------------------
361 restore the indices to their offsets
362 ------------------------------------
363 */
364 for ( ii = 0 ; ii < chvsize ; ii++ ) {
365 if ( chvind[ii] < 0 ) {
366 chvind[ii] = ichv - colind[iloc - chvind[ii]] ;
367 } else {
368 chvind[ii] = colind[chvind[ii] + iloc] - ichv ;
369 }
370 }
371 #if MYDEBUG > 0
372 fprintf(stdout, "\n restored indices") ;
373 IVfprintf(stdout, chvsize, chvind) ;
374 fflush(stdout) ;
375 #endif
376 }
377 return ; }
378
379 /*--------------------------------------------------------------------*/
380 /*
381 --------------------------------------------------------------
382 assemble Chv object chvI into Chv object chvJ.
383 note: the two objects must be of the same symmetry type,
384 the row indices of chvI must nest into those of chvJ,
385 the column indices of chvI must nest into those of chvJ.
386
387 created -- 98apr30, cca
388 --------------------------------------------------------------
389 */
390 void
Chv_assembleChv(Chv * chvJ,Chv * chvI)391 Chv_assembleChv (
392 Chv *chvJ,
393 Chv *chvI
394 ) {
395 double *diagI, *diagJ ;
396 int ii, ichvI, ichvJ, jj, ncolI, ncolJ, nDI, nDJ,
397 nLI, nLJ, nrowI, nrowJ, nUI, nUJ, offset ;
398 int *colindJ, *colindI, *rowindI, *rowindJ ;
399 /*
400 ---------------
401 check the input
402 ---------------
403 */
404 if ( chvJ == NULL || chvI == NULL ) {
405 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
406 "\n bad input\n", chvJ, chvI) ;
407 exit(-1) ;
408 }
409 if ( !(CHV_IS_SYMMETRIC(chvI)
410 || CHV_IS_HERMITIAN(chvI)
411 || CHV_IS_NONSYMMETRIC(chvI) ) ) {
412 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
413 "\n bad symflag %d\n", chvJ, chvI, chvI->symflag) ;
414 exit(-1) ;
415 }
416 if ( chvI->symflag != chvJ->symflag ) {
417 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
418 "\n chvI->symflag = %d, chvJ->symflag = %d\n",
419 chvJ, chvI, chvI->symflag, chvJ->symflag) ;
420 exit(-1) ;
421 }
422 /*
423 -------------------------------------
424 get the dimensions of the two objects
425 -------------------------------------
426 */
427 Chv_dimensions(chvJ, &nDJ, &nLJ, &nUJ) ;
428 Chv_dimensions(chvI, &nDI, &nLI, &nUI) ;
429 if ( nDI + nLI > nDJ + nLJ || nDI + nUI > nDJ + nUJ ) {
430 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
431 "\n bad dimensions"
432 "\n nDI = %d, nLI = %d, nUI = %d"
433 "\n nDI = %d, nLI = %d, nUI = %d",
434 chvJ, chvI, nDI, nLI, nUI, nDJ, nLJ, nUJ) ;
435 exit(-1) ;
436 }
437 /*
438 -----------------
439 get the local ids
440 -----------------
441 */
442 Chv_columnIndices(chvJ, &ncolJ, &colindJ) ;
443 Chv_columnIndices(chvI, &ncolI, &colindI) ;
444 #if MYDEBUG > 0
445 fprintf(stdout, "\n colindI") ;
446 IVfprintf(stdout, ncolI, colindI) ;
447 fprintf(stdout, "\n colindJ") ;
448 IVfprintf(stdout, ncolJ, colindJ) ;
449 #endif
450 for ( ii = 0, jj = 0 ; ii < ncolI ; ii++ ) {
451 while ( jj < ncolJ && colindI[ii] != colindJ[jj] ) {
452 #if MYDEBUG > 0
453 fprintf(stdout, "\n colindI[%d] = %d, colindJ[%d] = %d",
454 ii, colindI[ii], jj, colindJ[jj]) ;
455 #endif
456 jj++ ;
457 }
458 if ( jj == ncolJ ) {
459 break ;
460 }
461 colindI[ii] = jj ;
462 }
463 #if MYDEBUG > 0
464 fprintf(stdout, "\n local column indices") ;
465 IVfprintf(stdout, ncolI, colindI) ;
466 #endif
467 if ( jj == ncolJ ) {
468 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
469 "\n column indicesI do not nest in indicesJ\n", chvJ, chvI) ;
470 fprintf(stderr, "\n colindI") ;
471 IVfprintf(stderr, ncolI, colindI) ;
472 fprintf(stderr, "\n colindJ") ;
473 IVfprintf(stderr, ncolJ, colindJ) ;
474 exit(-1) ;
475 }
476 if ( CHV_IS_SYMMETRIC(chvJ) || CHV_IS_HERMITIAN(chvJ) ) {
477 /*
478 -------------------
479 symmetric structure
480 -------------------
481 */
482 nrowI = ncolI ;
483 rowindI = colindI ;
484 } else if ( CHV_IS_NONSYMMETRIC(chvJ) ) {
485 /*
486 ----------------------
487 nonsymmetric structure
488 ----------------------
489 */
490 Chv_rowIndices(chvJ, &nrowJ, &rowindJ) ;
491 Chv_rowIndices(chvI, &nrowI, &rowindI) ;
492 #if MYDEBUG > 0
493 fprintf(stdout, "\n rowindI") ;
494 IVfprintf(stdout, nrowI, rowindI) ;
495 fprintf(stdout, "\n rowindJ") ;
496 IVfprintf(stdout, nrowJ, rowindJ) ;
497 #endif
498 for ( ii = 0, jj = 0 ; ii < nrowI ; ii++ ) {
499 while ( jj < nrowJ && rowindI[ii] != rowindJ[jj] ) {
500 jj++ ;
501 }
502 if ( jj == nrowJ ) {
503 break ;
504 }
505 rowindI[ii] = jj ;
506 }
507 if ( jj == nrowJ ) {
508 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
509 "\n row indicesI do not nest in indicesJ\n", chvJ, chvI) ;
510 fprintf(stderr, "\n rowindI") ;
511 IVfprintf(stderr, nrowI, rowindI) ;
512 fprintf(stderr, "\n rowindJ") ;
513 IVfprintf(stderr, nrowJ, rowindJ) ;
514 exit(-1) ;
515 }
516 #if MYDEBUG > 0
517 fprintf(stdout, "\n local row indices") ;
518 IVfprintf(stdout, nrowI, rowindI) ;
519 #endif
520 }
521 #if MYDEBUG > 0
522 fprintf(stdout, "\n local column indices") ;
523 IVfprintf(stdout, ncolI, colindI) ;
524 fprintf(stdout, "\n local row indices") ;
525 IVfprintf(stdout, nrowI, rowindI) ;
526 #endif
527 /*
528 ---------------------------
529 loop over the chevrons in I
530 ---------------------------
531 */
532 for ( ichvI = 0 ; ichvI < nDI ; ichvI++ ) {
533 ichvJ = colindI[ichvI] ;
534 if ( ichvJ != rowindI[ichvI] ) {
535 fprintf(stderr, "\n fatal error in Chv_assembleChv(%p,%p)"
536 "\n ichvI = %d, ichvJ = %d, rowindI[ichvI] = %d",
537 chvJ, chvI, ichvI, ichvJ, rowindI[ichvI]) ;
538 exit(-1) ;
539 }
540 diagI = Chv_diagLocation(chvI, ichvI) ;
541 diagJ = Chv_diagLocation(chvJ, ichvJ) ;
542 #if MYDEBUG > 0
543 fprintf(stdout,
544 "\n ichvI = %d, diagI - entriesI = %d"
545 "\n ichvJ = %d, diagJ - entriesJ = %d",
546 ichvI, diagI - chvI->entries, ichvJ, diagJ - chvJ->entries) ;
547 #endif
548 if ( CHV_IS_REAL(chvJ) ) {
549 diagJ[0] += diagI[0] ;
550 } else if ( CHV_IS_COMPLEX(chvJ) ) {
551 diagJ[0] += diagI[0] ;
552 diagJ[1] += diagI[1] ;
553 }
554 if ( CHV_IS_SYMMETRIC(chvJ) || CHV_IS_HERMITIAN(chvJ) ) {
555 /*
556 ----------------------
557 symmetric or hermitian
558 ----------------------
559 */
560 if ( CHV_IS_REAL(chvJ) ) {
561 for ( ii = ichvI + 1, jj = 1 ; ii < ncolI ; ii++, jj++ ) {
562 offset = colindI[ii] - ichvJ ;
563 #if MYDEBUG > 0
564 fprintf(stdout,
565 "\n ii = %d, jj = %d, offset = %d", ii, jj, offset) ;
566 #endif
567 diagJ[offset] += diagI[jj] ;
568 }
569 } else if ( CHV_IS_COMPLEX(chvJ) ) {
570 for ( ii = ichvI + 1, jj = 1 ; ii < ncolI ; ii++, jj++ ) {
571 offset = colindI[ii] - ichvJ ;
572 #if MYDEBUG > 0
573 fprintf(stdout,
574 "\n ii = %d, jj = %d, offset = %d", ii, jj, offset) ;
575 #endif
576 diagJ[2*offset] += diagI[2*jj] ;
577 diagJ[2*offset+1] += diagI[2*jj+1] ;
578 }
579 }
580 } else if ( CHV_IS_NONSYMMETRIC(chvJ) ) {
581 /*
582 ------------
583 nonsymmetric
584 ------------
585 */
586 if ( CHV_IS_REAL(chvJ) ) {
587 for ( ii = ichvI + 1, jj = 1 ; ii < ncolI ; ii++, jj++ ) {
588 offset = colindI[ii] - ichvJ ;
589 #if MYDEBUG > 0
590 fprintf(stdout, "\n diagJ[%d] += diagI[%d]", offset, jj) ;
591 #endif
592 diagJ[offset] += diagI[jj] ;
593 }
594 for ( ii = ichvI + 1, jj = -1 ; ii < nrowI ; ii++, jj-- ) {
595 offset = rowindI[ii] - ichvJ ;
596 #if MYDEBUG > 0
597 fprintf(stdout, "\n diagJ[%d] += diagI[%d]", -offset, jj) ;
598 #endif
599 diagJ[-offset] += diagI[jj] ;
600 }
601 } else if ( CHV_IS_COMPLEX(chvJ) ) {
602 for ( ii = ichvI + 1, jj = 1 ; ii < ncolI ; ii++, jj++ ) {
603 offset = colindI[ii] - ichvJ ;
604 #if MYDEBUG > 0
605 fprintf(stdout, "\n diagJ[%d] += diagI[%d]", offset, jj) ;
606 #endif
607 diagJ[2*offset] += diagI[2*jj] ;
608 diagJ[2*offset+1] += diagI[2*jj+1] ;
609 }
610 for ( ii = ichvI + 1, jj = -1 ; ii < nrowI ; ii++, jj-- ) {
611 offset = rowindI[ii] - ichvJ ;
612 #if MYDEBUG > 0
613 fprintf(stdout, "\n diagJ[%d] += diagI[%d]", -offset, jj) ;
614 #endif
615 diagJ[-2*offset] += diagI[2*jj] ;
616 diagJ[-2*offset+1] += diagI[2*jj+1] ;
617 }
618 }
619 }
620 }
621 /*
622 -------------------
623 restore the indices
624 -------------------
625 */
626 for ( ii = 0 ; ii < ncolI ; ii++ ) {
627 colindI[ii] = colindJ[colindI[ii]] ;
628 }
629 if ( CHV_IS_NONSYMMETRIC(chvJ) ) {
630 for ( ii = 0 ; ii < nrowI ; ii++ ) {
631 rowindI[ii] = rowindJ[rowindI[ii]] ;
632 }
633 }
634 return ; }
635
636 /*--------------------------------------------------------------------*/
637 /*
638 ----------------------------------------------------------------
639 purpose -- assemble the postponed data from the children
640
641
642 newchv -- Chv object to contain fully assembled front
643 oldchv -- Chv object that contains former front
644 firstchild -- pointer to first child in the list of children
645 Chv objects to be merged into the new front
646
647 return value -- # of delayed rows and columns added to the front
648
649 created -- 98apr30, cca
650 ----------------------------------------------------------------
651 */
652 int
Chv_assemblePostponedData(Chv * newchv,Chv * oldchv,Chv * firstchild)653 Chv_assemblePostponedData (
654 Chv *newchv,
655 Chv *oldchv,
656 Chv *firstchild
657 ) {
658 Chv *child ;
659 int ierr, ndelay ;
660 /*
661 ---------------
662 check the input
663 ---------------
664 */
665 if ( newchv == NULL || oldchv == NULL || firstchild == NULL ) {
666 fprintf(stderr,
667 "\n fatal error in Chv_assemblePostponedData(%p,%p,%p)"
668 "\n bad input\n", newchv, oldchv, firstchild) ;
669 exit(-1) ;
670 }
671 Chv_zero(newchv) ;
672 #if MYDEBUG > 0
673 fprintf(stdout, "\n\n inside Chv_assemblePostponedData, J = %d",
674 oldchv->id) ;
675 fprintf(stdout, "\n\n old chevron %d", oldchv->id) ;
676 Chv_writeForHumanEye(oldchv, stdout) ;
677 for ( child = firstchild ; child != NULL ; child = child->next ) {
678 fprintf(stdout, "\n\n child %d", child->id) ;
679 Chv_writeForHumanEye(child, stdout) ;
680 }
681 fprintf(stdout, "\n\n new chevron %d", newchv->id) ;
682 Chv_writeForHumanEye(newchv, stdout) ;
683 fflush(stdout) ;
684 #endif
685 /*
686 ------------------------------------
687 copy the indices into the new object
688 ------------------------------------
689 */
690 ndelay = 0 ;
691 for ( child = firstchild ; child != NULL ; child = child->next ) {
692 IVcopy(child->nD, newchv->colind + ndelay, child->colind) ;
693 ndelay += child->nD ;
694 }
695 IVcopy(oldchv->nD + oldchv->nU,
696 newchv->colind + ndelay, oldchv->colind) ;
697 #if MYDEBUG > 0
698 fprintf(stdout, "\n ndelay = %d, new column indices", ndelay) ;
699 IVfp80(stdout, newchv->nD + newchv->nU, newchv->colind, 80, &ierr) ;
700 fflush(stdout) ;
701 #endif
702 if ( CHV_IS_NONSYMMETRIC(newchv) ) {
703 ndelay = 0 ;
704 for ( child = firstchild ; child != NULL ; child = child->next ) {
705 IVcopy(child->nD, newchv->rowind + ndelay, child->rowind) ;
706 ndelay += child->nD ;
707 }
708 IVcopy(oldchv->nD + oldchv->nU,
709 newchv->rowind + ndelay, oldchv->rowind) ;
710 #if MYDEBUG > 0
711 fprintf(stdout, "\n ndelay = %d, new row indices", ndelay) ;
712 IVfp80(stdout, newchv->nD + newchv->nL, newchv->rowind, 80, &ierr) ;
713 fflush(stdout) ;
714 #endif
715 }
716 /*
717 ------------------------
718 assemble the old chevron
719 ------------------------
720 */
721 Chv_assembleChv(newchv, oldchv) ;
722 #if MYDEBUG > 0
723 fprintf(stdout, "\n after merging old chevron") ;
724 Chv_writeForHumanEye(newchv, stdout) ;
725 #endif
726 /*
727 ------------------------------
728 assemble the children chevrons
729 ------------------------------
730 */
731 for ( child = firstchild ; child != NULL ; child = child->next ) {
732 Chv_assembleChv(newchv, child) ;
733 #if MYDEBUG > 0
734 fprintf(stdout, "\n after merging child %d", child->id) ;
735 Chv_writeForHumanEye(newchv, stdout) ;
736 #endif
737 }
738
739 return(ndelay) ; }
740
741 /*--------------------------------------------------------------------*/
742