1 /* solveupdH.c */
2
3 #include "../SubMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 static void
7 complex_updDenseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
8 static void
9 complex_updDenseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
10 static void
11 complex_updSparseRows ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
12 static void
13 complex_updSparseColumns ( SubMtx *mtxY, SubMtx *mtxA, SubMtx *mtxX ) ;
14 /*--------------------------------------------------------------------*/
15 /*
16 ------------------------------------------------------
17 purpose -- perform the matrix-matrix multiply
18 Y := Y - A^H * X used in the forward and backsolves
19 where
20 (1) rows(A) \subseteq rows(Y)
21 (2) rows(A) are local w.r.t. rows(Y)
22 (3) cols(A) \subseteq rows(X)
23 (4) cols(A) are local w.r.t. rows(X)
24 (5) cols(Y) = cols(X)
25 (6) Y and X have mode SUBMTX_DENSE_COLUMNS
26
27 created -- 98may02, cca
28 -----------------------------------------------------
29 */
30 void
SubMtx_solveupdH(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)31 SubMtx_solveupdH (
32 SubMtx *mtxY,
33 SubMtx *mtxA,
34 SubMtx *mtxX
35 ) {
36 /*
37 ---------------
38 check the input
39 ---------------
40 */
41 if ( mtxY == NULL || mtxA == NULL || mtxX == NULL ) {
42 fprintf(stderr, "\n fatal error in SubMtx_solveupdH(%p,%p,%p)"
43 "\n bad input\n", mtxY, mtxA, mtxX) ;
44 exit(-1) ;
45 }
46 if ( mtxA->type != SPOOLES_COMPLEX ) {
47 fprintf(stderr, "\n fatal error in SubMtx_solveupdH(%p,%p,%p)"
48 "\n Y has type %d, must be SPOOLES_COMPLEX\n",
49 mtxY, mtxA, mtxX, mtxA->type) ;
50 exit(-1) ;
51 }
52 if ( ! SUBMTX_IS_DENSE_COLUMNS(mtxY) ) {
53 fprintf(stderr, "\n fatal error in SubMtx_solveupdH(%p,%p,%p)"
54 "\n Y must have mode SUBMTX_DENSE_COLUMNS\n",
55 mtxY, mtxA, mtxX) ;
56 exit(-1) ;
57 }
58 if ( ! SUBMTX_IS_DENSE_COLUMNS(mtxX) ) {
59 fprintf(stderr, "\n fatal error in SubMtx_solveupdH(%p,%p,%p)"
60 "\n X must have mode SUBMTX_DENSE_COLUMNS\n",
61 mtxY, mtxA, mtxX) ;
62 exit(-1) ;
63 }
64 switch ( mtxA->mode ) {
65 case SUBMTX_DENSE_COLUMNS :
66 complex_updDenseColumns(mtxY, mtxA, mtxX) ;
67 break ;
68 case SUBMTX_DENSE_ROWS :
69 complex_updDenseRows(mtxY, mtxA, mtxX) ;
70 break ;
71 case SUBMTX_SPARSE_ROWS :
72 complex_updSparseRows(mtxY, mtxA, mtxX) ;
73 break ;
74 case SUBMTX_SPARSE_COLUMNS :
75 complex_updSparseColumns(mtxY, mtxA, mtxX) ;
76 break ;
77 default :
78 fprintf(stderr, "\n fatal error in SubMtx_solveupdH(%p,%p,%p)"
79 "\n unsupported mode %d for A\n",
80 mtxY, mtxA, mtxX, mtxA->mode) ;
81 SubMtx_writeForHumanEye(mtxA, stderr) ;
82 exit(-1) ;
83 break ;
84 }
85 return ; }
86
87 /*--------------------------------------------------------------------*/
88 /*
89 ----------------
90 A has dense rows
91 ----------------
92 */
93 static void
complex_updDenseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)94 complex_updDenseRows (
95 SubMtx *mtxY,
96 SubMtx *mtxA,
97 SubMtx *mtxX
98 ) {
99 double ai0, ai1, ai2, ar0, ar1, ar2,
100 xi00, xi01, xi02, xi10, xi11, xi12, xi20, xi21, xi22,
101 xr00, xr01, xr02, xr10, xr11, xr12, xr20, xr21, xr22 ;
102 double *colAT0, *colAT1, *colAT2, *colX0, *colX1, *colX2,
103 *colY0, *colY1, *colY2, *entA, *entX, *entY ;
104 int icolAT, ialoc, iloc, inc1, inc2, jcolX, krowAT,
105 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY, raloc, rloc ;
106 int *colindAT, *rowindAT ;
107
108 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
109 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
110 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
111 colX0 = entX ;
112 colY0 = entY ;
113 if ( ncolAT != nrowX ) {
114 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
115 } else {
116 colindAT = NULL ;
117 }
118 if ( nrowAT != nrowY ) {
119 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
120 } else {
121 rowindAT = NULL ;
122 }
123 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
124 colX1 = colX0 + 2*nrowX ;
125 colX2 = colX1 + 2*nrowX ;
126 colY1 = colY0 + 2*nrowY ;
127 colY2 = colY1 + 2*nrowY ;
128 colAT0 = entA ;
129 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
130 colAT1 = colAT0 + 2*nrowAT ;
131 colAT2 = colAT1 + 2*nrowAT ;
132 if ( ncolAT == nrowX ) {
133 rloc = 2*icolAT ; iloc = rloc + 1 ;
134 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
135 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
136 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
137 rloc += 2 ; iloc += 2 ;
138 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
139 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
140 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
141 rloc += 2 ; iloc += 2 ;
142 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
143 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
144 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
145 } else {
146 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
147 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
148 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
149 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
150 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
151 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
152 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
153 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
154 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
155 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
156 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
157 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
158 }
159 if ( nrowY == nrowAT ) {
160 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
161 rloc = 2*krowAT ; iloc = rloc + 1 ;
162 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
163 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
164 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
165 colY0[rloc] -= ar0*xr00 + ai0*xi00
166 + ar1*xr10 + ai1*xi10
167 + ar2*xr20 + ai2*xi20 ;
168 colY0[iloc] -= ar0*xi00 - ai0*xr00
169 + ar1*xi10 - ai1*xr10
170 + ar2*xi20 - ai2*xr20 ;
171 colY1[rloc] -= ar0*xr01 + ai0*xi01
172 + ar1*xr11 + ai1*xi11
173 + ar2*xr21 + ai2*xi21 ;
174 colY1[iloc] -= ar0*xi01 - ai0*xr01
175 + ar1*xi11 - ai1*xr11
176 + ar2*xi21 - ai2*xr21 ;
177 colY2[rloc] -= ar0*xr02 + ai0*xi02
178 + ar1*xr12 + ai1*xi12
179 + ar2*xr22 + ai2*xi22 ;
180 colY2[iloc] -= ar0*xi02 - ai0*xr02
181 + ar1*xi12 - ai1*xr12
182 + ar2*xi22 - ai2*xr22 ;
183 }
184 } else {
185 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
186 raloc = 2*krowAT ; ialoc = raloc + 1 ;
187 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
188 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
189 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
190 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
191 colY0[rloc] -= ar0*xr00 + ai0*xi00
192 + ar1*xr10 + ai1*xi10
193 + ar2*xr20 + ai2*xi20 ;
194 colY0[iloc] -= ar0*xi00 - ai0*xr00
195 + ar1*xi10 - ai1*xr10
196 + ar2*xi20 - ai2*xr20 ;
197 colY1[rloc] -= ar0*xr01 + ai0*xi01
198 + ar1*xr11 + ai1*xi11
199 + ar2*xr21 + ai2*xi21 ;
200 colY1[iloc] -= ar0*xi01 - ai0*xr01
201 + ar1*xi11 - ai1*xr11
202 + ar2*xi21 - ai2*xr21 ;
203 colY2[rloc] -= ar0*xr02 + ai0*xi02
204 + ar1*xr12 + ai1*xi12
205 + ar2*xr22 + ai2*xi22 ;
206 colY2[iloc] -= ar0*xi02 - ai0*xr02
207 + ar1*xi12 - ai1*xr12
208 + ar2*xi22 - ai2*xr22 ;
209 }
210 }
211 colAT0 = colAT2 + 2*nrowAT ;
212 }
213 if ( icolAT == ncolAT - 2 ) {
214 colAT1 = colAT0 + 2*nrowAT ;
215 if ( ncolAT == nrowX ) {
216 rloc = 2*icolAT ; iloc = rloc + 1 ;
217 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
218 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
219 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
220 rloc += 2 ; iloc += 2 ;
221 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
222 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
223 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
224 rloc += 2 ; iloc += 2 ;
225 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
226 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
227 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
228 } else {
229 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
230 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
231 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
232 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
233 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
234 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
235 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
236 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
237 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
238 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
239 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
240 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
241 }
242 if ( nrowY == nrowAT ) {
243 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
244 rloc = 2*krowAT ; iloc = rloc + 1 ;
245 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
246 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
247 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
248 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
249 colY1[rloc] -= ar0*xr01 + ai0*xi01 + ar1*xr11 + ai1*xi11 ;
250 colY1[iloc] -= ar0*xi01 - ai0*xr01 + ar1*xi11 - ai1*xr11 ;
251 colY2[rloc] -= ar0*xr02 + ai0*xi02 + ar1*xr12 + ai1*xi12 ;
252 colY2[iloc] -= ar0*xi02 - ai0*xr02 + ar1*xi12 - ai1*xr12 ;
253 }
254 } else {
255 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
256 raloc = 2*krowAT ; ialoc = raloc + 1 ;
257 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
258 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
259 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
260 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
261 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
262 colY1[rloc] -= ar0*xr01 + ai0*xi01 + ar1*xr11 + ai1*xi11 ;
263 colY1[iloc] -= ar0*xi01 - ai0*xr01 + ar1*xi11 - ai1*xr11 ;
264 colY2[rloc] -= ar0*xr02 + ai0*xi02 + ar1*xr12 + ai1*xi12 ;
265 colY2[iloc] -= ar0*xi02 - ai0*xr02 + ar1*xi12 - ai1*xr12 ;
266 }
267 }
268 } else if ( icolAT == ncolAT - 1 ) {
269 if ( ncolAT == nrowX ) {
270 rloc = 2*icolAT ; iloc = rloc + 1 ;
271 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
272 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
273 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
274 rloc += 2 ; iloc += 2 ;
275 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
276 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
277 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
278 rloc += 2 ; iloc += 2 ;
279 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
280 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
281 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
282 } else {
283 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
284 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
285 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
286 xr02 = colX2[rloc] ; xi02 = colX2[iloc] ;
287 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
288 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
289 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
290 xr12 = colX2[rloc] ; xi12 = colX2[iloc] ;
291 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
292 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
293 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
294 xr22 = colX2[rloc] ; xi22 = colX2[iloc] ;
295 }
296 if ( nrowY == nrowAT ) {
297 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
298 rloc = 2*krowAT ; iloc = rloc + 1 ;
299 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
300 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
301 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
302 colY1[rloc] -= ar0*xr01 + ai0*xi01 ;
303 colY1[iloc] -= ar0*xi01 - ai0*xr01 ;
304 colY2[rloc] -= ar0*xr02 + ai0*xi02 ;
305 colY2[iloc] -= ar0*xi02 - ai0*xr02 ;
306 }
307 } else {
308 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
309 raloc = 2*krowAT ; ialoc = raloc + 1 ;
310 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
311 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
312 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
313 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
314 colY1[rloc] -= ar0*xr01 + ai0*xi01 ;
315 colY1[iloc] -= ar0*xi01 - ai0*xr01 ;
316 colY2[rloc] -= ar0*xr02 + ai0*xi02 ;
317 colY2[iloc] -= ar0*xi02 - ai0*xr02 ;
318 }
319 }
320 }
321 colX0 = colX2 + 2*nrowX ;
322 colY0 = colY2 + 2*nrowY ;
323 }
324 /*
325 fprintf(stdout, "\n %% jcolX = %d", jcolX) ;
326 */
327 if ( jcolX == ncolX - 2 ) {
328 colX1 = colX0 + 2*nrowX ;
329 colY1 = colY0 + 2*nrowY ;
330 colAT0 = entA ;
331 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
332 colAT1 = colAT0 + 2*nrowAT ;
333 colAT2 = colAT1 + 2*nrowAT ;
334 if ( ncolAT == nrowX ) {
335 rloc = 2*icolAT ; iloc = rloc + 1 ;
336 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
337 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
338 rloc += 2 ; iloc += 2 ;
339 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
340 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
341 rloc += 2 ; iloc += 2 ;
342 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
343 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
344 } else {
345 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
346 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
347 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
348 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
349 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
350 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
351 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
352 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
353 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
354 }
355 if ( nrowY == nrowAT ) {
356 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
357 rloc = 2*krowAT ; iloc = rloc + 1 ;
358 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
359 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
360 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
361 colY0[rloc] -= ar0*xr00 + ai0*xi00
362 + ar1*xr10 + ai1*xi10
363 + ar2*xr20 + ai2*xi20 ;
364 colY0[iloc] -= ar0*xi00 - ai0*xr00
365 + ar1*xi10 - ai1*xr10
366 + ar2*xi20 - ai2*xr20 ;
367 colY1[rloc] -= ar0*xr01 + ai0*xi01
368 + ar1*xr11 + ai1*xi11
369 + ar2*xr21 + ai2*xi21 ;
370 colY1[iloc] -= ar0*xi01 - ai0*xr01
371 + ar1*xi11 - ai1*xr11
372 + ar2*xi21 - ai2*xr21 ;
373 }
374 } else {
375 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
376 raloc = 2*krowAT ; ialoc = raloc + 1 ;
377 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
378 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
379 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
380 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
381 colY0[rloc] -= ar0*xr00 + ai0*xi00
382 + ar1*xr10 + ai1*xi10
383 + ar2*xr20 + ai2*xi20 ;
384 colY0[iloc] -= ar0*xi00 - ai0*xr00
385 + ar1*xi10 - ai1*xr10
386 + ar2*xi20 - ai2*xr20 ;
387 colY1[rloc] -= ar0*xr01 + ai0*xi01
388 + ar1*xr11 + ai1*xi11
389 + ar2*xr21 + ai2*xi21 ;
390 colY1[iloc] -= ar0*xi01 - ai0*xr01
391 + ar1*xi11 - ai1*xr11
392 + ar2*xi21 - ai2*xr21 ;
393 }
394 }
395 colAT0 = colAT2 + 2*nrowAT ;
396 }
397 if ( icolAT == ncolAT - 2 ) {
398 colAT1 = colAT0 + 2*nrowAT ;
399 if ( ncolAT == nrowX ) {
400 rloc = 2*icolAT ; iloc = rloc + 1 ;
401 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
402 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
403 rloc += 2 ; iloc += 2 ;
404 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
405 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
406 rloc += 2 ; iloc += 2 ;
407 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
408 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
409 } else {
410 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
411 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
412 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
413 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
414 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
415 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
416 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
417 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
418 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
419 }
420 if ( nrowY == nrowAT ) {
421 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
422 rloc = 2*krowAT ; iloc = rloc + 1 ;
423 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
424 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
425 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
426 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
427 colY1[rloc] -= ar0*xr01 + ai0*xi01 + ar1*xr11 + ai1*xi11 ;
428 colY1[iloc] -= ar0*xi01 - ai0*xr01 + ar1*xi11 - ai1*xr11 ;
429 }
430 } else {
431 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
432 raloc = 2*krowAT ; ialoc = raloc + 1 ;
433 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
434 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
435 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
436 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
437 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
438 colY1[rloc] -= ar0*xr01 + ai0*xi01 + ar1*xr11 + ai1*xi11 ;
439 colY1[iloc] -= ar0*xi01 - ai0*xr01 + ar1*xi11 - ai1*xr11 ;
440 }
441 }
442 } else if ( icolAT == ncolAT - 1 ) {
443 if ( ncolAT == nrowX ) {
444 rloc = 2*icolAT ; iloc = rloc + 1 ;
445 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
446 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
447 rloc += 2 ; iloc += 2 ;
448 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
449 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
450 rloc += 2 ; iloc += 2 ;
451 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
452 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
453 } else {
454 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
455 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
456 xr01 = colX1[rloc] ; xi01 = colX1[iloc] ;
457 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
458 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
459 xr11 = colX1[rloc] ; xi11 = colX1[iloc] ;
460 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
461 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
462 xr21 = colX1[rloc] ; xi21 = colX1[iloc] ;
463 }
464 if ( nrowY == nrowAT ) {
465 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
466 rloc = 2*krowAT ; iloc = rloc + 1 ;
467 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
468 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
469 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
470 colY1[rloc] -= ar0*xr01 + ai0*xi01 ;
471 colY1[iloc] -= ar0*xi01 - ai0*xr01 ;
472 }
473 } else {
474 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
475 raloc = 2*krowAT ; ialoc = raloc + 1 ;
476 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
477 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
478 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
479 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
480 colY1[rloc] -= ar0*xr01 + ai0*xi01 ;
481 colY1[iloc] -= ar0*xi01 - ai0*xr01 ;
482 }
483 }
484 }
485 } else if ( jcolX == ncolX - 1 ) {
486 colAT0 = entA ;
487 for ( icolAT = 0 ; icolAT < ncolAT - 2 ; icolAT += 3 ) {
488 /*
489 fprintf(stdout, "\n %% icolAT = %d", icolAT) ;
490 */
491 colAT1 = colAT0 + 2*nrowAT ;
492 colAT2 = colAT1 + 2*nrowAT ;
493 if ( ncolAT == nrowX ) {
494 rloc = 2*icolAT ; iloc = rloc + 1 ;
495 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
496 rloc += 2 ; iloc += 2 ;
497 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
498 rloc += 2 ; iloc += 2 ;
499 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
500 } else {
501 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
502 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
503 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
504 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
505 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
506 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
507 }
508 /*
509 fprintf(stdout, "\n %% x00 = (%12.4e,%12.4e)", xr00, xi00) ;
510 fprintf(stdout, "\n %% x10 = (%12.4e,%12.4e)", xr10, xi10) ;
511 fprintf(stdout, "\n %% x20 = (%12.4e,%12.4e)", xr20, xi20) ;
512 */
513 if ( nrowY == nrowAT ) {
514 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
515 rloc = 2*krowAT ; iloc = rloc + 1 ;
516 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
517 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
518 ar2 = colAT2[rloc] ; ai2 = colAT2[iloc] ;
519 /*
520 fprintf(stdout, "\n %% rloc = %d, iloc = %d", rloc, iloc) ;
521 fprintf(stdout, "\n %% a0 = (%12.4e,%12.4e)", ar0, ai0) ;
522 fprintf(stdout, "\n %% a1 = (%12.4e,%12.4e)", ar1, ai1) ;
523 fprintf(stdout, "\n %% a2 = (%12.4e,%12.4e)", ar2, ai2) ;
524 */
525 colY0[rloc] -= ar0*xr00 + ai0*xi00
526 + ar1*xr10 + ai1*xi10
527 + ar2*xr20 + ai2*xi20 ;
528 colY0[iloc] -= ar0*xi00 - ai0*xr00
529 + ar1*xi10 - ai1*xr10
530 + ar2*xi20 - ai2*xr20 ;
531 }
532 } else {
533 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
534 raloc = 2*krowAT ; ialoc = raloc + 1 ;
535 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
536 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
537 ar2 = colAT2[raloc] ; ai2 = colAT2[ialoc] ;
538 /*
539 fprintf(stdout, "\n %% raloc = %d, ialoc = %d", raloc, ialoc) ;
540 fprintf(stdout, "\n %% a0 = (%12.4e,%12.4e)", ar0, ai0) ;
541 fprintf(stdout, "\n %% a1 = (%12.4e,%12.4e)", ar1, ai1) ;
542 fprintf(stdout, "\n %% a2 = (%12.4e,%12.4e)", ar2, ai2) ;
543 */
544 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
545 /*
546 fprintf(stdout, "\n %% rloc = %d, iloc = %d", rloc, iloc) ;
547 */
548 colY0[rloc] -= ar0*xr00 + ai0*xi00
549 + ar1*xr10 + ai1*xi10
550 + ar2*xr20 + ai2*xi20 ;
551 colY0[iloc] -= ar0*xi00 - ai0*xr00
552 + ar1*xi10 - ai1*xr10
553 + ar2*xi20 - ai2*xr20 ;
554 }
555 }
556 colAT0 = colAT2 + 2*nrowAT ;
557 }
558 if ( icolAT == ncolAT - 2 ) {
559 colAT1 = colAT0 + 2*nrowAT ;
560 if ( ncolAT == nrowX ) {
561 rloc = 2*icolAT ; iloc = rloc + 1 ;
562 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
563 rloc += 2 ; iloc += 2 ;
564 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
565 rloc += 2 ; iloc += 2 ;
566 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
567 } else {
568 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
569 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
570 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
571 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
572 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
573 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
574 }
575 if ( nrowY == nrowAT ) {
576 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
577 rloc = 2*krowAT ; iloc = rloc + 1 ;
578 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
579 ar1 = colAT1[rloc] ; ai1 = colAT1[iloc] ;
580 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
581 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
582 }
583 } else {
584 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
585 raloc = 2*krowAT ; ialoc = raloc + 1 ;
586 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
587 ar1 = colAT1[raloc] ; ai1 = colAT1[ialoc] ;
588 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
589 colY0[rloc] -= ar0*xr00 + ai0*xi00 + ar1*xr10 + ai1*xi10 ;
590 colY0[iloc] -= ar0*xi00 - ai0*xr00 + ar1*xi10 - ai1*xr10 ;
591 }
592 }
593 } else if ( icolAT == ncolAT - 1 ) {
594 if ( ncolAT == nrowX ) {
595 rloc = 2*icolAT ; iloc = rloc + 1 ;
596 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
597 rloc += 2 ; iloc += 2 ;
598 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
599 rloc += 2 ; iloc += 2 ;
600 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
601 } else {
602 rloc = 2*colindAT[icolAT] ; iloc = rloc + 1 ;
603 xr00 = colX0[rloc] ; xi00 = colX0[iloc] ;
604 rloc = 2*colindAT[icolAT+1] ; iloc = rloc + 1 ;
605 xr10 = colX0[rloc] ; xi10 = colX0[iloc] ;
606 rloc = 2*colindAT[icolAT+2] ; iloc = rloc + 1 ;
607 xr20 = colX0[rloc] ; xi20 = colX0[iloc] ;
608 }
609 if ( nrowY == nrowAT ) {
610 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
611 rloc = 2*krowAT ; iloc = rloc + 1 ;
612 ar0 = colAT0[rloc] ; ai0 = colAT0[iloc] ;
613 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
614 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
615 }
616 } else {
617 for ( krowAT = 0 ; krowAT < nrowAT ; krowAT++ ) {
618 raloc = 2*krowAT ; ialoc = raloc + 1 ;
619 ar0 = colAT0[raloc] ; ai0 = colAT0[ialoc] ;
620 rloc = 2*rowindAT[krowAT] ; iloc = rloc + 1 ;
621 colY0[rloc] -= ar0*xr00 + ai0*xi00 ;
622 colY0[iloc] -= ar0*xi00 - ai0*xr00 ;
623 }
624 }
625 }
626 }
627 return ; }
628
629 /*--------------------------------------------------------------------*/
630 /*
631 -------------------
632 A has dense columns
633 -------------------
634 */
635 static void
complex_updDenseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)636 complex_updDenseColumns (
637 SubMtx *mtxY,
638 SubMtx *mtxA,
639 SubMtx *mtxX
640 ) {
641 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
642 *rowAT0, *rowAT1, *rowAT2, *entA, *entX, *entY ;
643 int inc1, inc2, irowAT, jcolX, kcolAT,
644 ncolAT, ncolX, ncolY, nrowAT, nrowX, nrowY ;
645 int *colindAT, *rowindAT ;
646
647 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
648 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
649 SubMtx_denseInfo(mtxA, &ncolAT, &nrowAT, &inc1, &inc2, &entA) ;
650 if ( ncolAT != nrowX ) {
651 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
652 } else {
653 colindAT = NULL ;
654 }
655 if ( nrowAT != nrowY ) {
656 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
657 } else {
658 rowindAT = NULL ;
659 }
660 colX0 = entX ;
661 colY0 = entY ;
662 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
663 colX1 = colX0 + 2*nrowX ;
664 colX2 = colX1 + 2*nrowX ;
665 colY1 = colY0 + 2*nrowY ;
666 colY2 = colY1 + 2*nrowY ;
667 rowAT0 = entA ;
668 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
669 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum01, isum02,
670 isum10, isum11, isum12, isum20, isum21, isum22,
671 rsum00, rsum01, rsum02, rsum10, rsum11, rsum12,
672 rsum20, rsum21, rsum22, xi0, xi1, xi2, xr0, xr1, xr2 ;
673 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
674
675 isum00 = isum01 = isum02 = isum10 = isum11 = isum12
676 = isum20 = isum21 = isum22 = 0.0 ;
677 rsum00 = rsum01 = rsum02 = rsum10 = rsum11 = rsum12
678 = rsum20 = rsum21 = rsum22 = 0.0 ;
679 rowAT1 = rowAT0 + 2*ncolAT ;
680 rowAT2 = rowAT1 + 2*ncolAT ;
681 if ( ncolAT == nrowX ) {
682 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
683 rloc = 2*kcolAT ; iloc = rloc + 1 ;
684 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
685 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
686 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
687 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
688 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
689 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
690 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
691 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
692 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
693 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
694 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
695 rsum12 += ar1*xr2 + ai1*xi2 ; isum12 += ar1*xi2 - ai1*xr2 ;
696 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
697 rsum21 += ar2*xr1 + ai2*xi1 ; isum21 += ar2*xi1 - ai2*xr1 ;
698 rsum22 += ar2*xr2 + ai2*xi2 ; isum22 += ar2*xi2 - ai2*xr2 ;
699 }
700 } else {
701 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
702 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
703 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
704 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
705 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
706 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
707 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
708 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
709 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
710 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
711 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
712 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
713 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
714 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
715 rsum12 += ar1*xr2 + ai1*xi2 ; isum12 += ar1*xi2 - ai1*xr2 ;
716 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
717 rsum21 += ar2*xr1 + ai2*xi1 ; isum21 += ar2*xi1 - ai2*xr1 ;
718 rsum22 += ar2*xr2 + ai2*xi2 ; isum22 += ar2*xi2 - ai2*xr2 ;
719 }
720 }
721 if ( nrowY == nrowAT ) {
722 rloc = 2*irowAT ; iloc = rloc + 1 ;
723 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
724 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
725 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
726 rloc+= 2 ; iloc += 2 ;
727 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
728 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
729 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
730 rloc+= 2 ; iloc += 2 ;
731 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
732 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
733 colY2[rloc] -= rsum22 ; colY2[iloc] -= isum22 ;
734 } else {
735 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
736 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
737 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
738 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
739 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
740 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
741 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
742 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
743 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
744 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
745 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
746 colY2[rloc] -= rsum22 ; colY2[iloc] -= isum22 ;
747 }
748 rowAT0 = rowAT2 + 2*ncolAT ;
749 }
750 if ( irowAT == nrowAT - 2 ) {
751 double ai0, ai1, ar0, ar1, isum00, isum01, isum02,
752 isum10, isum11, isum12, rsum00, rsum01, rsum02,
753 rsum10, rsum11, rsum12, xi0, xi1, xi2, xr0, xr1, xr2 ;
754 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
755
756 isum00 = isum01 = isum02 = isum10 = isum11 = isum12 = 0.0 ;
757 rsum00 = rsum01 = rsum02 = rsum10 = rsum11 = rsum12 = 0.0 ;
758 rowAT1 = rowAT0 + 2*ncolAT ;
759 if ( ncolAT == nrowX ) {
760 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
761 rloc = 2*kcolAT ; iloc = rloc + 1 ;
762 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
763 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
764 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
765 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
766 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
767 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
768 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
769 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
770 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
771 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
772 rsum12 += ar1*xr2 + ai1*xi2 ; isum12 += ar1*xi2 - ai1*xr2 ;
773 }
774 } else {
775 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
776 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
777 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
778 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
779 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
780 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
781 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
782 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
783 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
784 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
785 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
786 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
787 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
788 rsum12 += ar1*xr2 + ai1*xi2 ; isum12 += ar1*xi2 - ai1*xr2 ;
789 }
790 }
791 if ( nrowY == nrowAT ) {
792 rloc = 2*irowAT ; iloc = rloc + 1 ;
793 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
794 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
795 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
796 rloc+= 2 ; iloc += 2 ;
797 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
798 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
799 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
800 } else {
801 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
802 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
803 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
804 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
805 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
806 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
807 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
808 colY2[rloc] -= rsum12 ; colY2[iloc] -= isum12 ;
809 }
810 } else if ( irowAT == nrowAT - 1 ) {
811 double ai0, ar0, isum00, isum01, isum02,
812 rsum00, rsum01, rsum02, xi0, xi1, xi2, xr0, xr1, xr2 ;
813 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
814
815 isum00 = isum01 = isum02 = 0.0 ;
816 rsum00 = rsum01 = rsum02 = 0.0 ;
817 if ( ncolAT == nrowX ) {
818 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
819 rloc = 2*kcolAT ; iloc = rloc + 1 ;
820 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
821 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
822 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
823 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
824 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
825 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
826 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
827 }
828 } else {
829 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
830 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
831 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
832 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
833 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
834 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
835 xr2 = colX2[rxloc] ; xi2 = colX2[ixloc] ;
836 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
837 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
838 rsum02 += ar0*xr2 + ai0*xi2 ; isum02 += ar0*xi2 - ai0*xr2 ;
839 }
840 }
841 if ( nrowY == nrowAT ) {
842 rloc = 2*irowAT ; iloc = rloc + 1 ;
843 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
844 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
845 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
846 } else {
847 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
848 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
849 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
850 colY2[rloc] -= rsum02 ; colY2[iloc] -= isum02 ;
851 }
852 }
853 colX0 = colX2 + 2*nrowX ;
854 colY0 = colY2 + 2*nrowY ;
855 }
856 if ( jcolX == ncolX - 2 ) {
857 colX1 = colX0 + 2*nrowX ;
858 colY1 = colY0 + 2*nrowY ;
859 rowAT0 = entA ;
860 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
861 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum01,
862 isum10, isum11, isum20, isum21, rsum00, rsum01, rsum10,
863 rsum11, rsum20, rsum21, xi0, xi1, xr0, xr1 ;
864 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
865
866 isum00 = isum01 = isum10 = isum11 = isum20 = isum21 = 0.0 ;
867 rsum00 = rsum01 = rsum10 = rsum11 = rsum20 = rsum21 = 0.0 ;
868 rowAT1 = rowAT0 + 2*ncolAT ;
869 rowAT2 = rowAT1 + 2*ncolAT ;
870 if ( ncolAT == nrowX ) {
871 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
872 rloc = 2*kcolAT ; iloc = rloc + 1 ;
873 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
874 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
875 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
876 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
877 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
878 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
879 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
880 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
881 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
882 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
883 rsum21 += ar2*xr1 + ai2*xi1 ; isum21 += ar2*xi1 - ai2*xr1 ;
884 }
885 } else {
886 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
887 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
888 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
889 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
890 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
891 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
892 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
893 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
894 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
895 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
896 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
897 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
898 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
899 rsum21 += ar2*xr1 + ai2*xi1 ; isum21 += ar2*xi1 - ai2*xr1 ;
900 }
901 }
902 if ( nrowY == nrowAT ) {
903 rloc = 2*irowAT ; iloc = rloc + 1 ;
904 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
905 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
906 rloc+= 2 ; iloc += 2 ;
907 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
908 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
909 rloc+= 2 ; iloc += 2 ;
910 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
911 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
912 } else {
913 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
914 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
915 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
916 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
917 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
918 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
919 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
920 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
921 colY1[rloc] -= rsum21 ; colY1[iloc] -= isum21 ;
922 }
923 rowAT0 = rowAT2 + 2*ncolAT ;
924 }
925 if ( irowAT == nrowAT - 2 ) {
926 double ai0, ai1, ar0, ar1, isum00, isum01, isum10, isum11,
927 rsum00, rsum01, rsum10, rsum11, xi0, xi1, xr0, xr1 ;
928 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
929
930 isum00 = isum01 = isum10 = isum11 = 0.0 ;
931 rsum00 = rsum01 = rsum10 = rsum11 = 0.0 ;
932 rowAT1 = rowAT0 + 2*ncolAT ;
933 if ( ncolAT == nrowX ) {
934 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
935 rloc = 2*kcolAT ; iloc = rloc + 1 ;
936 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
937 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
938 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
939 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
940 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
941 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
942 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
943 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
944 }
945 } else {
946 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
947 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
948 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
949 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
950 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
951 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
952 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
953 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
954 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
955 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
956 rsum11 += ar1*xr1 + ai1*xi1 ; isum11 += ar1*xi1 - ai1*xr1 ;
957 }
958 }
959 if ( nrowY == nrowAT ) {
960 rloc = 2*irowAT ; iloc = rloc + 1 ;
961 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
962 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
963 rloc+= 2 ; iloc += 2 ;
964 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
965 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
966 } else {
967 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
968 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
969 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
970 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
971 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
972 colY1[rloc] -= rsum11 ; colY1[iloc] -= isum11 ;
973 }
974 } else if ( irowAT == nrowAT - 1 ) {
975 double ai0, ar0, isum00, isum01,
976 rsum00, rsum01, xi0, xi1, xr0, xr1 ;
977 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
978
979 isum00 = isum01 = 0.0 ;
980 rsum00 = rsum01 = 0.0 ;
981 if ( ncolAT == nrowX ) {
982 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
983 rloc = 2*kcolAT ; iloc = rloc + 1 ;
984 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
985 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
986 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
987 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
988 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
989 }
990 } else {
991 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
992 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
993 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
994 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
995 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
996 xr1 = colX1[rxloc] ; xi1 = colX1[ixloc] ;
997 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
998 rsum01 += ar0*xr1 + ai0*xi1 ; isum01 += ar0*xi1 - ai0*xr1 ;
999 }
1000 }
1001 if ( nrowY == nrowAT ) {
1002 rloc = 2*irowAT ; iloc = rloc + 1 ;
1003 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1004 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1005 } else {
1006 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1007 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1008 colY1[rloc] -= rsum01 ; colY1[iloc] -= isum01 ;
1009 }
1010 }
1011 } else if ( jcolX == ncolX - 1 ) {
1012 rowAT0 = entA ;
1013 for ( irowAT = 0 ; irowAT < nrowAT - 2 ; irowAT += 3 ) {
1014 double ai0, ai1, ai2, ar0, ar1, ar2, isum00, isum10, isum20,
1015 rsum00, rsum10, rsum20, xi0, xr0 ;
1016 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1017
1018 isum00 = isum10 = isum20 = 0.0 ;
1019 rsum00 = rsum10 = rsum20 = 0.0 ;
1020 rowAT1 = rowAT0 + 2*ncolAT ;
1021 rowAT2 = rowAT1 + 2*ncolAT ;
1022 /*
1023 fprintf(stdout, "\n %% irowAT %d", irowAT) ;
1024 */
1025 if ( ncolAT == nrowX ) {
1026 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1027 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1028 /*
1029 fprintf(stdout, "\n %% rloc %d, iloc %d", rloc, iloc) ;
1030 */
1031 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1032 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
1033 ar2 = rowAT2[rloc] ; ai2 = rowAT2[iloc] ;
1034 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1035 /*
1036 fprintf(stdout, "\n %% a0 = (%12.4e,%12.4e)", ar0, ai0) ;
1037 fprintf(stdout, "\n %% a1 = (%12.4e,%12.4e)", ar1, ai1) ;
1038 fprintf(stdout, "\n %% a2 = (%12.4e,%12.4e)", ar2, ai2) ;
1039 fprintf(stdout, "\n %% x0 = (%12.4e,%12.4e)", xr0, xi0) ;
1040 */
1041 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1042 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
1043 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
1044 }
1045 } else {
1046 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1047 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1048 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1049 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
1050 ar2 = rowAT2[raloc] ; ai2 = rowAT2[ialoc] ;
1051 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1052 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1053 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1054 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
1055 rsum20 += ar2*xr0 + ai2*xi0 ; isum20 += ar2*xi0 - ai2*xr0 ;
1056 }
1057 }
1058 if ( nrowY == nrowAT ) {
1059 rloc = 2*irowAT ; iloc = rloc + 1 ;
1060 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1061 rloc+= 2 ; iloc += 2 ;
1062 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1063 rloc+= 2 ; iloc += 2 ;
1064 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
1065 } else {
1066 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1067 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1068 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
1069 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1070 rloc = 2*rowindAT[irowAT+2] ; iloc = rloc + 1 ;
1071 colY0[rloc] -= rsum20 ; colY0[iloc] -= isum20 ;
1072 }
1073 rowAT0 = rowAT2 + 2*ncolAT ;
1074 }
1075 if ( irowAT == nrowAT - 2 ) {
1076 double ai0, ai1, ar0, ar1, isum00, isum10, rsum00, rsum10,
1077 xi0, xr0 ;
1078 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1079
1080 isum00 = isum10 = 0.0 ;
1081 rsum00 = rsum10 = 0.0 ;
1082 rowAT1 = rowAT0 + 2*ncolAT ;
1083 if ( ncolAT == nrowX ) {
1084 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1085 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1086 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1087 ar1 = rowAT1[rloc] ; ai1 = rowAT1[iloc] ;
1088 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1089 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1090 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
1091 }
1092 } else {
1093 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1094 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1095 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1096 ar1 = rowAT1[raloc] ; ai1 = rowAT1[ialoc] ;
1097 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1098 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1099 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1100 rsum10 += ar1*xr0 + ai1*xi0 ; isum10 += ar1*xi0 - ai1*xr0 ;
1101 }
1102 }
1103 if ( nrowY == nrowAT ) {
1104 rloc = 2*irowAT ; iloc = rloc + 1 ;
1105 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1106 rloc+= 2 ; iloc += 2 ;
1107 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1108 } else {
1109 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1110 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1111 rloc = 2*rowindAT[irowAT+1] ; iloc = rloc + 1 ;
1112 colY0[rloc] -= rsum10 ; colY0[iloc] -= isum10 ;
1113 }
1114 } else if ( irowAT == nrowAT - 1 ) {
1115 double ai0, ar0, isum00, rsum00, xi0, xr0 ;
1116 int ialoc, iloc, ixloc, raloc, rloc, rxloc ;
1117
1118 isum00 = 0.0 ;
1119 rsum00 = 0.0 ;
1120 if ( ncolAT == nrowX ) {
1121 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1122 rloc = 2*kcolAT ; iloc = rloc + 1 ;
1123 ar0 = rowAT0[rloc] ; ai0 = rowAT0[iloc] ;
1124 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1125 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1126 }
1127 } else {
1128 for ( kcolAT = 0 ; kcolAT < ncolAT ; kcolAT++ ) {
1129 raloc = 2*kcolAT ; ialoc = raloc + 1 ;
1130 ar0 = rowAT0[raloc] ; ai0 = rowAT0[ialoc] ;
1131 rxloc = 2*colindAT[kcolAT] ; ixloc = rxloc + 1 ;
1132 xr0 = colX0[rxloc] ; xi0 = colX0[ixloc] ;
1133 rsum00 += ar0*xr0 + ai0*xi0 ; isum00 += ar0*xi0 - ai0*xr0 ;
1134 }
1135 }
1136 if ( nrowY == nrowAT ) {
1137 rloc = 2*irowAT ; iloc = rloc + 1 ;
1138 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1139 } else {
1140 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1141 colY0[rloc] -= rsum00 ; colY0[iloc] -= isum00 ;
1142 }
1143 }
1144 }
1145 return ; }
1146
1147 /*--------------------------------------------------------------------*/
1148 /*
1149 --------------------
1150 A has sparse columns
1151 --------------------
1152 */
1153 static void
complex_updSparseColumns(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)1154 complex_updSparseColumns (
1155 SubMtx *mtxY,
1156 SubMtx *mtxA,
1157 SubMtx *mtxX
1158 ) {
1159 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
1160 *entA, *entX, *entY ;
1161 int ii, iloc, inc1, inc2, irowAT, jcolX, kk, ncolAT, ncolX,
1162 ncolY, nentA, nrowAT, nrowX, nrowY, rloc, size ;
1163 int *colindAT, *indices, *rowindAT, *sizes ;
1164 /*
1165 fprintf(stdout, "\n UPDATE_SPARSE_ROWS(%d,%d)",
1166 mtxA->rowid, mtxA->colid) ;
1167 */
1168 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
1169 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
1170 SubMtx_sparseColumnsInfo(mtxA, &nrowAT, &nentA, &sizes, &indices, &entA) ;
1171 if ( (ncolAT = mtxA->nrow) != nrowX ) {
1172 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
1173 } else {
1174 colindAT = NULL ;
1175 }
1176 if ( (nrowAT = mtxA->ncol) != nrowY ) {
1177 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
1178 } else {
1179 rowindAT = NULL ;
1180 }
1181 colX0 = entX ;
1182 colY0 = entY ;
1183 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
1184 double ai, ar, isum0, isum1, isum2, rsum0, rsum1, rsum2,
1185 xi0, xi1, xi2, xr0, xr1, xr2 ;
1186
1187 colX1 = colX0 + 2*nrowX ;
1188 colX2 = colX1 + 2*nrowX ;
1189 colY1 = colY0 + 2*nrowY ;
1190 colY2 = colY1 + 2*nrowY ;
1191 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
1192 if ( (size = sizes[irowAT]) > 0 ) {
1193 isum0 = isum1 = isum2 = rsum0 = rsum1 = rsum2 = 0.0 ;
1194 if ( ncolAT == nrowX ) {
1195 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1196 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1197 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1198 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1199 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1200 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1201 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1202 rsum1 += ar*xr1 + ai*xi1 ; isum1 += ar*xi1 - ai*xr1 ;
1203 rsum2 += ar*xr2 + ai*xi2 ; isum2 += ar*xi2 - ai*xr2 ;
1204 }
1205 } else {
1206 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1207 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1208 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
1209 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1210 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1211 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1212 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1213 rsum1 += ar*xr1 + ai*xi1 ; isum1 += ar*xi1 - ai*xr1 ;
1214 rsum2 += ar*xr2 + ai*xi2 ; isum2 += ar*xi2 - ai*xr2 ;
1215 }
1216 }
1217 if ( nrowAT == nrowY ) {
1218 rloc = 2*irowAT ; iloc = rloc + 1 ;
1219 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1220 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
1221 colY2[rloc] -= rsum2 ; colY2[iloc] -= isum2 ;
1222 } else {
1223 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1224 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1225 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
1226 colY2[rloc] -= rsum2 ; colY2[iloc] -= isum2 ;
1227 }
1228 }
1229 }
1230 colX0 = colX2 + 2*nrowX ;
1231 colY0 = colY2 + 2*nrowY ;
1232 }
1233 if ( jcolX == ncolX - 2 ) {
1234 double ai, ar, isum0, isum1, rsum0, rsum1, xi0, xi1, xr0, xr1 ;
1235
1236 colX1 = colX0 + 2*nrowX ;
1237 colY1 = colY0 + 2*nrowY ;
1238 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
1239 if ( (size = sizes[irowAT]) > 0 ) {
1240 isum0 = isum1 = rsum0 = rsum1 = 0.0 ;
1241 if ( ncolAT == nrowX ) {
1242 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1243 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1244 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1245 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1246 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1247 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1248 rsum1 += ar*xr1 + ai*xi1 ; isum1 += ar*xi1 - ai*xr1 ;
1249 }
1250 } else {
1251 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1252 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1253 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
1254 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1255 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1256 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1257 rsum1 += ar*xr1 + ai*xi1 ; isum1 += ar*xi1 - ai*xr1 ;
1258 }
1259 }
1260 if ( nrowAT == nrowY ) {
1261 rloc = 2*irowAT ; iloc = rloc + 1 ;
1262 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1263 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
1264 } else {
1265 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1266 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1267 colY1[rloc] -= rsum1 ; colY1[iloc] -= isum1 ;
1268 }
1269 }
1270 }
1271 } else if ( jcolX == ncolX - 1 ) {
1272 double ai, ar, isum0, rsum0, xi0, xr0 ;
1273
1274 for ( irowAT = kk = 0 ; irowAT < nrowAT ; irowAT++ ) {
1275 if ( (size = sizes[irowAT]) > 0 ) {
1276 isum0 = rsum0 = 0.0 ;
1277 if ( ncolAT == nrowX ) {
1278 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1279 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1280 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1281 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1282 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1283 }
1284 } else {
1285 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1286 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1287 rloc = 2*colindAT[indices[kk]] ; iloc = rloc + 1 ;
1288 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1289 rsum0 += ar*xr0 + ai*xi0 ; isum0 += ar*xi0 - ai*xr0 ;
1290 }
1291 }
1292 if ( nrowAT == nrowY ) {
1293 rloc = 2*irowAT ; iloc = rloc + 1 ;
1294 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1295 } else {
1296 rloc = 2*rowindAT[irowAT] ; iloc = rloc + 1 ;
1297 colY0[rloc] -= rsum0 ; colY0[iloc] -= isum0 ;
1298 }
1299 }
1300 }
1301 }
1302 return ; }
1303
1304 /*--------------------------------------------------------------------*/
1305 /*
1306 -----------------
1307 A has sparse rows
1308 -----------------
1309 */
1310 static void
complex_updSparseRows(SubMtx * mtxY,SubMtx * mtxA,SubMtx * mtxX)1311 complex_updSparseRows (
1312 SubMtx *mtxY,
1313 SubMtx *mtxA,
1314 SubMtx *mtxX
1315 ) {
1316 double *colX0, *colX1, *colX2, *colY0, *colY1, *colY2,
1317 *entA, *entX, *entY ;
1318 int ii, inc1, inc2, jcolAT, jcolX, jrowX, kk,
1319 ncolAT, ncolX, ncolY, nentA, nrowAT, nrowX, nrowY, size ;
1320 int *colindAT, *indices, *rowindAT, *sizes ;
1321 /*
1322 fprintf(stdout, "\n UPDATE_SPARSE_COLUMNS(%d,%d)",
1323 mtxA->rowid, mtxA->colid) ;
1324 */
1325 SubMtx_denseInfo(mtxY, &nrowY, &ncolY, &inc1, &inc2, &entY) ;
1326 SubMtx_denseInfo(mtxX, &nrowX, &ncolX, &inc1, &inc2, &entX) ;
1327 SubMtx_sparseRowsInfo(mtxA, &ncolAT, &nentA, &sizes, &indices, &entA) ;
1328 if ( (ncolAT = mtxA->nrow) != nrowX ) {
1329 SubMtx_rowIndices(mtxA, &ncolAT, &colindAT) ;
1330 } else {
1331 colindAT = NULL ;
1332 }
1333 if ( (nrowAT = mtxA->ncol) != nrowY ) {
1334 SubMtx_columnIndices(mtxA, &nrowAT, &rowindAT) ;
1335 } else {
1336 rowindAT = NULL ;
1337 }
1338 colX0 = entX ;
1339 colY0 = entY ;
1340 for ( jcolX = 0 ; jcolX < ncolX - 2 ; jcolX += 3 ) {
1341 double ai, ar, xi0, xi1, xi2, xr0, xr1, xr2 ;
1342 int iloc, rloc ;
1343
1344 colX1 = colX0 + 2*nrowX ;
1345 colX2 = colX1 + 2*nrowX ;
1346 colY1 = colY0 + 2*nrowY ;
1347 colY2 = colY1 + 2*nrowY ;
1348 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1349 if ( (size = sizes[jcolAT]) > 0 ) {
1350 if ( ncolAT == nrowX ) {
1351 jrowX = jcolAT ;
1352 } else {
1353 jrowX = colindAT[jcolAT] ;
1354 }
1355 rloc = 2*jrowX ; iloc = rloc + 1 ;
1356 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1357 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1358 xr2 = colX2[rloc] ; xi2 = colX2[iloc] ;
1359 if ( nrowAT == nrowY ) {
1360 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1361 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1362 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1363 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1364 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1365 colY1[rloc] -= ar*xr1 + ai*xi1 ;
1366 colY1[iloc] -= ar*xi1 - ai*xr1 ;
1367 colY2[rloc] -= ar*xr2 + ai*xi2 ;
1368 colY2[iloc] -= ar*xi2 - ai*xr2 ;
1369 }
1370 } else {
1371 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1372 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1373 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
1374 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1375 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1376 colY1[rloc] -= ar*xr1 + ai*xi1 ;
1377 colY1[iloc] -= ar*xi1 - ai*xr1 ;
1378 colY2[rloc] -= ar*xr2 + ai*xi2 ;
1379 colY2[iloc] -= ar*xi2 - ai*xr2 ;
1380 }
1381 }
1382 }
1383 }
1384 colX0 = colX2 + 2*nrowX ;
1385 colY0 = colY2 + 2*nrowY ;
1386 }
1387 if ( jcolX == ncolX - 2 ) {
1388 double ai, ar, xi0, xi1, xr0, xr1 ;
1389 int iloc, rloc ;
1390
1391 colX1 = colX0 + 2*nrowX ;
1392 colY1 = colY0 + 2*nrowY ;
1393 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1394 if ( (size = sizes[jcolAT]) > 0 ) {
1395 if ( ncolAT == nrowX ) {
1396 jrowX = jcolAT ;
1397 } else {
1398 jrowX = colindAT[jcolAT] ;
1399 }
1400 rloc = 2*jrowX ; iloc = rloc + 1 ;
1401 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1402 xr1 = colX1[rloc] ; xi1 = colX1[iloc] ;
1403 if ( nrowAT == nrowY ) {
1404 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1405 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1406 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1407 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1408 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1409 colY1[rloc] -= ar*xr1 + ai*xi1 ;
1410 colY1[iloc] -= ar*xi1 - ai*xr1 ;
1411 }
1412 } else {
1413 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1414 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1415 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
1416 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1417 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1418 colY1[rloc] -= ar*xr1 + ai*xi1 ;
1419 colY1[iloc] -= ar*xi1 - ai*xr1 ;
1420 }
1421 }
1422 }
1423 }
1424 } else if ( jcolX == ncolX - 1 ) {
1425 double ai, ar, xi0, xr0 ;
1426 int iloc, rloc ;
1427
1428 for ( jcolAT = kk = 0 ; jcolAT < ncolAT ; jcolAT++ ) {
1429 if ( (size = sizes[jcolAT]) > 0 ) {
1430 if ( ncolAT == nrowX ) {
1431 jrowX = jcolAT ;
1432 } else {
1433 jrowX = colindAT[jcolAT] ;
1434 }
1435 rloc = 2*jrowX ; iloc = rloc + 1 ;
1436 xr0 = colX0[rloc] ; xi0 = colX0[iloc] ;
1437 if ( nrowAT == nrowY ) {
1438 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1439 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1440 rloc = 2*indices[kk] ; iloc = rloc + 1 ;
1441 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1442 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1443 }
1444 } else {
1445 for ( ii = 0 ; ii < size ; ii++, kk++ ) {
1446 ar = entA[2*kk] ; ai = entA[2*kk+1] ;
1447 rloc = 2*rowindAT[indices[kk]] ; iloc = rloc + 1 ;
1448 colY0[rloc] -= ar*xr0 + ai*xi0 ;
1449 colY0[iloc] -= ar*xi0 - ai*xr0 ;
1450 }
1451 }
1452 }
1453 }
1454 }
1455 return ; }
1456
1457 /*--------------------------------------------------------------------*/
1458