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