1 /* solve.c */
2
3 #include "../ILUMtx.h"
4
5 /*--------------------------------------------------------------------*/
6 /*
7 -------------------------------------------------------------------
8 purpose -- to solve the linear system
9 (L + I)D(I + U) X = B, (U^T + I)D(I + U) X = B or
10 (U^H + I)D(I + U) X = B. X and B are single vectors.
11
12 workDV is a temporary vector.
13 note, if workDV is different than B, then B is unchanged on return.
14 one can have X, B and workDV all point to the same object.
15
16 if pops is not NULL, then on return *pops has been incremented
17 by the number of operations performed in the solve.
18
19 return values ---
20 1 -- normal return
21 -1 -- mtx is NULL
22 -2 -- neqns <= 0
23 -3 -- bad type for mtx
24 -4 -- bad symmetryflag for mtx
25 -5 -- storage mode of L is invalid
26 -6 -- storage mode of U is invalid
27 -7 -- sizesL is NULL
28 -8 -- sizesU is NULL
29 -9 -- p_indL is NULL
30 -10 -- p_indU is NULL
31 -11 -- entD is NULL
32 -12 -- p_entL is NULL
33 -13 -- p_entU is NULL
34 -14 -- X is NULL
35 -15 -- size of X is incorrect
36 -16 -- entries of X are NULL
37 -17 -- B is NULL
38 -18 -- size of B is incorrect
39 -19 -- entries of B are NULL
40 -20 -- workDV is NULL
41 -21 -- size of workDV != neqns
42 -22 -- entries of workDV are NULL
43 -23 -- msglvl > 0 and msgFile is NULL
44
45 created -- 98oct03, cca
46 -------------------------------------------------------------------
47 */
48 int
ILUMtx_solveVector(ILUMtx * mtx,DV * X,DV * B,DV * workDV,double * pops,int msglvl,FILE * msgFile)49 ILUMtx_solveVector (
50 ILUMtx *mtx,
51 DV *X,
52 DV *B,
53 DV *workDV,
54 double *pops,
55 int msglvl,
56 FILE *msgFile
57 ) {
58 double ops ;
59 double *bent, *entD, *work, *xent ;
60 double **p_entL, **p_entU ;
61 int bsize, ieqn, neqns, wsize, xsize ;
62 int *sizesL, *sizesU ;
63 int **p_indL, **p_indU ;
64 /*
65 ---------------
66 check the input
67 ---------------
68 */
69 if ( mtx == NULL ) {
70 fprintf(stderr, "\n error in ILUM_solveVector(), mtx = NULL\n") ;
71 return(-1) ;
72 }
73 if ( (neqns = mtx->neqns) <= 0 ) {
74 fprintf(stderr, "\n error in ILUM_solveVector()"
75 "\n neqns = %d\n", neqns) ;
76 return(-2) ;
77 }
78 if ( !(ILUMTX_IS_REAL(mtx) || ILUMTX_IS_COMPLEX(mtx)) ) {
79 fprintf(stderr, "\n error in ILUM_solveVector()"
80 "\n type = %d\n", mtx->type) ;
81 return(-3) ;
82 }
83 if ( !(ILUMTX_IS_SYMMETRIC(mtx) || ILUMTX_IS_HERMITIAN(mtx)
84 || ILUMTX_IS_NONSYMMETRIC(mtx)) ) {
85 fprintf(stderr, "\n error in ILUMsolveVector()"
86 "\n mtx symmetry = %d\n", mtx->symmetryflag) ;
87 return(-4) ;
88 }
89 if ( !(ILUMTX_IS_L_BY_ROWS(mtx) || ILUMTX_IS_L_BY_COLUMNS(mtx)) ) {
90 fprintf(stderr, "\n error in ILUM_solveVector()"
91 "\n LstorageMode = %d\n", mtx->LstorageMode) ;
92 return(-5) ;
93 }
94 if ( !(ILUMTX_IS_U_BY_ROWS(mtx) || ILUMTX_IS_U_BY_COLUMNS(mtx)) ) {
95 fprintf(stderr, "\n error in ILUM_solveVector()"
96 "\n UstorageMode = %d\n", mtx->UstorageMode) ;
97 return(-6) ;
98 }
99 sizesU = mtx->sizesU ;
100 entD = mtx->entD ;
101 p_indU = mtx->p_indU ;
102 p_entU = mtx->p_entU ;
103 if ( ILUMTX_IS_SYMMETRIC(mtx) || ILUMTX_IS_HERMITIAN(mtx) ) {
104 sizesL = mtx->sizesU ;
105 p_indL = mtx->p_indU ;
106 p_entL = mtx->p_entU ;
107 } else {
108 sizesL = mtx->sizesL ;
109 p_indL = mtx->p_indL ;
110 p_entL = mtx->p_entL ;
111 }
112 if ( sizesL == NULL ) {
113 fprintf(stderr, "\n error in ILUM_solveVector(), sizesL = NULL\n") ;
114 return(-7) ;
115 }
116 if ( sizesU == NULL ) {
117 fprintf(stderr, "\n error in ILUM_solveVector(), sizesU = NULL\n") ;
118 return(-8) ;
119 }
120 if ( p_indL == NULL ) {
121 fprintf(stderr, "\n error in ILUM_solveVector(), p_indL = NULL\n") ;
122 return(-9) ;
123 }
124 if ( p_indU == NULL ) {
125 fprintf(stderr, "\n error in ILUM_solveVector(), p_indU = NULL\n") ;
126 return(-10) ;
127 }
128 if ( entD == NULL ) {
129 fprintf(stderr, "\n error in ILUM_solveVector(), entD = NULL\n") ;
130 return(-11) ;
131 }
132 if ( p_entL == NULL ) {
133 fprintf(stderr, "\n error in ILUM_solveVector(), p_entL = NULL\n") ;
134 return(-12) ;
135 }
136 if ( p_entU == NULL ) {
137 fprintf(stderr, "\n error in ILUM_solveVector(), p_entU = NULL\n") ;
138 return(-13) ;
139 }
140 if ( X == NULL ) {
141 fprintf(stderr, "\n error in ILUM_solveVector(), X = NULL\n") ;
142 return(-14) ;
143 }
144 DV_sizeAndEntries(X, &xsize, &xent) ;
145 if ( (ILUMTX_IS_REAL(mtx) && xsize != neqns)
146 || (ILUMTX_IS_COMPLEX(mtx) && xsize != 2*neqns) ) {
147 fprintf(stderr, "\n error in ILUM_solveVector()"
148 "\n neqns = %d, size of X = %d\n", neqns, xsize) ;
149 return(-15) ;
150 }
151 if ( xent == NULL ) {
152 fprintf(stderr, "\n error in ILUM_solveVector()"
153 "\n entries of X are NULL\n") ;
154 return(-16) ;
155 }
156 if ( B == NULL ) {
157 fprintf(stderr, "\n error in ILUM_solveVector(), B = NULL\n") ;
158 return(-17) ;
159 }
160 DV_sizeAndEntries(B, &bsize, &bent) ;
161 if ( (ILUMTX_IS_REAL(mtx) && bsize != neqns)
162 || (ILUMTX_IS_COMPLEX(mtx) && bsize != 2*neqns) ) {
163 fprintf(stderr, "\n error in ILUM_solveVector()"
164 "\n neqns = %d, size of B = %d\n", neqns, bsize) ;
165 return(-18) ;
166 }
167 if ( bent == NULL ) {
168 fprintf(stderr, "\n error in ILUM_solveVector()"
169 "\n entries of B are NULL\n") ;
170 return(-19) ;
171 }
172 if ( workDV == NULL ) {
173 fprintf(stderr, "\n error in ILUM_solveVector(), workDV = NULL\n") ;
174 return(-20) ;
175 }
176 DV_sizeAndEntries(workDV, &wsize, &work) ;
177 if ( (ILUMTX_IS_REAL(mtx) && wsize != neqns)
178 || (ILUMTX_IS_COMPLEX(mtx) && wsize != 2*neqns) ) {
179 fprintf(stderr, "\n error in ILUM_solveVector()"
180 "\n neqns = %d, size of workDV = %d\n", neqns, wsize) ;
181 return(-21) ;
182 }
183 if ( work == NULL ) {
184 fprintf(stderr, "\n error in ILUM_solveVector()"
185 "\n entries of workDV are NULL\n") ;
186 return(-22) ;
187 }
188 if ( msglvl > 0 && msgFile == NULL ) {
189 fprintf(stderr, "\n error in ILUM_solveVector()"
190 "\n msglvl = %d, msgFile is NULL\n", msglvl) ;
191 return(-23) ;
192 }
193 /*--------------------------------------------------------------------*/
194 /*
195 -------------------------
196 copy B vector into work[]
197 -------------------------
198 */
199 if ( msglvl > 4 ) {
200 fprintf(msgFile, "\n bent") ;
201 DVfprintf(msgFile, neqns, bent) ;
202 }
203 if ( bent != work ) {
204 if ( ILUMTX_IS_REAL(mtx) ) {
205 DVcopy(neqns, work, bent) ;
206 } else {
207 DVcopy(2*neqns, work, bent) ;
208 }
209 }
210 /*
211 -------------
212 solve L Y = B
213 -------------
214 */
215 ops = 0.0 ;
216 if ( ILUMTX_IS_REAL(mtx) ) {
217 if ( mtx->LstorageMode == SPOOLES_BY_COLUMNS ) {
218 int ii, Lsize, *ind ;
219 double yi, *ent ;
220
221 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
222 if ( (Lsize = sizesL[ieqn]) > 0 ) {
223 yi = work[ieqn] ;
224 ind = p_indL[ieqn] ;
225 ent = p_entL[ieqn] ;
226 for ( ii = 0 ; ii < Lsize ; ii++ ) {
227 work[ind[ii]] -= yi * ent[ii] ;
228 }
229 ops += 2*Lsize ;
230 }
231 }
232 } else {
233 int ii, Lsize, *ind ;
234 double sum, *ent ;
235
236 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
237 if ( (Lsize = sizesL[ieqn]) > 0 ) {
238 ind = p_indL[ieqn] ;
239 ent = p_entL[ieqn] ;
240 for ( ii = 0, sum = 0.0 ; ii < Lsize ; ii++ ) {
241 sum += work[ind[ii]] * ent[ii] ;
242 }
243 work[ieqn] -= sum ;
244 ops += 2*Lsize ;
245 }
246 }
247 }
248 } else {
249 if ( mtx->LstorageMode == SPOOLES_BY_COLUMNS ) {
250 int ii, jj, kk, Lsize, *ind ;
251 double LI, LR, yI, yR, *ent ;
252
253 if ( ILUMTX_IS_HERMITIAN(mtx) ) {
254 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
255 if ( (Lsize = sizesL[ieqn]) > 0 ) {
256 yR = work[2*ieqn] ;
257 yI = work[2*ieqn+1] ;
258 ind = p_indL[ieqn] ;
259 ent = p_entL[ieqn] ;
260 for ( ii = kk = 0 ; ii < Lsize ; ii++, kk += 2 ) {
261 LR = ent[kk] ; LI = -ent[kk+1] ;
262 jj = 2*ind[ii] ;
263 work[jj] -= yR*LR - yI*LI ;
264 work[jj+1] -= yR*LI + yI*LR ;
265 }
266 ops += 8*Lsize ;
267 }
268 }
269 } else {
270 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
271 if ( (Lsize = sizesL[ieqn]) > 0 ) {
272 yR = work[2*ieqn] ;
273 yI = work[2*ieqn+1] ;
274 ind = p_indL[ieqn] ;
275 ent = p_entL[ieqn] ;
276 for ( ii = kk = 0 ; ii < Lsize ; ii++, kk += 2 ) {
277 LR = ent[kk] ; LI = ent[kk+1] ;
278 jj = 2*ind[ii] ;
279 work[jj] -= yR*LR - yI*LI ;
280 work[jj+1] -= yR*LI + yI*LR ;
281 }
282 ops += 8*Lsize ;
283 }
284 }
285 }
286 } else {
287 int ii, jj, kk, Lsize, *ind ;
288 double LI, LR, sumI, sumR, yI, yR, *ent ;
289
290 if ( ILUMTX_IS_HERMITIAN(mtx) ) {
291 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
292 if ( (Lsize = sizesL[ieqn]) > 0 ) {
293 ind = p_indL[ieqn] ;
294 ent = p_entL[ieqn] ;
295 sumI = sumR = 0.0 ;
296 for ( ii = kk = 0 ; ii < Lsize ; ii++, kk += 2 ) {
297 jj = 2*ind[ii] ; yR = work[jj] ; yI = work[jj+1] ;
298 LR = ent[kk] ; LI = -ent[kk+1] ;
299 sumR += yR*LR - yI*LI ;
300 sumI += yR*LI + yI*LR ;
301 }
302 work[2*ieqn] -= sumR ;
303 work[2*ieqn+1] -= sumI ;
304 ops += 8*Lsize ;
305 }
306 }
307 } else {
308 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
309 if ( (Lsize = sizesL[ieqn]) > 0 ) {
310 ind = p_indL[ieqn] ;
311 ent = p_entL[ieqn] ;
312 sumI = sumR = 0.0 ;
313 for ( ii = kk = 0 ; ii < Lsize ; ii++, kk += 2 ) {
314 jj = 2*ind[ii] ; yR = work[jj] ; yI = work[jj+1] ;
315 LR = ent[kk] ; LI = ent[kk+1] ;
316 sumR += yR*LR - yI*LI ;
317 sumI += yR*LI + yI*LR ;
318 }
319 work[2*ieqn] -= sumR ;
320 work[2*ieqn+1] -= sumI ;
321 ops += 8*Lsize ;
322 }
323 }
324 }
325 }
326 }
327 if ( msglvl > 2 ) {
328 fprintf(msgFile, "\n %% after forward solve") ;
329 fprintf(msgFile, "\n Y = zeros(%d,1) ;", neqns) ;
330 if ( ILUMTX_IS_REAL(mtx) ) {
331 DV_writeForMatlab(workDV, "Y", msgFile) ;
332 } else {
333 workDV->size = neqns ;
334 ZV_writeForMatlab((ZV *) workDV, "Y", msgFile) ;
335 workDV->size = 2*neqns ;
336 }
337 fflush(msgFile) ;
338 }
339 /*
340 -------------
341 solve D Z = Y
342 -------------
343 */
344 if ( ILUMTX_IS_REAL(mtx) ) {
345 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
346 work[ieqn] = work[ieqn] / entD[ieqn] ;
347 }
348 ops += 2*neqns ;
349 } else {
350 double dI, dR, yI, yR ;
351 int rc ;
352
353 for ( ieqn = 0 ; ieqn < neqns ; ieqn++ ) {
354 rc = Zrecip(entD[2*ieqn], entD[2*ieqn+1], &dR, &dI) ;
355 if ( rc == 1 ) {
356 yR = work[2*ieqn] ; yI = work[2*ieqn+1] ;
357 work[2*ieqn] = yR*dR - yI*dI ;
358 work[2*ieqn+1] = yR*dI + yI*dR ;
359 }
360 }
361 ops += 11*neqns ;
362 }
363 if ( msglvl > 2 ) {
364 fprintf(msgFile, "\n %% after diagonal solve") ;
365 fprintf(msgFile, "\n Z = zeros(%d,1) ;", neqns) ;
366 if ( ILUMTX_IS_REAL(mtx) ) {
367 DV_writeForMatlab(workDV, "Z", msgFile) ;
368 } else {
369 workDV->size = neqns ;
370 ZV_writeForMatlab((ZV *) workDV, "Z", msgFile) ;
371 workDV->size = 2*neqns ;
372 }
373 fflush(msgFile) ;
374 }
375 /*
376 -------------
377 solve U X = Z
378 -------------
379 */
380 if ( ILUMTX_IS_REAL(mtx) ) {
381 if ( mtx->UstorageMode == SPOOLES_BY_COLUMNS ) {
382 int ii, Usize, *ind ;
383 double xi, *ent ;
384
385 for ( ieqn = neqns - 1 ; ieqn >= 0 ; ieqn-- ) {
386 if ( (Usize = sizesU[ieqn]) > 0 ) {
387 xi = work[ieqn] ;
388 ind = p_indU[ieqn] ;
389 ent = p_entU[ieqn] ;
390 for ( ii = 0 ; ii < Usize ; ii++ ) {
391 work[ind[ii]] -= xi * ent[ii] ;
392 }
393 ops += 2*Usize ;
394 }
395 }
396 } else {
397 int ii, Usize, *ind ;
398 double sum, *ent ;
399
400 for ( ieqn = neqns - 1 ; ieqn >= 0 ; ieqn-- ) {
401 if ( (Usize = sizesU[ieqn]) > 0 ) {
402 ind = p_indU[ieqn] ;
403 ent = p_entU[ieqn] ;
404 for ( ii = 0, sum = 0.0 ; ii < Usize ; ii++ ) {
405 sum += work[ind[ii]] * ent[ii] ;
406 }
407 work[ieqn] -= sum ;
408 ops += 2*Usize ;
409 }
410 }
411 }
412 } else {
413 if ( mtx->UstorageMode == SPOOLES_BY_COLUMNS ) {
414 int ii, jj, kk, Usize, *ind ;
415 double UI, UR, xI, xR, *ent ;
416
417 for ( ieqn = neqns - 1 ; ieqn >= 0 ; ieqn-- ) {
418 if ( (Usize = sizesU[ieqn]) > 0 ) {
419 xR = work[2*ieqn] ; xI = work[2*ieqn+1] ;
420 ind = p_indU[ieqn] ;
421 ent = p_entU[ieqn] ;
422 for ( ii = kk = 0 ; ii < Usize ; ii++, kk += 2 ) {
423 jj = 2*ind[ii] ;
424 UR = ent[kk] ; UI = ent[kk+1] ;
425 work[jj] -= xR*UR - xI*UI ;
426 work[jj+1] -= xR*UI + xI*UR ;
427 }
428 ops += 8*Usize ;
429 }
430 }
431 } else {
432 int ii, jj, kk, Usize, *ind ;
433 double sumI, sumR, UI, UR, xI, xR, *ent ;
434
435 for ( ieqn = neqns - 1 ; ieqn >= 0 ; ieqn-- ) {
436 if ( (Usize = sizesU[ieqn]) > 0 ) {
437 ind = p_indU[ieqn] ;
438 ent = p_entU[ieqn] ;
439 sumR = sumI = 0.0 ;
440 for ( ii = kk = 0 ; ii < Usize ; ii++, kk += 2 ) {
441 jj = 2*ind[ii] ; xR = work[jj] ; xI = work[jj+1] ;
442 UR = ent[kk] ; UI = ent[kk+1] ;
443 sumR += xR*UR - xI*UI ;
444 sumI += xR*UI + xI*UR ;
445 }
446 work[2*ieqn] -= sumR ;
447 work[2*ieqn+1] -= sumI ;
448 ops += 8*Usize ;
449 }
450 }
451 }
452 }
453 if ( msglvl > 2 ) {
454 fprintf(msgFile, "\n %% after backward solve") ;
455 fprintf(msgFile, "\n W = zeros(%d,1) ;", neqns) ;
456 if ( ILUMTX_IS_REAL(mtx) ) {
457 DV_writeForMatlab(workDV, "W", msgFile) ;
458 } else {
459 workDV->size = neqns ;
460 ZV_writeForMatlab((ZV *) workDV, "W", msgFile) ;
461 workDV->size = 2*neqns ;
462 }
463 fflush(msgFile) ;
464 }
465 /*
466 -----------------------
467 copy work vector into X
468 -----------------------
469 */
470 if ( work != xent ) {
471 if ( ILUMTX_IS_REAL(mtx) ) {
472 DVcopy(neqns, xent, work) ;
473 } else {
474 DVcopy(2*neqns, xent, work) ;
475 }
476 }
477 /*
478 -----------------------------------------
479 increment the operation count if not NULL
480 -----------------------------------------
481 */
482 if ( pops != NULL ) {
483 *pops += ops ;
484 }
485 return(1) ; }
486
487 /*--------------------------------------------------------------------*/
488