1 /* update.c */
2
3 #include "../Chv.h"
4
5 #define MYDEBUG 0
6
7 /*--------------------------------------------------------------------*/
8 /*
9 ---------------------------------------------------------------------
10 purpose -- perform the hermitian factor update
11 T_{\bnd{I} \cap J, \bnd{I} \cap J}
12 -= U_{I, \bnd{I} \cap J}^H D_{I, I} U_{I, \bnd{I} \cap J}
13 and
14 T_{\bnd{I} \cap J, \bnd{I} \cap \bnd{J}}
15 -= U_{I, \bnd{I} \cap J}^H D_{I, I} U_{I, \bnd{I} \cap \bnd{J}}
16
17 created -- 98apr17, cca
18 ---------------------------------------------------------------------
19 */
20 void
Chv_updateH(Chv * chvT,SubMtx * mtxD,SubMtx * mtxU,DV * tempDV)21 Chv_updateH (
22 Chv *chvT,
23 SubMtx *mtxD,
24 SubMtx *mtxU,
25 DV *tempDV
26 ) {
27 int firstInT, firstInU, jcolT, jcolU, lastInT, lastInU, ncolT, ncolU ;
28 int *colindT, *colindU ;
29 /*
30 ---------------
31 check the input
32 ---------------
33 */
34 if ( chvT == NULL || mtxD == NULL || mtxU == NULL || tempDV == NULL ) {
35 fprintf(stderr, "\n fatal error in Chv_updateH(%p,%p,%p,%p)"
36 "\n bad input\n", chvT, mtxD, mtxU, tempDV) ;
37 exit(-1) ;
38 }
39 if ( ! CHV_IS_COMPLEX(chvT) ) {
40 fprintf(stderr, "\n fatal error in Chv_updateH(%p,%p,%p,%p)"
41 "\n bad input, chvT is not complex\n",
42 chvT, mtxD, mtxU, tempDV) ;
43 exit(-1) ;
44 }
45 if ( ! SUBMTX_IS_COMPLEX(mtxD) ) {
46 fprintf(stderr, "\n fatal error in Chv_updateH(%p,%p,%p,%p)"
47 "\n bad input, mtxD is not complex\n",
48 chvT, mtxD, mtxU, tempDV) ;
49 exit(-1) ;
50 }
51 if ( ! SUBMTX_IS_COMPLEX(mtxU) ) {
52 fprintf(stderr, "\n fatal error in Chv_updateH(%p,%p,%p,%p)"
53 "\n bad input, mtxU is not complex\n",
54 chvT, mtxD, mtxU, tempDV) ;
55 exit(-1) ;
56 }
57 Chv_columnIndices(chvT, &ncolT, &colindT) ;
58 SubMtx_columnIndices(mtxU, &ncolU, &colindU) ;
59 /*
60 -----------------------------
61 locate first column of U in T
62 -----------------------------
63 */
64 firstInT = colindT[0] ;
65 lastInT = colindT[chvT->nD-1] ;
66 for ( jcolU = 0 ; jcolU < ncolU ; jcolU++ ) {
67 if ( firstInT <= colindU[jcolU] && colindU[jcolU] <= lastInT ) {
68 break ;
69 }
70 }
71 if ( (firstInU = jcolU) == ncolU ) {
72 return ;
73 }
74 /*
75 ----------------------------
76 locate last column of U in T
77 ----------------------------
78 */
79 lastInU = firstInU ;
80 for ( ; jcolU < ncolU ; jcolU++ ) {
81 if ( colindU[jcolU] <= lastInT ) {
82 lastInU = jcolU ;
83 } else {
84 break ;
85 }
86 }
87 #if MYDEBUG > 0
88 fprintf(stdout, "\n %% firstInU = %d, lastInU = %d",
89 firstInU, lastInU) ;
90 fflush(stdout) ;
91 #endif
92 /*
93 ----------------------------------------------------------
94 overwrite supported column indices with indices local to T
95 ----------------------------------------------------------
96 */
97 for ( jcolU = firstInU, jcolT = 0 ; jcolU < ncolU ; jcolU++ ) {
98 while ( colindU[jcolU] != colindT[jcolT] ) {
99 jcolT++ ;
100 }
101 colindU[jcolU] = jcolT ;
102 }
103 if ( SUBMTX_IS_DENSE_COLUMNS(mtxU) ) {
104 double isum, rsum ;
105 double sums[18] ;
106 double *base0, *base1, *base2, *colU0, *colU1, *colU2, *entU,
107 *rowUT0, *rowUT1, *rowUT2, *temp0, *temp1, *temp2 ;
108 int ichv0, ichv1, ichv2, ii, inc1, inc2, irowUT,
109 kloc0, kloc1, kloc2, nrowU ;
110
111 SubMtx_denseInfo(mtxU, &nrowU, &ncolU, &inc1, &inc2, &entU) ;
112 DV_setSize(tempDV, 6*nrowU) ;
113 temp0 = DV_entries(tempDV) ;
114 temp1 = temp0 + 2*nrowU ;
115 temp2 = temp1 + 2*nrowU ;
116 /*
117 --------------------------------------------
118 loop over the rows of U^H in groups of three
119 --------------------------------------------
120 */
121 rowUT0 = entU + 2*firstInU*nrowU ;
122 for ( irowUT = firstInU ; irowUT <= lastInU - 2 ; irowUT += 3 ) {
123 #if MYDEBUG > 0
124 fprintf(stdout, "\n %% working on rows %d, %d and %d",
125 colindU[irowUT], colindU[irowUT+1], colindU[irowUT+2]) ;
126 fflush(stdout) ;
127 #endif
128 rowUT1 = rowUT0 + 2*nrowU ;
129 rowUT2 = rowUT1 + 2*nrowU ;
130 /*
131 -----------------------------------------------------
132 get base locations for the chevron's diagonal entries
133 -----------------------------------------------------
134 */
135 ichv0 = colindU[irowUT] ;
136 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
137 ichv1 = colindU[irowUT+1] ;
138 base1 = Chv_diagLocation(chvT, ichv1) - 2*ichv1 ;
139 ichv2 = colindU[irowUT+2] ;
140 base2 = Chv_diagLocation(chvT, ichv2) - 2*ichv2 ;
141 /*
142 ------------------------------------
143 [ temp0 ] [ rowUT0 ]^H
144 compute [ temp1 ] = [ rowUT1 ] * D
145 [ temp2 ] [ rowUT2 ]
146 ------------------------------------
147 */
148 DVzero(6*nrowU, temp0) ;
149 SubMtx_scale3vec(mtxD, temp0, temp1, temp2,
150 rowUT0, rowUT1, rowUT2) ;
151 for ( ii = 0 ; ii < nrowU ; ii++ ) {
152 temp0[2*ii+1] = -temp0[2*ii+1] ;
153 temp1[2*ii+1] = -temp1[2*ii+1] ;
154 temp2[2*ii+1] = -temp2[2*ii+1] ;
155 }
156 /*
157 --------------------------------------------------
158 update the 3x3 upper triangle for these three rows
159 --------------------------------------------------
160 */
161 colU0 = rowUT0 ;
162 colU1 = rowUT1 ;
163 colU2 = rowUT2 ;
164 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
165 /*
166 if ( chvT->id != -1 ) {
167 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
168 chvT->id, mtxD->rowid, base0[2*ichv0+1]) ;
169 }
170 */
171 /*
172 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
173 */
174 base0[2*ichv0] -= rsum ;
175 /*
176 if ( chvT->id != -1 ) {
177 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
178 isum, base0[2*ichv0+1]) ;
179 }
180 */
181 ZVdotU(nrowU, temp0, colU1, &rsum, &isum) ;
182 base0[2*ichv1] -= rsum ; base0[2*ichv1+1] -= isum ;
183 ZVdotU(nrowU, temp0, colU2, &rsum, &isum) ;
184 base0[2*ichv2] -= rsum ; base0[2*ichv2+1] -= isum ;
185 ZVdotU(nrowU, temp1, colU1, &rsum, &isum) ;
186 /*
187 if ( chvT->id != -1 ) {
188 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
189 chvT->id, mtxD->rowid, base1[2*ichv1+1]) ;
190 }
191 */
192 /*
193 base1[2*ichv1] -= rsum ; base1[2*ichv1+1] -= isum ;
194 */
195 base1[2*ichv1] -= rsum ;
196 /*
197 if ( chvT->id != -1 ) {
198 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
199 isum, base1[2*ichv1+1]) ;
200 }
201 */
202 ZVdotU(nrowU, temp1, colU2, &rsum, &isum) ;
203 base1[2*ichv2] -= rsum ; base1[2*ichv2+1] -= isum ;
204 ZVdotU(nrowU, temp2, colU2, &rsum, &isum) ;
205 /*
206 if ( chvT->id != -1 ) {
207 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
208 chvT->id, mtxD->rowid, base2[2*ichv2+1]) ;
209 }
210 */
211 /*
212 base2[2*ichv2] -= rsum ; base2[2*ichv2+1] -= isum ;
213 */
214 base2[2*ichv2] -= rsum ;
215 /*
216 if ( chvT->id != -1 ) {
217 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
218 isum, base2[2*ichv2+1]) ;
219 }
220 */
221 colU0 = colU2 + 2*nrowU ;
222 /*
223 --------------------------------------
224 update the remainder of the three rows
225 --------------------------------------
226 */
227 for ( jcolU = irowUT + 3 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
228 colU1 = colU0 + 2*nrowU ;
229 colU2 = colU1 + 2*nrowU ;
230 ZVdotU33(nrowU, temp0, temp1, temp2,
231 colU0, colU1, colU2, sums) ;
232 kloc0 = 2*colindU[jcolU] ;
233 kloc1 = 2*colindU[jcolU+1] ;
234 kloc2 = 2*colindU[jcolU+2] ;
235 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
236 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
237 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
238 base1[kloc0] -= sums[ 6] ; base1[kloc0+1] -= sums[ 7] ;
239 base1[kloc1] -= sums[ 8] ; base1[kloc1+1] -= sums[ 9] ;
240 base1[kloc2] -= sums[10] ; base1[kloc2+1] -= sums[11] ;
241 base2[kloc0] -= sums[12] ; base2[kloc0+1] -= sums[13] ;
242 base2[kloc1] -= sums[14] ; base2[kloc1+1] -= sums[15] ;
243 base2[kloc2] -= sums[16] ; base2[kloc2+1] -= sums[17] ;
244 colU0 = colU2 + 2*nrowU ;
245 }
246 if ( jcolU == ncolU - 2 ) {
247 colU1 = colU0 + 2*nrowU ;
248 ZVdotU32(nrowU, temp0, temp1, temp2, colU0, colU1, sums) ;
249 kloc0 = 2*colindU[jcolU] ;
250 kloc1 = 2*colindU[jcolU+1] ;
251 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
252 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
253 base1[kloc0] -= sums[ 4] ; base1[kloc0+1] -= sums[ 5] ;
254 base1[kloc1] -= sums[ 6] ; base1[kloc1+1] -= sums[ 7] ;
255 base2[kloc0] -= sums[ 8] ; base2[kloc0+1] -= sums[ 9] ;
256 base2[kloc1] -= sums[10] ; base2[kloc1+1] -= sums[11] ;
257 } else if ( jcolU == ncolU - 1 ) {
258 ZVdotU31(nrowU, temp0, temp1, temp2, colU0, sums) ;
259 kloc0 = 2*colindU[jcolU] ;
260 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
261 base1[kloc0] -= sums[ 2] ; base1[kloc0+1] -= sums[ 3] ;
262 base2[kloc0] -= sums[ 4] ; base2[kloc0+1] -= sums[ 5] ;
263 }
264 rowUT0 = rowUT2 + 2*nrowU ;
265 }
266 if ( irowUT == lastInU - 1 ) {
267 #if MYDEBUG > 0
268 fprintf(stdout, "\n %% working on rows %d and %d",
269 colindU[irowUT], colindU[irowUT+1]) ;
270 fflush(stdout) ;
271 #endif
272 rowUT1 = rowUT0 + 2*nrowU ;
273 /*
274 -----------------------------------------------------
275 get base locations for the chevron's diagonal entries
276 -----------------------------------------------------
277 */
278 ichv0 = colindU[irowUT] ;
279 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
280 ichv1 = colindU[irowUT+1] ;
281 base1 = Chv_diagLocation(chvT, ichv1) - 2*ichv1 ;
282 /*
283 ------------------------------------
284 [ temp0 ] [ rowUT0 ]^H
285 compute [ temp1 ] = [ rowUT1 ] * D
286 ------------------------------------
287 */
288 DVzero(4*nrowU, temp0) ;
289 SubMtx_scale2vec(mtxD, temp0, temp1, rowUT0, rowUT1) ;
290 for ( ii = 0 ; ii < nrowU ; ii++ ) {
291 temp0[2*ii+1] = -temp0[2*ii+1] ;
292 temp1[2*ii+1] = -temp1[2*ii+1] ;
293 }
294 /*
295 ------------------------------------------------
296 update the 2x2 upper triangle for these two rows
297 ------------------------------------------------
298 */
299 colU0 = rowUT0 ;
300 colU1 = rowUT1 ;
301 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
302 /*
303 if ( chvT->id != -1 ) {
304 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
305 chvT->id, mtxD->rowid, base0[2*ichv0+1]) ;
306 }
307 */
308 /*
309 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
310 */
311 base0[2*ichv0] -= rsum ;
312 /*
313 if ( chvT->id != -1 ) {
314 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
315 isum, base0[2*ichv0+1]) ;
316 }
317 */
318 ZVdotU(nrowU, temp0, colU1, &rsum, &isum) ;
319 base0[2*ichv1] -= rsum ; base0[2*ichv1+1] -= isum ;
320 ZVdotU(nrowU, temp1, colU1, &rsum, &isum) ;
321 /*
322 if ( chvT->id != -1 ) {
323 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
324 chvT->id, mtxD->rowid, base1[2*ichv1+1]) ;
325 }
326 */
327 /*
328 base1[2*ichv1] -= rsum ; base1[2*ichv1+1] -= isum ;
329 */
330 base1[2*ichv1] -= rsum ;
331 /*
332 if ( chvT->id != -1 ) {
333 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
334 isum, base1[2*ichv1+1]) ;
335 }
336 */
337 colU0 = colU1 + 2*nrowU ;
338 /*
339 ------------------------------------
340 update the remainder of the two rows
341 ------------------------------------
342 */
343 for ( jcolU = irowUT + 2 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
344 colU1 = colU0 + 2*nrowU ;
345 colU2 = colU1 + 2*nrowU ;
346 ZVdotU23(nrowU, temp0, temp1, colU0, colU1, colU2, sums) ;
347 kloc0 = 2*colindU[jcolU] ;
348 kloc1 = 2*colindU[jcolU+1] ;
349 kloc2 = 2*colindU[jcolU+2] ;
350 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
351 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
352 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
353 base1[kloc0] -= sums[ 6] ; base1[kloc0+1] -= sums[ 7] ;
354 base1[kloc1] -= sums[ 8] ; base1[kloc1+1] -= sums[ 9] ;
355 base1[kloc2] -= sums[10] ; base1[kloc2+1] -= sums[11] ;
356 colU0 = colU2 + 2*nrowU ;
357 }
358 if ( jcolU == ncolU - 2 ) {
359 colU1 = colU0 + 2*nrowU ;
360 ZVdotU22(nrowU, temp0, temp1, colU0, colU1, sums) ;
361 kloc0 = 2*colindU[jcolU] ;
362 kloc1 = 2*colindU[jcolU+1] ;
363 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
364 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
365 base1[kloc0] -= sums[ 4] ; base1[kloc0+1] -= sums[ 5] ;
366 base1[kloc1] -= sums[ 6] ; base1[kloc1+1] -= sums[ 7] ;
367 } else if ( jcolU == ncolU - 1 ) {
368 ZVdotU21(nrowU, temp0, temp1, colU0, sums) ;
369 kloc0 = 2*colindU[jcolU] ;
370 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
371 base1[kloc0] -= sums[ 2] ; base1[kloc0+1] -= sums[ 3] ;
372 }
373 } else if ( irowUT == lastInU ) {
374 #if MYDEBUG > 0
375 fprintf(stdout, "\n %% working on row %d", colindU[irowUT]) ;
376 fflush(stdout) ;
377 #endif
378 /*
379 -----------------------------------------------------
380 get base locations for the chevron's diagonal entries
381 -----------------------------------------------------
382 */
383 ichv0 = colindU[irowUT] ;
384 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
385 /*
386 ------------------------------------
387 compute [ temp0 ] = [ rowUT0 ]^H * D
388 ------------------------------------
389 */
390 DVzero(2*nrowU, temp0) ;
391 SubMtx_scale1vec(mtxD, temp0, rowUT0) ;
392 for ( ii = 0 ; ii < nrowU ; ii++ ) {
393 temp0[2*ii+1] = -temp0[2*ii+1] ;
394 }
395 /*
396 ------------------------------------------
397 update the 1x1 upper triangle for this row
398 ------------------------------------------
399 */
400 colU0 = rowUT0 ;
401 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
402 /*
403 if ( chvT->id != -1 ) {
404 fprintf(stdout, "\n (%d,%d) : imag(diag) = %12.5e",
405 chvT->id, mtxD->rowid, base0[2*ichv0+1]) ;
406 }
407 */
408 /*
409 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
410 */
411 base0[2*ichv0] -= rsum ;
412 /*
413 if ( chvT->id != -1 ) {
414 fprintf(stdout, ", isum = %12.5e, imag(diag) = %12.5e",
415 isum, base0[2*ichv0+1]) ;
416 }
417 */
418 colU0 = colU0 + 2*nrowU ;
419 /*
420 -------------------------------
421 update the remainder of the row
422 -------------------------------
423 */
424 for ( jcolU = irowUT + 1 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
425 #if MYDEBUG > 0
426 fprintf(stdout, "\n %% working on columns %d, %d and %d",
427 jcolU, jcolU+1, jcolU+2) ;
428 fflush(stdout) ;
429 #endif
430 colU1 = colU0 + 2*nrowU ;
431 colU2 = colU1 + 2*nrowU ;
432 ZVdotU13(nrowU, temp0, colU0, colU1, colU2, sums) ;
433 kloc0 = 2*colindU[jcolU] ;
434 kloc1 = 2*colindU[jcolU+1] ;
435 kloc2 = 2*colindU[jcolU+2] ;
436 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
437 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
438 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
439 colU0 = colU2 + 2*nrowU ;
440 }
441 if ( jcolU == ncolU - 2 ) {
442 #if MYDEBUG > 0
443 fprintf(stdout, "\n %% working on columns %d and %d",
444 jcolU, jcolU+1) ;
445 fflush(stdout) ;
446 #endif
447 colU1 = colU0 + 2*nrowU ;
448 ZVdotU12(nrowU, temp0, colU0, colU1, sums) ;
449 kloc0 = 2*colindU[jcolU] ;
450 kloc1 = 2*colindU[jcolU+1] ;
451 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
452 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
453 } else if ( jcolU == ncolU - 1 ) {
454 #if MYDEBUG > 0
455 fprintf(stdout, "\n %% working on column %d", jcolU) ;
456 fflush(stdout) ;
457 #endif
458 ZVdotU11(nrowU, temp0, colU0, sums) ;
459 /*
460 fprintf(stdout, "\n SUMS[0] = %12.4e, SUMS[1] = %12.4e",
461 sums[0], sums[1]) ;
462 */
463 kloc0 = 2*colindU[jcolU] ;
464 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
465 /*
466 fprintf(stdout, "\n base0[%d] = %12.4e, base0[%d] = %12.4e",
467 kloc0, base0[kloc0], kloc0+1, base0[kloc0+1]) ;
468 */
469 }
470 }
471 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxU) ) {
472 double isum, rsum ;
473 double *base0, *colU0, *entU, *rowUT0, *temp0, *temp1 ;
474 int ichv0, ii, iloc, irowUT, kloc0, nentU, nrowU, offset,
475 rloc, sizeU, sizeUT ;
476 int *indU, *indU0, *indUT0, *sizes ;
477
478 SubMtx_sparseColumnsInfo(mtxU, &ncolU, &nentU, &sizes, &indU, &entU) ;
479 nrowU = mtxU->nrow ;
480 DV_setSize(tempDV, 4*nrowU) ;
481 temp0 = DV_entries(tempDV) ;
482 temp1 = temp0 + 2*nrowU ;
483 /*
484 -------------------------------------------
485 get the offset into the indices and entries
486 -------------------------------------------
487 */
488 for ( jcolU = offset = 0 ; jcolU < firstInU ; jcolU++ ) {
489 offset += sizes[jcolU] ;
490 }
491 /*
492 ------------------------------------
493 loop over the supporting rows of U^H
494 ------------------------------------
495 */
496 rowUT0 = entU + 2*offset ;
497 indUT0 = indU + offset ;
498 for ( irowUT = firstInU ; irowUT <= lastInU ; irowUT++ ) {
499 if ( (sizeUT = sizes[irowUT]) > 0 ) {
500 /*
501 -----------------------------------------------------
502 get base locations for the chevron's diagonal entries
503 -----------------------------------------------------
504 */
505 ichv0 = colindU[irowUT] ;
506 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
507 /*
508 ------------------------------------
509 compute [ temp0 ] = [ rowUT0 ]^H * D
510 ------------------------------------
511 */
512 DVzero(4*nrowU, temp0) ;
513 for ( ii = 0 ; ii < sizeUT ; ii++ ) {
514 rloc = 2*indUT0[ii] ; iloc = rloc + 1 ;
515 temp1[rloc] = rowUT0[2*ii] ;
516 temp1[iloc] = rowUT0[2*ii+1] ;
517 }
518 SubMtx_scale1vec(mtxD, temp0, temp1) ;
519 for ( ii = 0 ; ii < nrowU ; ii++ ) {
520 temp0[2*ii+1] = -temp0[2*ii+1] ;
521 }
522 /*
523 -------------------------------
524 loop over the following columns
525 -------------------------------
526 */
527 colU0 = rowUT0 ;
528 indU0 = indUT0 ;
529 for ( jcolU = irowUT ; jcolU < ncolU ; jcolU++ ) {
530 if ( (sizeU = sizes[jcolU]) > 0 ) {
531 ZVdotiU(sizeU, temp0, indU0, colU0, &rsum, &isum) ;
532 kloc0 = 2*colindU[jcolU] ;
533 base0[kloc0] -= rsum ; base0[kloc0+1] -= isum ;
534 colU0 += 2*sizeU ;
535 indU0 += sizeU ;
536 }
537 }
538 rowUT0 += 2*sizeUT ;
539 indUT0 += sizeUT ;
540 }
541 }
542 } else {
543 fprintf(stderr, "\n fatal error in Chv_updateH(%p,%p,%p,%p)"
544 "\n mtxU must have dense or sparse columns\n",
545 chvT, mtxD, mtxU, tempDV) ;
546 exit(-1) ;
547 }
548 /*
549 ---------------------------------------------------
550 overwrite the local indices with the global indices
551 ---------------------------------------------------
552 */
553 for ( jcolU = firstInU ; jcolU < ncolU ; jcolU++ ) {
554 colindU[jcolU] = colindT[colindU[jcolU]] ;
555 }
556 return ; }
557
558 /*--------------------------------------------------------------------*/
559 /*
560 ---------------------------------------------------------------------
561 purpose -- perform the symmetric factor update
562 T_{\bnd{I} \cap J, \bnd{I} \cap J}
563 -= U_{I, \bnd{I} \cap J}^T D_{I, I} U_{I, \bnd{I} \cap J}
564 and
565 T_{\bnd{I} \cap J, \bnd{I} \cap \bnd{J}}
566 -= U_{I, \bnd{I} \cap J}^T D_{I, I} U_{I, \bnd{I} \cap \bnd{J}}
567
568 created -- 98apr17, cca
569 ---------------------------------------------------------------------
570 */
571 void
Chv_updateS(Chv * chvT,SubMtx * mtxD,SubMtx * mtxU,DV * tempDV)572 Chv_updateS (
573 Chv *chvT,
574 SubMtx *mtxD,
575 SubMtx *mtxU,
576 DV *tempDV
577 ) {
578 int firstInT, firstInU, jcolT, jcolU, lastInT, lastInU, ncolT, ncolU ;
579 int *colindT, *colindU ;
580 /*
581 ---------------
582 check the input
583 ---------------
584 */
585 if ( chvT == NULL || mtxD == NULL || mtxU == NULL || tempDV == NULL ) {
586 fprintf(stderr, "\n fatal error in Chv_updateS(%p,%p,%p,%p)"
587 "\n bad input\n", chvT, mtxD, mtxU, tempDV) ;
588 exit(-1) ;
589 }
590 if ( CHV_IS_REAL(chvT) ) {
591 if ( ! SUBMTX_IS_REAL(mtxD) || ! SUBMTX_IS_REAL(mtxU) ) {
592 fprintf(stderr, "\n fatal error in Chv_updateT(%p,%p,%p,%p)"
593 "\n chvT is real, but mtxD and/or mtxU are not\n",
594 chvT, mtxD, mtxU, tempDV) ;
595 exit(-1) ;
596 }
597 } else if ( CHV_IS_COMPLEX(chvT) ) {
598 if ( ! SUBMTX_IS_COMPLEX(mtxD) || ! SUBMTX_IS_COMPLEX(mtxU) ) {
599 fprintf(stderr, "\n fatal error in Chv_updateT(%p,%p,%p,%p)"
600 "\n chvT is complex, but mtxD and/or mtxU are not\n",
601 chvT, mtxD, mtxU, tempDV) ;
602 exit(-1) ;
603 }
604 } else {
605 fprintf(stderr, "\n fatal error in Chv_updateT(%p,%p,%p,%p)"
606 "\n bad input, chvT is not real or complex\n",
607 chvT, mtxD, mtxU, tempDV) ;
608 exit(-1) ;
609 }
610 Chv_columnIndices(chvT, &ncolT, &colindT) ;
611 SubMtx_columnIndices(mtxU, &ncolU, &colindU) ;
612 #if MYDEBUG > 0
613 fprintf(stdout, "\n %% Chv column indices") ;
614 IVfprintf(stdout, ncolT, colindT) ;
615 fprintf(stdout, "\n %% submtx column indices") ;
616 IVfprintf(stdout, ncolU, colindU) ;
617 fflush(stdout) ;
618 #endif
619 /*
620 -----------------------------
621 locate first column of U in T
622 -----------------------------
623 */
624 firstInT = colindT[0] ;
625 lastInT = colindT[chvT->nD-1] ;
626 for ( jcolU = 0 ; jcolU < ncolU ; jcolU++ ) {
627 if ( firstInT <= colindU[jcolU] && colindU[jcolU] <= lastInT ) {
628 break ;
629 }
630 }
631 if ( (firstInU = jcolU) == ncolU ) {
632 return ;
633 }
634 /*
635 ----------------------------
636 locate last column of U in T
637 ----------------------------
638 */
639 lastInU = firstInU ;
640 for ( ; jcolU < ncolU ; jcolU++ ) {
641 if ( colindU[jcolU] <= lastInT ) {
642 lastInU = jcolU ;
643 } else {
644 break ;
645 }
646 }
647 #if MYDEBUG > 0
648 fprintf(stdout, "\n %% firstInU = %d, lastInU = %d",
649 firstInU, lastInU) ;
650 fflush(stdout) ;
651 #endif
652 /*
653 ----------------------------------------------------------
654 overwrite supported column indices with indices local to T
655 ----------------------------------------------------------
656 */
657 for ( jcolU = firstInU, jcolT = 0 ; jcolU < ncolU ; jcolU++ ) {
658 while ( colindU[jcolU] != colindT[jcolT] ) {
659 jcolT++ ;
660 }
661 colindU[jcolU] = jcolT ;
662 }
663 #if MYDEBUG > 0
664 fprintf(stdout, "\n %% submtx column indices") ;
665 IVfprintf(stdout, ncolU, colindU) ;
666 fflush(stdout) ;
667 #endif
668 if ( CHV_IS_REAL(chvT) ) {
669 if ( SUBMTX_IS_DENSE_COLUMNS(mtxU) ) {
670 double sums[9] ;
671 double *base0, *base1, *base2, *colU0, *colU1, *colU2, *entU,
672 *rowUT0, *rowUT1, *rowUT2, *temp0, *temp1, *temp2 ;
673 int ichv0, ichv1, ichv2, inc1, inc2, irowUT,
674 kloc0, kloc1, kloc2, nrowU ;
675
676 SubMtx_denseInfo(mtxU, &nrowU, &ncolU, &inc1, &inc2, &entU) ;
677 DV_setSize(tempDV, 3*nrowU) ;
678 temp0 = DV_entries(tempDV) ;
679 temp1 = temp0 + nrowU ;
680 temp2 = temp1 + nrowU ;
681 /*
682 --------------------------------------------
683 loop over the rows of U^T in groups of three
684 --------------------------------------------
685 */
686 rowUT0 = entU + firstInU*nrowU ;
687 for ( irowUT = firstInU ; irowUT <= lastInU - 2 ; irowUT += 3 ) {
688 #if MYDEBUG > 0
689 fprintf(stdout, "\n %% working on rows %d, %d and %d",
690 colindU[irowUT], colindU[irowUT+1], colindU[irowUT+2]) ;
691 fflush(stdout) ;
692 #endif
693 rowUT1 = rowUT0 + nrowU ;
694 rowUT2 = rowUT1 + nrowU ;
695 /*
696 -----------------------------------------------------
697 get base locations for the chevron's diagonal entries
698 -----------------------------------------------------
699 */
700 ichv0 = colindU[irowUT] ;
701 base0 = Chv_diagLocation(chvT, ichv0) - ichv0 ;
702 ichv1 = colindU[irowUT+1] ;
703 base1 = Chv_diagLocation(chvT, ichv1) - ichv1 ;
704 ichv2 = colindU[irowUT+2] ;
705 base2 = Chv_diagLocation(chvT, ichv2) - ichv2 ;
706 /*
707 ------------------------------------
708 [ temp0 ] [ rowUT0 ]^T
709 compute [ temp1 ] = [ rowUT1 ] * D
710 [ temp2 ] [ rowUT2 ]
711 ------------------------------------
712 */
713 DVzero(3*nrowU, temp0) ;
714 SubMtx_scale3vec(mtxD, temp0, temp1, temp2,
715 rowUT0, rowUT1, rowUT2) ;
716 /*
717 --------------------------------------------------
718 update the 3x3 upper triangle for these three rows
719 --------------------------------------------------
720 */
721 colU0 = rowUT0 ;
722 colU1 = rowUT1 ;
723 colU2 = rowUT2 ;
724 base0[ichv0] -= DVdot(nrowU, temp0, colU0) ;
725 base0[ichv1] -= DVdot(nrowU, temp0, colU1) ;
726 base0[ichv2] -= DVdot(nrowU, temp0, colU2) ;
727 base1[ichv1] -= DVdot(nrowU, temp1, colU1) ;
728 base1[ichv2] -= DVdot(nrowU, temp1, colU2) ;
729 base2[ichv2] -= DVdot(nrowU, temp2, colU2) ;
730 colU0 = colU2 + nrowU ;
731 /*
732 --------------------------------------
733 update the remainder of the three rows
734 --------------------------------------
735 */
736 for ( jcolU = irowUT + 3 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
737 colU1 = colU0 + nrowU ;
738 colU2 = colU1 + nrowU ;
739 DVdot33(nrowU, temp0, temp1, temp2,
740 colU0, colU1, colU2, sums) ;
741 kloc0 = colindU[jcolU] ;
742 kloc1 = colindU[jcolU+1] ;
743 kloc2 = colindU[jcolU+2] ;
744 base0[kloc0] -= sums[0] ;
745 base0[kloc1] -= sums[1] ;
746 base0[kloc2] -= sums[2] ;
747 base1[kloc0] -= sums[3] ;
748 base1[kloc1] -= sums[4] ;
749 base1[kloc2] -= sums[5] ;
750 base2[kloc0] -= sums[6] ;
751 base2[kloc1] -= sums[7] ;
752 base2[kloc2] -= sums[8] ;
753 colU0 = colU2 + nrowU ;
754 }
755 if ( jcolU == ncolU - 2 ) {
756 colU1 = colU0 + nrowU ;
757 DVdot32(nrowU, temp0, temp1, temp2, colU0, colU1, sums) ;
758 kloc0 = colindU[jcolU] ;
759 kloc1 = colindU[jcolU+1] ;
760 base0[kloc0] -= sums[0] ;
761 base0[kloc1] -= sums[1] ;
762 base1[kloc0] -= sums[2] ;
763 base1[kloc1] -= sums[3] ;
764 base2[kloc0] -= sums[4] ;
765 base2[kloc1] -= sums[5] ;
766 } else if ( jcolU == ncolU - 1 ) {
767 DVdot31(nrowU, temp0, temp1, temp2, colU0, sums) ;
768 kloc0 = colindU[jcolU] ;
769 base0[kloc0] -= sums[0] ;
770 base1[kloc0] -= sums[1] ;
771 base2[kloc0] -= sums[2] ;
772 }
773 rowUT0 = rowUT2 + nrowU ;
774 }
775 if ( irowUT == lastInU - 1 ) {
776 #if MYDEBUG > 0
777 fprintf(stdout, "\n %% working on rows %d and %d",
778 colindU[irowUT], colindU[irowUT+1]) ;
779 fflush(stdout) ;
780 #endif
781 rowUT1 = rowUT0 + nrowU ;
782 /*
783 -----------------------------------------------------
784 get base locations for the chevron's diagonal entries
785 -----------------------------------------------------
786 */
787 ichv0 = colindU[irowUT] ;
788 base0 = Chv_diagLocation(chvT, ichv0) - ichv0 ;
789 ichv1 = colindU[irowUT+1] ;
790 base1 = Chv_diagLocation(chvT, ichv1) - ichv1 ;
791 /*
792 ------------------------------------
793 [ temp0 ] [ rowUT0 ]^T
794 compute [ temp1 ] = [ rowUT1 ] * D
795 ------------------------------------
796 */
797 DVzero(2*nrowU, temp0) ;
798 SubMtx_scale2vec(mtxD, temp0, temp1, rowUT0, rowUT1) ;
799 /*
800 ------------------------------------------------
801 update the 2x2 upper triangle for these two rows
802 ------------------------------------------------
803 */
804 colU0 = rowUT0 ;
805 colU1 = rowUT1 ;
806 base0[ichv0] -= DVdot(nrowU, temp0, colU0) ;
807 base0[ichv1] -= DVdot(nrowU, temp0, colU1) ;
808 base1[ichv1] -= DVdot(nrowU, temp1, colU1) ;
809 colU0 = colU1 + nrowU ;
810 /*
811 ------------------------------------
812 update the remainder of the two rows
813 ------------------------------------
814 */
815 for ( jcolU = irowUT + 2 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
816 colU1 = colU0 + nrowU ;
817 colU2 = colU1 + nrowU ;
818 DVdot23(nrowU, temp0, temp1, colU0, colU1, colU2, sums) ;
819 kloc0 = colindU[jcolU] ;
820 kloc1 = colindU[jcolU+1] ;
821 kloc2 = colindU[jcolU+2] ;
822 base0[kloc0] -= sums[0] ;
823 base0[kloc1] -= sums[1] ;
824 base0[kloc2] -= sums[2] ;
825 base1[kloc0] -= sums[3] ;
826 base1[kloc1] -= sums[4] ;
827 base1[kloc2] -= sums[5] ;
828 colU0 = colU2 + nrowU ;
829 }
830 if ( jcolU == ncolU - 2 ) {
831 colU1 = colU0 + nrowU ;
832 DVdot22(nrowU, temp0, temp1, colU0, colU1, sums) ;
833 kloc0 = colindU[jcolU] ;
834 kloc1 = colindU[jcolU+1] ;
835 base0[kloc0] -= sums[0] ;
836 base0[kloc1] -= sums[1] ;
837 base1[kloc0] -= sums[2] ;
838 base1[kloc1] -= sums[3] ;
839 } else if ( jcolU == ncolU - 1 ) {
840 DVdot21(nrowU, temp0, temp1, colU0, sums) ;
841 kloc0 = colindU[jcolU] ;
842 base0[kloc0] -= sums[0] ;
843 base1[kloc0] -= sums[1] ;
844 }
845 } else if ( irowUT == lastInU ) {
846 #if MYDEBUG > 0
847 fprintf(stdout, "\n %% working on row %d", colindU[irowUT]) ;
848 fflush(stdout) ;
849 #endif
850 /*
851 -----------------------------------------------------
852 get base locations for the chevron's diagonal entries
853 -----------------------------------------------------
854 */
855 ichv0 = colindU[irowUT] ;
856 base0 = Chv_diagLocation(chvT, ichv0) - ichv0 ;
857 /*
858 ------------------------------------
859 compute [ temp0 ] = [ rowUT0 ]^T * D
860 ------------------------------------
861 */
862 DVzero(nrowU, temp0) ;
863 SubMtx_scale1vec(mtxD, temp0, rowUT0) ;
864 /*
865 ------------------------------------------
866 update the 1x1 upper triangle for this row
867 ------------------------------------------
868 */
869 colU0 = rowUT0 ;
870 base0[ichv0] -= DVdot(nrowU, temp0, colU0) ;
871 colU0 = colU0 + nrowU ;
872 /*
873 -------------------------------
874 update the remainder of the row
875 -------------------------------
876 */
877 for ( jcolU = irowUT + 1 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
878 #if MYDEBUG > 0
879 fprintf(stdout, "\n %% working on columns %d, %d and %d",
880 jcolU, jcolU+1, jcolU+2) ;
881 fflush(stdout) ;
882 #endif
883 colU1 = colU0 + nrowU ;
884 colU2 = colU1 + nrowU ;
885 DVdot13(nrowU, temp0, colU0, colU1, colU2, sums) ;
886 kloc0 = colindU[jcolU] ;
887 kloc1 = colindU[jcolU+1] ;
888 kloc2 = colindU[jcolU+2] ;
889 base0[kloc0] -= sums[0] ;
890 base0[kloc1] -= sums[1] ;
891 base0[kloc2] -= sums[2] ;
892 colU0 = colU2 + nrowU ;
893 }
894 if ( jcolU == ncolU - 2 ) {
895 #if MYDEBUG > 0
896 fprintf(stdout, "\n %% working on columns %d and %d",
897 jcolU, jcolU+1) ;
898 fflush(stdout) ;
899 #endif
900 colU1 = colU0 + nrowU ;
901 DVdot12(nrowU, temp0, colU0, colU1, sums) ;
902 kloc0 = colindU[jcolU] ;
903 kloc1 = colindU[jcolU+1] ;
904 base0[kloc0] -= sums[0] ;
905 base0[kloc1] -= sums[1] ;
906 } else if ( jcolU == ncolU - 1 ) {
907 #if MYDEBUG > 0
908 fprintf(stdout, "\n %% working on column %d", jcolU) ;
909 fflush(stdout) ;
910 #endif
911 DVdot11(nrowU, temp0, colU0, sums) ;
912 kloc0 = colindU[jcolU] ;
913 base0[kloc0] -= sums[0] ;
914 }
915 }
916 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxU) ) {
917 double sum ;
918 double *base0, *colU0, *entU, *rowUT0, *temp0, *temp1 ;
919 int ichv0, ii, irowUT, kloc0, loc, nentU, nrowU, offset,
920 sizeU, sizeUT ;
921 int *indU, *indU0, *indUT0, *sizes ;
922
923 SubMtx_sparseColumnsInfo(mtxU,
924 &ncolU, &nentU, &sizes, &indU, &entU) ;
925 nrowU = mtxU->nrow ;
926 DV_setSize(tempDV, 2*nrowU) ;
927 temp0 = DV_entries(tempDV) ;
928 temp1 = temp0 + nrowU ;
929 /*
930 -------------------------------------------
931 get the offset into the indices and entries
932 -------------------------------------------
933 */
934 for ( jcolU = offset = 0 ; jcolU < firstInU ; jcolU++ ) {
935 offset += sizes[jcolU] ;
936 }
937 /*
938 ------------------------------------
939 loop over the supporting rows of U^T
940 ------------------------------------
941 */
942 rowUT0 = entU + offset ;
943 indUT0 = indU + offset ;
944 for ( irowUT = firstInU ; irowUT <= lastInU ; irowUT++ ) {
945 if ( (sizeUT = sizes[irowUT]) > 0 ) {
946 /*
947 -----------------------------------------------------
948 get base locations for the chevron's diagonal entries
949 -----------------------------------------------------
950 */
951 ichv0 = colindU[irowUT] ;
952 base0 = Chv_diagLocation(chvT, ichv0) - ichv0 ;
953 /*
954 ------------------------------------
955 compute [ temp0 ] = [ rowUT0 ]^T * D
956 ------------------------------------
957 */
958 DVzero(2*nrowU, temp0) ;
959 for ( ii = 0 ; ii < sizeUT ; ii++ ) {
960 loc = indUT0[ii] ;
961 temp1[loc] = rowUT0[ii] ;
962 }
963 SubMtx_scale1vec(mtxD, temp0, temp1) ;
964 /*
965 -------------------------------
966 loop over the following columns
967 -------------------------------
968 */
969 colU0 = rowUT0 ;
970 indU0 = indUT0 ;
971 for ( jcolU = irowUT ; jcolU < ncolU ; jcolU++ ) {
972 if ( (sizeU = sizes[jcolU]) > 0 ) {
973 sum = DVdoti(sizeU, temp0, indU0, colU0) ;
974 kloc0 = colindU[jcolU] ;
975 base0[kloc0] -= sum ;
976 colU0 += sizeU ;
977 indU0 += sizeU ;
978 }
979 }
980 rowUT0 += sizeUT ;
981 indUT0 += sizeUT ;
982 }
983 }
984 } else {
985 fprintf(stderr, "\n fatal error in Chv_updateT(%p,%p,%p,%p)"
986 "\n bad input, mtxU must have dense or sparse columns\n",
987 chvT, mtxD, mtxU, tempDV) ;
988 exit(-1) ;
989 }
990 } else if ( CHV_IS_COMPLEX(chvT) ) {
991 if ( SUBMTX_IS_DENSE_COLUMNS(mtxU) ) {
992 double isum, rsum ;
993 double sums[18] ;
994 double *base0, *base1, *base2, *colU0, *colU1, *colU2, *entU,
995 *rowUT0, *rowUT1, *rowUT2, *temp0, *temp1, *temp2 ;
996 int ichv0, ichv1, ichv2, inc1, inc2, irowUT,
997 kloc0, kloc1, kloc2, nrowU ;
998
999 SubMtx_denseInfo(mtxU, &nrowU, &ncolU, &inc1, &inc2, &entU) ;
1000 DV_setSize(tempDV, 6*nrowU) ;
1001 temp0 = DV_entries(tempDV) ;
1002 temp1 = temp0 + 2*nrowU ;
1003 temp2 = temp1 + 2*nrowU ;
1004 /*
1005 --------------------------------------------
1006 loop over the rows of U^T in groups of three
1007 --------------------------------------------
1008 */
1009 rowUT0 = entU + 2*firstInU*nrowU ;
1010 for ( irowUT = firstInU ; irowUT <= lastInU - 2 ; irowUT += 3 ) {
1011 #if MYDEBUG > 0
1012 fprintf(stdout, "\n %% working on rows %d, %d and %d",
1013 colindU[irowUT], colindU[irowUT+1], colindU[irowUT+2]) ;
1014 fflush(stdout) ;
1015 #endif
1016 rowUT1 = rowUT0 + 2*nrowU ;
1017 rowUT2 = rowUT1 + 2*nrowU ;
1018 /*
1019 -----------------------------------------------------
1020 get base locations for the chevron's diagonal entries
1021 -----------------------------------------------------
1022 */
1023 ichv0 = colindU[irowUT] ;
1024 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
1025 ichv1 = colindU[irowUT+1] ;
1026 base1 = Chv_diagLocation(chvT, ichv1) - 2*ichv1 ;
1027 ichv2 = colindU[irowUT+2] ;
1028 base2 = Chv_diagLocation(chvT, ichv2) - 2*ichv2 ;
1029 /*
1030 ------------------------------------
1031 [ temp0 ] [ rowUT0 ]^T
1032 compute [ temp1 ] = [ rowUT1 ] * D
1033 [ temp2 ] [ rowUT2 ]
1034 ------------------------------------
1035 */
1036 DVzero(6*nrowU, temp0) ;
1037 SubMtx_scale3vec(mtxD, temp0, temp1, temp2,
1038 rowUT0, rowUT1, rowUT2) ;
1039 /*
1040 --------------------------------------------------
1041 update the 3x3 upper triangle for these three rows
1042 --------------------------------------------------
1043 */
1044 colU0 = rowUT0 ;
1045 colU1 = rowUT1 ;
1046 colU2 = rowUT2 ;
1047 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
1048 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
1049 ZVdotU(nrowU, temp0, colU1, &rsum, &isum) ;
1050 base0[2*ichv1] -= rsum ; base0[2*ichv1+1] -= isum ;
1051 ZVdotU(nrowU, temp0, colU2, &rsum, &isum) ;
1052 base0[2*ichv2] -= rsum ; base0[2*ichv2+1] -= isum ;
1053 ZVdotU(nrowU, temp1, colU1, &rsum, &isum) ;
1054 base1[2*ichv1] -= rsum ; base1[2*ichv1+1] -= isum ;
1055 ZVdotU(nrowU, temp1, colU2, &rsum, &isum) ;
1056 base1[2*ichv2] -= rsum ; base1[2*ichv2+1] -= isum ;
1057 ZVdotU(nrowU, temp2, colU2, &rsum, &isum) ;
1058 base2[2*ichv2] -= rsum ; base2[2*ichv2+1] -= isum ;
1059 colU0 = colU2 + 2*nrowU ;
1060 /*
1061 --------------------------------------
1062 update the remainder of the three rows
1063 --------------------------------------
1064 */
1065 for ( jcolU = irowUT + 3 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1066 colU1 = colU0 + 2*nrowU ;
1067 colU2 = colU1 + 2*nrowU ;
1068 ZVdotU33(nrowU, temp0, temp1, temp2,
1069 colU0, colU1, colU2, sums) ;
1070 kloc0 = 2*colindU[jcolU] ;
1071 kloc1 = 2*colindU[jcolU+1] ;
1072 kloc2 = 2*colindU[jcolU+2] ;
1073 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1074 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1075 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
1076 base1[kloc0] -= sums[ 6] ; base1[kloc0+1] -= sums[ 7] ;
1077 base1[kloc1] -= sums[ 8] ; base1[kloc1+1] -= sums[ 9] ;
1078 base1[kloc2] -= sums[10] ; base1[kloc2+1] -= sums[11] ;
1079 base2[kloc0] -= sums[12] ; base2[kloc0+1] -= sums[13] ;
1080 base2[kloc1] -= sums[14] ; base2[kloc1+1] -= sums[15] ;
1081 base2[kloc2] -= sums[16] ; base2[kloc2+1] -= sums[17] ;
1082 colU0 = colU2 + 2*nrowU ;
1083 }
1084 if ( jcolU == ncolU - 2 ) {
1085 colU1 = colU0 + 2*nrowU ;
1086 ZVdotU32(nrowU, temp0, temp1, temp2, colU0, colU1, sums) ;
1087 kloc0 = 2*colindU[jcolU] ;
1088 kloc1 = 2*colindU[jcolU+1] ;
1089 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1090 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1091 base1[kloc0] -= sums[ 4] ; base1[kloc0+1] -= sums[ 5] ;
1092 base1[kloc1] -= sums[ 6] ; base1[kloc1+1] -= sums[ 7] ;
1093 base2[kloc0] -= sums[ 8] ; base2[kloc0+1] -= sums[ 9] ;
1094 base2[kloc1] -= sums[10] ; base2[kloc1+1] -= sums[11] ;
1095 } else if ( jcolU == ncolU - 1 ) {
1096 ZVdotU31(nrowU, temp0, temp1, temp2, colU0, sums) ;
1097 kloc0 = 2*colindU[jcolU] ;
1098 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1099 base1[kloc0] -= sums[ 2] ; base1[kloc0+1] -= sums[ 3] ;
1100 base2[kloc0] -= sums[ 4] ; base2[kloc0+1] -= sums[ 5] ;
1101 }
1102 rowUT0 = rowUT2 + 2*nrowU ;
1103 }
1104 if ( irowUT == lastInU - 1 ) {
1105 #if MYDEBUG > 0
1106 fprintf(stdout, "\n %% working on rows %d and %d",
1107 colindU[irowUT], colindU[irowUT+1]) ;
1108 fflush(stdout) ;
1109 #endif
1110 rowUT1 = rowUT0 + 2*nrowU ;
1111 /*
1112 -----------------------------------------------------
1113 get base locations for the chevron's diagonal entries
1114 -----------------------------------------------------
1115 */
1116 ichv0 = colindU[irowUT] ;
1117 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
1118 ichv1 = colindU[irowUT+1] ;
1119 base1 = Chv_diagLocation(chvT, ichv1) - 2*ichv1 ;
1120 /*
1121 ------------------------------------
1122 [ temp0 ] [ rowUT0 ]^T
1123 compute [ temp1 ] = [ rowUT1 ] * D
1124 ------------------------------------
1125 */
1126 DVzero(4*nrowU, temp0) ;
1127 SubMtx_scale2vec(mtxD, temp0, temp1, rowUT0, rowUT1) ;
1128 /*
1129 ------------------------------------------------
1130 update the 2x2 upper triangle for these two rows
1131 ------------------------------------------------
1132 */
1133 colU0 = rowUT0 ;
1134 colU1 = rowUT1 ;
1135 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
1136 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
1137 ZVdotU(nrowU, temp0, colU1, &rsum, &isum) ;
1138 base0[2*ichv1] -= rsum ; base0[2*ichv1+1] -= isum ;
1139 ZVdotU(nrowU, temp1, colU1, &rsum, &isum) ;
1140 base1[2*ichv1] -= rsum ; base1[2*ichv1+1] -= isum ;
1141 colU0 = colU1 + 2*nrowU ;
1142 /*
1143 ------------------------------------
1144 update the remainder of the two rows
1145 ------------------------------------
1146 */
1147 for ( jcolU = irowUT + 2 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1148 colU1 = colU0 + 2*nrowU ;
1149 colU2 = colU1 + 2*nrowU ;
1150 ZVdotU23(nrowU, temp0, temp1, colU0, colU1, colU2, sums) ;
1151 kloc0 = 2*colindU[jcolU] ;
1152 kloc1 = 2*colindU[jcolU+1] ;
1153 kloc2 = 2*colindU[jcolU+2] ;
1154 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1155 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1156 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
1157 base1[kloc0] -= sums[ 6] ; base1[kloc0+1] -= sums[ 7] ;
1158 base1[kloc1] -= sums[ 8] ; base1[kloc1+1] -= sums[ 9] ;
1159 base1[kloc2] -= sums[10] ; base1[kloc2+1] -= sums[11] ;
1160 colU0 = colU2 + 2*nrowU ;
1161 }
1162 if ( jcolU == ncolU - 2 ) {
1163 colU1 = colU0 + 2*nrowU ;
1164 ZVdotU22(nrowU, temp0, temp1, colU0, colU1, sums) ;
1165 kloc0 = 2*colindU[jcolU] ;
1166 kloc1 = 2*colindU[jcolU+1] ;
1167 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1168 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1169 base1[kloc0] -= sums[ 4] ; base1[kloc0+1] -= sums[ 5] ;
1170 base1[kloc1] -= sums[ 6] ; base1[kloc1+1] -= sums[ 7] ;
1171 } else if ( jcolU == ncolU - 1 ) {
1172 ZVdotU21(nrowU, temp0, temp1, colU0, sums) ;
1173 kloc0 = 2*colindU[jcolU] ;
1174 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1175 base1[kloc0] -= sums[ 2] ; base1[kloc0+1] -= sums[ 3] ;
1176 }
1177 } else if ( irowUT == lastInU ) {
1178 #if MYDEBUG > 0
1179 fprintf(stdout, "\n %% working on row %d", colindU[irowUT]) ;
1180 fflush(stdout) ;
1181 #endif
1182 /*
1183 -----------------------------------------------------
1184 get base locations for the chevron's diagonal entries
1185 -----------------------------------------------------
1186 */
1187 ichv0 = colindU[irowUT] ;
1188 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
1189 /*
1190 ------------------------------------
1191 compute [ temp0 ] = [ rowUT0 ]^T * D
1192 ------------------------------------
1193 */
1194 DVzero(2*nrowU, temp0) ;
1195 SubMtx_scale1vec(mtxD, temp0, rowUT0) ;
1196 /*
1197 ------------------------------------------
1198 update the 1x1 upper triangle for this row
1199 ------------------------------------------
1200 */
1201 colU0 = rowUT0 ;
1202 ZVdotU(nrowU, temp0, colU0, &rsum, &isum) ;
1203 base0[2*ichv0] -= rsum ; base0[2*ichv0+1] -= isum ;
1204 colU0 = colU0 + 2*nrowU ;
1205 /*
1206 -------------------------------
1207 update the remainder of the row
1208 -------------------------------
1209 */
1210 for ( jcolU = irowUT + 1 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1211 #if MYDEBUG > 0
1212 fprintf(stdout, "\n %% working on columns %d, %d and %d",
1213 jcolU, jcolU+1, jcolU+2) ;
1214 fflush(stdout) ;
1215 #endif
1216 colU1 = colU0 + 2*nrowU ;
1217 colU2 = colU1 + 2*nrowU ;
1218 ZVdotU13(nrowU, temp0, colU0, colU1, colU2, sums) ;
1219 kloc0 = 2*colindU[jcolU] ;
1220 kloc1 = 2*colindU[jcolU+1] ;
1221 kloc2 = 2*colindU[jcolU+2] ;
1222 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1223 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1224 base0[kloc2] -= sums[ 4] ; base0[kloc2+1] -= sums[ 5] ;
1225 colU0 = colU2 + 2*nrowU ;
1226 }
1227 if ( jcolU == ncolU - 2 ) {
1228 #if MYDEBUG > 0
1229 fprintf(stdout, "\n %% working on columns %d and %d",
1230 jcolU, jcolU+1) ;
1231 fflush(stdout) ;
1232 #endif
1233 colU1 = colU0 + 2*nrowU ;
1234 ZVdotU12(nrowU, temp0, colU0, colU1, sums) ;
1235 kloc0 = 2*colindU[jcolU] ;
1236 kloc1 = 2*colindU[jcolU+1] ;
1237 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1238 base0[kloc1] -= sums[ 2] ; base0[kloc1+1] -= sums[ 3] ;
1239 } else if ( jcolU == ncolU - 1 ) {
1240 #if MYDEBUG > 0
1241 fprintf(stdout, "\n %% working on column %d", jcolU) ;
1242 fflush(stdout) ;
1243 #endif
1244 ZVdotU11(nrowU, temp0, colU0, sums) ;
1245 kloc0 = 2*colindU[jcolU] ;
1246 base0[kloc0] -= sums[ 0] ; base0[kloc0+1] -= sums[ 1] ;
1247 }
1248 }
1249 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxU) ) {
1250 double isum, rsum ;
1251 double *base0, *colU0, *entU, *rowUT0, *temp0, *temp1 ;
1252 int ichv0, ii, iloc, irowUT, kloc0, nentU, nrowU, offset,
1253 rloc, sizeU, sizeUT ;
1254 int *indU, *indU0, *indUT0, *sizes ;
1255
1256 SubMtx_sparseColumnsInfo(mtxU,
1257 &ncolU, &nentU, &sizes, &indU, &entU) ;
1258 nrowU = mtxU->nrow ;
1259 DV_setSize(tempDV, 4*nrowU) ;
1260 temp0 = DV_entries(tempDV) ;
1261 temp1 = temp0 + 2*nrowU ;
1262 /*
1263 -------------------------------------------
1264 get the offset into the indices and entries
1265 -------------------------------------------
1266 */
1267 for ( jcolU = offset = 0 ; jcolU < firstInU ; jcolU++ ) {
1268 offset += sizes[jcolU] ;
1269 }
1270 /*
1271 ------------------------------------
1272 loop over the supporting rows of U^T
1273 ------------------------------------
1274 */
1275 rowUT0 = entU + 2*offset ;
1276 indUT0 = indU + offset ;
1277 for ( irowUT = firstInU ; irowUT <= lastInU ; irowUT++ ) {
1278 if ( (sizeUT = sizes[irowUT]) > 0 ) {
1279 /*
1280 -----------------------------------------------------
1281 get base locations for the chevron's diagonal entries
1282 -----------------------------------------------------
1283 */
1284 ichv0 = colindU[irowUT] ;
1285 base0 = Chv_diagLocation(chvT, ichv0) - 2*ichv0 ;
1286 /*
1287 ------------------------------------
1288 compute [ temp0 ] = [ rowUT0 ]^T * D
1289 ------------------------------------
1290 */
1291 DVzero(4*nrowU, temp0) ;
1292 for ( ii = 0 ; ii < sizeUT ; ii++ ) {
1293 rloc = 2*indUT0[ii] ; iloc = rloc + 1 ;
1294 temp1[rloc] = rowUT0[2*ii] ;
1295 temp1[iloc] = rowUT0[2*ii+1] ;
1296 }
1297 SubMtx_scale1vec(mtxD, temp0, temp1) ;
1298 /*
1299 -------------------------------
1300 loop over the following columns
1301 -------------------------------
1302 */
1303 colU0 = rowUT0 ;
1304 indU0 = indUT0 ;
1305 for ( jcolU = irowUT ; jcolU < ncolU ; jcolU++ ) {
1306 if ( (sizeU = sizes[jcolU]) > 0 ) {
1307 ZVdotiU(sizeU, temp0, indU0, colU0, &rsum, &isum) ;
1308 kloc0 = 2*colindU[jcolU] ;
1309 base0[kloc0] -= rsum ; base0[kloc0+1] -= isum ;
1310 colU0 += 2*sizeU ;
1311 indU0 += sizeU ;
1312 }
1313 }
1314 rowUT0 += 2*sizeUT ;
1315 indUT0 += sizeUT ;
1316 }
1317 }
1318 } else {
1319 fprintf(stderr, "\n fatal error in Chv_updateT(%p,%p,%p,%p)"
1320 "\n bad input, mtxU must have dense or sparse columns\n",
1321 chvT, mtxD, mtxU, tempDV) ;
1322 exit(-1) ;
1323 }
1324 }
1325 /*
1326 ---------------------------------------------------
1327 overwrite the local indices with the global indices
1328 ---------------------------------------------------
1329 */
1330 for ( jcolU = firstInU ; jcolU < ncolU ; jcolU++ ) {
1331 colindU[jcolU] = colindT[colindU[jcolU]] ;
1332 }
1333 return ; }
1334
1335 /*--------------------------------------------------------------------*/
1336 /*
1337 ---------------------------------------------------------------------
1338 purpose -- perform the nonsymmetric factor update
1339 T_{\bnd{I} \cap J, \bnd{I} \cap J}
1340 -= L_{\bnd{I} \cap J, I} D_{I, I} U_{I, \bnd{I} \cap J}
1341 and
1342 T_{\bnd{I} \cap J, \bnd{I} \cap \bnd{J}}
1343 -= L_{\bnd{I} \cap J, I} D_{I, I} U_{I, \bnd{I} \cap \bnd{J}}
1344 and
1345 T_{\bnd{I} \cap \bnd{J}, \bnd{I} \cap J}
1346 -= L_{\bnd{I} \cap \bnd{J}, I} D_{I, I} U_{I, \bnd{I} \cap J}
1347
1348 created -- 98feb27, cca
1349 ---------------------------------------------------------------------
1350 */
1351 void
Chv_updateN(Chv * chvT,SubMtx * mtxL,SubMtx * mtxD,SubMtx * mtxU,DV * tempDV)1352 Chv_updateN (
1353 Chv *chvT,
1354 SubMtx *mtxL,
1355 SubMtx *mtxD,
1356 SubMtx *mtxU,
1357 DV *tempDV
1358 ) {
1359 int firstInT, firstInU, jcolT, jcolU, lastInT, lastInU, ncolT, ncolU ;
1360 int *colindT, *colindU ;
1361 /*
1362 ---------------
1363 check the input
1364 ---------------
1365 */
1366 if ( chvT == NULL || mtxD == NULL || mtxU == NULL || tempDV == NULL ) {
1367 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p,%p)"
1368 "\n bad input\n", chvT, mtxL, mtxD, mtxU, tempDV) ;
1369 exit(-1) ;
1370 }
1371 if ( CHV_IS_REAL(chvT) ) {
1372 if ( ! SUBMTX_IS_REAL(mtxD) || ! SUBMTX_IS_REAL(mtxL)
1373 || ! SUBMTX_IS_REAL(mtxU) ) {
1374 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p)"
1375 "\n chvT is real, but mtxD, mtxL and/or mtxU are not\n",
1376 chvT, mtxD, mtxU, tempDV) ;
1377 exit(-1) ;
1378 }
1379 } else if ( CHV_IS_COMPLEX(chvT) ) {
1380 if ( ! SUBMTX_IS_COMPLEX(mtxD) || ! SUBMTX_IS_COMPLEX(mtxL)
1381 || ! SUBMTX_IS_COMPLEX(mtxU) ) {
1382 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p)"
1383 "\n chvT is complex, but mtxD, mtxL and/or mtxU are not\n",
1384 chvT, mtxD, mtxU, tempDV) ;
1385 exit(-1) ;
1386 }
1387 } else {
1388 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p)"
1389 "\n bad input, chvT is not real or complex\n",
1390 chvT, mtxD, mtxU, tempDV) ;
1391 exit(-1) ;
1392 }
1393 Chv_columnIndices(chvT, &ncolT, &colindT) ;
1394 SubMtx_columnIndices(mtxU, &ncolU, &colindU) ;
1395 /*
1396 -----------------------------
1397 locate first column of U in T
1398 -----------------------------
1399 */
1400 firstInT = colindT[0] ;
1401 lastInT = colindT[chvT->nD-1] ;
1402 for ( jcolU = 0 ; jcolU < ncolU ; jcolU++ ) {
1403 if ( firstInT <= colindU[jcolU] && colindU[jcolU] <= lastInT ) {
1404 break ;
1405 }
1406 }
1407 if ( (firstInU = jcolU) == ncolU ) {
1408 return ;
1409 }
1410 /*
1411 ----------------------------
1412 locate last column of U in T
1413 ----------------------------
1414 */
1415 lastInU = firstInU ;
1416 for ( ; jcolU < ncolU ; jcolU++ ) {
1417 if ( colindU[jcolU] <= lastInT ) {
1418 lastInU = jcolU ;
1419 } else {
1420 break ;
1421 }
1422 }
1423 #if MYDEBUG > 0
1424 fprintf(stdout, "\n %% firstInU = %d, lastInU = %d",
1425 firstInU, lastInU) ;
1426 fflush(stdout) ;
1427 #endif
1428 /*
1429 ----------------------------------------------------------
1430 overwrite supported column indices with indices local to T
1431 ----------------------------------------------------------
1432 */
1433 for ( jcolU = firstInU, jcolT = 0 ; jcolU < ncolU ; jcolU++ ) {
1434 while ( colindU[jcolU] != colindT[jcolT] ) {
1435 jcolT++ ;
1436 }
1437 colindU[jcolU] = jcolT ;
1438 }
1439 if ( CHV_IS_REAL(chvT) ) {
1440 if ( SUBMTX_IS_DENSE_COLUMNS(mtxU) && SUBMTX_IS_DENSE_ROWS(mtxL) ) {
1441 double sums[9] ;
1442 double *base0, *base1, *base2, *colU0, *colU1, *colU2, *entL,
1443 *entU, *Ltemp0, *Ltemp1, *Ltemp2, *rowL0, *rowL1, *rowL2,
1444 *Utemp0, *Utemp1, *Utemp2 ;
1445 int ichv0, ichv1, ichv2, inc1, inc2, irowL, jcolU,
1446 loc, loc0, loc1, loc2, ncolL, nrowL, nrowU, offset ;
1447
1448 SubMtx_denseInfo(mtxL, &nrowL, &ncolL, &inc1, &inc2, &entL) ;
1449 SubMtx_denseInfo(mtxU, &nrowU, &ncolU, &inc1, &inc2, &entU) ;
1450 DV_setSize(tempDV, 6*nrowU) ;
1451 Ltemp0 = DV_entries(tempDV) ;
1452 Ltemp1 = Ltemp0 + nrowU ;
1453 Ltemp2 = Ltemp1 + nrowU ;
1454 Utemp0 = Ltemp2 + nrowU ;
1455 Utemp1 = Utemp0 + nrowU ;
1456 Utemp2 = Utemp1 + nrowU ;
1457 /*
1458 -----------------------------------------------------
1459 loop over the rows of L in groups of three
1460 updating the diagonal and upper blocks of the chevron
1461 -----------------------------------------------------
1462 */
1463 offset = firstInU*nrowU ;
1464 for ( irowL = firstInU ; irowL <= lastInU - 2 ; irowL += 3 ) {
1465 #if MYDEBUG > 0
1466 fprintf(stdout, "\n %% working on rows %d, %d and %d",
1467 colindU[irowL], colindU[irowL+1], colindU[irowL+2]) ;
1468 fflush(stdout) ;
1469 #endif
1470 rowL0 = entL + offset ;
1471 rowL1 = rowL0 + nrowU ;
1472 rowL2 = rowL1 + nrowU ;
1473 colU0 = entU + offset ;
1474 colU1 = colU0 + nrowU ;
1475 colU2 = colU1 + nrowU ;
1476 /*
1477 -----------------------------------------------------
1478 get base locations for the chevron's diagonal entries
1479 -----------------------------------------------------
1480 */
1481 ichv0 = colindU[irowL] ;
1482 base0 = Chv_diagLocation(chvT, ichv0) ;
1483 ichv1 = colindU[irowL+1] ;
1484 base1 = Chv_diagLocation(chvT, ichv1) ;
1485 ichv2 = colindU[irowL+2] ;
1486 base2 = Chv_diagLocation(chvT, ichv2) ;
1487 /*
1488 -----------------------------------
1489 [ Ltemp0 ] [ rowL0 ]
1490 compute [ Ltemp1 ] = [ rowL1 ] * D
1491 [ Ltemp2 ] [ rowL2 ]
1492 -----------------------------------
1493 */
1494 DVzero(3*nrowU, Ltemp0) ;
1495 SubMtx_scale3vec(mtxD, Ltemp0, Ltemp1, Ltemp2,
1496 rowL0, rowL1, rowL2) ;
1497 /*
1498 -----------------------------------
1499 [ Utemp0 ] [ colU0 ]
1500 compute [ Utemp1 ] = D * [ colU0 ]
1501 [ Utemp2 ] [ colU0 ]
1502 -----------------------------------
1503 */
1504 DVzero(3*nrowU, Utemp0) ;
1505 SubMtx_scale3vec(mtxD, Utemp0, Utemp1, Utemp2,
1506 colU0, colU1, colU2) ;
1507 /*
1508 --------------------------------------------------------------
1509 update the 3x3 diagonal block for these three rows and columns
1510 --------------------------------------------------------------
1511 */
1512 DVdot33(nrowU, Ltemp0, Ltemp1, Ltemp2,
1513 colU0, colU1, colU2, sums) ;
1514 /*
1515 -------------------------------------
1516 inject the nine sums into the chevron
1517 -------------------------------------
1518 */
1519 base0[0] -= sums[0] ;
1520 loc = ichv1 - ichv0 ;
1521 base0[loc] -= sums[1] ;
1522 base0[-loc] -= sums[3] ;
1523 loc = ichv2 - ichv0 ;
1524 base0[loc] -= sums[2] ;
1525 base0[-loc] -= sums[6] ;
1526 base1[0] -= sums[4] ;
1527 loc = ichv2 - ichv1 ;
1528 base1[loc] -= sums[5] ;
1529 base1[-loc] -= sums[7] ;
1530 base2[0] -= sums[8] ;
1531 /*
1532 ------------------------------------------
1533 update the lower and upper blocks together
1534 ------------------------------------------
1535 */
1536 rowL0 = rowL2 + nrowU ;
1537 colU0 = colU2 + nrowU ;
1538 for ( jcolU = irowL + 3 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1539 rowL1 = rowL0 + nrowU ;
1540 rowL2 = rowL1 + nrowU ;
1541 colU1 = colU0 + nrowU ;
1542 colU2 = colU1 + nrowU ;
1543 /*
1544 ------------------------------------------------------
1545 get the local indices for these three rows and columns
1546 ------------------------------------------------------
1547 */
1548 loc0 = colindU[jcolU] ;
1549 loc1 = colindU[jcolU+1] ;
1550 loc2 = colindU[jcolU+2] ;
1551 /*
1552 ---------------------
1553 upper 3x3 block first
1554 ---------------------
1555 */
1556 DVdot33(nrowU, Ltemp0, Ltemp1, Ltemp2,
1557 colU0, colU1, colU2, sums) ;
1558 /*
1559 -------------------------------------
1560 inject the nine sums into the chevron
1561 -------------------------------------
1562 */
1563 base0 -= ichv0 ;
1564 base1 -= ichv1 ;
1565 base2 -= ichv2 ;
1566 base0[loc0] -= sums[0] ;
1567 base0[loc1] -= sums[1] ;
1568 base0[loc2] -= sums[2] ;
1569 base1[loc0] -= sums[3] ;
1570 base1[loc1] -= sums[4] ;
1571 base1[loc2] -= sums[5] ;
1572 base2[loc0] -= sums[6] ;
1573 base2[loc1] -= sums[7] ;
1574 base2[loc2] -= sums[8] ;
1575 base0 += ichv0 ;
1576 base1 += ichv1 ;
1577 base2 += ichv2 ;
1578 /*
1579 ----------------------
1580 lower 3x3 block second
1581 ----------------------
1582 */
1583 DVdot33(nrowU, rowL0, rowL1, rowL2,
1584 Utemp0, Utemp1, Utemp2, sums) ;
1585 /*
1586 -------------------------------------
1587 inject the nine sums into the chevron
1588 -------------------------------------
1589 */
1590 base0 += ichv0 ;
1591 base1 += ichv1 ;
1592 base2 += ichv2 ;
1593 base0[-loc0] -= sums[0] ;
1594 base1[-loc0] -= sums[1] ;
1595 base2[-loc0] -= sums[2] ;
1596 base0[-loc1] -= sums[3] ;
1597 base1[-loc1] -= sums[4] ;
1598 base2[-loc1] -= sums[5] ;
1599 base0[-loc2] -= sums[6] ;
1600 base1[-loc2] -= sums[7] ;
1601 base2[-loc2] -= sums[8] ;
1602 base0 -= ichv0 ;
1603 base1 -= ichv1 ;
1604 base2 -= ichv2 ;
1605 /*
1606 -------------------------------------
1607 increment the row and column pointers
1608 -------------------------------------
1609 */
1610 rowL0 = rowL2 + nrowU ;
1611 colU0 = colU2 + nrowU ;
1612 }
1613 if ( jcolU == ncolU - 2 ) {
1614 rowL1 = rowL0 + nrowU ;
1615 colU1 = colU0 + nrowU ;
1616 /*
1617 ----------------------------------------------------
1618 get the local indices for these two rows and columns
1619 ----------------------------------------------------
1620 */
1621 loc0 = colindU[jcolU] ;
1622 loc1 = colindU[jcolU+1] ;
1623 /*
1624 ---------------------
1625 upper 3x2 block first
1626 ---------------------
1627 */
1628 DVdot32(nrowU, Ltemp0, Ltemp1, Ltemp2, colU0, colU1, sums) ;
1629 /*
1630 ------------------------------------
1631 inject the six sums into the chevron
1632 ------------------------------------
1633 */
1634 base0 -= ichv0 ;
1635 base1 -= ichv1 ;
1636 base2 -= ichv2 ;
1637 base0[loc0] -= sums[0] ;
1638 base0[loc1] -= sums[1] ;
1639 base1[loc0] -= sums[2] ;
1640 base1[loc1] -= sums[3] ;
1641 base2[loc0] -= sums[4] ;
1642 base2[loc1] -= sums[5] ;
1643 base0 += ichv0 ;
1644 base1 += ichv1 ;
1645 base2 += ichv2 ;
1646 /*
1647 ----------------------
1648 lower 2x3 block second
1649 ----------------------
1650 */
1651 DVdot23(nrowU, rowL0, rowL1, Utemp0, Utemp1, Utemp2, sums) ;
1652 /*
1653 ------------------------------------
1654 inject the six sums into the chevron
1655 ------------------------------------
1656 */
1657 base0 += ichv0 ;
1658 base1 += ichv1 ;
1659 base2 += ichv2 ;
1660 base0[-loc0] -= sums[0] ;
1661 base1[-loc0] -= sums[1] ;
1662 base2[-loc0] -= sums[2] ;
1663 base0[-loc1] -= sums[3] ;
1664 base1[-loc1] -= sums[4] ;
1665 base2[-loc1] -= sums[5] ;
1666 base0 -= ichv0 ;
1667 base1 -= ichv1 ;
1668 base2 -= ichv2 ;
1669 } else if ( jcolU == ncolU - 1 ) {
1670 /*
1671 ---------------------------------------------
1672 get the local indices for this row and column
1673 ---------------------------------------------
1674 */
1675 loc0 = colindU[jcolU] ;
1676 /*
1677 ---------------------
1678 upper 3x1 block first
1679 ---------------------
1680 */
1681 DVdot31(nrowU, Ltemp0, Ltemp1, Ltemp2, colU0, sums) ;
1682 /*
1683 --------------------------------------
1684 inject the three sums into the chevron
1685 --------------------------------------
1686 */
1687 base0 -= ichv0 ;
1688 base1 -= ichv1 ;
1689 base2 -= ichv2 ;
1690 base0[loc0] -= sums[0] ;
1691 base1[loc0] -= sums[1] ;
1692 base2[loc0] -= sums[2] ;
1693 base0 += ichv0 ;
1694 base1 += ichv1 ;
1695 base2 += ichv2 ;
1696 /*
1697 ----------------------
1698 lower 1x3 block second
1699 ----------------------
1700 */
1701 DVdot13(nrowU, rowL0, Utemp0, Utemp1, Utemp2, sums) ;
1702 /*
1703 --------------------------------------
1704 inject the three sums into the chevron
1705 --------------------------------------
1706 */
1707 base0 += ichv0 ;
1708 base1 += ichv1 ;
1709 base2 += ichv2 ;
1710 base0[-loc0] -= sums[0] ;
1711 base1[-loc0] -= sums[1] ;
1712 base2[-loc0] -= sums[2] ;
1713 base0 -= ichv0 ;
1714 base1 -= ichv1 ;
1715 base2 -= ichv2 ;
1716 }
1717 offset += 3*nrowU ;
1718 }
1719 if ( irowL == lastInU - 1 ) {
1720 #if MYDEBUG > 0
1721 fprintf(stdout, "\n %% working on rows %d and %d",
1722 colindU[irowL], colindU[irowL+1]) ;
1723 fflush(stdout) ;
1724 #endif
1725 rowL0 = entL + offset ;
1726 rowL1 = rowL0 + nrowU ;
1727 colU0 = entU + offset ;
1728 colU1 = colU0 + nrowU ;
1729 /*
1730 -----------------------------------------------------
1731 get base locations for the chevron's diagonal entries
1732 -----------------------------------------------------
1733 */
1734 ichv0 = colindU[irowL] ;
1735 base0 = Chv_diagLocation(chvT, ichv0) ;
1736 ichv1 = colindU[irowL+1] ;
1737 base1 = Chv_diagLocation(chvT, ichv1) ;
1738 /*
1739 -----------------------------------
1740 [ Ltemp0 ] [ rowL0 ]
1741 compute [ Ltemp1 ] = [ rowL1 ] * D
1742 -----------------------------------
1743 */
1744 DVzero(2*nrowU, Ltemp0) ;
1745 SubMtx_scale2vec(mtxD, Ltemp0, Ltemp1, rowL0, rowL1) ;
1746 /*
1747 -----------------------------------
1748 [ Utemp0 ] [ colU0 ]
1749 compute [ Utemp1 ] = D * [ colU0 ]
1750 -----------------------------------
1751 */
1752 DVzero(2*nrowU, Utemp0) ;
1753 SubMtx_scale2vec(mtxD, Utemp0, Utemp1, colU0, colU1) ;
1754 /*
1755 ------------------------------------------------------------
1756 update the 2x2 diagonal block for these two rows and columns
1757 ------------------------------------------------------------
1758 */
1759 DVdot22(nrowU, Ltemp0, Ltemp1, colU0, colU1, sums) ;
1760 /*
1761 -------------------------------------
1762 inject the four sums into the chevron
1763 -------------------------------------
1764 */
1765 base0[0] -= sums[ 0] ;
1766 loc = ichv1 - ichv0 ;
1767 base0[loc] -= sums[ 1] ;
1768 base0[-loc] -= sums[ 2] ;
1769 base1[0] -= sums[ 3] ;
1770 /*
1771 ------------------------------------------
1772 update the lower and upper blocks together
1773 ------------------------------------------
1774 */
1775 rowL0 = rowL1 + nrowU ;
1776 colU0 = colU1 + nrowU ;
1777 for ( jcolU = irowL + 2 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1778 rowL1 = rowL0 + nrowU ;
1779 rowL2 = rowL1 + nrowU ;
1780 colU1 = colU0 + nrowU ;
1781 colU2 = colU1 + nrowU ;
1782 /*
1783 ------------------------------------------------------
1784 get the local indices for these three rows and columns
1785 ------------------------------------------------------
1786 */
1787 loc0 = colindU[jcolU] ;
1788 loc1 = colindU[jcolU+1] ;
1789 loc2 = colindU[jcolU+2] ;
1790 /*
1791 ---------------------
1792 upper 2x3 block first
1793 ---------------------
1794 */
1795 DVdot23(nrowU, Ltemp0, Ltemp1, colU0, colU1, colU2, sums) ;
1796 /*
1797 ------------------------------------
1798 inject the six sums into the chevron
1799 ------------------------------------
1800 */
1801 base0 -= ichv0 ;
1802 base1 -= ichv1 ;
1803 base0[loc0] -= sums[0] ;
1804 base0[loc1] -= sums[1] ;
1805 base0[loc2] -= sums[2] ;
1806 base1[loc0] -= sums[3] ;
1807 base1[loc1] -= sums[4] ;
1808 base1[loc2] -= sums[5] ;
1809 base0 += ichv0 ;
1810 base1 += ichv1 ;
1811 /*
1812 ----------------------
1813 lower 3x2 block second
1814 ----------------------
1815 */
1816 DVdot32(nrowU, rowL0, rowL1, rowL2, Utemp0, Utemp1, sums) ;
1817 /*
1818 ------------------------------------
1819 inject the six sums into the chevron
1820 ------------------------------------
1821 */
1822 base0 += ichv0 ;
1823 base1 += ichv1 ;
1824 base0[-loc0] -= sums[0] ;
1825 base1[-loc0] -= sums[1] ;
1826 base0[-loc1] -= sums[2] ;
1827 base1[-loc1] -= sums[3] ;
1828 base0[-loc2] -= sums[4] ;
1829 base1[-loc2] -= sums[5] ;
1830 base0 -= ichv0 ;
1831 base1 -= ichv1 ;
1832 /*
1833 -------------------------------------
1834 increment the row and column pointers
1835 -------------------------------------
1836 */
1837 rowL0 = rowL2 + nrowU ;
1838 colU0 = colU2 + nrowU ;
1839 }
1840 if ( jcolU == ncolU - 2 ) {
1841 rowL1 = rowL0 + nrowU ;
1842 colU1 = colU0 + nrowU ;
1843 /*
1844 ----------------------------------------------------
1845 get the local indices for these two rows and columns
1846 ----------------------------------------------------
1847 */
1848 loc0 = colindU[jcolU] ;
1849 loc1 = colindU[jcolU+1] ;
1850 /*
1851 ---------------------
1852 upper 2x2 block first
1853 ---------------------
1854 */
1855 DVdot22(nrowU, Ltemp0, Ltemp1, colU0, colU1, sums) ;
1856 /*
1857 -------------------------------------
1858 inject the four sums into the chevron
1859 -------------------------------------
1860 */
1861 base0 -= ichv0 ;
1862 base1 -= ichv1 ;
1863 base0[loc0] -= sums[0] ;
1864 base0[loc1] -= sums[1] ;
1865 base1[loc0] -= sums[2] ;
1866 base1[loc1] -= sums[3] ;
1867 base0 += ichv0 ;
1868 base1 += ichv1 ;
1869 /*
1870 ----------------------
1871 lower 2x2 block second
1872 ----------------------
1873 */
1874 DVdot22(nrowU, rowL0, rowL1, Utemp0, Utemp1, sums) ;
1875 /*
1876 -------------------------------------
1877 inject the four sums into the chevron
1878 -------------------------------------
1879 */
1880 base0 += ichv0 ;
1881 base1 += ichv1 ;
1882 base0[-loc0] -= sums[0] ;
1883 base1[-loc0] -= sums[1] ;
1884 base0[-loc1] -= sums[2] ;
1885 base1[-loc1] -= sums[3] ;
1886 base0 -= ichv0 ;
1887 base1 -= ichv1 ;
1888 } else if ( jcolU == ncolU - 1 ) {
1889 /*
1890 ---------------------------------------------
1891 get the local indices for this row and column
1892 ---------------------------------------------
1893 */
1894 loc0 = colindU[jcolU] ;
1895 /*
1896 ---------------------
1897 upper 2x1 block first
1898 ---------------------
1899 */
1900 DVdot21(nrowU, Ltemp0, Ltemp1, colU0, sums) ;
1901 /*
1902 ------------------------------------
1903 inject the two sums into the chevron
1904 ------------------------------------
1905 */
1906 base0 -= ichv0 ;
1907 base1 -= ichv1 ;
1908 base0[loc0] -= sums[0] ;
1909 base1[loc0] -= sums[1] ;
1910 base0 += ichv0 ;
1911 base1 += ichv1 ;
1912 /*
1913 ----------------------
1914 lower 1x2 block second
1915 ----------------------
1916 */
1917 DVdot12(nrowU, rowL0, Utemp0, Utemp1, sums) ;
1918 /*
1919 ------------------------------------
1920 inject the two sums into the chevron
1921 ------------------------------------
1922 */
1923 base0 += ichv0 ;
1924 base1 += ichv1 ;
1925 base0[-loc0] -= sums[0] ;
1926 base1[-loc0] -= sums[1] ;
1927 base0 -= ichv0 ;
1928 base1 -= ichv1 ;
1929 }
1930 } else if ( irowL == lastInU ) {
1931 #if MYDEBUG > 0
1932 fprintf(stdout, "\n %% working on row %d", colindU[irowL]) ;
1933 fflush(stdout) ;
1934 #endif
1935 rowL0 = entL + offset ;
1936 colU0 = entU + offset ;
1937 /*
1938 -----------------------------------------------------
1939 get base locations for the chevron's diagonal entries
1940 -----------------------------------------------------
1941 */
1942 ichv0 = colindU[irowL] ;
1943 base0 = Chv_diagLocation(chvT, ichv0) ;
1944 /*
1945 -----------------------------------
1946 compute [ Ltemp0 ] = [ rowL0 ] * D
1947 -----------------------------------
1948 */
1949 DVzero(nrowU, Ltemp0) ;
1950 SubMtx_scale1vec(mtxD, Ltemp0, rowL0) ;
1951 /*
1952 -----------------------------------
1953 compute [ Utemp0 ] = D * [ colU0 ]
1954 -----------------------------------
1955 */
1956 DVzero(nrowU, Utemp0) ;
1957 SubMtx_scale1vec(mtxD, Utemp0, colU0) ;
1958 /*
1959 ------------------------------------------------------------
1960 update the 1x1 diagonal block for these two rows and columns
1961 ------------------------------------------------------------
1962 */
1963 DVdot11(nrowU, Ltemp0, colU0, sums) ;
1964 /*
1965 -------------------------------
1966 inject the sum into the chevron
1967 -------------------------------
1968 */
1969 base0[0] -= sums[0] ;
1970 /*
1971 ------------------------------------------
1972 update the lower and upper blocks together
1973 ------------------------------------------
1974 */
1975 rowL0 = rowL0 + nrowU ;
1976 colU0 = colU0 + nrowU ;
1977 for ( jcolU = irowL + 1 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
1978 rowL1 = rowL0 + nrowU ;
1979 rowL2 = rowL1 + nrowU ;
1980 colU1 = colU0 + nrowU ;
1981 colU2 = colU1 + nrowU ;
1982 /*
1983 ------------------------------------------------------
1984 get the local indices for these three rows and columns
1985 ------------------------------------------------------
1986 */
1987 loc0 = colindU[jcolU] ;
1988 loc1 = colindU[jcolU+1] ;
1989 loc2 = colindU[jcolU+2] ;
1990 /*
1991 ---------------------
1992 upper 1x3 block first
1993 ---------------------
1994 */
1995 DVdot13(nrowU, Ltemp0, colU0, colU1, colU2, sums) ;
1996 /*
1997 --------------------------------------
1998 inject the three sums into the chevron
1999 --------------------------------------
2000 */
2001 base0 -= ichv0 ;
2002 base0[loc0] -= sums[0] ;
2003 base0[loc1] -= sums[1] ;
2004 base0[loc2] -= sums[2] ;
2005 base0 += ichv0 ;
2006 /*
2007 ----------------------
2008 lower 3x1 block second
2009 ----------------------
2010 */
2011 DVdot31(nrowU, rowL0, rowL1, rowL2, Utemp0, sums) ;
2012 /*
2013 --------------------------------------
2014 inject the three sums into the chevron
2015 --------------------------------------
2016 */
2017 base0 += ichv0 ;
2018 base0[-loc0] -= sums[0] ;
2019 base0[-loc1] -= sums[1] ;
2020 base0[-loc2] -= sums[2] ;
2021 base0 -= ichv0 ;
2022 /*
2023 -------------------------------------
2024 increment the row and column pointers
2025 -------------------------------------
2026 */
2027 rowL0 = rowL2 + nrowU ;
2028 colU0 = colU2 + nrowU ;
2029 }
2030 if ( jcolU == ncolU - 2 ) {
2031 rowL1 = rowL0 + nrowU ;
2032 colU1 = colU0 + nrowU ;
2033 /*
2034 ----------------------------------------------------
2035 get the local indices for these two rows and columns
2036 ----------------------------------------------------
2037 */
2038 loc0 = colindU[jcolU] ;
2039 loc1 = colindU[jcolU+1] ;
2040 /*
2041 ---------------------
2042 upper 1x2 block first
2043 ---------------------
2044 */
2045 DVdot12(nrowU, Ltemp0, colU0, colU1, sums) ;
2046 /*
2047 ------------------------------------
2048 inject the two sums into the chevron
2049 ------------------------------------
2050 */
2051 base0 -= ichv0 ;
2052 base0[loc0] -= sums[0] ;
2053 base0[loc1] -= sums[1] ;
2054 base0 += ichv0 ;
2055 /*
2056 ----------------------
2057 lower 2x1 block second
2058 ----------------------
2059 */
2060 DVdot21(nrowU, rowL0, rowL1, Utemp0, sums) ;
2061 /*
2062 ------------------------------------
2063 inject the two sums into the chevron
2064 ------------------------------------
2065 */
2066 base0 += ichv0 ;
2067 base0[-loc0] -= sums[0] ;
2068 base0[-loc1] -= sums[1] ;
2069 base0 -= ichv0 ;
2070 } else if ( jcolU == ncolU - 1 ) {
2071 /*
2072 ---------------------------------------------
2073 get the local indices for this row and column
2074 ---------------------------------------------
2075 */
2076 loc0 = colindU[jcolU] ;
2077 /*
2078 ---------------------
2079 upper 1x1 block first
2080 ---------------------
2081 */
2082 DVdot11(nrowU, Ltemp0, colU0, sums) ;
2083 /*
2084 -------------------------------
2085 inject the sum into the chevron
2086 -------------------------------
2087 */
2088 base0 -= ichv0 ;
2089 base0[loc0] -= sums[0] ;
2090 base0 += ichv0 ;
2091 /*
2092 ----------------------
2093 lower 1x1 block second
2094 ----------------------
2095 */
2096 DVdot11(nrowU, rowL0, Utemp0, sums) ;
2097 /*
2098 -------------------------------
2099 inject the sum into the chevron
2100 -------------------------------
2101 */
2102 base0 += ichv0 ;
2103 base0[-loc0] -= sums[0] ;
2104 base0 -= ichv0 ;
2105 }
2106 }
2107 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxU)
2108 && SUBMTX_IS_SPARSE_ROWS(mtxL) ) {
2109 double *base, *colU0, *colU1, *entL, *entU, *rowL0, *rowL1,
2110 *temp0, *temp1, *temp2 ;
2111 int ichv, ii, irow0, irow1, jj, loc, ncolL, ncolU, nentL,
2112 nentU, nrowL, nrowU, offsetL, offsetU, sizeL0, sizeL1,
2113 sizeU0, sizeU1 ;
2114 int *indL, *indL0, *indL1, *indU, *indU0, *indU1,
2115 *sizesL, *sizesU ;
2116
2117 SubMtx_sparseColumnsInfo(mtxU,
2118 &ncolU, &nentU, &sizesU, &indU, &entU) ;
2119 SubMtx_sparseRowsInfo(mtxL,
2120 &nrowL, &nentL, &sizesL, &indL, &entL) ;
2121 nrowU = mtxU->nrow ;
2122 ncolL = mtxL->ncol ;
2123 DV_setSize(tempDV, 3*nrowU) ;
2124 temp0 = DV_entries(tempDV) ;
2125 temp1 = temp0 + nrowU ;
2126 temp2 = temp1 + nrowU ;
2127 /*
2128 --------------------------------------------
2129 get the offsets into the indices and entries
2130 --------------------------------------------
2131 */
2132 for ( jcolU = offsetL = offsetU = 0 ;
2133 jcolU < firstInU ;
2134 jcolU++ ) {
2135 offsetL += sizesL[jcolU] ;
2136 offsetU += sizesU[jcolU] ;
2137 }
2138 #if MYDEBUG > 0
2139 fprintf(stdout, "\n\n %% offsetL %d, offsetU %d",
2140 offsetL, offsetU) ;
2141 fflush(stdout) ;
2142 #endif
2143 /*
2144 ---------------------------------------------------
2145 loop over the supporting rows of L and columns of U
2146 ---------------------------------------------------
2147 */
2148 for ( irow0 = firstInU ; irow0 <= lastInU ; irow0++ ) {
2149 rowL0 = entL + offsetL ;
2150 indL0 = indL + offsetL ;
2151 colU0 = entU + offsetU ;
2152 indU0 = indU + offsetU ;
2153 sizeL0 = sizesL[irow0] ;
2154 sizeU0 = sizesU[irow0] ;
2155 #if MYDEBUG > 0
2156 fprintf(stdout,
2157 "\n\n %% irow0 %d, offsetL %d, offsetU %d, sizeL0 %d, sizeU0 %d",
2158 irow0, offsetL, offsetU, sizeL0, sizeU0) ;
2159 fflush(stdout) ;
2160 #endif
2161 if ( sizeL0 > 0 || sizeU0 > 0 ) {
2162 /*
2163 -----------------------------------------------------
2164 get base locations for the chevron's diagonal entries
2165 -----------------------------------------------------
2166 */
2167 ichv = colindU[irow0] ;
2168 base = Chv_diagLocation(chvT, ichv) ;
2169 #if MYDEBUG > 0
2170 fprintf(stdout, "\n\n %% ichv = %d, base = %p", ichv, base) ;
2171 fflush(stdout) ;
2172 #endif
2173 if ( sizeL0 > 0 ) {
2174 /*
2175 ----------------------------------
2176 compute [ temp0 ] = D * [ rowL0 ]
2177 ----------------------------------
2178 */
2179 #if MYDEBUG > 0
2180 fprintf(stdout, "\n\n %% loading temp1") ;
2181 fflush(stdout) ;
2182 #endif
2183 DVzero(2*nrowU, temp0) ;
2184 for ( ii = 0 ; ii < sizeL0 ; ii++ ) {
2185 jj = indL0[ii] ;
2186 #if MYDEBUG > 0
2187 fprintf(stdout, "\n\n %% ii = %d, jj = %d", ii, jj) ;
2188 fflush(stdout) ;
2189 #endif
2190 temp1[jj] = rowL0[ii] ;
2191 }
2192 SubMtx_scale1vec(mtxD, temp0, temp1) ;
2193 #if MYDEBUG > 0
2194 fprintf(stdout, "\n\n %% temp0 = L * D") ;
2195 DVfprintf(stdout, 2*nrowU, temp0) ;
2196 fflush(stdout) ;
2197 #endif
2198 }
2199 if ( sizeU0 > 0 ) {
2200 /*
2201 ----------------------------------
2202 compute [ temp1 ] = D * [ colU0 ]
2203 ----------------------------------
2204 */
2205 #if MYDEBUG > 0
2206 fprintf(stdout, "\n\n %% loading temp2") ;
2207 fflush(stdout) ;
2208 #endif
2209 DVzero(2*nrowU, temp1) ;
2210 for ( ii = 0 ; ii < sizeU0 ; ii++ ) {
2211 jj = indU0[ii] ;
2212 #if MYDEBUG > 0
2213 fprintf(stdout, "\n\n %% ii = %d, jj = %d", ii, jj) ;
2214 fflush(stdout) ;
2215 #endif
2216 temp2[jj] = colU0[ii] ;
2217 }
2218 SubMtx_scale1vec(mtxD, temp1, temp2) ;
2219 #if MYDEBUG > 0
2220 fprintf(stdout, "\n\n %% temp1 = D * U") ;
2221 DVfprintf(stdout, nrowU, temp1) ;
2222 fflush(stdout) ;
2223 #endif
2224 }
2225 if ( sizeL0 > 0 && sizeU0 > 0 ) {
2226 /*
2227 -------------------------
2228 update the diagonal entry
2229 -------------------------
2230 */
2231 base[0] -= DVdoti(sizeU0, temp0, indU0, colU0) ;
2232 }
2233 /*
2234 ----------------------------------------
2235 loop over the following rows and columns
2236 ----------------------------------------
2237 */
2238 if ( sizeU0 > 0 ) {
2239 /*
2240 ------------------------
2241 update the lower entries
2242 ------------------------
2243 */
2244 base += ichv ;
2245 rowL1 = rowL0 + sizeL0 ;
2246 indL1 = indL0 + sizeL0 ;
2247 for ( irow1 = irow0 + 1 ; irow1 < ncolU ; irow1++ ) {
2248 #if MYDEBUG > 0
2249 fprintf(stdout, "\n\n %% irow1 %d, sizeL1 %d",
2250 irow1, sizesL[irow1]) ;
2251 fflush(stdout) ;
2252 #endif
2253 if ( (sizeL1 = sizesL[irow1]) > 0 ) {
2254 loc = colindU[irow1] ;
2255 #if MYDEBUG > 0
2256 fprintf(stdout,
2257 "\n\n %% base[%d] = %12.4e, base[%d] = %12.4e",
2258 -loc, base[-loc], -loc + 1, base[-loc+1]) ;
2259 fflush(stdout) ;
2260 #endif
2261 base[-loc] -= DVdoti(sizeL1, temp1, indL1, rowL1) ;
2262 rowL1 += sizeL1 ;
2263 indL1 += sizeL1 ;
2264 }
2265 }
2266 base -= ichv ;
2267 }
2268 if ( sizeL0 > 0 ) {
2269 /*
2270 ------------------------
2271 update the upper entries
2272 ------------------------
2273 */
2274 base -= ichv ;
2275 colU1 = colU0 + sizeU0 ;
2276 indU1 = indU0 + sizeU0 ;
2277 for ( irow1 = irow0 + 1 ; irow1 < ncolU ; irow1++ ) {
2278 #if MYDEBUG > 0
2279 fprintf(stdout, "\n\n %% irow1 %d, sizeU1 %d",
2280 irow1, sizesU[irow1]) ;
2281 fflush(stdout) ;
2282 #endif
2283 if ( (sizeU1 = sizesU[irow1]) > 0 ) {
2284 loc = colindU[irow1] ;
2285 base[loc] -= DVdoti(sizeU1, temp0, indU1, colU1) ;
2286 colU1 += sizeU1 ;
2287 indU1 += sizeU1 ;
2288 }
2289 }
2290 base -= ichv ;
2291 }
2292 }
2293 offsetL += sizeL0 ;
2294 offsetU += sizeU0 ;
2295 }
2296 } else {
2297 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p)"
2298 "\n bad input, mtxU must have dense or sparse columns"
2299 "\n and mtxL must have dense or sparse rows\n",
2300 chvT, mtxD, mtxU, tempDV) ;
2301 exit(-1) ;
2302 }
2303 } else if ( CHV_IS_COMPLEX(chvT) ) {
2304 if ( SUBMTX_IS_DENSE_COLUMNS(mtxU) && SUBMTX_IS_DENSE_ROWS(mtxL) ) {
2305 double sums[18] ;
2306 double *base0, *base1, *base2, *colU0, *colU1, *colU2, *entL,
2307 *entU, *Ltemp0, *Ltemp1, *Ltemp2, *rowL0, *rowL1, *rowL2,
2308 *Utemp0, *Utemp1, *Utemp2 ;
2309 int ichv0, ichv1, ichv2, inc1, inc2, irowL, jcolU,
2310 loc, loc0, loc1, loc2, ncolL, nrowL, nrowU, offset ;
2311
2312 SubMtx_denseInfo(mtxL, &nrowL, &ncolL, &inc1, &inc2, &entL) ;
2313 SubMtx_denseInfo(mtxU, &nrowU, &ncolU, &inc1, &inc2, &entU) ;
2314 DV_setSize(tempDV, 12*nrowU) ;
2315 Ltemp0 = DV_entries(tempDV) ;
2316 Ltemp1 = Ltemp0 + 2*nrowU ;
2317 Ltemp2 = Ltemp1 + 2*nrowU ;
2318 Utemp0 = Ltemp2 + 2*nrowU ;
2319 Utemp1 = Utemp0 + 2*nrowU ;
2320 Utemp2 = Utemp1 + 2*nrowU ;
2321 /*
2322 -----------------------------------------------------
2323 loop over the rows of L in groups of three
2324 updating the diagonal and upper blocks of the chevron
2325 -----------------------------------------------------
2326 */
2327 offset = firstInU*nrowU ;
2328 for ( irowL = firstInU ; irowL <= lastInU - 2 ; irowL += 3 ) {
2329 #if MYDEBUG > 0
2330 fprintf(stdout, "\n %% working on rows %d, %d and %d",
2331 colindU[irowL], colindU[irowL+1], colindU[irowL+2]) ;
2332 fflush(stdout) ;
2333 #endif
2334 rowL0 = entL + 2*offset ;
2335 rowL1 = rowL0 + 2*nrowU ;
2336 rowL2 = rowL1 + 2*nrowU ;
2337 colU0 = entU + 2*offset ;
2338 colU1 = colU0 + 2*nrowU ;
2339 colU2 = colU1 + 2*nrowU ;
2340 /*
2341 -----------------------------------------------------
2342 get base locations for the chevron's diagonal entries
2343 -----------------------------------------------------
2344 */
2345 ichv0 = colindU[irowL] ;
2346 base0 = Chv_diagLocation(chvT, ichv0) ;
2347 ichv1 = colindU[irowL+1] ;
2348 base1 = Chv_diagLocation(chvT, ichv1) ;
2349 ichv2 = colindU[irowL+2] ;
2350 base2 = Chv_diagLocation(chvT, ichv2) ;
2351 /*
2352 -----------------------------------
2353 [ Ltemp0 ] [ rowL0 ]
2354 compute [ Ltemp1 ] = [ rowL1 ] * D
2355 [ Ltemp2 ] [ rowL2 ]
2356 -----------------------------------
2357 */
2358 DVzero(6*nrowU, Ltemp0) ;
2359 SubMtx_scale3vec(mtxD, Ltemp0, Ltemp1, Ltemp2,
2360 rowL0, rowL1, rowL2) ;
2361 /*
2362 -----------------------------------
2363 [ Utemp0 ] [ colU0 ]
2364 compute [ Utemp1 ] = D * [ colU0 ]
2365 [ Utemp2 ] [ colU0 ]
2366 -----------------------------------
2367 */
2368 DVzero(6*nrowU, Utemp0) ;
2369 SubMtx_scale3vec(mtxD, Utemp0, Utemp1, Utemp2,
2370 colU0, colU1, colU2) ;
2371 /*
2372 --------------------------------------------------------------
2373 update the 3x3 diagonal block for these three rows and columns
2374 --------------------------------------------------------------
2375 */
2376 ZVdotU33(nrowU, Ltemp0, Ltemp1, Ltemp2,
2377 colU0, colU1, colU2, sums) ;
2378 /*
2379 -------------------------------------
2380 inject the nine sums into the chevron
2381 -------------------------------------
2382 */
2383 base0[0] -= sums[ 0] ; base0[1] -= sums[ 1] ;
2384 loc = 2*(ichv1 - ichv0) ;
2385 base0[loc] -= sums[ 2] ; base0[loc+1] -= sums[ 3] ;
2386 base0[-loc] -= sums[ 6] ; base0[-loc+1] -= sums[ 7] ;
2387 loc = 2*(ichv2 - ichv0) ;
2388 base0[loc] -= sums[ 4] ; base0[loc+1] -= sums[ 5] ;
2389 base0[-loc] -= sums[12] ; base0[-loc+1] -= sums[13] ;
2390 base1[0] -= sums[ 8] ; base1[1] -= sums[ 9] ;
2391 loc = 2*(ichv2 - ichv1) ;
2392 base1[loc] -= sums[10] ; base1[loc+1] -= sums[11] ;
2393 base1[-loc] -= sums[14] ; base1[-loc+1] -= sums[15] ;
2394 base2[0] -= sums[16] ; base2[1] -= sums[17] ;
2395 /*
2396 ------------------------------------------
2397 update the lower and upper blocks together
2398 ------------------------------------------
2399 */
2400 rowL0 = rowL2 + 2*nrowU ;
2401 colU0 = colU2 + 2*nrowU ;
2402 for ( jcolU = irowL + 3 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
2403 rowL1 = rowL0 + 2*nrowU ;
2404 rowL2 = rowL1 + 2*nrowU ;
2405 colU1 = colU0 + 2*nrowU ;
2406 colU2 = colU1 + 2*nrowU ;
2407 /*
2408 ------------------------------------------------------
2409 get the local indices for these three rows and columns
2410 ------------------------------------------------------
2411 */
2412 loc0 = 2*colindU[jcolU] ;
2413 loc1 = 2*colindU[jcolU+1] ;
2414 loc2 = 2*colindU[jcolU+2] ;
2415 /*
2416 ---------------------
2417 upper 3x3 block first
2418 ---------------------
2419 */
2420 ZVdotU33(nrowU, Ltemp0, Ltemp1, Ltemp2,
2421 colU0, colU1, colU2, sums) ;
2422 /*
2423 -------------------------------------
2424 inject the nine sums into the chevron
2425 -------------------------------------
2426 */
2427 base0 -= 2*ichv0 ;
2428 base1 -= 2*ichv1 ;
2429 base2 -= 2*ichv2 ;
2430 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2431 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2432 base0[loc2] -= sums[ 4] ; base0[loc2+1] -= sums[ 5] ;
2433 base1[loc0] -= sums[ 6] ; base1[loc0+1] -= sums[ 7] ;
2434 base1[loc1] -= sums[ 8] ; base1[loc1+1] -= sums[ 9] ;
2435 base1[loc2] -= sums[10] ; base1[loc2+1] -= sums[11] ;
2436 base2[loc0] -= sums[12] ; base2[loc0+1] -= sums[13] ;
2437 base2[loc1] -= sums[14] ; base2[loc1+1] -= sums[15] ;
2438 base2[loc2] -= sums[16] ; base2[loc2+1] -= sums[17] ;
2439 base0 += 2*ichv0 ;
2440 base1 += 2*ichv1 ;
2441 base2 += 2*ichv2 ;
2442 /*
2443 ----------------------
2444 lower 3x3 block second
2445 ----------------------
2446 */
2447 ZVdotU33(nrowU, rowL0, rowL1, rowL2,
2448 Utemp0, Utemp1, Utemp2, sums) ;
2449 /*
2450 -------------------------------------
2451 inject the nine sums into the chevron
2452 -------------------------------------
2453 */
2454 base0 += 2*ichv0 ;
2455 base1 += 2*ichv1 ;
2456 base2 += 2*ichv2 ;
2457 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2458 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2459 base2[-loc0] -= sums[ 4] ; base2[-loc0+1] -= sums[ 5] ;
2460 base0[-loc1] -= sums[ 6] ; base0[-loc1+1] -= sums[ 7] ;
2461 base1[-loc1] -= sums[ 8] ; base1[-loc1+1] -= sums[ 9] ;
2462 base2[-loc1] -= sums[10] ; base2[-loc1+1] -= sums[11] ;
2463 base0[-loc2] -= sums[12] ; base0[-loc2+1] -= sums[13] ;
2464 base1[-loc2] -= sums[14] ; base1[-loc2+1] -= sums[15] ;
2465 base2[-loc2] -= sums[16] ; base2[-loc2+1] -= sums[17] ;
2466 base0 -= 2*ichv0 ;
2467 base1 -= 2*ichv1 ;
2468 base2 -= 2*ichv2 ;
2469 /*
2470 -------------------------------------
2471 increment the row and column pointers
2472 -------------------------------------
2473 */
2474 rowL0 = rowL2 + 2*nrowU ;
2475 colU0 = colU2 + 2*nrowU ;
2476 }
2477 if ( jcolU == ncolU - 2 ) {
2478 rowL1 = rowL0 + 2*nrowU ;
2479 colU1 = colU0 + 2*nrowU ;
2480 /*
2481 ----------------------------------------------------
2482 get the local indices for these two rows and columns
2483 ----------------------------------------------------
2484 */
2485 loc0 = 2*colindU[jcolU] ;
2486 loc1 = 2*colindU[jcolU+1] ;
2487 /*
2488 ---------------------
2489 upper 3x2 block first
2490 ---------------------
2491 */
2492 ZVdotU32(nrowU, Ltemp0, Ltemp1, Ltemp2,
2493 colU0, colU1, sums) ;
2494 /*
2495 ------------------------------------
2496 inject the six sums into the chevron
2497 ------------------------------------
2498 */
2499 base0 -= 2*ichv0 ;
2500 base1 -= 2*ichv1 ;
2501 base2 -= 2*ichv2 ;
2502 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2503 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2504 base1[loc0] -= sums[ 4] ; base1[loc0+1] -= sums[ 5] ;
2505 base1[loc1] -= sums[ 6] ; base1[loc1+1] -= sums[ 7] ;
2506 base2[loc0] -= sums[ 8] ; base2[loc0+1] -= sums[ 9] ;
2507 base2[loc1] -= sums[10] ; base2[loc1+1] -= sums[11] ;
2508 base0 += 2*ichv0 ;
2509 base1 += 2*ichv1 ;
2510 base2 += 2*ichv2 ;
2511 /*
2512 ----------------------
2513 lower 2x3 block second
2514 ----------------------
2515 */
2516 ZVdotU23(nrowU, rowL0, rowL1,
2517 Utemp0, Utemp1, Utemp2, sums) ;
2518 /*
2519 ------------------------------------
2520 inject the six sums into the chevron
2521 ------------------------------------
2522 */
2523 base0 += 2*ichv0 ;
2524 base1 += 2*ichv1 ;
2525 base2 += 2*ichv2 ;
2526 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2527 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2528 base2[-loc0] -= sums[ 4] ; base2[-loc0+1] -= sums[ 5] ;
2529 base0[-loc1] -= sums[ 6] ; base0[-loc1+1] -= sums[ 7] ;
2530 base1[-loc1] -= sums[ 8] ; base1[-loc1+1] -= sums[ 9] ;
2531 base2[-loc1] -= sums[10] ; base2[-loc1+1] -= sums[11] ;
2532 base0 -= 2*ichv0 ;
2533 base1 -= 2*ichv1 ;
2534 base2 -= 2*ichv2 ;
2535 } else if ( jcolU == ncolU - 1 ) {
2536 /*
2537 ---------------------------------------------
2538 get the local indices for this row and column
2539 ---------------------------------------------
2540 */
2541 loc0 = 2*colindU[jcolU] ;
2542 /*
2543 ---------------------
2544 upper 3x1 block first
2545 ---------------------
2546 */
2547 ZVdotU31(nrowU, Ltemp0, Ltemp1, Ltemp2, colU0, sums) ;
2548 /*
2549 --------------------------------------
2550 inject the three sums into the chevron
2551 --------------------------------------
2552 */
2553 base0 -= 2*ichv0 ;
2554 base1 -= 2*ichv1 ;
2555 base2 -= 2*ichv2 ;
2556 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2557 base1[loc0] -= sums[ 2] ; base1[loc0+1] -= sums[ 3] ;
2558 base2[loc0] -= sums[ 4] ; base2[loc0+1] -= sums[ 5] ;
2559 base0 += 2*ichv0 ;
2560 base1 += 2*ichv1 ;
2561 base2 += 2*ichv2 ;
2562 /*
2563 ----------------------
2564 lower 1x3 block second
2565 ----------------------
2566 */
2567 ZVdotU13(nrowU, rowL0, Utemp0, Utemp1, Utemp2, sums) ;
2568 /*
2569 --------------------------------------
2570 inject the three sums into the chevron
2571 --------------------------------------
2572 */
2573 base0 += 2*ichv0 ;
2574 base1 += 2*ichv1 ;
2575 base2 += 2*ichv2 ;
2576 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2577 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2578 base2[-loc0] -= sums[ 4] ; base2[-loc0+1] -= sums[ 5] ;
2579 base0 -= 2*ichv0 ;
2580 base1 -= 2*ichv1 ;
2581 base2 -= 2*ichv2 ;
2582 }
2583 offset += 3*nrowU ;
2584 }
2585 if ( irowL == lastInU - 1 ) {
2586 #if MYDEBUG > 0
2587 fprintf(stdout, "\n %% working on rows %d and %d",
2588 colindU[irowL], colindU[irowL+1]) ;
2589 fflush(stdout) ;
2590 #endif
2591 rowL0 = entL + 2*offset ;
2592 rowL1 = rowL0 + 2*nrowU ;
2593 colU0 = entU + 2*offset ;
2594 colU1 = colU0 + 2*nrowU ;
2595 /*
2596 -----------------------------------------------------
2597 get base locations for the chevron's diagonal entries
2598 -----------------------------------------------------
2599 */
2600 ichv0 = colindU[irowL] ;
2601 base0 = Chv_diagLocation(chvT, ichv0) ;
2602 ichv1 = colindU[irowL+1] ;
2603 base1 = Chv_diagLocation(chvT, ichv1) ;
2604 /*
2605 -----------------------------------
2606 [ Ltemp0 ] [ rowL0 ]
2607 compute [ Ltemp1 ] = [ rowL1 ] * D
2608 -----------------------------------
2609 */
2610 DVzero(4*nrowU, Ltemp0) ;
2611 SubMtx_scale2vec(mtxD, Ltemp0, Ltemp1, rowL0, rowL1) ;
2612 /*
2613 -----------------------------------
2614 [ Utemp0 ] [ colU0 ]
2615 compute [ Utemp1 ] = D * [ colU0 ]
2616 -----------------------------------
2617 */
2618 DVzero(4*nrowU, Utemp0) ;
2619 SubMtx_scale2vec(mtxD, Utemp0, Utemp1, colU0, colU1) ;
2620 /*
2621 ------------------------------------------------------------
2622 update the 2x2 diagonal block for these two rows and columns
2623 ------------------------------------------------------------
2624 */
2625 ZVdotU22(nrowU, Ltemp0, Ltemp1, colU0, colU1, sums) ;
2626 /*
2627 -------------------------------------
2628 inject the four sums into the chevron
2629 -------------------------------------
2630 */
2631 base0[0] -= sums[ 0] ; base0[1] -= sums[ 1] ;
2632 loc = 2*(ichv1 - ichv0) ;
2633 base0[loc] -= sums[ 2] ; base0[loc+1] -= sums[ 3] ;
2634 base0[-loc] -= sums[ 4] ; base0[-loc+1] -= sums[ 5] ;
2635 base1[0] -= sums[ 6] ; base1[1] -= sums[ 7] ;
2636 /*
2637 ------------------------------------------
2638 update the lower and upper blocks together
2639 ------------------------------------------
2640 */
2641 rowL0 = rowL1 + 2*nrowU ;
2642 colU0 = colU1 + 2*nrowU ;
2643 for ( jcolU = irowL + 2 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
2644 rowL1 = rowL0 + 2*nrowU ;
2645 rowL2 = rowL1 + 2*nrowU ;
2646 colU1 = colU0 + 2*nrowU ;
2647 colU2 = colU1 + 2*nrowU ;
2648 /*
2649 ------------------------------------------------------
2650 get the local indices for these three rows and columns
2651 ------------------------------------------------------
2652 */
2653 loc0 = 2*colindU[jcolU] ;
2654 loc1 = 2*colindU[jcolU+1] ;
2655 loc2 = 2*colindU[jcolU+2] ;
2656 /*
2657 ---------------------
2658 upper 2x3 block first
2659 ---------------------
2660 */
2661 ZVdotU23(nrowU, Ltemp0, Ltemp1, colU0, colU1, colU2, sums) ;
2662 /*
2663 ------------------------------------
2664 inject the six sums into the chevron
2665 ------------------------------------
2666 */
2667 base0 -= 2*ichv0 ;
2668 base1 -= 2*ichv1 ;
2669 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2670 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2671 base0[loc2] -= sums[ 4] ; base0[loc2+1] -= sums[ 5] ;
2672 base1[loc0] -= sums[ 6] ; base1[loc0+1] -= sums[ 7] ;
2673 base1[loc1] -= sums[ 8] ; base1[loc1+1] -= sums[ 9] ;
2674 base1[loc2] -= sums[10] ; base1[loc2+1] -= sums[11] ;
2675 base0 += 2*ichv0 ;
2676 base1 += 2*ichv1 ;
2677 /*
2678 ----------------------
2679 lower 3x2 block second
2680 ----------------------
2681 */
2682 ZVdotU32(nrowU, rowL0, rowL1, rowL2, Utemp0, Utemp1, sums) ;
2683 /*
2684 ------------------------------------
2685 inject the six sums into the chevron
2686 ------------------------------------
2687 */
2688 base0 += 2*ichv0 ;
2689 base1 += 2*ichv1 ;
2690 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2691 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2692 base0[-loc1] -= sums[ 4] ; base0[-loc1+1] -= sums[ 5] ;
2693 base1[-loc1] -= sums[ 6] ; base1[-loc1+1] -= sums[ 7] ;
2694 base0[-loc2] -= sums[ 8] ; base0[-loc2+1] -= sums[ 9] ;
2695 base1[-loc2] -= sums[10] ; base1[-loc2+1] -= sums[11] ;
2696 base0 -= 2*ichv0 ;
2697 base1 -= 2*ichv1 ;
2698 /*
2699 -------------------------------------
2700 increment the row and column pointers
2701 -------------------------------------
2702 */
2703 rowL0 = rowL2 + 2*nrowU ;
2704 colU0 = colU2 + 2*nrowU ;
2705 }
2706 if ( jcolU == ncolU - 2 ) {
2707 rowL1 = rowL0 + 2*nrowU ;
2708 colU1 = colU0 + 2*nrowU ;
2709 /*
2710 ----------------------------------------------------
2711 get the local indices for these two rows and columns
2712 ----------------------------------------------------
2713 */
2714 loc0 = 2*colindU[jcolU] ;
2715 loc1 = 2*colindU[jcolU+1] ;
2716 /*
2717 ---------------------
2718 upper 2x2 block first
2719 ---------------------
2720 */
2721 ZVdotU22(nrowU, Ltemp0, Ltemp1, colU0, colU1, sums) ;
2722 /*
2723 -------------------------------------
2724 inject the four sums into the chevron
2725 -------------------------------------
2726 */
2727 base0 -= 2*ichv0 ;
2728 base1 -= 2*ichv1 ;
2729 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2730 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2731 base1[loc0] -= sums[ 4] ; base1[loc0+1] -= sums[ 5] ;
2732 base1[loc1] -= sums[ 6] ; base1[loc1+1] -= sums[ 7] ;
2733 base0 += 2*ichv0 ;
2734 base1 += 2*ichv1 ;
2735 /*
2736 ----------------------
2737 lower 2x2 block second
2738 ----------------------
2739 */
2740 ZVdotU22(nrowU, rowL0, rowL1, Utemp0, Utemp1, sums) ;
2741 /*
2742 -------------------------------------
2743 inject the four sums into the chevron
2744 -------------------------------------
2745 */
2746 base0 += 2*ichv0 ;
2747 base1 += 2*ichv1 ;
2748 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2749 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2750 base0[-loc1] -= sums[ 4] ; base0[-loc1+1] -= sums[ 5] ;
2751 base1[-loc1] -= sums[ 6] ; base1[-loc1+1] -= sums[ 7] ;
2752 base0 -= 2*ichv0 ;
2753 base1 -= 2*ichv1 ;
2754 } else if ( jcolU == ncolU - 1 ) {
2755 /*
2756 ---------------------------------------------
2757 get the local indices for this row and column
2758 ---------------------------------------------
2759 */
2760 loc0 = 2*colindU[jcolU] ;
2761 /*
2762 ---------------------
2763 upper 2x1 block first
2764 ---------------------
2765 */
2766 ZVdotU21(nrowU, Ltemp0, Ltemp1, colU0, sums) ;
2767 /*
2768 ------------------------------------
2769 inject the two sums into the chevron
2770 ------------------------------------
2771 */
2772 base0 -= 2*ichv0 ;
2773 base1 -= 2*ichv1 ;
2774 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2775 base1[loc0] -= sums[ 2] ; base1[loc0+1] -= sums[ 3] ;
2776 base0 += 2*ichv0 ;
2777 base1 += 2*ichv1 ;
2778 /*
2779 ----------------------
2780 lower 1x2 block second
2781 ----------------------
2782 */
2783 ZVdotU12(nrowU, rowL0, Utemp0, Utemp1, sums) ;
2784 /*
2785 ------------------------------------
2786 inject the two sums into the chevron
2787 ------------------------------------
2788 */
2789 base0 += 2*ichv0 ;
2790 base1 += 2*ichv1 ;
2791 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2792 base1[-loc0] -= sums[ 2] ; base1[-loc0+1] -= sums[ 3] ;
2793 base0 -= 2*ichv0 ;
2794 base1 -= 2*ichv1 ;
2795 }
2796 } else if ( irowL == lastInU ) {
2797 #if MYDEBUG > 0
2798 fprintf(stdout, "\n %% working on row %d", colindU[irowL]) ;
2799 fflush(stdout) ;
2800 #endif
2801 rowL0 = entL + 2*offset ;
2802 colU0 = entU + 2*offset ;
2803 /*
2804 -----------------------------------------------------
2805 get base locations for the chevron's diagonal entries
2806 -----------------------------------------------------
2807 */
2808 ichv0 = colindU[irowL] ;
2809 base0 = Chv_diagLocation(chvT, ichv0) ;
2810 /*
2811 -----------------------------------
2812 compute [ Ltemp0 ] = [ rowL0 ] * D
2813 -----------------------------------
2814 */
2815 DVzero(2*nrowU, Ltemp0) ;
2816 SubMtx_scale1vec(mtxD, Ltemp0, rowL0) ;
2817 /*
2818 -----------------------------------
2819 compute [ Utemp0 ] = D * [ colU0 ]
2820 -----------------------------------
2821 */
2822 DVzero(2*nrowU, Utemp0) ;
2823 SubMtx_scale1vec(mtxD, Utemp0, colU0) ;
2824 /*
2825 ------------------------------------------------------------
2826 update the 1x1 diagonal block for these two rows and columns
2827 ------------------------------------------------------------
2828 */
2829 ZVdotU11(nrowU, Ltemp0, colU0, sums) ;
2830 /*
2831 -------------------------------
2832 inject the sum into the chevron
2833 -------------------------------
2834 */
2835 base0[0] -= sums[ 0] ; base0[1] -= sums[ 1] ;
2836 /*
2837 ------------------------------------------
2838 update the lower and upper blocks together
2839 ------------------------------------------
2840 */
2841 rowL0 = rowL0 + 2*nrowU ;
2842 colU0 = colU0 + 2*nrowU ;
2843 for ( jcolU = irowL + 1 ; jcolU < ncolU - 2 ; jcolU += 3 ) {
2844 rowL1 = rowL0 + 2*nrowU ;
2845 rowL2 = rowL1 + 2*nrowU ;
2846 colU1 = colU0 + 2*nrowU ;
2847 colU2 = colU1 + 2*nrowU ;
2848 /*
2849 ------------------------------------------------------
2850 get the local indices for these three rows and columns
2851 ------------------------------------------------------
2852 */
2853 loc0 = 2*colindU[jcolU] ;
2854 loc1 = 2*colindU[jcolU+1] ;
2855 loc2 = 2*colindU[jcolU+2] ;
2856 /*
2857 ---------------------
2858 upper 1x3 block first
2859 ---------------------
2860 */
2861 ZVdotU13(nrowU, Ltemp0, colU0, colU1, colU2, sums) ;
2862 /*
2863 --------------------------------------
2864 inject the three sums into the chevron
2865 --------------------------------------
2866 */
2867 base0 -= 2*ichv0 ;
2868 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2869 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2870 base0[loc2] -= sums[ 4] ; base0[loc2+1] -= sums[ 5] ;
2871 base0 += 2*ichv0 ;
2872 /*
2873 ----------------------
2874 lower 3x1 block second
2875 ----------------------
2876 */
2877 ZVdotU31(nrowU, rowL0, rowL1, rowL2, Utemp0, sums) ;
2878 /*
2879 --------------------------------------
2880 inject the three sums into the chevron
2881 --------------------------------------
2882 */
2883 base0 += 2*ichv0 ;
2884 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2885 base0[-loc1] -= sums[ 2] ; base0[-loc1+1] -= sums[ 3] ;
2886 base0[-loc2] -= sums[ 4] ; base0[-loc2+1] -= sums[ 5] ;
2887 base0 -= 2*ichv0 ;
2888 /*
2889 -------------------------------------
2890 increment the row and column pointers
2891 -------------------------------------
2892 */
2893 rowL0 = rowL2 + 2*nrowU ;
2894 colU0 = colU2 + 2*nrowU ;
2895 }
2896 if ( jcolU == ncolU - 2 ) {
2897 rowL1 = rowL0 + 2*nrowU ;
2898 colU1 = colU0 + 2*nrowU ;
2899 /*
2900 ----------------------------------------------------
2901 get the local indices for these two rows and columns
2902 ----------------------------------------------------
2903 */
2904 loc0 = 2*colindU[jcolU] ;
2905 loc1 = 2*colindU[jcolU+1] ;
2906 /*
2907 ---------------------
2908 upper 1x2 block first
2909 ---------------------
2910 */
2911 ZVdotU12(nrowU, Ltemp0, colU0, colU1, sums) ;
2912 /*
2913 ------------------------------------
2914 inject the two sums into the chevron
2915 ------------------------------------
2916 */
2917 base0 -= 2*ichv0 ;
2918 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2919 base0[loc1] -= sums[ 2] ; base0[loc1+1] -= sums[ 3] ;
2920 base0 += 2*ichv0 ;
2921 /*
2922 ----------------------
2923 lower 2x1 block second
2924 ----------------------
2925 */
2926 ZVdotU21(nrowU, rowL0, rowL1, Utemp0, sums) ;
2927 /*
2928 ------------------------------------
2929 inject the two sums into the chevron
2930 ------------------------------------
2931 */
2932 base0 += 2*ichv0 ;
2933 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2934 base0[-loc1] -= sums[ 2] ; base0[-loc1+1] -= sums[ 3] ;
2935 base0 -= 2*ichv0 ;
2936 } else if ( jcolU == ncolU - 1 ) {
2937 /*
2938 ---------------------------------------------
2939 get the local indices for this row and column
2940 ---------------------------------------------
2941 */
2942 loc0 = 2*colindU[jcolU] ;
2943 /*
2944 ---------------------
2945 upper 1x1 block first
2946 ---------------------
2947 */
2948 ZVdotU11(nrowU, Ltemp0, colU0, sums) ;
2949 /*
2950 -------------------------------
2951 inject the sum into the chevron
2952 -------------------------------
2953 */
2954 base0 -= 2*ichv0 ;
2955 base0[loc0] -= sums[ 0] ; base0[loc0+1] -= sums[ 1] ;
2956 base0 += 2*ichv0 ;
2957 /*
2958 ----------------------
2959 lower 1x1 block second
2960 ----------------------
2961 */
2962 ZVdotU11(nrowU, rowL0, Utemp0, sums) ;
2963 /*
2964 -------------------------------
2965 inject the sum into the chevron
2966 -------------------------------
2967 */
2968 base0 += 2*ichv0 ;
2969 base0[-loc0] -= sums[ 0] ; base0[-loc0+1] -= sums[ 1] ;
2970 base0 -= 2*ichv0 ;
2971 }
2972 }
2973 } else if ( SUBMTX_IS_SPARSE_COLUMNS(mtxU)
2974 && SUBMTX_IS_SPARSE_ROWS(mtxL) ) {
2975 double imag, real ;
2976 double *base, *colU0, *colU1, *entL, *entU, *rowL0, *rowL1,
2977 *temp0, *temp1, *temp2 ;
2978 int ichv, ii, irow0, irow1, jj, loc, ncolL, ncolU, nentL,
2979 nentU, nrowL, nrowU, offsetL, offsetU, sizeL0, sizeL1,
2980 sizeU0, sizeU1 ;
2981 int *indL, *indL0, *indL1, *indU, *indU0, *indU1,
2982 *sizesL, *sizesU ;
2983
2984 SubMtx_sparseColumnsInfo(mtxU,
2985 &ncolU, &nentU, &sizesU, &indU, &entU) ;
2986 SubMtx_sparseRowsInfo(mtxL,
2987 &nrowL, &nentL, &sizesL, &indL, &entL) ;
2988 nrowU = mtxU->nrow ;
2989 ncolL = mtxL->ncol ;
2990 DV_setSize(tempDV, 6*nrowU) ;
2991 temp0 = DV_entries(tempDV) ;
2992 temp1 = temp0 + 2*nrowU ;
2993 temp2 = temp1 + 2*nrowU ;
2994 /*
2995 --------------------------------------------
2996 get the offsets into the indices and entries
2997 --------------------------------------------
2998 */
2999 for ( jcolU = offsetL = offsetU = 0 ;
3000 jcolU < firstInU ;
3001 jcolU++ ) {
3002 offsetL += sizesL[jcolU] ;
3003 offsetU += sizesU[jcolU] ;
3004 }
3005 #if MYDEBUG > 0
3006 fprintf(stdout, "\n\n %% offsetL %d, offsetU %d",
3007 offsetL, offsetU) ;
3008 fflush(stdout) ;
3009 #endif
3010 /*
3011 ---------------------------------------------------
3012 loop over the supporting rows of L and columns of U
3013 ---------------------------------------------------
3014 */
3015 for ( irow0 = firstInU ; irow0 <= lastInU ; irow0++ ) {
3016 rowL0 = entL + 2*offsetL ;
3017 indL0 = indL + offsetL ;
3018 colU0 = entU + 2*offsetU ;
3019 indU0 = indU + offsetU ;
3020 sizeL0 = sizesL[irow0] ;
3021 sizeU0 = sizesU[irow0] ;
3022 #if MYDEBUG > 0
3023 fprintf(stdout,
3024 "\n\n %% irow0 %d, offsetL %d, offsetU %d, sizeL0 %d, sizeU0 %d",
3025 irow0, offsetL, offsetU, sizeL0, sizeU0) ;
3026 fflush(stdout) ;
3027 #endif
3028 if ( sizeL0 > 0 || sizeU0 > 0 ) {
3029 /*
3030 -----------------------------------------------------
3031 get base locations for the chevron's diagonal entries
3032 -----------------------------------------------------
3033 */
3034 ichv = colindU[irow0] ;
3035 base = Chv_diagLocation(chvT, ichv) ;
3036 #if MYDEBUG > 0
3037 fprintf(stdout, "\n\n %% ichv = %d, base = %p", ichv, base) ;
3038 fflush(stdout) ;
3039 #endif
3040 if ( sizeL0 > 0 ) {
3041 /*
3042 ----------------------------------
3043 compute [ temp0 ] = D * [ rowL0 ]
3044 ----------------------------------
3045 */
3046 #if MYDEBUG > 0
3047 fprintf(stdout, "\n\n %% loading temp1") ;
3048 fflush(stdout) ;
3049 #endif
3050 DVzero(4*nrowU, temp0) ;
3051 for ( ii = 0 ; ii < sizeL0 ; ii++ ) {
3052 jj = indL0[ii] ;
3053 #if MYDEBUG > 0
3054 fprintf(stdout, "\n\n %% ii = %d, jj = %d", ii, jj) ;
3055 fflush(stdout) ;
3056 #endif
3057 temp1[2*jj] = rowL0[2*ii] ;
3058 temp1[2*jj+1] = rowL0[2*ii+1] ;
3059 #if MYDEBUG > 0
3060 fprintf(stdout, ", (%12.4e,%12.4e)",
3061 rowL0[2*ii], rowL0[2*ii+1]) ;
3062 fflush(stdout) ;
3063 #endif
3064 }
3065 SubMtx_scale1vec(mtxD, temp0, temp1) ;
3066 #if MYDEBUG > 0
3067 fprintf(stdout, "\n\n %% temp0 = L * D") ;
3068 DVfprintf(stdout, 2*nrowU, temp0) ;
3069 fflush(stdout) ;
3070 #endif
3071 }
3072 if ( sizeU0 > 0 ) {
3073 /*
3074 ----------------------------------
3075 compute [ temp1 ] = D * [ colU0 ]
3076 ----------------------------------
3077 */
3078 #if MYDEBUG > 0
3079 fprintf(stdout, "\n\n %% loading temp2") ;
3080 fflush(stdout) ;
3081 #endif
3082 DVzero(4*nrowU, temp1) ;
3083 for ( ii = 0 ; ii < sizeU0 ; ii++ ) {
3084 jj = indU0[ii] ;
3085 #if MYDEBUG > 0
3086 fprintf(stdout, "\n\n %% ii = %d, jj = %d", ii, jj) ;
3087 fflush(stdout) ;
3088 #endif
3089 temp2[2*jj] = colU0[2*ii] ;
3090 temp2[2*jj+1] = colU0[2*ii+1] ;
3091 #if MYDEBUG > 0
3092 fprintf(stdout, ", (%12.4e,%12.4e)",
3093 colU0[2*ii], colU0[2*ii+1]) ;
3094 fflush(stdout) ;
3095 #endif
3096 }
3097 SubMtx_scale1vec(mtxD, temp1, temp2) ;
3098 #if MYDEBUG > 0
3099 fprintf(stdout, "\n\n %% temp1 = D * U") ;
3100 DVfprintf(stdout, 2*nrowU, temp1) ;
3101 fflush(stdout) ;
3102 #endif
3103 }
3104 if ( sizeL0 > 0 && sizeU0 > 0 ) {
3105 /*
3106 -------------------------
3107 update the diagonal entry
3108 -------------------------
3109 */
3110 ZVdotiU(sizeU0, temp0, indU0, colU0, &real, &imag) ;
3111 base[0] -= real ; base[1] -= imag ;
3112 #if MYDEBUG > 0
3113 fprintf(stdout,
3114 "\n\n %% sum = %12.4e + %12.4e*i, "
3115 "\n base[%d] = %12.4e, base[%d] = %12.4e ",
3116 real, imag, 0, base[0], 1, base[1]) ;
3117 fflush(stdout) ;
3118 #endif
3119 }
3120 /*
3121 ----------------------------------------
3122 loop over the following rows and columns
3123 ----------------------------------------
3124 */
3125 if ( sizeU0 > 0 ) {
3126 /*
3127 ------------------------
3128 update the lower entries
3129 ------------------------
3130 */
3131 base += 2*ichv ;
3132 rowL1 = rowL0 + 2*sizeL0 ;
3133 indL1 = indL0 + sizeL0 ;
3134 for ( irow1 = irow0 + 1 ; irow1 < ncolU ; irow1++ ) {
3135 #if MYDEBUG > 0
3136 fprintf(stdout, "\n\n %% irow1 %d, sizeL1 %d",
3137 irow1, sizesL[irow1]) ;
3138 fflush(stdout) ;
3139 #endif
3140 if ( (sizeL1 = sizesL[irow1]) > 0 ) {
3141 loc = 2*colindU[irow1] ;
3142 #if MYDEBUG > 0
3143 fprintf(stdout,
3144 "\n\n %% base[%d] = %12.4e, base[%d] = %12.4e",
3145 -loc, base[-loc], -loc + 1, base[-loc+1]) ;
3146 fflush(stdout) ;
3147 #endif
3148 ZVdotiU(sizeL1, temp1, indL1, rowL1, &real, &imag);
3149 base[-loc] -= real ; base[-loc+1] -= imag ;
3150 #if MYDEBUG > 0
3151 fprintf(stdout,
3152 "\n\n %% sum = %12.4e + %12.4e*i, "
3153 "\n base[%d] = %12.4e, base[%d] = %12.4e ",
3154 real, imag, -loc, base[-loc],
3155 -loc + 1, base[-loc+1]) ;
3156 fflush(stdout) ;
3157 #endif
3158 rowL1 += 2*sizeL1 ;
3159 indL1 += sizeL1 ;
3160 }
3161 }
3162 base -= 2*ichv ;
3163 }
3164 if ( sizeL0 > 0 ) {
3165 /*
3166 ------------------------
3167 update the upper entries
3168 ------------------------
3169 */
3170 base -= 2*ichv ;
3171 colU1 = colU0 + 2*sizeU0 ;
3172 indU1 = indU0 + sizeU0 ;
3173 for ( irow1 = irow0 + 1 ; irow1 < ncolU ; irow1++ ) {
3174 #if MYDEBUG > 0
3175 fprintf(stdout, "\n\n %% irow1 %d, sizeU1 %d",
3176 irow1, sizesU[irow1]) ;
3177 fflush(stdout) ;
3178 #endif
3179 if ( (sizeU1 = sizesU[irow1]) > 0 ) {
3180 loc = 2*colindU[irow1] ;
3181 #if MYDEBUG > 0
3182 fprintf(stdout,
3183 "\n\n %% base[%d] = %12.4e, base[%d] = %12.4e",
3184 loc, base[loc], loc + 1, base[loc+1]) ;
3185 fflush(stdout) ;
3186 #endif
3187 ZVdotiU(sizeU1, temp0, indU1, colU1, &real, &imag);
3188 base[loc] -= real ; base[loc+1] -= imag ;
3189 #if MYDEBUG > 0
3190 fprintf(stdout,
3191 "\n\n %% sum = %12.4e + %12.4e*i, "
3192 "\n base[%d] = %12.4e, base[%d] = %12.4e ",
3193 real, imag, loc, base[loc], loc + 1, base[loc+1]) ;
3194 fflush(stdout) ;
3195 #endif
3196 colU1 += 2*sizeU1 ;
3197 indU1 += sizeU1 ;
3198 }
3199 }
3200 base -= ichv ;
3201 }
3202 }
3203 offsetL += sizeL0 ;
3204 offsetU += sizeU0 ;
3205 }
3206 } else {
3207 fprintf(stderr, "\n fatal error in Chv_updateN(%p,%p,%p,%p)"
3208 "\n bad input, mtxU must have dense or sparse columns"
3209 "\n and mtxL must have dense or sparse rows\n",
3210 chvT, mtxD, mtxU, tempDV) ;
3211 exit(-1) ;
3212 }
3213 }
3214 /*
3215 ---------------------------------------------------
3216 overwrite the local indices with the global indices
3217 ---------------------------------------------------
3218 */
3219 for ( jcolU = firstInU ; jcolU < ncolU ; jcolU++ ) {
3220 colindU[jcolU] = colindT[colindU[jcolU]] ;
3221 }
3222 return ; }
3223
3224 /*--------------------------------------------------------------------*/
3225