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