1 /* solveupdT.c */
2
3 #include "../SubMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 static void
7 real_updDenseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
8 static void
9 real_updDenseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
10 static void
11 real_updSparseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
12 static void
13 real_updSparseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
14 static void
15 complex_updDenseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
16 static void
17 complex_updDenseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
18 static void
19 complex_updSparseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
20 static void
21 complex_updSparseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
22 /*--------------------------------------------------------------------*/
23 /*
24 ------------------------------------------------------
25 purpose -- perform the matrix-matrix multiply
26 Y := Y - A^T * X used in the forward and backsolves
27 where
28 (1) rows(A) \subseteq rows(Y)
29 (2) rows(A) are local w.r.t. rows(Y)
30 (3) cols(A) \subseteq rows(X)
31 (4) cols(A) are local w.r.t. rows(X)
32 (5) cols(Y) = cols(X)
33 (6) Y and X have type SUBMTX_DENSE_COLUMNS
34
35 created -- 98may02, cca
36 -----------------------------------------------------
37 */
38 void
SubMtx_solveupdT(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)39 SubMtx_solveupdT (
40 SubMtx *mtxY,
41 SubMtx *mtxA,
42 SubMtx *mtxX
43 ) {
44 /*
45 ---------------
46 check the input
47 ---------------
48 */
49 if ( mtxY == NULL || mtxA == NULL || mtxX == NULL ) {
50 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
51 "\n bad input\n", mtxY, mtxA, mtxX) ;
52 exit(-1) ;
53 }
54 if ( ! SUBMTX_IS_DENSE_COLUMNS(mtxY) ) {
55 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
56 "\n Y must have mode SUBMTX_DENSE_COLUMNS\n",
57 mtxY, mtxA, mtxX) ;
58 exit(-1) ;
59 }
60 if ( ! SUBMTX_IS_DENSE_COLUMNS(mtxX) ) {
61 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
62 "\n X must have mode SUBMTX_DENSE_COLUMNS\n",
63 mtxY, mtxA, mtxX) ;
64 exit(-1) ;
65 }
66 switch ( mtxA->type ) {
67 case SPOOLES_REAL :
68 switch ( mtxA->mode ) {
69 case SUBMTX_DENSE_COLUMNS :
70 real_updDenseColumns(mtxY, mtxA, mtxX) ;
71 break ;
72 case SUBMTX_DENSE_ROWS :
73 real_updDenseRows(mtxY, mtxA, mtxX) ;
74 break ;
75 case SUBMTX_SPARSE_ROWS :
76 real_updSparseRows(mtxY, mtxA, mtxX) ;
77 break ;
78 case SUBMTX_SPARSE_COLUMNS :
79 real_updSparseColumns(mtxY, mtxA, mtxX) ;
80 break ;
81 default :
82 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
83 "\n unsupported mode %d for A\n",
84 mtxY, mtxA, mtxX, mtxA->mode) ;
85 exit(-1) ;
86 break ;
87 }
88 break ;
89 case SPOOLES_COMPLEX :
90 switch ( mtxA->mode ) {
91 case SUBMTX_DENSE_COLUMNS :
92 complex_updDenseColumns(mtxY, mtxA, mtxX) ;
93 break ;
94 case SUBMTX_DENSE_ROWS :
95 complex_updDenseRows(mtxY, mtxA, mtxX) ;
96 break ;
97 case SUBMTX_SPARSE_ROWS :
98 complex_updSparseRows(mtxY, mtxA, mtxX) ;
99 break ;
100 case SUBMTX_SPARSE_COLUMNS :
101 complex_updSparseColumns(mtxY, mtxA, mtxX) ;
102 break ;
103 default :
104 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
105 "\n unsupported mode %d for A\n",
106 mtxY, mtxA, mtxX, mtxA->mode) ;
107 exit(-1) ;
108 break ;
109 }
110 break ;
111 default :
112 fprintf(stderr, "\n fatal error in SubMtx_solveupdT(%p,%p,%p)"
113 "\n unsupported type %d for A\n",
114 mtxY, mtxA, mtxX, mtxA->type) ;
115 exit(-1) ;
116 break ;
117 }
118 return ; }
119
120 /*--------------------------------------------------------------------*/
121 /*
122 ----------------
123 A has dense rows
124 ----------------
125 */
126 static void
real_updDenseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)127 real_updDenseRows (
128 SubMtx *mtxY,
129 SubMtx *mtxA,
130 SubMtx *mtxX
131 ) {
132 double Ak0, Ak1, Ak2, x00, x01, x02, x10, x11, x12,
133 x20, x21, x22 ;
134 double *colAT0, *colAT1, *colAT2, *colX0, *colX1, *colX2,
135 *colY0, *colY1, *colY2, *entA, *entX, *entY ;
136 int icolAT, inc1, inc2, irowX, jcolX, krowAT, krowY,
137 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY ;
138 int *colindAT, *rowindAT ;
139
140 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
141 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
142 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
143 colX0 = entX ;
144 colY0 = entY ;
145 if ( ncolAT != nrowX ) {
146 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
147 } else {
148 colindAT = NULL ;
149 }
150 if ( nrowAT != nrowY ) {
151 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
152 } else {
153 rowindAT = NULL ;
154 }
155 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
156 colX1 = colX0 + nrowX ;
157 colX2 = colX1 + nrowX ;
158 colY1 = colY0 + nrowY ;
159 colY2 = colY1 + nrowY ;
160 colAT0 = entA ;
161 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
162 colAT1 = colAT0 + nrowAT ;
163 colAT2 = colAT1 + nrowAT ;
164 if ( ncolAT == nrowX ) {
165 x00 = colX0[icolAT] ;
166 x01 = colX1[icolAT] ;
167 x02 = colX2[icolAT] ;
168 x10 = colX0[icolAT+1] ;
169 x11 = colX1[icolAT+1] ;
170 x12 = colX2[icolAT+1] ;
171 x20 = colX0[icolAT+2] ;
172 x21 = colX1[icolAT+2] ;
173 x22 = colX2[icolAT+2] ;
174 } else {
175 irowX = colindAT[icolAT] ;
176 x00 = colX0[irowX] ;
177 x01 = colX1[irowX] ;
178 x02 = colX2[irowX] ;
179 irowX = colindAT[icolAT+1] ;
180 x10 = colX0[irowX] ;
181 x11 = colX1[irowX] ;
182 x12 = colX2[irowX] ;
183 irowX = colindAT[icolAT+2] ;
184 x20 = colX0[irowX] ;
185 x21 = colX1[irowX] ;
186 x22 = colX2[irowX] ;
187 }
188 if ( nrowY == nrowAT ) {
189 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
190 Ak0 = colAT0[krowAT] ;
191 Ak1 = colAT1[krowAT] ;
192 Ak2 = colAT2[krowAT] ;
193 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
194 colY1[krowAT] -= Ak0 * x01 + Ak1 * x11 + Ak2 * x21 ;
195 colY2[krowAT] -= Ak0 * x02 + Ak1 * x12 + Ak2 * x22 ;
196 }
197 } else {
198 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
199 Ak0 = colAT0[krowAT] ;
200 Ak1 = colAT1[krowAT] ;
201 Ak2 = colAT2[krowAT] ;
202 krowY = rowindAT[krowAT] ;
203 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
204 colY1[krowY] -= Ak0 * x01 + Ak1 * x11 + Ak2 * x21 ;
205 colY2[krowY] -= Ak0 * x02 + Ak1 * x12 + Ak2 * x22 ;
206 }
207 }
208 colAT0 = colAT2 + nrowAT ;
209 }
210 if ( icolAT == ncolAT - 2 ) {
211 colAT1 = colAT0 + nrowAT ;
212 if ( ncolAT == nrowX ) {
213 x00 = colX0[icolAT] ;
214 x01 = colX1[icolAT] ;
215 x02 = colX2[icolAT] ;
216 x10 = colX0[icolAT+1] ;
217 x11 = colX1[icolAT+1] ;
218 x12 = colX2[icolAT+1] ;
219 } else {
220 irowX = colindAT[icolAT] ;
221 x00 = colX0[irowX] ;
222 x01 = colX1[irowX] ;
223 x02 = colX2[irowX] ;
224 irowX = colindAT[icolAT+1] ;
225 x10 = colX0[irowX] ;
226 x11 = colX1[irowX] ;
227 x12 = colX2[irowX] ;
228 }
229 if ( nrowY == nrowAT ) {
230 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
231 Ak0 = colAT0[krowAT] ;
232 Ak1 = colAT1[krowAT] ;
233 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 ;
234 colY1[krowAT] -= Ak0 * x01 + Ak1 * x11 ;
235 colY2[krowAT] -= Ak0 * x02 + Ak1 * x12 ;
236 }
237 } else {
238 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
239 Ak0 = colAT0[krowAT] ;
240 Ak1 = colAT1[krowAT] ;
241 krowY = rowindAT[krowAT] ;
242 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 ;
243 colY1[krowY] -= Ak0 * x01 + Ak1 * x11 ;
244 colY2[krowY] -= Ak0 * x02 + Ak1 * x12 ;
245 }
246 }
247 } else if ( icolAT == ncolAT - 1 ) {
248 if ( ncolAT == nrowX ) {
249 x00 = colX0[icolAT] ;
250 x01 = colX1[icolAT] ;
251 x02 = colX2[icolAT] ;
252 } else {
253 irowX = colindAT[icolAT] ;
254 x00 = colX0[irowX] ;
255 x01 = colX1[irowX] ;
256 x02 = colX2[irowX] ;
257 }
258 if ( nrowY == nrowAT ) {
259 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
260 Ak0 = colAT0[krowAT] ;
261 colY0[krowAT] -= Ak0 * x00 ;
262 colY1[krowAT] -= Ak0 * x01 ;
263 colY2[krowAT] -= Ak0 * x02 ;
264 }
265 } else {
266 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
267 Ak0 = colAT0[krowAT] ;
268 krowY = rowindAT[krowAT] ;
269 colY0[krowY] -= Ak0 * x00 ;
270 colY1[krowY] -= Ak0 * x01 ;
271 colY2[krowY] -= Ak0 * x02 ;
272 }
273 }
274 }
275 colX0 = colX2 + nrowX ;
276 colY0 = colY2 + nrowY ;
277 }
278 if ( jcolX == ncolX - 2 ) {
279 colX1 = colX0 + nrowX ;
280 colY1 = colY0 + nrowY ;
281 colAT0 = entA ;
282 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
283 colAT1 = colAT0 + nrowAT ;
284 colAT2 = colAT1 + nrowAT ;
285 if ( ncolAT == nrowX ) {
286 x00 = colX0[icolAT] ;
287 x01 = colX1[icolAT] ;
288 x10 = colX0[icolAT+1] ;
289 x11 = colX1[icolAT+1] ;
290 x20 = colX0[icolAT+2] ;
291 x21 = colX1[icolAT+2] ;
292 } else {
293 irowX = colindAT[icolAT] ;
294 x00 = colX0[irowX] ;
295 x01 = colX1[irowX] ;
296 irowX = colindAT[icolAT+1] ;
297 x10 = colX0[irowX] ;
298 x11 = colX1[irowX] ;
299 irowX = colindAT[icolAT+2] ;
300 x20 = colX0[irowX] ;
301 x21 = colX1[irowX] ;
302 }
303 if ( nrowY == nrowAT ) {
304 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
305 Ak0 = colAT0[krowAT] ;
306 Ak1 = colAT1[krowAT] ;
307 Ak2 = colAT2[krowAT] ;
308 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
309 colY1[krowAT] -= Ak0 * x01 + Ak1 * x11 + Ak2 * x21 ;
310 }
311 } else {
312 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
313 Ak0 = colAT0[krowAT] ;
314 Ak1 = colAT1[krowAT] ;
315 Ak2 = colAT2[krowAT] ;
316 krowY = rowindAT[krowAT] ;
317 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
318 colY1[krowY] -= Ak0 * x01 + Ak1 * x11 + Ak2 * x21 ;
319 }
320 }
321 colAT0 = colAT2 + nrowAT ;
322 }
323 if ( icolAT == ncolAT - 2 ) {
324 colAT1 = colAT0 + nrowAT ;
325 if ( ncolAT == nrowX ) {
326 x00 = colX0[icolAT] ;
327 x01 = colX1[icolAT] ;
328 x10 = colX0[icolAT+1] ;
329 x11 = colX1[icolAT+1] ;
330 } else {
331 irowX = colindAT[icolAT] ;
332 x00 = colX0[irowX] ;
333 x01 = colX1[irowX] ;
334 irowX = colindAT[icolAT+1] ;
335 x10 = colX0[irowX] ;
336 x11 = colX1[irowX] ;
337 }
338 if ( nrowY == nrowAT ) {
339 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
340 Ak0 = colAT0[krowAT] ;
341 Ak1 = colAT1[krowAT] ;
342 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 ;
343 colY1[krowAT] -= Ak0 * x01 + Ak1 * x11 ;
344 }
345 } else {
346 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
347 Ak0 = colAT0[krowAT] ;
348 Ak1 = colAT1[krowAT] ;
349 krowY = rowindAT[krowAT] ;
350 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 ;
351 colY1[krowY] -= Ak0 * x01 + Ak1 * x11 ;
352 }
353 }
354 } else if ( icolAT == ncolAT - 1 ) {
355 if ( ncolAT == nrowX ) {
356 x00 = colX0[icolAT] ;
357 x01 = colX1[icolAT] ;
358 } else {
359 irowX = colindAT[icolAT] ;
360 x00 = colX0[irowX] ;
361 x01 = colX1[irowX] ;
362 }
363 if ( nrowY == nrowAT ) {
364 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
365 Ak0 = colAT0[krowAT] ;
366 colY0[krowAT] -= Ak0 * x00 ;
367 colY1[krowAT] -= Ak0 * x01 ;
368 }
369 } else {
370 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
371 Ak0 = colAT0[krowAT] ;
372 krowY = rowindAT[krowAT] ;
373 colY0[krowY] -= Ak0 * x00 ;
374 colY1[krowY] -= Ak0 * x01 ;
375 }
376 }
377 }
378 } else if ( jcolX == ncolX - 1 ) {
379 colAT0 = entA ;
380 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
381 colAT1 = colAT0 + nrowAT ;
382 colAT2 = colAT1 + nrowAT ;
383 if ( ncolAT == nrowX ) {
384 x00 = colX0[icolAT] ;
385 x10 = colX0[icolAT+1] ;
386 x20 = colX0[icolAT+2] ;
387 } else {
388 irowX = colindAT[icolAT] ;
389 x00 = colX0[irowX] ;
390 irowX = colindAT[icolAT+1] ;
391 x10 = colX0[irowX] ;
392 irowX = colindAT[icolAT+2] ;
393 x20 = colX0[irowX] ;
394 }
395 if ( nrowY == nrowAT ) {
396 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
397 Ak0 = colAT0[krowAT] ;
398 Ak1 = colAT1[krowAT] ;
399 Ak2 = colAT2[krowAT] ;
400 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
401 }
402 } else {
403 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
404 Ak0 = colAT0[krowAT] ;
405 Ak1 = colAT1[krowAT] ;
406 Ak2 = colAT2[krowAT] ;
407 krowY = rowindAT[krowAT] ;
408 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 + Ak2 * x20 ;
409 }
410 }
411 colAT0 = colAT2 + nrowAT ;
412 }
413 if ( icolAT == ncolAT - 2 ) {
414 colAT1 = colAT0 + nrowAT ;
415 if ( ncolAT == nrowX ) {
416 x00 = colX0[icolAT] ;
417 x10 = colX0[icolAT+1] ;
418 } else {
419 irowX = colindAT[icolAT] ;
420 x00 = colX0[irowX] ;
421 irowX = colindAT[icolAT+1] ;
422 x10 = colX0[irowX] ;
423 }
424 if ( nrowY == nrowAT ) {
425 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
426 Ak0 = colAT0[krowAT] ;
427 Ak1 = colAT1[krowAT] ;
428 colY0[krowAT] -= Ak0 * x00 + Ak1 * x10 ;
429 }
430 } else {
431 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
432 Ak0 = colAT0[krowAT] ;
433 Ak1 = colAT1[krowAT] ;
434 krowY = rowindAT[krowAT] ;
435 colY0[krowY] -= Ak0 * x00 + Ak1 * x10 ;
436 }
437 }
438 } else if ( icolAT == ncolAT - 1 ) {
439 if ( ncolAT == nrowX ) {
440 x00 = colX0[icolAT] ;
441 } else {
442 irowX = colindAT[icolAT] ;
443 x00 = colX0[irowX] ;
444 }
445 if ( nrowY == nrowAT ) {
446 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
447 Ak0 = colAT0[krowAT] ;
448 colY0[krowAT] -= Ak0 * x00 ;
449 }
450 } else {
451 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
452 Ak0 = colAT0[krowAT] ;
453 krowY = rowindAT[krowAT] ;
454 colY0[krowY] -= Ak0 * x00 ;
455 }
456 }
457 }
458 }
459 return ; }
460
461 /*--------------------------------------------------------------------*/
462 /*
463 ----------------
464 A has dense rows
465 ----------------
466 */
467 static void
real_updDenseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)468 real_updDenseColumns (
469 SubMtx *mtxY,
470 SubMtx *mtxA,
471 SubMtx *mtxX
472 ) {
473 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
474 *rowAT0, *rowAT1, *rowAT2, *entA, *entX, *entY ;
475 int inc1, inc2, irowAT, irowY, jcolX, kcolAT, krowX,
476 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY ;
477 int *colindAT, *rowindAT ;
478
479 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
480 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
481 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
482 if ( ncolAT != nrowX ) {
483 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
484 } else {
485 colindAT = NULL ;
486 }
487 if ( nrowAT != nrowY ) {
488 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
489 } else {
490 rowindAT = NULL ;
491 }
492 colX0 = entX ;
493 colY0 = entY ;
494 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
495 colX1 = colX0 + nrowX ;
496 colX2 = colX1 + nrowX ;
497 colY1 = colY0 + nrowY ;
498 colY2 = colY1 + nrowY ;
499 rowAT0 = entA ;
500 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
501 double A0k, A1k, A2k, Xk0, Xk1, Xk2 ;
502 double sum00, sum01, sum02,
503 sum10, sum11, sum12,
504 sum20, sum21, sum22 ;
505
506 sum00 = sum01 = sum02
507 = sum10 = sum11 = sum12 = sum20 = sum21 = sum22 = 0.0 ;
508 rowAT1 = rowAT0 + ncolAT ;
509 rowAT2 = rowAT1 + ncolAT ;
510 if ( ncolAT == nrowX ) {
511 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
512 A0k = rowAT0[kcolAT] ;
513 A1k = rowAT1[kcolAT] ;
514 A2k = rowAT2[kcolAT] ;
515 Xk0 = colX0[kcolAT] ;
516 Xk1 = colX1[kcolAT] ;
517 Xk2 = colX2[kcolAT] ;
518 sum00 += A0k * Xk0 ;
519 sum01 += A0k * Xk1 ;
520 sum02 += A0k * Xk2 ;
521 sum10 += A1k * Xk0 ;
522 sum11 += A1k * Xk1 ;
523 sum12 += A1k * Xk2 ;
524 sum20 += A2k * Xk0 ;
525 sum21 += A2k * Xk1 ;
526 sum22 += A2k * Xk2 ;
527 }
528 } else {
529 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
530 A0k = rowAT0[kcolAT] ;
531 A1k = rowAT1[kcolAT] ;
532 A2k = rowAT2[kcolAT] ;
533 krowX = colindAT[kcolAT] ;
534 Xk0 = colX0[krowX] ;
535 Xk1 = colX1[krowX] ;
536 Xk2 = colX2[krowX] ;
537 sum00 += A0k * Xk0 ;
538 sum01 += A0k * Xk1 ;
539 sum02 += A0k * Xk2 ;
540 sum10 += A1k * Xk0 ;
541 sum11 += A1k * Xk1 ;
542 sum12 += A1k * Xk2 ;
543 sum20 += A2k * Xk0 ;
544 sum21 += A2k * Xk1 ;
545 sum22 += A2k * Xk2 ;
546 }
547 }
548 if ( nrowY == nrowAT ) {
549 colY0[irowAT] -= sum00 ;
550 colY1[irowAT] -= sum01 ;
551 colY2[irowAT] -= sum02 ;
552 colY0[irowAT+1] -= sum10 ;
553 colY1[irowAT+1] -= sum11 ;
554 colY2[irowAT+1] -= sum12 ;
555 colY0[irowAT+2] -= sum20 ;
556 colY1[irowAT+2] -= sum21 ;
557 colY2[irowAT+2] -= sum22 ;
558 } else {
559 irowY = rowindAT[irowAT] ;
560 colY0[irowY] -= sum00 ;
561 colY1[irowY] -= sum01 ;
562 colY2[irowY] -= sum02 ;
563 irowY = rowindAT[irowAT+1] ;
564 colY0[irowY] -= sum10 ;
565 colY1[irowY] -= sum11 ;
566 colY2[irowY] -= sum12 ;
567 irowY = rowindAT[irowAT+2] ;
568 colY0[irowY] -= sum20 ;
569 colY1[irowY] -= sum21 ;
570 colY2[irowY] -= sum22 ;
571 }
572 rowAT0 = rowAT2 + ncolAT ;
573 }
574 if ( irowAT == nrowAT - 2 ) {
575 double A0k, A1k, Xk0, Xk1, Xk2 ;
576 double sum00, sum01, sum02, sum10, sum11, sum12 ;
577
578 sum00 = sum01 = sum02 = sum10 = sum11 = sum12 = 0.0 ;
579 rowAT1 = rowAT0 + ncolAT ;
580 if ( ncolAT == nrowX ) {
581 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
582 A0k = rowAT0[kcolAT] ;
583 A1k = rowAT1[kcolAT] ;
584 Xk0 = colX0[kcolAT] ;
585 Xk1 = colX1[kcolAT] ;
586 Xk2 = colX2[kcolAT] ;
587 sum00 += A0k * Xk0 ;
588 sum01 += A0k * Xk1 ;
589 sum02 += A0k * Xk2 ;
590 sum10 += A1k * Xk0 ;
591 sum11 += A1k * Xk1 ;
592 sum12 += A1k * Xk2 ;
593 }
594 } else {
595 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
596 A0k = rowAT0[kcolAT] ;
597 A1k = rowAT1[kcolAT] ;
598 krowX = colindAT[kcolAT] ;
599 Xk0 = colX0[krowX] ;
600 Xk1 = colX1[krowX] ;
601 Xk2 = colX2[krowX] ;
602 sum00 += A0k * Xk0 ;
603 sum01 += A0k * Xk1 ;
604 sum02 += A0k * Xk2 ;
605 sum10 += A1k * Xk0 ;
606 sum11 += A1k * Xk1 ;
607 sum12 += A1k * Xk2 ;
608 }
609 }
610 if ( nrowY == nrowAT ) {
611 colY0[irowAT] -= sum00 ;
612 colY1[irowAT] -= sum01 ;
613 colY2[irowAT] -= sum02 ;
614 colY0[irowAT+1] -= sum10 ;
615 colY1[irowAT+1] -= sum11 ;
616 colY2[irowAT+1] -= sum12 ;
617 } else {
618 irowY = rowindAT[irowAT] ;
619 colY0[irowY] -= sum00 ;
620 colY1[irowY] -= sum01 ;
621 colY2[irowY] -= sum02 ;
622 irowY = rowindAT[irowAT+1] ;
623 colY0[irowY] -= sum10 ;
624 colY1[irowY] -= sum11 ;
625 colY2[irowY] -= sum12 ;
626 }
627 } else if ( irowAT == nrowAT - 1 ) {
628 double A0k, Xk0, Xk1, Xk2 ;
629 double sum00, sum01, sum02 ;
630
631 sum00 = sum01 = sum02 = 0.0 ;
632 if ( ncolAT == nrowX ) {
633 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
634 A0k = rowAT0[kcolAT] ;
635 Xk0 = colX0[kcolAT] ;
636 Xk1 = colX1[kcolAT] ;
637 Xk2 = colX2[kcolAT] ;
638 sum00 += A0k * Xk0 ;
639 sum01 += A0k * Xk1 ;
640 sum02 += A0k * Xk2 ;
641 }
642 } else {
643 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
644 A0k = rowAT0[kcolAT] ;
645 krowX = colindAT[kcolAT] ;
646 Xk0 = colX0[krowX] ;
647 Xk1 = colX1[krowX] ;
648 Xk2 = colX2[krowX] ;
649 sum00 += A0k * Xk0 ;
650 sum01 += A0k * Xk1 ;
651 sum02 += A0k * Xk2 ;
652 }
653 }
654 if ( nrowY == nrowAT ) {
655 colY0[irowAT] -= sum00 ;
656 colY1[irowAT] -= sum01 ;
657 colY2[irowAT] -= sum02 ;
658 } else {
659 irowY = rowindAT[irowAT] ;
660 colY0[irowY] -= sum00 ;
661 colY1[irowY] -= sum01 ;
662 colY2[irowY] -= sum02 ;
663 }
664 }
665 colX0 = colX2 + nrowX ;
666 colY0 = colY2 + nrowY ;
667 }
668 if ( jcolX == ncolX - 2 ) {
669 colX1 = colX0 + nrowX ;
670 colY1 = colY0 + nrowY ;
671 rowAT0 = entA ;
672 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
673 double A0k, A1k, A2k, Xk0, Xk1 ;
674 double sum00, sum01, sum10, sum11, sum20, sum21 ;
675
676 sum00 = sum01 = sum10 = sum11 = sum20 = sum21 = 0.0 ;
677 rowAT1 = rowAT0 + ncolAT ;
678 rowAT2 = rowAT1 + ncolAT ;
679 if ( ncolAT == nrowX ) {
680 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
681 A0k = rowAT0[kcolAT] ;
682 A1k = rowAT1[kcolAT] ;
683 A2k = rowAT2[kcolAT] ;
684 Xk0 = colX0[kcolAT] ;
685 Xk1 = colX1[kcolAT] ;
686 sum00 += A0k * Xk0 ;
687 sum01 += A0k * Xk1 ;
688 sum10 += A1k * Xk0 ;
689 sum11 += A1k * Xk1 ;
690 sum20 += A2k * Xk0 ;
691 sum21 += A2k * Xk1 ;
692 }
693 } else {
694 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
695 A0k = rowAT0[kcolAT] ;
696 A1k = rowAT1[kcolAT] ;
697 A2k = rowAT2[kcolAT] ;
698 krowX = colindAT[kcolAT] ;
699 Xk0 = colX0[krowX] ;
700 Xk1 = colX1[krowX] ;
701 sum00 += A0k * Xk0 ;
702 sum01 += A0k * Xk1 ;
703 sum10 += A1k * Xk0 ;
704 sum11 += A1k * Xk1 ;
705 sum20 += A2k * Xk0 ;
706 sum21 += A2k * Xk1 ;
707 }
708 }
709 if ( nrowY == nrowAT ) {
710 colY0[irowAT] -= sum00 ;
711 colY1[irowAT] -= sum01 ;
712 colY0[irowAT+1] -= sum10 ;
713 colY1[irowAT+1] -= sum11 ;
714 colY0[irowAT+2] -= sum20 ;
715 colY1[irowAT+2] -= sum21 ;
716 } else {
717 irowY = rowindAT[irowAT] ;
718 colY0[irowY] -= sum00 ;
719 colY1[irowY] -= sum01 ;
720 irowY = rowindAT[irowAT+1] ;
721 colY0[irowY] -= sum10 ;
722 colY1[irowY] -= sum11 ;
723 irowY = rowindAT[irowAT+2] ;
724 colY0[irowY] -= sum20 ;
725 colY1[irowY] -= sum21 ;
726 }
727 rowAT0 = rowAT2 + ncolAT ;
728 }
729 if ( irowAT == nrowAT - 2 ) {
730 double A0k, A1k, Xk0, Xk1 ;
731 double sum00, sum01, sum10, sum11 ;
732
733 sum00 = sum01 = sum10 = sum11 = 0.0 ;
734 rowAT1 = rowAT0 + ncolAT ;
735 if ( ncolAT == nrowX ) {
736 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
737 A0k = rowAT0[kcolAT] ;
738 A1k = rowAT1[kcolAT] ;
739 Xk0 = colX0[kcolAT] ;
740 Xk1 = colX1[kcolAT] ;
741 sum00 += A0k * Xk0 ;
742 sum01 += A0k * Xk1 ;
743 sum10 += A1k * Xk0 ;
744 sum11 += A1k * Xk1 ;
745 }
746 } else {
747 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
748 A0k = rowAT0[kcolAT] ;
749 A1k = rowAT1[kcolAT] ;
750 krowX = colindAT[kcolAT] ;
751 Xk0 = colX0[krowX] ;
752 Xk1 = colX1[krowX] ;
753 sum00 += A0k * Xk0 ;
754 sum01 += A0k * Xk1 ;
755 sum10 += A1k * Xk0 ;
756 sum11 += A1k * Xk1 ;
757 }
758 }
759 if ( nrowY == nrowAT ) {
760 colY0[irowAT] -= sum00 ;
761 colY1[irowAT] -= sum01 ;
762 colY0[irowAT+1] -= sum10 ;
763 colY1[irowAT+1] -= sum11 ;
764 } else {
765 irowY = rowindAT[irowAT] ;
766 colY0[irowY] -= sum00 ;
767 colY1[irowY] -= sum01 ;
768 irowY = rowindAT[irowAT+1] ;
769 colY0[irowY] -= sum10 ;
770 colY1[irowY] -= sum11 ;
771 }
772 } else if ( irowAT == nrowAT - 1 ) {
773 double A0k, Xk0, Xk1 ;
774 double sum00, sum01 ;
775
776 sum00 = sum01 = 0.0 ;
777 if ( ncolAT == nrowX ) {
778 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
779 A0k = rowAT0[kcolAT] ;
780 Xk0 = colX0[kcolAT] ;
781 Xk1 = colX1[kcolAT] ;
782 sum00 += A0k * Xk0 ;
783 sum01 += A0k * Xk1 ;
784 }
785 } else {
786 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
787 A0k = rowAT0[kcolAT] ;
788 krowX = colindAT[kcolAT] ;
789 Xk0 = colX0[krowX] ;
790 Xk1 = colX1[krowX] ;
791 sum00 += A0k * Xk0 ;
792 sum01 += A0k * Xk1 ;
793 }
794 }
795 if ( nrowY == nrowAT ) {
796 colY0[irowAT] -= sum00 ;
797 colY1[irowAT] -= sum01 ;
798 } else {
799 irowY = rowindAT[irowAT] ;
800 colY0[irowY] -= sum00 ;
801 colY1[irowY] -= sum01 ;
802 }
803 }
804 } else if ( jcolX == ncolX - 1 ) {
805 rowAT0 = entA ;
806 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
807 double A0k, A1k, A2k, Xk0 ;
808 double sum00, sum10, sum20 ;
809
810 sum00 = sum10 = sum20 = 0.0 ;
811 rowAT1 = rowAT0 + ncolAT ;
812 rowAT2 = rowAT1 + ncolAT ;
813 if ( ncolAT == nrowX ) {
814 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
815 A0k = rowAT0[kcolAT] ;
816 A1k = rowAT1[kcolAT] ;
817 A2k = rowAT2[kcolAT] ;
818 Xk0 = colX0[kcolAT] ;
819 sum00 += A0k * Xk0 ;
820 sum10 += A1k * Xk0 ;
821 sum20 += A2k * Xk0 ;
822 }
823 } else {
824 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
825 A0k = rowAT0[kcolAT] ;
826 A1k = rowAT1[kcolAT] ;
827 A2k = rowAT2[kcolAT] ;
828 krowX = colindAT[kcolAT] ;
829 Xk0 = colX0[krowX] ;
830 sum00 += A0k * Xk0 ;
831 sum10 += A1k * Xk0 ;
832 sum20 += A2k * Xk0 ;
833 }
834 }
835 if ( nrowY == nrowAT ) {
836 colY0[irowAT] -= sum00 ;
837 colY0[irowAT+1] -= sum10 ;
838 colY0[irowAT+2] -= sum20 ;
839 } else {
840 irowY = rowindAT[irowAT] ;
841 colY0[irowY] -= sum00 ;
842 irowY = rowindAT[irowAT+1] ;
843 colY0[irowY] -= sum10 ;
844 irowY = rowindAT[irowAT+2] ;
845 colY0[irowY] -= sum20 ;
846 }
847 rowAT0 = rowAT2 + ncolAT ;
848 }
849 if ( irowAT == nrowAT - 2 ) {
850 double A0k, A1k, Xk0 ;
851 double sum00, sum10 ;
852
853 sum00 = sum10 = 0.0 ;
854 rowAT1 = rowAT0 + ncolAT ;
855 if ( ncolAT == nrowX ) {
856 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
857 A0k = rowAT0[kcolAT] ;
858 A1k = rowAT1[kcolAT] ;
859 Xk0 = colX0[kcolAT] ;
860 sum00 += A0k * Xk0 ;
861 sum10 += A1k * Xk0 ;
862 }
863 } else {
864 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
865 A0k = rowAT0[kcolAT] ;
866 A1k = rowAT1[kcolAT] ;
867 krowX = colindAT[kcolAT] ;
868 Xk0 = colX0[krowX] ;
869 sum00 += A0k * Xk0 ;
870 sum10 += A1k * Xk0 ;
871 }
872 }
873 if ( nrowY == nrowAT ) {
874 colY0[irowAT] -= sum00 ;
875 colY0[irowAT+1] -= sum10 ;
876 } else {
877 irowY = rowindAT[irowAT] ;
878 colY0[irowY] -= sum00 ;
879 irowY = rowindAT[irowAT+1] ;
880 colY0[irowY] -= sum10 ;
881 }
882 } else if ( irowAT == nrowAT - 1 ) {
883 double A0k, Xk0 ;
884 double sum00 ;
885
886 sum00 = 0.0 ;
887 if ( ncolAT == nrowX ) {
888 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
889 A0k = rowAT0[kcolAT] ;
890 Xk0 = colX0[kcolAT] ;
891 sum00 += A0k * Xk0 ;
892 }
893 } else {
894 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
895 A0k = rowAT0[kcolAT] ;
896 krowX = colindAT[kcolAT] ;
897 Xk0 = colX0[krowX] ;
898 sum00 += A0k * Xk0 ;
899 }
900 }
901 if ( nrowY == nrowAT ) {
902 colY0[irowAT] -= sum00 ;
903 } else {
904 irowY = rowindAT[irowAT] ;
905 colY0[irowY] -= sum00 ;
906 }
907 }
908 }
909 return ; }
910
911 /*--------------------------------------------------------------------*/
912 /*
913 --------------------
914 A has sparse columns
915 --------------------
916 */
917 static void
real_updSparseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)918 real_updSparseColumns (
919 SubMtx *mtxY,
920 SubMtx *mtxA,
921 SubMtx *mtxX
922 ) {
923 double Aik, sum0, sum1, sum2 ;
924 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
925 *entA, *entX, *entY ;
926 int ii, inc1, inc2, irowAT, irowY, jcolX, kk, krowX,
927 ncolAT, ncolX, ncolY, nentA, nrowAT, nrowX, nrowY, size ;
928 int *colindAT, *indices, *rowindAT, *sizes ;
929 /*
930 fprintf(stdout, "\n UPDATE_SPARSE_ROWS(%d,%d)",
931 mtxA->rowid, mtxA->colid) ;
932 */
933 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
934 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
935 SubMtx_sparseColumnsInfo(mtxA, &nrowAT, &nentA, &sizes, &indices, &entA) ;
936 if ( (ncolAT = mtxA->nrow) != nrowX ) {
937 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
938 } else {
939 colindAT = NULL ;
940 }
941 if ( (nrowAT = mtxA->ncol) != nrowY ) {
942 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
943 } else {
944 rowindAT = NULL ;
945 }
946 colX0 = entX ;
947 colY0 = entY ;
948 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
949 colX1 = colX0 + nrowX ;
950 colX2 = colX1 + nrowX ;
951 colY1 = colY0 + nrowY ;
952 colY2 = colY1 + nrowY ;
953 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
954 if ( (size = sizes[irowAT]) > 0 ) {
955 sum0 = sum1 = sum2 = 0.0 ;
956 if ( ncolAT == nrowX ) {
957 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
958 Aik = entA[kk] ;
959 krowX = indices[kk] ;
960 sum0 += Aik * colX0[krowX] ;
961 sum1 += Aik * colX1[krowX] ;
962 sum2 += Aik * colX2[krowX] ;
963 }
964 } else {
965 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
966 Aik = entA[kk] ;
967 krowX = colindAT[indices[kk]] ;
968 sum0 += Aik * colX0[krowX] ;
969 sum1 += Aik * colX1[krowX] ;
970 sum2 += Aik * colX2[krowX] ;
971 }
972 }
973 if ( nrowAT == nrowY ) {
974 colY0[irowAT] -= sum0 ;
975 colY1[irowAT] -= sum1 ;
976 colY2[irowAT] -= sum2 ;
977 } else {
978 irowY = rowindAT[irowAT] ;
979 colY0[irowY] -= sum0 ;
980 colY1[irowY] -= sum1 ;
981 colY2[irowY] -= sum2 ;
982 }
983 }
984 }
985 colX0 = colX2 + nrowX ;
986 colY0 = colY2 + nrowY ;
987 }
988 if ( jcolX == ncolX - 2 ) {
989 colX1 = colX0 + nrowX ;
990 colY1 = colY0 + nrowY ;
991 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
992 if ( (size = sizes[irowAT]) > 0 ) {
993 sum0 = sum1 = 0.0 ;
994 if ( ncolAT == nrowX ) {
995 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
996 Aik = entA[kk] ;
997 krowX = indices[kk] ;
998 sum0 += Aik * colX0[krowX] ;
999 sum1 += Aik * colX1[krowX] ;
1000 }
1001 } else {
1002 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1003 Aik = entA[kk] ;
1004 krowX = colindAT[indices[kk]] ;
1005 sum0 += Aik * colX0[krowX] ;
1006 sum1 += Aik * colX1[krowX] ;
1007 }
1008 }
1009 if ( nrowAT == nrowY ) {
1010 colY0[irowAT] -= sum0 ;
1011 colY1[irowAT] -= sum1 ;
1012 } else {
1013 irowY = rowindAT[irowAT] ;
1014 colY0[irowY] -= sum0 ;
1015 colY1[irowY] -= sum1 ;
1016 }
1017 }
1018 }
1019 } else if ( jcolX == ncolX - 1 ) {
1020 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
1021 if ( (size = sizes[irowAT]) > 0 ) {
1022 sum0 = 0.0 ;
1023 if ( ncolAT == nrowX ) {
1024 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1025 Aik = entA[kk] ;
1026 krowX = indices[kk] ;
1027 sum0 += Aik * colX0[krowX] ;
1028 }
1029 } else {
1030 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1031 Aik = entA[kk] ;
1032 krowX = colindAT[indices[kk]] ;
1033 sum0 += Aik * colX0[krowX] ;
1034 }
1035 }
1036 if ( nrowAT == nrowY ) {
1037 colY0[irowAT] -= sum0 ;
1038 } else {
1039 irowY = rowindAT[irowAT] ;
1040 colY0[irowY] -= sum0 ;
1041 }
1042 }
1043 }
1044 }
1045 return ; }
1046
1047 /*--------------------------------------------------------------------*/
1048 /*
1049 -----------------
1050 A has sparse rows
1051 -----------------
1052 */
1053 static void
real_updSparseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)1054 real_updSparseRows (
1055 SubMtx *mtxY,
1056 SubMtx *mtxA,
1057 SubMtx *mtxX
1058 ) {
1059 double Aij, Xj0, Xj1, Xj2 ;
1060 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
1061 *entA, *entX, *entY ;
1062 int ii, inc1, inc2, irowY, jcolAT, jcolX, jrowX, kk,
1063 ncolAT, ncolX, ncolY, nentA, nrowAT, nrowX, nrowY, size ;
1064 int *colindAT, *indices, *rowindAT, *sizes ;
1065 /*
1066 fprintf(stdout, "\n UPDATE_SPARSE_COLUMNS(%d,%d)",
1067 mtxA->rowid, mtxA->colid) ;
1068 */
1069 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
1070 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
1071 SubMtx_sparseRowsInfo(mtxA, &ncolAT, &nentA, &sizes, &indices, &entA) ;
1072 if ( (ncolAT = mtxA->nrow) != nrowX ) {
1073 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
1074 } else {
1075 colindAT = NULL ;
1076 }
1077 if ( (nrowAT = mtxA->ncol) != nrowY ) {
1078 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
1079 } else {
1080 rowindAT = NULL ;
1081 }
1082 colX0 = entX ;
1083 colY0 = entY ;
1084 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
1085 colX1 = colX0 + nrowX ;
1086 colX2 = colX1 + nrowX ;
1087 colY1 = colY0 + nrowY ;
1088 colY2 = colY1 + nrowY ;
1089 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1090 if ( (size = sizes[jcolAT]) > 0 ) {
1091 if ( ncolAT == nrowX ) {
1092 jrowX = jcolAT ;
1093 } else {
1094 jrowX = colindAT[jcolAT] ;
1095 }
1096 Xj0 = colX0[jrowX] ;
1097 Xj1 = colX1[jrowX] ;
1098 Xj2 = colX2[jrowX] ;
1099 if ( nrowAT == nrowY ) {
1100 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1101 Aij = entA[kk] ;
1102 irowY = indices[kk] ;
1103 colY0[irowY] -= Aij * Xj0 ;
1104 colY1[irowY] -= Aij * Xj1 ;
1105 colY2[irowY] -= Aij * Xj2 ;
1106 }
1107 } else {
1108 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1109 Aij = entA[kk] ;
1110 irowY = rowindAT[indices[kk]] ;
1111 colY0[irowY] -= Aij * Xj0 ;
1112 colY1[irowY] -= Aij * Xj1 ;
1113 colY2[irowY] -= Aij * Xj2 ;
1114 }
1115 }
1116 }
1117 }
1118 colX0 = colX2 + nrowX ;
1119 colY0 = colY2 + nrowY ;
1120 }
1121 if ( jcolX == ncolX - 2 ) {
1122 colX1 = colX0 + nrowX ;
1123 colY1 = colY0 + nrowY ;
1124 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1125 if ( (size = sizes[jcolAT]) > 0 ) {
1126 if ( ncolAT == nrowX ) {
1127 jrowX = jcolAT ;
1128 } else {
1129 jrowX = colindAT[jcolAT] ;
1130 }
1131 Xj0 = colX0[jrowX] ;
1132 Xj1 = colX1[jrowX] ;
1133 if ( nrowAT == nrowY ) {
1134 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1135 Aij = entA[kk] ;
1136 irowY = indices[kk] ;
1137 colY0[irowY] -= Aij * Xj0 ;
1138 colY1[irowY] -= Aij * Xj1 ;
1139 }
1140 } else {
1141 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1142 Aij = entA[kk] ;
1143 irowY = rowindAT[indices[kk]] ;
1144 colY0[irowY] -= Aij * Xj0 ;
1145 colY1[irowY] -= Aij * Xj1 ;
1146 }
1147 }
1148 }
1149 }
1150 } else if ( jcolX == ncolX - 1 ) {
1151 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1152 if ( (size = sizes[jcolAT]) > 0 ) {
1153 if ( ncolAT == nrowX ) {
1154 jrowX = jcolAT ;
1155 } else {
1156 jrowX = colindAT[jcolAT] ;
1157 }
1158 Xj0 = colX0[jrowX] ;
1159 if ( nrowAT == nrowY ) {
1160 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1161 Aij = entA[kk] ;
1162 irowY = indices[kk] ;
1163 colY0[irowY] -= Aij * Xj0 ;
1164 }
1165 } else {
1166 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1167 Aij = entA[kk] ;
1168 irowY = rowindAT[indices[kk]] ;
1169 colY0[irowY] -= Aij * Xj0 ;
1170 }
1171 }
1172 }
1173 }
1174 }
1175 return ; }
1176
1177 /*--------------------------------------------------------------------*/
1178 /*
1179 ----------------
1180 A has dense rows
1181 ----------------
1182 */
1183 static void
complex_updDenseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)1184 complex_updDenseRows (
1185 SubMtx *mtxY,
1186 SubMtx *mtxA,
1187 SubMtx *mtxX
1188 ) {
1189 double ai0, ai1, ai2, ar0, ar1, ar2,
1190 xi00, xi01, xi02, xi10, xi11, xi12, xi20, xi21, xi22,
1191 xr00, xr01, xr02, xr10, xr11, xr12, xr20, xr21, xr22 ;
1192 double *colAT0, *colAT1, *colAT2, *colX0, *colX1, *colX2,
1193 *colY0, *colY1, *colY2, *entA, *entX, *entY ;
1194 int icolAT, ialoc, iloc, inc1, inc2, jcolX, krowAT,
1195 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY, raloc, rloc ;
1196 int *colindAT, *rowindAT ;
1197
1198 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
1199 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
1200 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
1201 colX0 = entX ;
1202 colY0 = entY ;
1203 if ( ncolAT != nrowX ) {
1204 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
1205 } else {
1206 colindAT = NULL ;
1207 }
1208 if ( nrowAT != nrowY ) {
1209 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
1210 } else {
1211 rowindAT = NULL ;
1212 }
1213 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
1214 colX1 = colX0 + 2*nrowX ;
1215 colX2 = colX1 + 2*nrowX ;
1216 colY1 = colY0 + 2*nrowY ;
1217 colY2 = colY1 + 2*nrowY ;
1218 colAT0 = entA ;
1219 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
1220 colAT1 = colAT0 + 2*nrowAT ;
1221 colAT2 = colAT1 + 2*nrowAT ;
1222 if ( ncolAT == nrowX ) {
1223 rloc = 2*icolAT ; iloc = rloc + 1 ;
1224 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1225 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1226 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1227 rloc += 2 ; iloc += 2 ;
1228 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1229 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1230 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1231 rloc += 2 ; iloc += 2 ;
1232 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1233 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1234 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1235 } else {
1236 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1237 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1238 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1239 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1240 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1241 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1242 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1243 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1244 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1245 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1246 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1247 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1248 }
1249 if ( nrowY == nrowAT ) {
1250 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1251 rloc = 2*krowAT ; iloc = rloc + 1 ;
1252 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1253 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1254 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
1255 colY0[rloc] -= ar0*xr00 - ai0*xi00
1256 + ar1*xr10 - ai1*xi10
1257 + ar2*xr20 - ai2*xi20 ;
1258 colY0[iloc] -= ar0*xi00 + ai0*xr00
1259 + ar1*xi10 + ai1*xr10
1260 + ar2*xi20 + ai2*xr20 ;
1261 colY1[rloc] -= ar0*xr01 - ai0*xi01
1262 + ar1*xr11 - ai1*xi11
1263 + ar2*xr21 - ai2*xi21 ;
1264 colY1[iloc] -= ar0*xi01 + ai0*xr01
1265 + ar1*xi11 + ai1*xr11
1266 + ar2*xi21 + ai2*xr21 ;
1267 colY2[rloc] -= ar0*xr02 - ai0*xi02
1268 + ar1*xr12 - ai1*xi12
1269 + ar2*xr22 - ai2*xi22 ;
1270 colY2[iloc] -= ar0*xi02 + ai0*xr02
1271 + ar1*xi12 + ai1*xr12
1272 + ar2*xi22 + ai2*xr22 ;
1273 }
1274 } else {
1275 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1276 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1277 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1278 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1279 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
1280 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1281 colY0[rloc] -= ar0*xr00 - ai0*xi00
1282 + ar1*xr10 - ai1*xi10
1283 + ar2*xr20 - ai2*xi20 ;
1284 colY0[iloc] -= ar0*xi00 + ai0*xr00
1285 + ar1*xi10 + ai1*xr10
1286 + ar2*xi20 + ai2*xr20 ;
1287 colY1[rloc] -= ar0*xr01 - ai0*xi01
1288 + ar1*xr11 - ai1*xi11
1289 + ar2*xr21 - ai2*xi21 ;
1290 colY1[iloc] -= ar0*xi01 + ai0*xr01
1291 + ar1*xi11 + ai1*xr11
1292 + ar2*xi21 + ai2*xr21 ;
1293 colY2[rloc] -= ar0*xr02 - ai0*xi02
1294 + ar1*xr12 - ai1*xi12
1295 + ar2*xr22 - ai2*xi22 ;
1296 colY2[iloc] -= ar0*xi02 + ai0*xr02
1297 + ar1*xi12 + ai1*xr12
1298 + ar2*xi22 + ai2*xr22 ;
1299 }
1300 }
1301 colAT0 = colAT2 + 2*nrowAT ;
1302 }
1303 if ( icolAT == ncolAT - 2 ) {
1304 colAT1 = colAT0 + 2*nrowAT ;
1305 if ( ncolAT == nrowX ) {
1306 rloc = 2*icolAT ; iloc = rloc + 1 ;
1307 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1308 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1309 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1310 rloc += 2 ; iloc += 2 ;
1311 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1312 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1313 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1314 rloc += 2 ; iloc += 2 ;
1315 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1316 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1317 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1318 } else {
1319 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1320 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1321 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1322 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1323 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1324 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1325 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1326 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1327 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1328 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1329 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1330 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1331 }
1332 if ( nrowY == nrowAT ) {
1333 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1334 rloc = 2*krowAT ; iloc = rloc + 1 ;
1335 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1336 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1337 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1338 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1339 colY1[rloc] -= ar0*xr01 - ai0*xi01 + ar1*xr11 - ai1*xi11 ;
1340 colY1[iloc] -= ar0*xi01 + ai0*xr01 + ar1*xi11 + ai1*xr11 ;
1341 colY2[rloc] -= ar0*xr02 - ai0*xi02 + ar1*xr12 - ai1*xi12 ;
1342 colY2[iloc] -= ar0*xi02 + ai0*xr02 + ar1*xi12 + ai1*xr12 ;
1343 }
1344 } else {
1345 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1346 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1347 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1348 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1349 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1350 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1351 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1352 colY1[rloc] -= ar0*xr01 - ai0*xi01 + ar1*xr11 - ai1*xi11 ;
1353 colY1[iloc] -= ar0*xi01 + ai0*xr01 + ar1*xi11 + ai1*xr11 ;
1354 colY2[rloc] -= ar0*xr02 - ai0*xi02 + ar1*xr12 - ai1*xi12 ;
1355 colY2[iloc] -= ar0*xi02 + ai0*xr02 + ar1*xi12 + ai1*xr12 ;
1356 }
1357 }
1358 } else if ( icolAT == ncolAT - 1 ) {
1359 if ( ncolAT == nrowX ) {
1360 rloc = 2*icolAT ; iloc = rloc + 1 ;
1361 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1362 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1363 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1364 rloc += 2 ; iloc += 2 ;
1365 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1366 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1367 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1368 rloc += 2 ; iloc += 2 ;
1369 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1370 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1371 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1372 } else {
1373 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1374 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1375 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1376 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
1377 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1378 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1379 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1380 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
1381 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1382 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1383 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1384 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
1385 }
1386 if ( nrowY == nrowAT ) {
1387 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1388 rloc = 2*krowAT ; iloc = rloc + 1 ;
1389 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1390 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1391 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1392 colY1[rloc] -= ar0*xr01 - ai0*xi01 ;
1393 colY1[iloc] -= ar0*xi01 + ai0*xr01 ;
1394 colY2[rloc] -= ar0*xr02 - ai0*xi02 ;
1395 colY2[iloc] -= ar0*xi02 + ai0*xr02 ;
1396 }
1397 } else {
1398 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1399 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1400 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1401 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1402 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1403 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1404 colY1[rloc] -= ar0*xr01 - ai0*xi01 ;
1405 colY1[iloc] -= ar0*xi01 + ai0*xr01 ;
1406 colY2[rloc] -= ar0*xr02 - ai0*xi02 ;
1407 colY2[iloc] -= ar0*xi02 + ai0*xr02 ;
1408 }
1409 }
1410 }
1411 colX0 = colX2 + 2*nrowX ;
1412 colY0 = colY2 + 2*nrowY ;
1413 }
1414 /*
1415 fprintf(stdout, "\n %% jcolX = %d", jcolX) ;
1416 */
1417 if ( jcolX == ncolX - 2 ) {
1418 colX1 = colX0 + 2*nrowX ;
1419 colY1 = colY0 + 2*nrowY ;
1420 colAT0 = entA ;
1421 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
1422 colAT1 = colAT0 + 2*nrowAT ;
1423 colAT2 = colAT1 + 2*nrowAT ;
1424 if ( ncolAT == nrowX ) {
1425 rloc = 2*icolAT ; iloc = rloc + 1 ;
1426 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1427 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1428 rloc += 2 ; iloc += 2 ;
1429 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1430 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1431 rloc += 2 ; iloc += 2 ;
1432 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1433 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1434 } else {
1435 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1436 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1437 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1438 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1439 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1440 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1441 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1442 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1443 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1444 }
1445 if ( nrowY == nrowAT ) {
1446 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1447 rloc = 2*krowAT ; iloc = rloc + 1 ;
1448 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1449 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1450 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
1451 colY0[rloc] -= ar0*xr00 - ai0*xi00
1452 + ar1*xr10 - ai1*xi10
1453 + ar2*xr20 - ai2*xi20 ;
1454 colY0[iloc] -= ar0*xi00 + ai0*xr00
1455 + ar1*xi10 + ai1*xr10
1456 + ar2*xi20 + ai2*xr20 ;
1457 colY1[rloc] -= ar0*xr01 - ai0*xi01
1458 + ar1*xr11 - ai1*xi11
1459 + ar2*xr21 - ai2*xi21 ;
1460 colY1[iloc] -= ar0*xi01 + ai0*xr01
1461 + ar1*xi11 + ai1*xr11
1462 + ar2*xi21 + ai2*xr21 ;
1463 }
1464 } else {
1465 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1466 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1467 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1468 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1469 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
1470 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1471 colY0[rloc] -= ar0*xr00 - ai0*xi00
1472 + ar1*xr10 - ai1*xi10
1473 + ar2*xr20 - ai2*xi20 ;
1474 colY0[iloc] -= ar0*xi00 + ai0*xr00
1475 + ar1*xi10 + ai1*xr10
1476 + ar2*xi20 + ai2*xr20 ;
1477 colY1[rloc] -= ar0*xr01 - ai0*xi01
1478 + ar1*xr11 - ai1*xi11
1479 + ar2*xr21 - ai2*xi21 ;
1480 colY1[iloc] -= ar0*xi01 + ai0*xr01
1481 + ar1*xi11 + ai1*xr11
1482 + ar2*xi21 + ai2*xr21 ;
1483 }
1484 }
1485 colAT0 = colAT2 + 2*nrowAT ;
1486 }
1487 if ( icolAT == ncolAT - 2 ) {
1488 colAT1 = colAT0 + 2*nrowAT ;
1489 if ( ncolAT == nrowX ) {
1490 rloc = 2*icolAT ; iloc = rloc + 1 ;
1491 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1492 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1493 rloc += 2 ; iloc += 2 ;
1494 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1495 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1496 rloc += 2 ; iloc += 2 ;
1497 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1498 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1499 } else {
1500 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1501 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1502 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1503 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1504 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1505 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1506 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1507 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1508 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1509 }
1510 if ( nrowY == nrowAT ) {
1511 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1512 rloc = 2*krowAT ; iloc = rloc + 1 ;
1513 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1514 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1515 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1516 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1517 colY1[rloc] -= ar0*xr01 - ai0*xi01 + ar1*xr11 - ai1*xi11 ;
1518 colY1[iloc] -= ar0*xi01 + ai0*xr01 + ar1*xi11 + ai1*xr11 ;
1519 }
1520 } else {
1521 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1522 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1523 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1524 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1525 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1526 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1527 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1528 colY1[rloc] -= ar0*xr01 - ai0*xi01 + ar1*xr11 - ai1*xi11 ;
1529 colY1[iloc] -= ar0*xi01 + ai0*xr01 + ar1*xi11 + ai1*xr11 ;
1530 }
1531 }
1532 } else if ( icolAT == ncolAT - 1 ) {
1533 if ( ncolAT == nrowX ) {
1534 rloc = 2*icolAT ; iloc = rloc + 1 ;
1535 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1536 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1537 rloc += 2 ; iloc += 2 ;
1538 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1539 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1540 rloc += 2 ; iloc += 2 ;
1541 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1542 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1543 } else {
1544 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1545 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1546 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
1547 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1548 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1549 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
1550 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1551 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1552 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
1553 }
1554 if ( nrowY == nrowAT ) {
1555 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1556 rloc = 2*krowAT ; iloc = rloc + 1 ;
1557 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1558 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1559 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1560 colY1[rloc] -= ar0*xr01 - ai0*xi01 ;
1561 colY1[iloc] -= ar0*xi01 + ai0*xr01 ;
1562 }
1563 } else {
1564 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1565 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1566 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1567 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1568 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1569 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1570 colY1[rloc] -= ar0*xr01 - ai0*xi01 ;
1571 colY1[iloc] -= ar0*xi01 + ai0*xr01 ;
1572 }
1573 }
1574 }
1575 } else if ( jcolX == ncolX - 1 ) {
1576 colAT0 = entA ;
1577 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
1578 /*
1579 fprintf(stdout, "\n %% icolAT = %d", icolAT) ;
1580 */
1581 colAT1 = colAT0 + 2*nrowAT ;
1582 colAT2 = colAT1 + 2*nrowAT ;
1583 if ( ncolAT == nrowX ) {
1584 rloc = 2*icolAT ; iloc = rloc + 1 ;
1585 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1586 rloc += 2 ; iloc += 2 ;
1587 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1588 rloc += 2 ; iloc += 2 ;
1589 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1590 } else {
1591 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1592 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1593 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1594 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1595 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1596 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1597 }
1598 /*
1599 fprintf(stdout, "\n %% x00 = (%12.4e,%12.4e)", xr00, xi00) ;
1600 fprintf(stdout, "\n %% x10 = (%12.4e,%12.4e)", xr10, xi10) ;
1601 fprintf(stdout, "\n %% x20 = (%12.4e,%12.4e)", xr20, xi20) ;
1602 */
1603 if ( nrowY == nrowAT ) {
1604 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1605 rloc = 2*krowAT ; iloc = rloc + 1 ;
1606 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1607 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1608 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
1609 colY0[rloc] -= ar0*xr00 - ai0*xi00
1610 + ar1*xr10 - ai1*xi10
1611 + ar2*xr20 - ai2*xi20 ;
1612 colY0[iloc] -= ar0*xi00 + ai0*xr00
1613 + ar1*xi10 + ai1*xr10
1614 + ar2*xi20 + ai2*xr20 ;
1615 }
1616 } else {
1617 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1618 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1619 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1620 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1621 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
1622 /*
1623 fprintf(stdout, "\n %% raloc = %d, ialoc = %d", raloc, ialoc) ;
1624 fprintf(stdout, "\n %% a0 = (%12.4e,%12.4e)", ar0, ai0) ;
1625 fprintf(stdout, "\n %% a1 = (%12.4e,%12.4e)", ar1, ai1) ;
1626 fprintf(stdout, "\n %% a2 = (%12.4e,%12.4e)", ar2, ai2) ;
1627 */
1628 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1629 /*
1630 fprintf(stdout, "\n %% rloc = %d, iloc = %d", rloc, iloc) ;
1631 */
1632 colY0[rloc] -= ar0*xr00 - ai0*xi00
1633 + ar1*xr10 - ai1*xi10
1634 + ar2*xr20 - ai2*xi20 ;
1635 colY0[iloc] -= ar0*xi00 + ai0*xr00
1636 + ar1*xi10 + ai1*xr10
1637 + ar2*xi20 + ai2*xr20 ;
1638 }
1639 }
1640 colAT0 = colAT2 + 2*nrowAT ;
1641 }
1642 if ( icolAT == ncolAT - 2 ) {
1643 colAT1 = colAT0 + 2*nrowAT ;
1644 if ( ncolAT == nrowX ) {
1645 rloc = 2*icolAT ; iloc = rloc + 1 ;
1646 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1647 rloc += 2 ; iloc += 2 ;
1648 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1649 rloc += 2 ; iloc += 2 ;
1650 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1651 } else {
1652 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1653 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1654 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1655 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1656 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1657 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1658 }
1659 if ( nrowY == nrowAT ) {
1660 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1661 rloc = 2*krowAT ; iloc = rloc + 1 ;
1662 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1663 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
1664 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1665 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1666 }
1667 } else {
1668 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1669 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1670 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1671 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
1672 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1673 colY0[rloc] -= ar0*xr00 - ai0*xi00 + ar1*xr10 - ai1*xi10 ;
1674 colY0[iloc] -= ar0*xi00 + ai0*xr00 + ar1*xi10 + ai1*xr10 ;
1675 }
1676 }
1677 } else if ( icolAT == ncolAT - 1 ) {
1678 if ( ncolAT == nrowX ) {
1679 rloc = 2*icolAT ; iloc = rloc + 1 ;
1680 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1681 rloc += 2 ; iloc += 2 ;
1682 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1683 rloc += 2 ; iloc += 2 ;
1684 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1685 } else {
1686 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
1687 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
1688 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
1689 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
1690 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
1691 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
1692 }
1693 if ( nrowY == nrowAT ) {
1694 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1695 rloc = 2*krowAT ; iloc = rloc + 1 ;
1696 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
1697 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1698 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1699 }
1700 } else {
1701 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
1702 raloc = 2*krowAT ; ialoc = raloc + 1 ;
1703 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
1704 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
1705 colY0[rloc] -= ar0*xr00 - ai0*xi00 ;
1706 colY0[iloc] -= ar0*xi00 + ai0*xr00 ;
1707 }
1708 }
1709 }
1710 }
1711 return ; }
1712
1713 /*--------------------------------------------------------------------*/
1714 /*
1715 ----------------
1716 A has dense rows
1717 ----------------
1718 */
1719 static void
complex_updDenseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)1720 complex_updDenseColumns (
1721 SubMtx *mtxY,
1722 SubMtx *mtxA,
1723 SubMtx *mtxX
1724 ) {
1725 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
1726 *rowAT0, *rowAT1, *rowAT2, *entA, *entX, *entY ;
1727 int inc1, inc2, irowAT, jcolX, kcolAT,
1728 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY ;
1729 int *colindAT, *rowindAT ;
1730
1731 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
1732 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
1733 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
1734 if ( ncolAT != nrowX ) {
1735 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
1736 } else {
1737 colindAT = NULL ;
1738 }
1739 if ( nrowAT != nrowY ) {
1740 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
1741 } else {
1742 rowindAT = NULL ;
1743 }
1744 colX0 = entX ;
1745 colY0 = entY ;
1746 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
1747 colX1 = colX0 + 2*nrowX ;
1748 colX2 = colX1 + 2*nrowX ;
1749 colY1 = colY0 + 2*nrowY ;
1750 colY2 = colY1 + 2*nrowY ;
1751 rowAT0 = entA ;
1752 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
1753 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum01, isum02,
1754 isum10, isum11, isum12, isum20, isum21, isum22,
1755 rsum00, rsum01, rsum02, rsum10, rsum11, rsum12,
1756 rsum20, rsum21, rsum22, xi0, xi1, xi2, xr0, xr1, xr2 ;
1757 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1758
1759 isum00 = isum01 = isum02 = isum10 = isum11 = isum12
1760 = isum20 = isum21 = isum22 = 0.0 ;
1761 rsum00 = rsum01 = rsum02 = rsum10 = rsum11 = rsum12
1762 = rsum20 = rsum21 = rsum22 = 0.0 ;
1763 rowAT1 = rowAT0 + 2*ncolAT ;
1764 rowAT2 = rowAT1 + 2*ncolAT ;
1765 if ( ncolAT == nrowX ) {
1766 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1767 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1768 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1769 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
1770 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
1771 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1772 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1773 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1774 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1775 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1776 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1777 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1778 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1779 rsum12 += ar1*xr2 - ai1*xi2 ; isum12 += ar1*xi2 + ai1*xr2 ;
1780 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
1781 rsum21 += ar2*xr1 - ai2*xi1 ; isum21 += ar2*xi1 + ai2*xr1 ;
1782 rsum22 += ar2*xr2 - ai2*xi2 ; isum22 += ar2*xi2 + ai2*xr2 ;
1783 }
1784 } else {
1785 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1786 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1787 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1788 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
1789 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
1790 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1791 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1792 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
1793 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
1794 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1795 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1796 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1797 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1798 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1799 rsum12 += ar1*xr2 - ai1*xi2 ; isum12 += ar1*xi2 + ai1*xr2 ;
1800 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
1801 rsum21 += ar2*xr1 - ai2*xi1 ; isum21 += ar2*xi1 + ai2*xr1 ;
1802 rsum22 += ar2*xr2 - ai2*xi2 ; isum22 += ar2*xi2 + ai2*xr2 ;
1803 }
1804 }
1805 if ( nrowY == nrowAT ) {
1806 rloc = 2*irowAT ; iloc = rloc + 1 ;
1807 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1808 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1809 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1810 rloc+= 2 ; iloc += 2 ;
1811 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1812 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
1813 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
1814 rloc+= 2 ; iloc += 2 ;
1815 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
1816 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
1817 colY2[rloc] -= rsum22 ; colY2[iloc] -= isum22 ;
1818 } else {
1819 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1820 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1821 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1822 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1823 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
1824 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1825 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
1826 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
1827 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
1828 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
1829 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
1830 colY2[rloc] -= rsum22 ; colY2[iloc] -= isum22 ;
1831 }
1832 rowAT0 = rowAT2 + 2*ncolAT ;
1833 }
1834 if ( irowAT == nrowAT - 2 ) {
1835 double ai0, ai1, ar0, ar1, isum00, isum01, isum02,
1836 isum10, isum11, isum12, rsum00, rsum01, rsum02,
1837 rsum10, rsum11, rsum12, xi0, xi1, xi2, xr0, xr1, xr2 ;
1838 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1839
1840 isum00 = isum01 = isum02 = isum10 = isum11 = isum12 = 0.0 ;
1841 rsum00 = rsum01 = rsum02 = rsum10 = rsum11 = rsum12 = 0.0 ;
1842 rowAT1 = rowAT0 + 2*ncolAT ;
1843 if ( ncolAT == nrowX ) {
1844 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1845 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1846 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1847 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
1848 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1849 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1850 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1851 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1852 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1853 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1854 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1855 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1856 rsum12 += ar1*xr2 - ai1*xi2 ; isum12 += ar1*xi2 + ai1*xr2 ;
1857 }
1858 } else {
1859 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1860 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1861 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1862 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
1863 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1864 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1865 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
1866 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
1867 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1868 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1869 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1870 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1871 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1872 rsum12 += ar1*xr2 - ai1*xi2 ; isum12 += ar1*xi2 + ai1*xr2 ;
1873 }
1874 }
1875 if ( nrowY == nrowAT ) {
1876 rloc = 2*irowAT ; iloc = rloc + 1 ;
1877 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1878 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1879 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1880 rloc+= 2 ; iloc += 2 ;
1881 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1882 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
1883 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
1884 } else {
1885 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1886 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1887 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1888 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1889 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
1890 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1891 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
1892 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
1893 }
1894 } else if ( irowAT == nrowAT - 1 ) {
1895 double ai0, ar0, isum00, isum01, isum02,
1896 rsum00, rsum01, rsum02, xi0, xi1, xi2, xr0, xr1, xr2 ;
1897 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1898
1899 isum00 = isum01 = isum02 = 0.0 ;
1900 rsum00 = rsum01 = rsum02 = 0.0 ;
1901 if ( ncolAT == nrowX ) {
1902 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1903 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1904 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1905 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1906 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1907 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1908 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1909 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1910 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1911 }
1912 } else {
1913 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1914 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1915 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1916 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1917 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1918 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
1919 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
1920 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1921 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1922 rsum02 += ar0*xr2 - ai0*xi2 ; isum02 += ar0*xi2 + ai0*xr2 ;
1923 }
1924 }
1925 if ( nrowY == nrowAT ) {
1926 rloc = 2*irowAT ; iloc = rloc + 1 ;
1927 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1928 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1929 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1930 } else {
1931 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1932 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1933 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1934 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
1935 }
1936 }
1937 colX0 = colX2 + 2*nrowX ;
1938 colY0 = colY2 + 2*nrowY ;
1939 }
1940 if ( jcolX == ncolX - 2 ) {
1941 colX1 = colX0 + 2*nrowX ;
1942 colY1 = colY0 + 2*nrowY ;
1943 rowAT0 = entA ;
1944 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
1945 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum01,
1946 isum10, isum11, isum20, isum21, rsum00, rsum01, rsum10,
1947 rsum11, rsum20, rsum21, xi0, xi1, xr0, xr1 ;
1948 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1949
1950 isum00 = isum01 = isum10 = isum11 = isum20 = isum21 = 0.0 ;
1951 rsum00 = rsum01 = rsum10 = rsum11 = rsum20 = rsum21 = 0.0 ;
1952 rowAT1 = rowAT0 + 2*ncolAT ;
1953 rowAT2 = rowAT1 + 2*ncolAT ;
1954 if ( ncolAT == nrowX ) {
1955 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1956 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1957 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1958 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
1959 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
1960 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1961 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1962 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1963 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1964 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1965 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1966 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
1967 rsum21 += ar2*xr1 - ai2*xi1 ; isum21 += ar2*xi1 + ai2*xr1 ;
1968 }
1969 } else {
1970 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1971 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1972 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1973 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
1974 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
1975 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1976 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1977 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
1978 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
1979 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
1980 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
1981 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
1982 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
1983 rsum21 += ar2*xr1 - ai2*xi1 ; isum21 += ar2*xi1 + ai2*xr1 ;
1984 }
1985 }
1986 if ( nrowY == nrowAT ) {
1987 rloc = 2*irowAT ; iloc = rloc + 1 ;
1988 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1989 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1990 rloc+= 2 ; iloc += 2 ;
1991 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1992 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
1993 rloc+= 2 ; iloc += 2 ;
1994 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
1995 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
1996 } else {
1997 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1998 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1999 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
2000 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
2001 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2002 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
2003 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
2004 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
2005 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
2006 }
2007 rowAT0 = rowAT2 + 2*ncolAT ;
2008 }
2009 if ( irowAT == nrowAT - 2 ) {
2010 double ai0, ai1, ar0, ar1, isum00, isum01, isum10, isum11,
2011 rsum00, rsum01, rsum10, rsum11, xi0, xi1, xr0, xr1 ;
2012 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
2013
2014 isum00 = isum01 = isum10 = isum11 = 0.0 ;
2015 rsum00 = rsum01 = rsum10 = rsum11 = 0.0 ;
2016 rowAT1 = rowAT0 + 2*ncolAT ;
2017 if ( ncolAT == nrowX ) {
2018 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2019 rloc = 2*kcolAT ; iloc = rloc + 1 ;
2020 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
2021 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
2022 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2023 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2024 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2025 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
2026 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2027 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
2028 }
2029 } else {
2030 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2031 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
2032 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
2033 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
2034 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
2035 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
2036 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
2037 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2038 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
2039 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2040 rsum11 += ar1*xr1 - ai1*xi1 ; isum11 += ar1*xi1 + ai1*xr1 ;
2041 }
2042 }
2043 if ( nrowY == nrowAT ) {
2044 rloc = 2*irowAT ; iloc = rloc + 1 ;
2045 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2046 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
2047 rloc+= 2 ; iloc += 2 ;
2048 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2049 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
2050 } else {
2051 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2052 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2053 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
2054 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
2055 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2056 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
2057 }
2058 } else if ( irowAT == nrowAT - 1 ) {
2059 double ai0, ar0, isum00, isum01,
2060 rsum00, rsum01, xi0, xi1, xr0, xr1 ;
2061 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
2062
2063 isum00 = isum01 = 0.0 ;
2064 rsum00 = rsum01 = 0.0 ;
2065 if ( ncolAT == nrowX ) {
2066 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2067 rloc = 2*kcolAT ; iloc = rloc + 1 ;
2068 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
2069 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2070 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2071 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2072 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
2073 }
2074 } else {
2075 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2076 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
2077 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
2078 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
2079 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
2080 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
2081 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2082 rsum01 += ar0*xr1 - ai0*xi1 ; isum01 += ar0*xi1 + ai0*xr1 ;
2083 }
2084 }
2085 if ( nrowY == nrowAT ) {
2086 rloc = 2*irowAT ; iloc = rloc + 1 ;
2087 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2088 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
2089 } else {
2090 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2091 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2092 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
2093 }
2094 }
2095 } else if ( jcolX == ncolX - 1 ) {
2096 rowAT0 = entA ;
2097 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
2098 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum10, isum20,
2099 rsum00, rsum10, rsum20, xi0, xr0 ;
2100 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
2101
2102 isum00 = isum10 = isum20 = 0.0 ;
2103 rsum00 = rsum10 = rsum20 = 0.0 ;
2104 rowAT1 = rowAT0 + 2*ncolAT ;
2105 rowAT2 = rowAT1 + 2*ncolAT ;
2106 /*
2107 fprintf(stdout, "\n %% irowAT %d", irowAT) ;
2108 */
2109 if ( ncolAT == nrowX ) {
2110 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2111 rloc = 2*kcolAT ; iloc = rloc + 1 ;
2112 /*
2113 fprintf(stdout, "\n %% rloc %d, iloc %d", rloc, iloc) ;
2114 */
2115 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
2116 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
2117 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
2118 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2119 /*
2120 fprintf(stdout, "\n %% a0 = (%12.4e,%12.4e)", ar0, ai0) ;
2121 fprintf(stdout, "\n %% a1 = (%12.4e,%12.4e)", ar1, ai1) ;
2122 fprintf(stdout, "\n %% a2 = (%12.4e,%12.4e)", ar2, ai2) ;
2123 fprintf(stdout, "\n %% x0 = (%12.4e,%12.4e)", xr0, xi0) ;
2124 */
2125 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2126 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2127 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
2128 }
2129 } else {
2130 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2131 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
2132 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
2133 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
2134 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
2135 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
2136 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
2137 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2138 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2139 rsum20 += ar2*xr0 - ai2*xi0 ; isum20 += ar2*xi0 + ai2*xr0 ;
2140 }
2141 }
2142 if ( nrowY == nrowAT ) {
2143 rloc = 2*irowAT ; iloc = rloc + 1 ;
2144 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2145 rloc+= 2 ; iloc += 2 ;
2146 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2147 rloc+= 2 ; iloc += 2 ;
2148 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
2149 } else {
2150 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2151 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2152 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
2153 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2154 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
2155 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
2156 }
2157 rowAT0 = rowAT2 + 2*ncolAT ;
2158 }
2159 if ( irowAT == nrowAT - 2 ) {
2160 double ai0, ai1, ar0, ar1, isum00, isum10, rsum00, rsum10,
2161 xi0, xr0 ;
2162 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
2163
2164 isum00 = isum10 = 0.0 ;
2165 rsum00 = rsum10 = 0.0 ;
2166 rowAT1 = rowAT0 + 2*ncolAT ;
2167 if ( ncolAT == nrowX ) {
2168 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2169 rloc = 2*kcolAT ; iloc = rloc + 1 ;
2170 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
2171 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
2172 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2173 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2174 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2175 }
2176 } else {
2177 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2178 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
2179 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
2180 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
2181 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
2182 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
2183 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2184 rsum10 += ar1*xr0 - ai1*xi0 ; isum10 += ar1*xi0 + ai1*xr0 ;
2185 }
2186 }
2187 if ( nrowY == nrowAT ) {
2188 rloc = 2*irowAT ; iloc = rloc + 1 ;
2189 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2190 rloc+= 2 ; iloc += 2 ;
2191 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2192 } else {
2193 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2194 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2195 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
2196 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
2197 }
2198 } else if ( irowAT == nrowAT - 1 ) {
2199 double ai0, ar0, isum00, rsum00, xi0, xr0 ;
2200 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
2201
2202 isum00 = 0.0 ;
2203 rsum00 = 0.0 ;
2204 if ( ncolAT == nrowX ) {
2205 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2206 rloc = 2*kcolAT ; iloc = rloc + 1 ;
2207 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
2208 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2209 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2210 }
2211 } else {
2212 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
2213 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
2214 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
2215 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
2216 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
2217 rsum00 += ar0*xr0 - ai0*xi0 ; isum00 += ar0*xi0 + ai0*xr0 ;
2218 }
2219 }
2220 if ( nrowY == nrowAT ) {
2221 rloc = 2*irowAT ; iloc = rloc + 1 ;
2222 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2223 } else {
2224 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2225 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
2226 }
2227 }
2228 }
2229 return ; }
2230
2231 /*--------------------------------------------------------------------*/
2232 /*
2233 --------------------
2234 A has sparse columns
2235 --------------------
2236 */
2237 static void
complex_updSparseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)2238 complex_updSparseColumns (
2239 SubMtx *mtxY,
2240 SubMtx *mtxA,
2241 SubMtx *mtxX
2242 ) {
2243 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
2244 *entA, *entX, *entY ;
2245 int ii, iloc, inc1, inc2, irowAT, jcolX, kk, ncolAT, ncolX,
2246 ncolY, nentA, nrowAT, nrowX, nrowY, rloc, size ;
2247 int *colindAT, *indices, *rowindAT, *sizes ;
2248 /*
2249 fprintf(stdout, "\n UPDATE_SPARSE_ROWS(%d,%d)",
2250 mtxA->rowid, mtxA->colid) ;
2251 */
2252 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
2253 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
2254 SubMtx_sparseColumnsInfo(mtxA, &nrowAT, &nentA, &sizes, &indices, &entA) ;
2255 if ( (ncolAT = mtxA->nrow) != nrowX ) {
2256 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
2257 } else {
2258 colindAT = NULL ;
2259 }
2260 if ( (nrowAT = mtxA->ncol) != nrowY ) {
2261 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
2262 } else {
2263 rowindAT = NULL ;
2264 }
2265 colX0 = entX ;
2266 colY0 = entY ;
2267 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
2268 double ai, ar, isum0, isum1, isum2, rsum0, rsum1, rsum2,
2269 xi0, xi1, xi2, xr0, xr1, xr2 ;
2270
2271 colX1 = colX0 + 2*nrowX ;
2272 colX2 = colX1 + 2*nrowX ;
2273 colY1 = colY0 + 2*nrowY ;
2274 colY2 = colY1 + 2*nrowY ;
2275 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
2276 if ( (size = sizes[irowAT]) > 0 ) {
2277 isum0 = isum1 = isum2 = rsum0 = rsum1 = rsum2 = 0.0 ;
2278 if ( ncolAT == nrowX ) {
2279 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2280 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2281 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2282 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2283 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2284 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
2285 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2286 rsum1 += ar*xr1 - ai*xi1 ; isum1 += ar*xi1 + ai*xr1 ;
2287 rsum2 += ar*xr2 - ai*xi2 ; isum2 += ar*xi2 + ai*xr2 ;
2288 }
2289 } else {
2290 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2291 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2292 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
2293 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2294 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2295 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
2296 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2297 rsum1 += ar*xr1 - ai*xi1 ; isum1 += ar*xi1 + ai*xr1 ;
2298 rsum2 += ar*xr2 - ai*xi2 ; isum2 += ar*xi2 + ai*xr2 ;
2299 }
2300 }
2301 if ( nrowAT == nrowY ) {
2302 rloc = 2*irowAT ; iloc = rloc + 1 ;
2303 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2304 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
2305 colY2[rloc] -= rsum2 ; colY2[iloc] -= isum2 ;
2306 } else {
2307 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2308 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2309 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
2310 colY2[rloc] -= rsum2 ; colY2[iloc] -= isum2 ;
2311 }
2312 }
2313 }
2314 colX0 = colX2 + 2*nrowX ;
2315 colY0 = colY2 + 2*nrowY ;
2316 }
2317 if ( jcolX == ncolX - 2 ) {
2318 double ai, ar, isum0, isum1, rsum0, rsum1, xi0, xi1, xr0, xr1 ;
2319
2320 colX1 = colX0 + 2*nrowX ;
2321 colY1 = colY0 + 2*nrowY ;
2322 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
2323 if ( (size = sizes[irowAT]) > 0 ) {
2324 isum0 = isum1 = rsum0 = rsum1 = 0.0 ;
2325 if ( ncolAT == nrowX ) {
2326 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2327 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2328 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2329 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2330 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2331 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2332 rsum1 += ar*xr1 - ai*xi1 ; isum1 += ar*xi1 + ai*xr1 ;
2333 }
2334 } else {
2335 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2336 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2337 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
2338 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2339 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2340 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2341 rsum1 += ar*xr1 - ai*xi1 ; isum1 += ar*xi1 + ai*xr1 ;
2342 }
2343 }
2344 if ( nrowAT == nrowY ) {
2345 rloc = 2*irowAT ; iloc = rloc + 1 ;
2346 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2347 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
2348 } else {
2349 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2350 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2351 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
2352 }
2353 }
2354 }
2355 } else if ( jcolX == ncolX - 1 ) {
2356 double ai, ar, isum0, rsum0, xi0, xr0 ;
2357
2358 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
2359 if ( (size = sizes[irowAT]) > 0 ) {
2360 isum0 = rsum0 = 0.0 ;
2361 if ( ncolAT == nrowX ) {
2362 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2363 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2364 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2365 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2366 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2367 }
2368 } else {
2369 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2370 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2371 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
2372 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2373 rsum0 += ar*xr0 - ai*xi0 ; isum0 += ar*xi0 + ai*xr0 ;
2374 }
2375 }
2376 if ( nrowAT == nrowY ) {
2377 rloc = 2*irowAT ; iloc = rloc + 1 ;
2378 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2379 } else {
2380 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
2381 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
2382 }
2383 }
2384 }
2385 }
2386 return ; }
2387
2388 /*--------------------------------------------------------------------*/
2389 /*
2390 -----------------
2391 A has sparse rows
2392 -----------------
2393 */
2394 static void
complex_updSparseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)2395 complex_updSparseRows (
2396 SubMtx *mtxY,
2397 SubMtx *mtxA,
2398 SubMtx *mtxX
2399 ) {
2400 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
2401 *entA, *entX, *entY ;
2402 int ii, inc1, inc2, jcolAT, jcolX, jrowX, kk,
2403 ncolAT, ncolX, ncolY, nentA, nrowAT, nrowX, nrowY, size ;
2404 int *colindAT, *indices, *rowindAT, *sizes ;
2405 /*
2406 fprintf(stdout, "\n UPDATE_SPARSE_COLUMNS(%d,%d)",
2407 mtxA->rowid, mtxA->colid) ;
2408 */
2409 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
2410 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
2411 SubMtx_sparseRowsInfo(mtxA, &ncolAT, &nentA, &sizes, &indices, &entA) ;
2412 if ( (ncolAT = mtxA->nrow) != nrowX ) {
2413 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
2414 } else {
2415 colindAT = NULL ;
2416 }
2417 if ( (nrowAT = mtxA->ncol) != nrowY ) {
2418 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
2419 } else {
2420 rowindAT = NULL ;
2421 }
2422 colX0 = entX ;
2423 colY0 = entY ;
2424 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
2425 double ai, ar, xi0, xi1, xi2, xr0, xr1, xr2 ;
2426 int iloc, rloc ;
2427
2428 colX1 = colX0 + 2*nrowX ;
2429 colX2 = colX1 + 2*nrowX ;
2430 colY1 = colY0 + 2*nrowY ;
2431 colY2 = colY1 + 2*nrowY ;
2432 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
2433 if ( (size = sizes[jcolAT]) > 0 ) {
2434 if ( ncolAT == nrowX ) {
2435 jrowX = jcolAT ;
2436 } else {
2437 jrowX = colindAT[jcolAT] ;
2438 }
2439 rloc = 2*jrowX ; iloc = rloc + 1 ;
2440 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2441 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2442 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
2443 if ( nrowAT == nrowY ) {
2444 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2445 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2446 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2447 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2448 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2449 colY1[rloc] -= ar*xr1 - ai*xi1 ;
2450 colY1[iloc] -= ar*xi1 + ai*xr1 ;
2451 colY2[rloc] -= ar*xr2 - ai*xi2 ;
2452 colY2[iloc] -= ar*xi2 + ai*xr2 ;
2453 }
2454 } else {
2455 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2456 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2457 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
2458 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2459 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2460 colY1[rloc] -= ar*xr1 - ai*xi1 ;
2461 colY1[iloc] -= ar*xi1 + ai*xr1 ;
2462 colY2[rloc] -= ar*xr2 - ai*xi2 ;
2463 colY2[iloc] -= ar*xi2 + ai*xr2 ;
2464 }
2465 }
2466 }
2467 }
2468 colX0 = colX2 + 2*nrowX ;
2469 colY0 = colY2 + 2*nrowY ;
2470 }
2471 if ( jcolX == ncolX - 2 ) {
2472 double ai, ar, xi0, xi1, xr0, xr1 ;
2473 int iloc, rloc ;
2474
2475 colX1 = colX0 + 2*nrowX ;
2476 colY1 = colY0 + 2*nrowY ;
2477 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
2478 if ( (size = sizes[jcolAT]) > 0 ) {
2479 if ( ncolAT == nrowX ) {
2480 jrowX = jcolAT ;
2481 } else {
2482 jrowX = colindAT[jcolAT] ;
2483 }
2484 rloc = 2*jrowX ; iloc = rloc + 1 ;
2485 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2486 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
2487 if ( nrowAT == nrowY ) {
2488 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2489 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2490 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2491 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2492 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2493 colY1[rloc] -= ar*xr1 - ai*xi1 ;
2494 colY1[iloc] -= ar*xi1 + ai*xr1 ;
2495 }
2496 } else {
2497 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2498 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2499 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
2500 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2501 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2502 colY1[rloc] -= ar*xr1 - ai*xi1 ;
2503 colY1[iloc] -= ar*xi1 + ai*xr1 ;
2504 }
2505 }
2506 }
2507 }
2508 } else if ( jcolX == ncolX - 1 ) {
2509 double ai, ar, xi0, xr0 ;
2510 int iloc, rloc ;
2511
2512 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
2513 if ( (size = sizes[jcolAT]) > 0 ) {
2514 if ( ncolAT == nrowX ) {
2515 jrowX = jcolAT ;
2516 } else {
2517 jrowX = colindAT[jcolAT] ;
2518 }
2519 rloc = 2*jrowX ; iloc = rloc + 1 ;
2520 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
2521 if ( nrowAT == nrowY ) {
2522 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2523 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2524 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
2525 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2526 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2527 }
2528 } else {
2529 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
2530 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
2531 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
2532 colY0[rloc] -= ar*xr0 - ai*xi0 ;
2533 colY0[iloc] -= ar*xi0 + ai*xr0 ;
2534 }
2535 }
2536 }
2537 }
2538 }
2539 return ; }
2540
2541 /*--------------------------------------------------------------------*/
2542