1 /*  solveH.c  */
2 
3 #include "../SubMtx.h"
4 
5 #define MYDEBUG 0
6 
7 /*--------------------------------------------------------------------*/
8 static void solveDenseSubrows ( SubMtx *mtxA, SubMtx *mtxB ) ;
9 static void solveDenseSubcolumns ( SubMtx *mtxA, SubMtx *mtxB ) ;
10 static void solveSparseRows ( SubMtx *mtxA, SubMtx *mtxB ) ;
11 static void solveSparseColumns ( SubMtx *mtxA, SubMtx *mtxB ) ;
12 /*--------------------------------------------------------------------*/
13 /*
14    -------------------------------------------------------------
15    purpose -- solve (A^H + I) X = B,
16       where
17         (1) X overwrites B
18         (2) A must be strict lower or upper triangular
19         (2) A, B and X are complex
20         (4) columns(A) = rows(X)
21         (5) rows(A)    = rows(B)
22         (6) B has mode SUBMTX_DENSE_COLUMNS
23         (7) if A is SUBMTX_DENSE_SUBROWS or SUBMTX_SPARSE_ROWS
24             then A must be strict lower triangular
25         (8) if A is SUBMTX_DENSE_SUBCOLUMNS or SUBMTX_SPARSE_COLUMNS
26             then A must be strict upper triangular
27 
28    created -- 98may01, cca
29    -------------------------------------------------------------
30 */
31 void
SubMtx_solveH(SubMtx * mtxA,SubMtx * mtxB)32 SubMtx_solveH (
33    SubMtx   *mtxA,
34    SubMtx   *mtxB
35 ) {
36 /*
37    ---------------
38    check the input
39    ---------------
40 */
41 if ( mtxA == NULL || mtxB == NULL ) {
42    fprintf(stderr, "\n fatal error in SubMtx_solveH(%p,%p)"
43            "\n bad input\n", mtxA, mtxB) ;
44    exit(-1) ;
45 }
46 if ( ! SUBMTX_IS_COMPLEX(mtxB) ) {
47    fprintf(stderr, "\n fatal error in SubMtx_solveH(%p,%p)"
48            "\n mtxB has bad type %d\n", mtxA, mtxB, mtxB->type) ;
49    exit(-1) ;
50 }
51 if ( ! SUBMTX_IS_DENSE_COLUMNS(mtxB) ) {
52    fprintf(stderr, "\n fatal error in SubMtx_solveH(%p,%p)"
53            "\n mtxB has bad mode %d\n", mtxA, mtxB, mtxB->mode) ;
54    exit(-1) ;
55 }
56 if ( mtxA->nrow != mtxB->nrow ) {
57    fprintf(stderr, "\n fatal error in SubMtx_solveH(%p,%p)"
58            "\n mtxA->nrow = %d, mtxB->nrwo = %d\n",
59            mtxA, mtxB, mtxA->nrow, mtxB->nrow) ;
60    exit(-1) ;
61 }
62 /*
63    -------------------------
64    switch over the mode of A
65    -------------------------
66 */
67 switch ( mtxA->mode ) {
68 case SUBMTX_DENSE_SUBROWS :
69    solveDenseSubrows(mtxA, mtxB) ;
70    break ;
71 case SUBMTX_SPARSE_ROWS :
72    solveSparseRows(mtxA, mtxB) ;
73    break ;
74 case SUBMTX_DENSE_SUBCOLUMNS :
75    solveDenseSubcolumns(mtxA, mtxB) ;
76    break ;
77 case SUBMTX_SPARSE_COLUMNS :
78    solveSparseColumns(mtxA, mtxB) ;
79    break ;
80 default :
81    fprintf(stderr, "\n fatal error in SubMtx_solveH(%p,%p)"
82            "\n bad mode %d\n", mtxA, mtxB, mtxA->mode) ;
83    exit(-1) ;
84    break ;
85 }
86 return ; }
87 
88 /*--------------------------------------------------------------------*/
89 /*
90    ---------------------------------------
91    purpose -- solve (A^T + I) X = B, where
92      (1) A is strictly upper triangular
93      (2) X overwrites B
94      (B) B has mode SUBMTX_DENSE_COLUMNS
95 
96    created -- 98may01, cca
97    ---------------------------------------
98 */
99 static void
solveDenseSubcolumns(SubMtx * mtxA,SubMtx * mtxB)100 solveDenseSubcolumns (
101    SubMtx   *mtxA,
102    SubMtx   *mtxB
103 ) {
104 double   ai, ar, bi0, bi1, bi2, br0, br1, br2,
105          isum0, isum1, isum2, rsum0, rsum1, rsum2 ;
106 double   *colB0, *colB1, *colB2, *entriesA, *entriesB ;
107 int      first, ii, iloc, inc1, inc2, irowA, jcolB, kk, last,
108          ncolB, nentA, nrowA, nrowB, rloc ;
109 int      *firstlocsA, *sizesA ;
110 /*
111    ----------------------------------------------------
112    extract the pointer and dimensions from two matrices
113    ----------------------------------------------------
114 */
115 SubMtx_denseSubcolumnsInfo(mtxA, &nrowA, &nentA,
116                          &firstlocsA, &sizesA, &entriesA) ;
117 SubMtx_denseInfo(mtxB, &nrowB, &ncolB, &inc1, &inc2, &entriesB) ;
118 #if MYDEBUG > 0
119    fprintf(stdout, "\n nentA = %d", nentA) ;
120    fflush(stdout) ;
121 #endif
122 colB0 = entriesB ;
123 for ( jcolB = 0 ; jcolB < ncolB - 2 ; jcolB += 3 ) {
124    colB1 = colB0 + 2*nrowB ;
125    colB2 = colB1 + 2*nrowB ;
126 #if MYDEBUG > 0
127    fprintf(stdout, "\n %% jcolB = %d", jcolB) ;
128    fflush(stdout) ;
129 #endif
130    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
131 #if MYDEBUG > 0
132       fprintf(stdout, "\n %% irowA %d, size %d", irowA, sizesA[irowA]) ;
133       fflush(stdout) ;
134 #endif
135       if ( sizesA[irowA] > 0 ) {
136          first = firstlocsA[irowA] ;
137          last  = first + sizesA[irowA] - 1 ;
138 #if MYDEBUG > 0
139          fprintf(stdout, ", first %d, last %d", first, last) ;
140          fflush(stdout) ;
141 #endif
142          rsum0 = isum0 = 0.0 ;
143          rsum1 = isum1 = 0.0 ;
144          rsum2 = isum2 = 0.0 ;
145          for ( ii = first ; ii <= last ; ii++, kk++ ) {
146             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
147 #if MYDEBUG > 0
148             fprintf(stdout, "\n %%   A(%d,%d) = (%12.4e,%12.4e)",
149                     irowA+1, ii+1, ar, ai) ;
150             fflush(stdout) ;
151 #endif
152             rloc = 2*ii ; iloc = rloc + 1 ;
153             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
154             br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
155             br2 = colB2[rloc] ; bi2 = colB2[iloc] ;
156             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
157             rsum1 += ar*br1 + ai*bi1 ; isum1 += ar*bi1 - ai*br1 ;
158             rsum2 += ar*br2 + ai*bi2 ; isum2 += ar*bi2 - ai*br2 ;
159          }
160          rloc = 2*irowA ; iloc = rloc + 1 ;
161          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
162          colB1[rloc] -= rsum1 ; colB1[iloc] -= isum1 ;
163          colB2[rloc] -= rsum2 ; colB2[iloc] -= isum2 ;
164       }
165    }
166 #if MYDEBUG > 0
167    fprintf(stdout, "\n %% kk = %d", kk) ;
168    fflush(stdout) ;
169 #endif
170    colB0 = colB2 + 2*nrowB ;
171 }
172 if ( jcolB == ncolB - 2 ) {
173    colB1 = colB0 + 2*nrowB ;
174    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
175 #if MYDEBUG > 0
176       fprintf(stdout, "\n %% irowA %d, size %d", irowA, sizesA[irowA]) ;
177       fflush(stdout) ;
178 #endif
179       if ( sizesA[irowA] > 0 ) {
180          first = firstlocsA[irowA] ;
181          last  = first + sizesA[irowA] - 1 ;
182 #if MYDEBUG > 0
183          fprintf(stdout, ", first %d, last %d", first, last) ;
184          fflush(stdout) ;
185 #endif
186          rsum0 = isum0 = 0.0 ;
187          rsum1 = isum1 = 0.0 ;
188          for ( ii = first ; ii <= last ; ii++, kk++ ) {
189             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
190 #if MYDEBUG > 0
191             fprintf(stdout, "\n %%   A(%d,%d) = (%12.4e,%12.4e)",
192                     irowA+1, ii+1, ar, ai) ;
193             fflush(stdout) ;
194 #endif
195             rloc = 2*ii ; iloc = rloc + 1 ;
196             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
197             br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
198             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
199             rsum1 += ar*br1 + ai*bi1 ; isum1 += ar*bi1 - ai*br1 ;
200          }
201          rloc = 2*irowA ; iloc = rloc + 1 ;
202          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
203          colB1[rloc] -= rsum1 ; colB1[iloc] -= isum1 ;
204       }
205 #if MYDEBUG > 0
206       fprintf(stdout, "\n %% kk = %d", kk) ;
207       fflush(stdout) ;
208 #endif
209    }
210 } else if ( jcolB == ncolB - 1 ) {
211    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
212 #if MYDEBUG > 0
213       fprintf(stdout, "\n %% irowA %d, size %d", irowA, sizesA[irowA]) ;
214       fflush(stdout) ;
215 #endif
216       if ( sizesA[irowA] > 0 ) {
217          first = firstlocsA[irowA] ;
218          last  = first + sizesA[irowA] - 1 ;
219 #if MYDEBUG > 0
220          fprintf(stdout, ", first %d, last %d", first, last) ;
221          fflush(stdout) ;
222 #endif
223          rsum0 = isum0 = 0.0 ;
224          for ( ii = first ; ii <= last ; ii++, kk++ ) {
225             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
226 #if MYDEBUG > 0
227             fprintf(stdout, "\n %%   A(%d,%d) = (%12.4e,%12.4e)",
228                     irowA+1, ii+1, ar, ai) ;
229             fflush(stdout) ;
230 #endif
231             rloc = 2*ii ; iloc = rloc + 1 ;
232             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
233             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
234          }
235          rloc = 2*irowA ; iloc = rloc + 1 ;
236          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
237       }
238 #if MYDEBUG > 0
239       fprintf(stdout, "\n %% kk = %d", kk) ;
240       fflush(stdout) ;
241 #endif
242    }
243 }
244 return ; }
245 
246 /*--------------------------------------------------------------------*/
247 /*
248    ---------------------------------------
249    purpose -- solve (A^T + I) X = B, where
250      (1) A is strictly upper triangular
251      (2) X overwrites B
252      (B) B has mode SUBMTX_DENSE_COLUMNS
253 
254    created -- 98may01, cca
255    ---------------------------------------
256 */
257 static void
solveSparseColumns(SubMtx * mtxA,SubMtx * mtxB)258 solveSparseColumns (
259    SubMtx   *mtxA,
260    SubMtx   *mtxB
261 ) {
262 double   ai, ar, bi0, bi1, bi2, br0, br1, br2,
263          isum0, isum1, isum2, rsum0, rsum1, rsum2 ;
264 double   *colB0, *colB1, *colB2, *entriesA, *entriesB ;
265 int      ii, iloc, inc1, inc2, irowA, jcolB, jj, kk,
266          ncolB, nentA, nrowA, nrowB, rloc, size ;
267 int      *indicesA, *sizesA ;
268 /*
269    ----------------------------------------------------
270    extract the pointer and dimensions from two matrices
271    ----------------------------------------------------
272 */
273 SubMtx_sparseColumnsInfo(mtxA, &nrowA, &nentA,
274                        &sizesA, &indicesA, &entriesA) ;
275 SubMtx_denseInfo(mtxB, &nrowB, &ncolB, &inc1, &inc2, &entriesB) ;
276 colB0 = entriesB ;
277 for ( jcolB = 0 ; jcolB < ncolB - 2 ; jcolB += 3 ) {
278    colB1 = colB0 + 2*nrowB ;
279    colB2 = colB1 + 2*nrowB ;
280 #if MYDEBUG > 0
281    fprintf(stdout, "\n jcolB = %d", jcolB) ;
282    fflush(stdout) ;
283 #endif
284    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
285       if ( (size = sizesA[irowA]) > 0 ) {
286          rsum0 = isum0 = 0.0 ;
287          rsum1 = isum1 = 0.0 ;
288          rsum2 = isum2 = 0.0 ;
289          for ( ii = 0 ; ii < size ; ii++, kk++ ) {
290             ar = entriesA[2*kk] ;
291             ai = entriesA[2*kk+1] ;
292             jj  = indicesA[kk] ;
293             if ( jj < 0 || jj >= irowA ) {
294                fprintf(stderr,
295             "\n fatal error, irowA = %d, kk =%d, ii = %d, jj = %d",
296             irowA, kk, ii, jj) ;
297                exit(-1) ;
298             }
299             rloc = 2*jj ;
300             iloc = rloc + 1 ;
301             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
302             br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
303             br2 = colB2[rloc] ; bi2 = colB2[iloc] ;
304             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
305             rsum1 += ar*br1 + ai*bi1 ; isum1 += ar*bi1 - ai*br1 ;
306             rsum2 += ar*br2 + ai*bi2 ; isum2 += ar*bi2 - ai*br2 ;
307          }
308          rloc = 2*irowA ;
309          iloc = rloc + 1 ;
310          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
311          colB1[rloc] -= rsum1 ; colB1[iloc] -= isum1 ;
312          colB2[rloc] -= rsum2 ; colB2[iloc] -= isum2 ;
313       }
314    }
315    colB0 = colB2 + 2*nrowB ;
316 }
317 if ( jcolB == ncolB - 2 ) {
318    colB1 = colB0 + 2*nrowB ;
319 #if MYDEBUG > 0
320    fprintf(stdout, "\n jcolB = %d", jcolB) ;
321    fflush(stdout) ;
322 #endif
323    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
324       if ( (size = sizesA[irowA]) > 0 ) {
325          rsum0 = isum0 = 0.0 ;
326          rsum1 = isum1 = 0.0 ;
327          for ( ii = 0 ; ii < size ; ii++, kk++ ) {
328             ar = entriesA[2*kk] ;
329             ai = entriesA[2*kk+1] ;
330             jj  = indicesA[kk] ;
331             if ( jj < 0 || jj >= irowA ) {
332                fprintf(stderr,
333             "\n fatal error, irowA = %d, kk =%d, ii = %d, jj = %d",
334             irowA, kk, ii, jj) ;
335                exit(-1) ;
336             }
337             rloc = 2*jj ;
338             iloc = rloc + 1 ;
339             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
340             br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
341             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
342             rsum1 += ar*br1 + ai*bi1 ; isum1 += ar*bi1 - ai*br1 ;
343          }
344          rloc = 2*irowA ;
345          iloc = rloc + 1 ;
346          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
347          colB1[rloc] -= rsum1 ; colB1[iloc] -= isum1 ;
348       }
349    }
350 } else if ( jcolB == ncolB - 1 ) {
351    for ( irowA = kk = 0 ; irowA < nrowA ; irowA++ ) {
352       if ( (size = sizesA[irowA]) > 0 ) {
353          rsum0 = isum0 = 0.0 ;
354          for ( ii = 0 ; ii < size ; ii++, kk++ ) {
355             ar = entriesA[2*kk] ;
356             ai = entriesA[2*kk+1] ;
357             jj  = indicesA[kk] ;
358             if ( jj < 0 || jj >= irowA ) {
359                fprintf(stderr,
360             "\n fatal error, irowA = %d, kk =%d, ii = %d, jj = %d",
361             irowA, kk, ii, jj) ;
362                exit(-1) ;
363             }
364             rloc = 2*jj ;
365             iloc = rloc + 1 ;
366             br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
367             rsum0 += ar*br0 + ai*bi0 ; isum0 += ar*bi0 - ai*br0 ;
368          }
369          rloc = 2*irowA ;
370          iloc = rloc + 1 ;
371          colB0[rloc] -= rsum0 ; colB0[iloc] -= isum0 ;
372       }
373    }
374 }
375 return ; }
376 
377 /*--------------------------------------------------------------------*/
378 /*
379    ---------------------------------------
380    purpose -- solve (I + A^T) X = B, where
381      (1) A is strictly lower triangular
382      (2) X overwrites B
383      (B) B has mode SUBMTX_DENSE_COLUMNS
384 
385    created -- 98may01, cca
386    ---------------------------------------
387 */
388 static void
solveDenseSubrows(SubMtx * mtxA,SubMtx * mtxB)389 solveDenseSubrows (
390    SubMtx   *mtxA,
391    SubMtx   *mtxB
392 ) {
393 double   ai, ar, bi0, bi1, bi2, br0, br1, br2 ;
394 double   *colB0, *colB1, *colB2, *entriesA, *entriesB ;
395 int      colstart, first, iloc, inc1, inc2, irowA, jcolB,
396          jj, kk, last, ncolB, nentA, nrowA, nrowB, rloc ;
397 int      *firstlocsA, *sizesA ;
398 /*
399    ----------------------------------------------------
400    extract the pointer and dimensions from two matrices
401    ----------------------------------------------------
402 */
403 SubMtx_denseSubrowsInfo(mtxA, &nrowA, &nentA,
404                       &firstlocsA, &sizesA, &entriesA) ;
405 SubMtx_denseInfo(mtxB, &nrowB, &ncolB, &inc1, &inc2, &entriesB) ;
406 #if MYDEBUG > 0
407    fprintf(stdout, "\n nrowA = %d, ncolA = %d", nrowA, nentA) ;
408    fflush(stdout) ;
409 #endif
410 colB0 = entriesB ;
411 for ( jcolB = 0 ; jcolB < ncolB - 2 ; jcolB += 3 ) {
412    colB1 = colB0 + 2*nrowB ;
413    colB2 = colB1 + 2*nrowB ;
414 #if MYDEBUG > 0
415    fprintf(stdout, "\n jcolB = %d", jcolB) ;
416    fflush(stdout) ;
417 #endif
418    for ( irowA = nrowA - 1, colstart = nentA ;
419          irowA >= 0 ;
420          irowA-- ) {
421       if ( sizesA[irowA] > 0 ) {
422          first = firstlocsA[irowA] ;
423          last  = first + sizesA[irowA] - 1 ;
424          colstart -= last - first + 1 ;
425          rloc = 2*irowA ;
426          iloc = rloc + 1 ;
427          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
428          br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
429          br2 = colB2[rloc] ; bi2 = colB2[iloc] ;
430          for ( jj = first, kk = colstart ; jj <= last ; jj++, kk++ ) {
431             ar = entriesA[2*kk] ;
432             ai = entriesA[2*kk+1] ;
433             rloc = 2*jj ;
434             iloc = rloc + 1 ;
435             colB0[rloc] -= ar*br0 + ai*bi0 ;
436             colB0[iloc] -= ar*bi0 - ai*br0 ;
437             colB1[rloc] -= ar*br1 + ai*bi1 ;
438             colB1[iloc] -= ar*bi1 - ai*br1 ;
439             colB2[rloc] -= ar*br2 + ai*bi2 ;
440             colB2[iloc] -= ar*bi2 - ai*br2 ;
441          }
442       }
443    }
444    colB0 = colB2 + 2*nrowB ;
445 }
446 if ( jcolB == ncolB - 2 ) {
447    colB1 = colB0 + 2*nrowB ;
448 #if MYDEBUG > 0
449    fprintf(stdout, "\n jcolB = %d", jcolB) ;
450    fflush(stdout) ;
451 #endif
452    for ( irowA = nrowA - 1, colstart = nentA ;
453          irowA >= 0 ;
454          irowA-- ) {
455       if ( sizesA[irowA] > 0 ) {
456          first = firstlocsA[irowA] ;
457          last  = first + sizesA[irowA] - 1 ;
458          colstart -= last - first + 1 ;
459          rloc = 2*irowA ;
460          iloc = rloc + 1 ;
461          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
462          br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
463          for ( jj = first, kk = colstart ; jj <= last ; jj++, kk++ ) {
464             ar = entriesA[2*kk] ;
465             ai = entriesA[2*kk+1] ;
466             rloc = 2*jj ;
467             iloc = rloc + 1 ;
468             colB0[rloc] -= ar*br0 + ai*bi0 ;
469             colB0[iloc] -= ar*bi0 - ai*br0 ;
470             colB1[rloc] -= ar*br1 + ai*bi1 ;
471             colB1[iloc] -= ar*bi1 - ai*br1 ;
472          }
473       }
474    }
475 } else if ( jcolB == ncolB - 1 ) {
476 #if MYDEBUG > 0
477    fprintf(stdout, "\n jcolB = %d", jcolB) ;
478    fflush(stdout) ;
479 #endif
480    for ( irowA = nrowA - 1, colstart = nentA ;
481          irowA >= 0 ;
482          irowA-- ) {
483       if ( sizesA[irowA] > 0 ) {
484          first = firstlocsA[irowA] ;
485          last  = first + sizesA[irowA] - 1 ;
486          colstart -= last - first + 1 ;
487          rloc = 2*irowA ;
488          iloc = rloc + 1 ;
489          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
490          for ( jj = first, kk = colstart ; jj <= last ; jj++, kk++ ) {
491             ar = entriesA[2*kk] ;
492             ai = entriesA[2*kk+1] ;
493             rloc = 2*jj ;
494             iloc = rloc + 1 ;
495             colB0[rloc] -= ar*br0 + ai*bi0 ;
496             colB0[iloc] -= ar*bi0 - ai*br0 ;
497          }
498       }
499    }
500 }
501 return ; }
502 
503 /*--------------------------------------------------------------------*/
504 /*
505    ---------------------------------------
506    purpose -- solve (I + A^T) X = B, where
507      (1) A is strictly lower triangular
508      (2) X overwrites B
509      (B) B has mode SUBMTX_DENSE_COLUMNS
510 
511    created -- 98may01, cca
512    ---------------------------------------
513 */
514 static void
solveSparseRows(SubMtx * mtxA,SubMtx * mtxB)515 solveSparseRows (
516    SubMtx   *mtxA,
517    SubMtx   *mtxB
518 ) {
519 double   ai, ar, bi0, bi1, bi2, br0, br1, br2 ;
520 double   *colB0, *colB1, *colB2, *entriesA, *entriesB ;
521 int      colstart, ii, iloc, inc1, inc2, jcolA, jcolB,
522          jj, kk, ncolB, nentA, nrowA, nrowB, rloc, size ;
523 int      *indicesA, *sizesA ;
524 /*
525    ----------------------------------------------------
526    extract the pointer and dimensions from two matrices
527    ----------------------------------------------------
528 */
529 SubMtx_sparseRowsInfo(mtxA, &nrowA, &nentA,
530                     &sizesA, &indicesA, &entriesA) ;
531 SubMtx_denseInfo(mtxB, &nrowB, &ncolB, &inc1, &inc2, &entriesB) ;
532 #if MYDEBUG > 0
533    fprintf(stdout, "\n nrowA = %d, ncolA = %d", nrowA, nentA) ;
534    fflush(stdout) ;
535 #endif
536 colB0 = entriesB ;
537 for ( jcolB = 0 ; jcolB < ncolB - 2 ; jcolB += 3 ) {
538    colB1 = colB0 + 2*nrowB ;
539    colB2 = colB1 + 2*nrowB ;
540 #if MYDEBUG > 0
541    fprintf(stdout, "\n jcolB = %d", jcolB) ;
542    fflush(stdout) ;
543 #endif
544    for ( jcolA = nrowA - 1, colstart = nentA ;
545          jcolA >= 0 ;
546          jcolA-- ) {
547       if ( (size = sizesA[jcolA]) > 0 ) {
548          colstart -= size ;
549          rloc = 2*jcolA ;
550          iloc = rloc + 1 ;
551          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
552          br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
553          br2 = colB2[rloc] ; bi2 = colB2[iloc] ;
554          for ( ii = 0, kk = colstart ; ii < size ; ii++, kk++ ) {
555             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
556             jj  = indicesA[kk] ;
557             rloc = 2*jj ;
558             iloc = rloc + 1 ;
559             colB0[rloc] -= ar*br0 + ai*bi0 ;
560             colB0[iloc] -= ar*bi0 - ai*br0 ;
561             colB1[rloc] -= ar*br1 + ai*bi1 ;
562             colB1[iloc] -= ar*bi1 - ai*br1 ;
563             colB2[rloc] -= ar*br2 + ai*bi2 ;
564             colB2[iloc] -= ar*bi2 - ai*br2 ;
565          }
566       }
567    }
568    colB0 = colB2 + 2*nrowB ;
569 }
570 if ( jcolB == ncolB - 2 ) {
571    colB1 = colB0 + 2*nrowB ;
572 #if MYDEBUG > 0
573    fprintf(stdout, "\n jcolB = %d", jcolB) ;
574    fflush(stdout) ;
575 #endif
576    for ( jcolA = nrowA - 1, colstart = nentA ;
577          jcolA >= 0 ;
578          jcolA-- ) {
579       if ( (size = sizesA[jcolA]) > 0 ) {
580          colstart -= size ;
581          rloc = 2*jcolA ;
582          iloc = rloc + 1 ;
583          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
584          br1 = colB1[rloc] ; bi1 = colB1[iloc] ;
585          for ( ii = 0, kk = colstart ; ii < size ; ii++, kk++ ) {
586             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
587             jj  = indicesA[kk] ;
588             rloc = 2*jj ;
589             iloc = rloc + 1 ;
590             colB0[rloc] -= ar*br0 + ai*bi0 ;
591             colB0[iloc] -= ar*bi0 - ai*br0 ;
592             colB1[rloc] -= ar*br1 + ai*bi1 ;
593             colB1[iloc] -= ar*bi1 - ai*br1 ;
594          }
595       }
596    }
597 } else if ( jcolB == ncolB - 1 ) {
598 #if MYDEBUG > 0
599    fprintf(stdout, "\n jcolB = %d", jcolB) ;
600    fflush(stdout) ;
601 #endif
602    for ( jcolA = nrowA - 1, colstart = nentA ;
603          jcolA >= 0 ;
604          jcolA-- ) {
605       if ( (size = sizesA[jcolA]) > 0 ) {
606          colstart -= size ;
607          rloc = 2*jcolA ;
608          iloc = rloc + 1 ;
609          br0 = colB0[rloc] ; bi0 = colB0[iloc] ;
610          for ( ii = 0, kk = colstart ; ii < size ; ii++, kk++ ) {
611             ar = entriesA[2*kk] ; ai = entriesA[2*kk+1] ;
612             jj  = indicesA[kk] ;
613             rloc = 2*jj ;
614             iloc = rloc + 1 ;
615             colB0[rloc] -= ar*br0 + ai*bi0 ;
616             colB0[iloc] -= ar*bi0 - ai*br0 ;
617          }
618       }
619    }
620 }
621 return ; }
622 
623 /*--------------------------------------------------------------------*/
624