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