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