"Fossies" - the Fresh Open Source Software Archive 
Member "laspack/itersolv.c" (6 Apr 1995, 43357 Bytes) of package /linux/privat/old/laspack.tgz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
1 /****************************************************************************/
2 /* itersolv.c */
3 /****************************************************************************/
4 /* */
5 /* ITERative SOLVers for systems of linear equations */
6 /* */
7 /* Copyright (C) 1992-1995 Tomas Skalicky. All rights reserved. */
8 /* */
9 /****************************************************************************/
10 /* */
11 /* ANY USE OF THIS CODE CONSTITUTES ACCEPTANCE OF THE TERMS */
12 /* OF THE COPYRIGHT NOTICE (SEE FILE COPYRGHT.H) */
13 /* */
14 /****************************************************************************/
15
16 #include <stddef.h>
17 #include <math.h>
18
19 #include "laspack/itersolv.h"
20 #include "laspack/elcmp.h"
21 #include "laspack/errhandl.h"
22 #include "laspack/operats.h"
23 #include "laspack/rtc.h"
24 #include "laspack/copyrght.h"
25
26 /* number of GMRES steps bevore restart */
27 static int GMRESSteps = 10;
28
29 Vector *JacobiIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
30 PrecondProcType Dummy, double Omega)
31 {
32 int Iter;
33 double bNorm;
34 size_t Dim;
35 Vector r;
36
37 Q_Lock(A);
38 V_Lock(x);
39 V_Lock(b);
40
41 Dim = Q_GetDim(A);
42 V_Constr(&r, "r", Dim, Normal, True);
43
44 if (LASResult() == LASOK) {
45 bNorm = l2Norm_V(b);
46
47 Iter = 0;
48 /* r = b - A * x(i) */
49 if (!IsZero(l1Norm_V(x) / Dim)) {
50 if (Q_KerDefined(A))
51 OrthoRightKer_VQ(x, A);
52 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
53 } else {
54 Asgn_VV(&r, b);
55 }
56 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, JacobiIterId)
57 && Iter < MaxIter) {
58 Iter++;
59 /* x(i+1) = x(i) + Omega * D^(-1) * r */
60 AddAsgn_VV(x, Mul_SV(Omega, MulInv_QV(Diag_Q(A), &r)));
61 if (Q_KerDefined(A))
62 OrthoRightKer_VQ(x, A);
63 /* r = b - A * x(i) */
64 if (Iter < MaxIter)
65 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
66 }
67 }
68
69 V_Destr(&r);
70
71 Q_Unlock(A);
72 V_Unlock(x);
73 V_Unlock(b);
74
75 return(x);
76 }
77
78 Vector *SORForwIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
79 PrecondProcType Dummy, double Omega)
80 {
81 int Iter;
82 double bNorm;
83 size_t Dim;
84 Vector r;
85
86 Q_Lock(A);
87 V_Lock(x);
88 V_Lock(b);
89
90 Dim = Q_GetDim(A);
91 V_Constr(&r, "r", Dim, Normal, True);
92
93 if (LASResult() == LASOK) {
94 bNorm = l2Norm_V(b);
95
96 Iter = 0;
97 /* r = b - A * x(i) */
98 if (!IsZero(l1Norm_V(x) / Dim)) {
99 if (Q_KerDefined(A))
100 OrthoRightKer_VQ(x, A);
101 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
102 } else {
103 Asgn_VV(&r, b);
104 }
105 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SORForwIterId)
106 && Iter < MaxIter) {
107 Iter++;
108 /* x(i+1) = x(i) + (D / Omega + L)^(-1) r */
109 AddAsgn_VV(x, MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A)), &r));
110 if (Q_KerDefined(A))
111 OrthoRightKer_VQ(x, A);
112 /* r = b - A * x(i) */
113 if (Iter < MaxIter)
114 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
115 }
116 }
117
118 V_Destr(&r);
119
120 Q_Unlock(A);
121 V_Unlock(x);
122 V_Unlock(b);
123
124 return(x);
125 }
126
127 Vector *SORBackwIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
128 PrecondProcType Dummy, double Omega)
129 {
130 int Iter;
131 double bNorm;
132 size_t Dim;
133 Vector r;
134
135 Q_Lock(A);
136 V_Lock(x);
137 V_Lock(b);
138
139 Dim = Q_GetDim(A);
140 V_Constr(&r, "r", Dim, Normal, True);
141
142 if (LASResult() == LASOK) {
143 bNorm = l2Norm_V(b);
144
145 Iter = 0;
146 /* r = b - A * x(i) */
147 if (!IsZero(l1Norm_V(x) / Dim)) {
148 if (Q_KerDefined(A))
149 OrthoRightKer_VQ(x, A);
150 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
151 } else {
152 Asgn_VV(&r, b);
153 }
154 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SORBackwIterId)
155 && Iter < MaxIter) {
156 Iter++;
157 /* x(i+1) = x(i) + (D / Omega + U)^(-1) r */
158 AddAsgn_VV(x, MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Upper_Q(A)), &r));
159 if (Q_KerDefined(A))
160 OrthoRightKer_VQ(x, A);
161 /* r = b - A * x(i) */
162 if (Iter < MaxIter)
163 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
164 }
165 }
166
167 V_Destr(&r);
168
169 Q_Unlock(A);
170 V_Unlock(x);
171 V_Unlock(b);
172
173 return(x);
174 }
175
176 Vector *SSORIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
177 PrecondProcType Dummy, double Omega)
178 {
179 int Iter;
180 double bNorm;
181 size_t Dim;
182 Vector r;
183
184 Q_Lock(A);
185 V_Lock(x);
186 V_Lock(b);
187
188 Dim = Q_GetDim(A);
189 V_Constr(&r, "r", Dim, Normal, True);
190
191 if (LASResult() == LASOK) {
192 bNorm = l2Norm_V(b);
193
194 Iter = 0;
195 /* r = b - A * x(i) */
196 if (!IsZero(l1Norm_V(x) / Dim)) {
197 if (Q_KerDefined(A))
198 OrthoRightKer_VQ(x, A);
199 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
200 } else {
201 Asgn_VV(&r, b);
202 }
203 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, SSORIterId)
204 && Iter < MaxIter) {
205 Iter++;
206 /* x(i+1) = x(i) + (2 - Omega) * (Diag(A) / Omega + Upper(A))^(-1)
207 * Diag(A) / Omega * (Diag(A) / Omega + Lower(A))^(-1) r */
208 AddAsgn_VV(x, Mul_SV((2.0 - Omega) / Omega,
209 MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Upper_Q(A)),
210 Mul_QV(Diag_Q(A),
211 MulInv_QV(Add_QQ(Mul_SQ(1.0 / Omega, Diag_Q(A)), Lower_Q(A)), &r)))));
212 if (Q_KerDefined(A))
213 OrthoRightKer_VQ(x, A);
214 /* r = b - A * x(i) */
215 if (Iter < MaxIter)
216 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
217 }
218 }
219
220 V_Destr(&r);
221
222 Q_Unlock(A);
223 V_Unlock(x);
224 V_Unlock(b);
225
226 return(x);
227 }
228
229 Vector *ChebyshevIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
230 PrecondProcType PrecondProc, double OmegaPrecond)
231 {
232 /*
233 * for details to the algorithm see
234 *
235 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
236 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
237 * Templates for the Solution of Linear Systems: Building Blocks
238 * for Iterative Solvers;
239 * SIAM, Philadelphia, 1994
240 *
241 * The original algorithm doesn't seem to work correctly.
242 * The three term recursion was therefore in the curent version
243 * adopted after
244 *
245 * W. Hackbusch:
246 * Iterative Solution of Large Sparse Systems of Equations,
247 * Springer-Verlag, Berlin, 1994
248 *
249 */
250
251 int Iter;
252 double MaxEigenval, MinEigenval;
253 double c, d, Alpha, Beta;
254 double T = 0.0, TOld = 0.0, TNew; /* values of Chebyshev polynomials */
255 double bNorm;
256 size_t Dim;
257 Vector r, p, z;
258
259
260 Q_Lock(A);
261 V_Lock(x);
262 V_Lock(b);
263
264 Dim = Q_GetDim(A);
265 V_Constr(&r, "r", Dim, Normal, True);
266 V_Constr(&p, "p", Dim, Normal, True);
267 if (PrecondProc != NULL || Q_KerDefined(A))
268 V_Constr(&z, "z", Dim, Normal, True);
269
270 if (LASResult() == LASOK) {
271 bNorm = l2Norm_V(b);
272
273 MaxEigenval = GetMaxEigenval(A, PrecondProc, OmegaPrecond);
274 MinEigenval = GetMinEigenval(A, PrecondProc, OmegaPrecond);
275
276 c = (MaxEigenval - MinEigenval) / 2.0;
277 d = (MaxEigenval + MinEigenval) / 2.0;
278
279 Iter = 0;
280 /* r = b - A * x(i) */
281 if (!IsZero(l1Norm_V(x) / Dim)) {
282 if (Q_KerDefined(A))
283 OrthoRightKer_VQ(x, A);
284 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
285 } else {
286 Asgn_VV(&r, b);
287 }
288 if (PrecondProc != NULL || Q_KerDefined(A)) {
289 /* preconditioned Chebyshev method */
290 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, ChebyshevIterId)
291 && Iter < MaxIter) {
292 Iter++;
293 if (PrecondProc != NULL)
294 (*PrecondProc)(A, &z, &r, OmegaPrecond);
295 else
296 Asgn_VV(&z, &r);
297 if (Q_KerDefined(A))
298 OrthoRightKer_VQ(&z, A);
299 if (Iter == 1) {
300 TOld = 1.0;
301 T = d / c;
302 Alpha = 1.0 / d;
303 Asgn_VV(&p, Mul_SV(Alpha, &z));
304 } else {
305 TNew = 2.0 * d / c * T - TOld;
306 TOld = T;
307 T = TNew;
308 Alpha = 2.0 / c * TOld / T;
309 Beta = 2.0 * d / c * TOld / T - 1.0;
310 Asgn_VV(&p, Add_VV(Mul_SV(Alpha, &z), Mul_SV(Beta, &p)));
311 }
312 AddAsgn_VV(x, &p);
313 if (Iter < MaxIter)
314 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
315 }
316 } else {
317 /* plain Chebyshev method (z = r) */
318 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, ChebyshevIterId)
319 && Iter < MaxIter) {
320 Iter++;
321 if (Iter == 1) {
322 TOld = 1.0;
323 T = d / c;
324 Alpha = 1.0 / d;
325 Asgn_VV(&p, Mul_SV(Alpha, &r));
326 } else {
327 TNew = 2.0 * d / c * T - TOld;
328 TOld = T;
329 T = TNew;
330 Alpha = 2.0 / c * TOld / T;
331 Beta = 2.0 * d / c * TOld / T - 1.0;
332 Asgn_VV(&p, Add_VV(Mul_SV(Alpha, &r), Mul_SV(Beta, &p)));
333 }
334 AddAsgn_VV(x, &p);
335 if (Iter < MaxIter)
336 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
337 }
338 }
339 if (Q_KerDefined(A))
340 OrthoRightKer_VQ(x, A);
341 }
342
343 V_Destr(&r);
344 V_Destr(&p);
345 if (PrecondProc != NULL || Q_KerDefined(A))
346 V_Destr(&z);
347
348 Q_Unlock(A);
349 V_Unlock(x);
350 V_Unlock(b);
351
352 return(x);
353 }
354
355 Vector *CGIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
356 PrecondProcType PrecondProc, double OmegaPrecond)
357 {
358 /*
359 * for details to the algorithm see
360 *
361 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
362 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
363 * Templates for the Solution of Linear Systems: Building Blocks
364 * for Iterative Solvers;
365 * SIAM, Philadelphia, 1994
366 *
367 */
368
369 int Iter;
370 double Alpha, Beta, Rho, RhoOld = 0.0;
371 double bNorm;
372 size_t Dim;
373 Vector r, p, q, z;
374
375 Q_Lock(A);
376 V_Lock(x);
377 V_Lock(b);
378
379 Dim = Q_GetDim(A);
380 V_Constr(&r, "r", Dim, Normal, True);
381 V_Constr(&p, "p", Dim, Normal, True);
382 V_Constr(&q, "q", Dim, Normal, True);
383 if (PrecondProc != NULL || Q_KerDefined(A))
384 V_Constr(&z, "z", Dim, Normal, True);
385
386 if (LASResult() == LASOK) {
387 bNorm = l2Norm_V(b);
388
389 Iter = 0;
390 /* r = b - A * x(i) */
391 if (!IsZero(l1Norm_V(x) / Dim)) {
392 if (Q_KerDefined(A))
393 OrthoRightKer_VQ(x, A);
394 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
395 } else {
396 Asgn_VV(&r, b);
397 }
398 if (PrecondProc != NULL || Q_KerDefined(A)) {
399 /* preconditioned CG */
400 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGIterId)
401 && Iter < MaxIter) {
402 Iter++;
403 if (PrecondProc != NULL)
404 (*PrecondProc)(A, &z, &r, OmegaPrecond);
405 else
406 Asgn_VV(&z, &r);
407 if (Q_KerDefined(A))
408 OrthoRightKer_VQ(&z, A);
409 Rho = Mul_VV(&r, &z);
410 if (Iter == 1) {
411 Asgn_VV(&p, &z);
412 } else {
413 Beta = Rho / RhoOld;
414 Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
415 }
416 Asgn_VV(&q, Mul_QV(A, &p));
417 Alpha = Rho / Mul_VV(&p, &q);
418 AddAsgn_VV(x, Mul_SV(Alpha, &p));
419 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
420 RhoOld = Rho;
421 }
422 } else {
423 /* plain CG (z = r) */
424 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGIterId)
425 && Iter < MaxIter) {
426 Iter++;
427 Rho = pow(l2Norm_V(&r), 2.0);
428 if (Iter == 1) {
429 Asgn_VV(&p, &r);
430 } else {
431 Beta = Rho / RhoOld;
432 Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, &p)));
433 }
434 Asgn_VV(&q, Mul_QV(A, &p));
435 Alpha = Rho / Mul_VV(&p, &q);
436 AddAsgn_VV(x, Mul_SV(Alpha, &p));
437 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
438 RhoOld = Rho;
439 }
440 }
441 if (Q_KerDefined(A))
442 OrthoRightKer_VQ(x, A);
443 }
444
445 V_Destr(&r);
446 V_Destr(&p);
447 V_Destr(&q);
448 if (PrecondProc != NULL || Q_KerDefined(A))
449 V_Destr(&z);
450
451 Q_Unlock(A);
452 V_Unlock(x);
453 V_Unlock(b);
454
455 return(x);
456 }
457
458 Vector *CGNIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
459 PrecondProcType PrecondProc, double OmegaPrecond)
460 {
461 int Iter;
462 double Alpha, Beta, Rho, RhoOld = 0.0;
463 double bNorm;
464 size_t Dim;
465 Vector r, p, q, s, z;
466
467 Q_Lock(A);
468 V_Lock(x);
469 V_Lock(b);
470
471 Dim = Q_GetDim(A);
472 V_Constr(&r, "r", Dim, Normal, True);
473 V_Constr(&p, "p", Dim, Normal, True);
474 V_Constr(&q, "q", Dim, Normal, True);
475 V_Constr(&z, "z", Dim, Normal, True);
476 if (PrecondProc != NULL || Q_KerDefined(A))
477 V_Constr(&s, "s", Dim, Normal, True);
478
479 if (LASResult() == LASOK) {
480 bNorm = l2Norm_V(b);
481
482 Iter = 0;
483 /* r = b - A * x(i) */
484 if (!IsZero(l1Norm_V(x) / Dim)) {
485 if (Q_KerDefined(A))
486 OrthoRightKer_VQ(x, A);
487 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
488 } else {
489 Asgn_VV(&r, b);
490 }
491 if (PrecondProc != NULL || Q_KerDefined(A)) {
492 /* preconditioned CGN */
493 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGNIterId)
494 && Iter < MaxIter) {
495 Iter++;
496 if (PrecondProc != NULL) {
497 (*PrecondProc)(A, &z, &r, OmegaPrecond);
498 (*PrecondProc)(Transp_Q(A), &q, &z, OmegaPrecond);
499 Asgn_VV(&z, Mul_QV(Transp_Q(A), &q));
500 } else {
501 Asgn_VV(&z, &r);
502 }
503 if (Q_KerDefined(A))
504 OrthoRightKer_VQ(&z, A);
505 Rho = pow(l2Norm_V(&z), 2.0);
506 if (Iter == 1) {
507 Asgn_VV(&p, &z);
508 } else {
509 Beta = Rho / RhoOld;
510 Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
511 }
512 Asgn_VV(&q, Mul_QV(A, &p));
513 if (PrecondProc != NULL)
514 (*PrecondProc)(A, &s, &q, OmegaPrecond);
515 else
516 Asgn_VV(&s, &q);
517 if (Q_KerDefined(A))
518 OrthoRightKer_VQ(&s, A);
519 Alpha = Rho / pow(l2Norm_V(&s), 2.0);
520 AddAsgn_VV(x, Mul_SV(Alpha, &p));
521 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
522 RhoOld = Rho;
523 }
524 } else {
525 /* plain CGN */
526 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGNIterId)
527 && Iter < MaxIter) {
528 Iter++;
529 Asgn_VV(&z, Mul_QV(Transp_Q(A), &r));
530 Rho = pow(l2Norm_V(&z), 2.0);
531 if (Iter == 1) {
532 Asgn_VV(&p, &z);
533 } else {
534 Beta = Rho / RhoOld;
535 Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
536 }
537 Asgn_VV(&q, Mul_QV(A, &p));
538 Alpha = Rho / pow(l2Norm_V(&q), 2.0);
539 AddAsgn_VV(x, Mul_SV(Alpha, &p));
540 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
541 RhoOld = Rho;
542 }
543 }
544 if (Q_KerDefined(A))
545 OrthoRightKer_VQ(x, A);
546 }
547
548 V_Destr(&r);
549 V_Destr(&p);
550 V_Destr(&q);
551 V_Destr(&z);
552 if (PrecondProc != NULL || Q_KerDefined(A))
553 V_Destr(&s);
554
555 Q_Unlock(A);
556 V_Unlock(x);
557 V_Unlock(b);
558
559 return(x);
560 }
561
562 Vector *GMRESIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
563 PrecondProcType PrecondProc, double OmegaPrecond)
564 {
565 /*
566 * for details to the algorithm see
567 *
568 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
569 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
570 * Templates for the Solution of Linear Systems: Building Blocks
571 * for Iterative Solvers;
572 * SIAM, Philadelphia, 1994
573 *
574 */
575
576 char vName[10];
577 int Iter, i, j, k;
578 double h1, h2, r;
579 double bNorm;
580 double **h, *y, *s, *c1, *c2;
581 size_t Dim;
582 Boolean AllocOK;
583 Vector *v;
584
585 Q_Lock(A);
586 V_Lock(x);
587 V_Lock(b);
588
589 /* allocation of matrix H and vectors y, s, c1 and c2 */
590 AllocOK = True;
591 h = (double **)malloc((GMRESSteps + 1) * sizeof(double *));
592 if (h == NULL) {
593 AllocOK = False;
594 } else {
595 for (i = 1; i <= GMRESSteps; i++) {
596 h[i] = (double *)malloc((GMRESSteps + 2) * sizeof(double));
597 if (h[i] == NULL)
598 AllocOK = False;
599 }
600 }
601 y = (double *)malloc((GMRESSteps + 1) * sizeof(double));
602 if (y == NULL)
603 AllocOK = False;
604 s = (double *)malloc((GMRESSteps + 2) * sizeof(double));
605 if (s == NULL)
606 AllocOK = False;
607 c1 = (double *)malloc((GMRESSteps + 1) * sizeof(double));
608 if (c1 == NULL)
609 AllocOK = False;
610 c2 = (double *)malloc((GMRESSteps + 1) * sizeof(double));
611 if (c2 == NULL)
612 AllocOK = False;
613
614 /* ... and vectors u */
615 Dim = Q_GetDim(A);
616 v = (Vector *)malloc((GMRESSteps + 2) * sizeof(Vector));
617 if (v == NULL)
618 AllocOK = False;
619 else
620 for (i = 1; i <= GMRESSteps + 1; i++) {
621 sprintf(vName, "v[%d]", i % 1000);
622 V_Constr(&v[i], vName, Dim, Normal, True);
623 }
624
625 if (!AllocOK)
626 LASError(LASMemAllocErr, "GMRESIter", NULL, NULL, NULL);
627
628 if (LASResult() == LASOK) {
629 bNorm = l2Norm_V(b);
630
631 /* loop for 'MaxIter' GMRES cycles */
632 Iter = 0;
633 /* v[1] = r = b - A * x(i) */
634 if (!IsZero(l1Norm_V(x) / Dim)) {
635 if (Q_KerDefined(A))
636 OrthoRightKer_VQ(x, A);
637 Asgn_VV(&v[1], Sub_VV(b, Mul_QV(A, x)));
638 } else {
639 Asgn_VV(&v[1], b);
640 }
641 while (!RTCResult(Iter, l2Norm_V(&v[1]), bNorm, GMRESIterId)
642 && Iter < MaxIter) {
643 if (PrecondProc != NULL)
644 (*PrecondProc)(A, &v[1], &v[1], OmegaPrecond);
645 s[1] = l2Norm_V(&v[1]);
646 MulAsgn_VS(&v[1], 1.0 / s[1]);
647
648 /* GMRES iteration */
649 i = 0;
650 while ((PrecondProc != NULL ? True : !RTCResult(Iter, fabs(s[i+1]),
651 bNorm, GMRESIterId)) && i < GMRESSteps && Iter < MaxIter) {
652 i++;
653 Iter++;
654 /* w = v[i+1] */
655 if (PrecondProc != NULL)
656 (*PrecondProc)(A, &v[i+1], Mul_QV(A, &v[i]), OmegaPrecond);
657 else
658 Asgn_VV(&v[i+1], Mul_QV(A, &v[i]));
659
660 /* modified Gram-Schmidt orthogonalization */
661 for (k = 1; k <= i; k++) {
662 h[i][k] = Mul_VV(&v[i+1], &v[k]);
663 SubAsgn_VV(&v[i+1], Mul_SV(h[i][k], &v[k]));
664 }
665
666 h[i][i+1] = l2Norm_V(&v[i+1]);
667 MulAsgn_VS(&v[i+1], 1.0 / h[i][i+1]);
668
669 /* Q-R algorithm */
670 for (k = 1; k < i; k++) {
671 h1 = c1[k] * h[i][k] + c2[k] * h[i][k+1];
672 h2 = - c2[k] * h[i][k] + c1[k] * h[i][k+1];
673 h[i][k] = h1;
674 h[i][k+1] = h2;
675 }
676 r = sqrt(h[i][i] * h[i][i] + h[i][i+1] * h[i][i+1]);
677 c1[i] = h[i][i] / r;
678 c2[i] = h[i][i+1] / r;
679 h[i][i] = r;
680 h[i][i+1] = 0.0;
681 s[i+1] = - c2[i] * s[i];
682 s[i] = c1[i] * s[i];
683 }
684
685 /* Solving of the system of equations : H y = s */
686 for (j = i; j > 0; j--) {
687 y[j] = s[j] / h[j][j];
688 for (k = j - 1; k > 0; k--)
689 s[k] -= h[j][k] * y[j];
690 }
691
692 /* updating solution */
693 for (j = i; j > 0; j--)
694 AddAsgn_VV(x, Mul_SV(y[j], &v[j]));
695 if (Q_KerDefined(A))
696 OrthoRightKer_VQ(x, A);
697
698 /* computing new residual */
699 Asgn_VV(&v[1], Sub_VV(b, Mul_QV(A, x)));
700 }
701 }
702
703 /* release of vectors u, matrix H and vectors y, s, c1 and c2 */
704 if (v != NULL) {
705 for (i = 1; i <= GMRESSteps + 1; i++)
706 V_Destr(&v[i]);
707 free(v);
708 }
709
710 if (h != NULL) {
711 for (i = 1; i <= GMRESSteps; i++)
712 if (h[i] != NULL)
713 free(h[i]);
714 free(h);
715 }
716 if (y != NULL)
717 free(y);
718 if (s != NULL)
719 free(s);
720 if (c1 != NULL)
721 free(c1);
722 if (c2 != NULL)
723 free(c2);
724
725 Q_Unlock(A);
726 V_Unlock(x);
727 V_Unlock(b);
728
729 return(x);
730 }
731
732 Vector *BiCGIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
733 PrecondProcType PrecondProc, double OmegaPrecond)
734 {
735 /*
736 * for details to the algorithm see
737 *
738 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
739 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
740 * Templates for the Solution of Linear Systems: Building Blocks
741 * for Iterative Solvers;
742 * SIAM, Philadelphia, 1994
743 *
744 */
745
746 int Iter;
747 double Alpha, Beta, Rho, RhoOld = 0.0;
748 double bNorm;
749 size_t Dim;
750 Vector r, r_, p, p_, q, z, z_;
751
752 Q_Lock(A);
753 V_Lock(x);
754 V_Lock(b);
755
756 Dim = Q_GetDim(A);
757 V_Constr(&r, "r", Dim, Normal, True);
758 V_Constr(&r_, "r_", Dim, Normal, True);
759 V_Constr(&p, "p", Dim, Normal, True);
760 V_Constr(&p_, "p_", Dim, Normal, True);
761 V_Constr(&q, "q", Dim, Normal, True);
762 if (PrecondProc != NULL || Q_KerDefined(A)) {
763 V_Constr(&z, "z", Dim, Normal, True);
764 V_Constr(&z_, "z_", Dim, Normal, True);
765 }
766
767 if (LASResult() == LASOK) {
768 bNorm = l2Norm_V(b);
769
770 Iter = 0;
771 /* r = b - A * x(i) */
772 if (!IsZero(l1Norm_V(x) / Dim)) {
773 if (Q_KerDefined(A))
774 OrthoRightKer_VQ(x, A);
775 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
776 } else {
777 Asgn_VV(&r, b);
778 }
779 if (PrecondProc != NULL || Q_KerDefined(A)) {
780 /* preconditioned BiCG */
781 Asgn_VV(&r_, &r);
782 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGIterId)
783 && Iter < MaxIter) {
784 Iter++;
785 if (PrecondProc != NULL)
786 (*PrecondProc)(A, &z, &r, OmegaPrecond);
787 else
788 Asgn_VV(&z, &r);
789 if (Q_KerDefined(A))
790 OrthoRightKer_VQ(&z, A);
791 if (PrecondProc != NULL)
792 (*PrecondProc)(Transp_Q(A), &z_, &r_, OmegaPrecond);
793 else
794 Asgn_VV(&z_, &r_);
795 if (Q_KerDefined(A))
796 OrthoRightKer_VQ(&z_, Transp_Q(A));
797 Rho = Mul_VV(&z, &r_);
798 if (IsZero(Rho)){
799 LASError(LASBreakdownErr, "BiCGIter", "Rho", NULL, NULL);
800 break;
801 }
802 if (Iter == 1) {
803 Asgn_VV(&p, &z);
804 Asgn_VV(&p_, &z_);
805 } else {
806 Beta = Rho / RhoOld;
807 Asgn_VV(&p, Add_VV(&z, Mul_SV(Beta, &p)));
808 Asgn_VV(&p_, Add_VV(&z_, Mul_SV(Beta, &p_)));
809 }
810 Asgn_VV(&q, Mul_QV(A, &p));
811 Alpha = Rho / Mul_VV(&p_, &q);
812 AddAsgn_VV(x, Mul_SV(Alpha, &p));
813 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
814 SubAsgn_VV(&r_, Mul_SV(Alpha, Mul_QV(Transp_Q(A), &p_)));
815 RhoOld = Rho;
816 }
817 } else {
818 /* plain BiCG (z = r, z_ = r_) */
819 Asgn_VV(&r_, &r);
820 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGIterId)
821 && Iter < MaxIter) {
822 Iter++;
823 Rho = Mul_VV(&r, &r_);
824 if (IsZero(Rho)) {
825 LASError(LASBreakdownErr, "BiCGIter", "Rho", NULL, NULL);
826 break;
827 }
828 if (Iter == 1) {
829 Asgn_VV(&p, &r);
830 Asgn_VV(&p_, &r_);
831 } else {
832 Beta = Rho / RhoOld;
833 Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, &p)));
834 Asgn_VV(&p_, Add_VV(&r_, Mul_SV(Beta, &p_)));
835 }
836 Asgn_VV(&q, Mul_QV(A, &p));
837 Alpha = Rho / Mul_VV(&p_, &q);
838 AddAsgn_VV(x, Mul_SV(Alpha, &p));
839 SubAsgn_VV(&r, Mul_SV(Alpha, &q));
840 SubAsgn_VV(&r_, Mul_SV(Alpha, Mul_QV(Transp_Q(A), &p_)));
841 RhoOld = Rho;
842 }
843 }
844 if (Q_KerDefined(A))
845 OrthoRightKer_VQ(x, A);
846 }
847
848 V_Destr(&r);
849 V_Destr(&r_);
850 V_Destr(&p);
851 V_Destr(&p_);
852 V_Destr(&q);
853 if (PrecondProc != NULL || Q_KerDefined(A)) {
854 V_Destr(&z);
855 V_Destr(&z_);
856 }
857
858 Q_Unlock(A);
859 V_Unlock(x);
860 V_Unlock(b);
861
862 return(x);
863 }
864
865 Vector *QMRIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
866 PrecondProcType PrecondProc, double OmegaPrecond)
867 {
868 /*
869 * for details to the algorithm see
870 *
871 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
872 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
873 * Templates for the Solution of Linear Systems: Building Blocks
874 * for Iterative Solvers;
875 * SIAM, Philadelphia, 1994
876 *
877 */
878
879 int Iter;
880 double Rho, RhoNew, Xi, Gamma, GammaOld, Eta, Delta, Beta, Epsilon = 0.0,
881 Theta, ThetaOld = 0.0;
882 double bNorm;
883 size_t Dim;
884 Vector r, p, p_, q, y, y_, v, v_, w, w_, z, z_, s, d;
885
886 Q_Lock(A);
887 V_Lock(x);
888 V_Lock(b);
889
890 Dim = Q_GetDim(A);
891 V_Constr(&r, "r", Dim, Normal, True);
892 V_Constr(&p, "p", Dim, Normal, True);
893 V_Constr(&p_, "p_", Dim, Normal, True);
894 V_Constr(&q, "q", Dim, Normal, True);
895 V_Constr(&v, "v", Dim, Normal, True);
896 V_Constr(&v_, "v_", Dim, Normal, True);
897 V_Constr(&w, "w", Dim, Normal, True);
898 V_Constr(&w_, "w_", Dim, Normal, True);
899 V_Constr(&s, "s", Dim, Normal, True);
900 V_Constr(&d, "d", Dim, Normal, True);
901 if (PrecondProc != NULL || Q_KerDefined(A)) {
902 V_Constr(&y, "y", Dim, Normal, True);
903 V_Constr(&y_, "y_", Dim, Normal, True);
904 V_Constr(&z, "z", Dim, Normal, True);
905 V_Constr(&z_, "z_", Dim, Normal, True);
906 }
907
908 if (LASResult() == LASOK) {
909 bNorm = l2Norm_V(b);
910
911 Iter = 0;
912 /* r = b - A * x(i) */
913 if (!IsZero(l1Norm_V(x) / Dim)) {
914 if (Q_KerDefined(A))
915 OrthoRightKer_VQ(x, A);
916 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
917 } else {
918 Asgn_VV(&r, b);
919 }
920 if (PrecondProc != NULL || Q_KerDefined(A)) {
921 /* preconditioned QMR (M1 ~ A, M2 = I) */
922 Asgn_VV(&v_, &r);
923 if (PrecondProc != NULL)
924 (*PrecondProc)(A, &y, &v_, OmegaPrecond);
925 else
926 Asgn_VV(&y, &v_);
927 if (Q_KerDefined(A))
928 OrthoRightKer_VQ(&y, A);
929 RhoNew = l2Norm_V(&y);
930 Asgn_VV(&w_, &r);
931 Asgn_VV(&z, &w_);
932 Xi = l2Norm_V(&z);
933 Gamma = 1.0;
934 Eta = - 1.0;
935 GammaOld = Gamma;
936 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, QMRIterId)
937 && Iter < MaxIter) {
938 Iter++;
939 Rho = RhoNew;
940 if (IsZero(Rho) || IsZero(Xi)) {
941 LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
942 break;
943 }
944 Asgn_VV(&v, Mul_SV(1.0 / Rho, &v_));
945 MulAsgn_VS(&y, 1.0 / Rho);
946 Asgn_VV(&w, Mul_SV(1.0 / Xi, &w_));
947 MulAsgn_VS(&z, 1.0 / Xi);
948 Delta = Mul_VV(&z, &y);
949 if (IsZero(Delta)) {
950 LASError(LASBreakdownErr, "QMRIter", "Delta", NULL, NULL);
951 break;
952 }
953 Asgn_VV(&y_, &y);
954 if (PrecondProc != NULL)
955 (*PrecondProc)(Transp_Q(A), &z_, &z, OmegaPrecond);
956 else
957 Asgn_VV(&z_, &z);
958 if (Q_KerDefined(A))
959 OrthoRightKer_VQ(&z_, Transp_Q(A));
960 if (Iter == 1) {
961 Asgn_VV(&p, &y_);
962 Asgn_VV(&q, &z_);
963 } else {
964 Asgn_VV(&p, Sub_VV(&y_, Mul_SV(Xi * Delta / Epsilon, &p)));
965 Asgn_VV(&q, Sub_VV(&z_, Mul_SV(Rho * Delta / Epsilon, &q)));
966 }
967 Asgn_VV(&p_, Mul_QV(A, &p));
968 Epsilon = Mul_VV(&q, &p_);
969 if (IsZero(Epsilon)) {
970 LASError(LASBreakdownErr, "QMRIter", "Epsilon", NULL, NULL);
971 break;
972 }
973 Beta = Epsilon / Delta;
974 if (IsZero(Beta)) {
975 LASError(LASBreakdownErr, "QMRIter", "Beta", NULL, NULL);
976 break;
977 }
978 Asgn_VV(&v_, Sub_VV(&p_, Mul_SV(Beta, &v)));
979 if (PrecondProc != NULL)
980 (*PrecondProc)(A, &y, &v_, OmegaPrecond);
981 else
982 Asgn_VV(&y, &v_);
983 if (Q_KerDefined(A))
984 OrthoRightKer_VQ(&y, A);
985 RhoNew = l2Norm_V(&y);
986 Asgn_VV(&w_, Sub_VV(Mul_QV(Transp_Q(A), &q), Mul_SV(Beta, &w)));
987 Asgn_VV(&z, &w_);
988 Xi = l2Norm_V(&z);
989 Theta = RhoNew / (GammaOld * fabs(Beta));
990 Gamma = 1.0 / sqrt(1.0 + Theta * Theta);
991 if (IsZero(Gamma)) {
992 LASError(LASBreakdownErr, "QMRIter", "Gamma", NULL, NULL);
993 break;
994 }
995 Eta = - Eta * Rho * Gamma * Gamma / (Beta * GammaOld * GammaOld);
996 if (Iter == 1) {
997 Asgn_VV(&d, Mul_SV(Eta, &p));
998 Asgn_VV(&s, Mul_SV(Eta, &p_));
999 } else {
1000 Asgn_VV(&d, Add_VV(Mul_SV(Eta, &p),
1001 Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &d)));
1002 Asgn_VV(&s, Add_VV(Mul_SV(Eta, &p_),
1003 Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &s)));
1004 }
1005 AddAsgn_VV(x, &d);
1006 SubAsgn_VV(&r, &s);
1007 GammaOld = Gamma;
1008 ThetaOld = Theta;
1009 }
1010 } else {
1011 /* plain QMR (M1 ~ A, M2 = I, y = y_ = v_, z = z_ = w_) */
1012 Asgn_VV(&v_, &r);
1013 RhoNew = l2Norm_V(&v_);
1014 Asgn_VV(&w_, &r);
1015 Xi = l2Norm_V(&w_);
1016 Gamma = 1.0;
1017 Eta = - 1.0;
1018 GammaOld = Gamma;
1019 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, QMRIterId)
1020 && Iter < MaxIter) {
1021 Iter++;
1022 Rho = RhoNew;
1023 if (IsZero(Rho) || IsZero(Xi)) {
1024 LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
1025 break;
1026 }
1027 Asgn_VV(&v, Mul_SV(1.0 / Rho, &v_));
1028 MulAsgn_VS(&v_, 1.0 / Rho);
1029 Asgn_VV(&w, Mul_SV(1.0 / Xi, &w_));
1030 MulAsgn_VS(&w_, 1.0 / Xi);
1031 Delta = Mul_VV(&w_, &v_);
1032 if (IsZero(Delta)) {
1033 LASError(LASBreakdownErr, "QMRIter", "Rho", "Xi", NULL);
1034 break;
1035 }
1036 if (Iter == 1) {
1037 Asgn_VV(&p, &v_);
1038 Asgn_VV(&q, &w_);
1039 } else {
1040 Asgn_VV(&p, Sub_VV(&v_, Mul_SV(Xi * Delta / Epsilon, &p)));
1041 Asgn_VV(&q, Sub_VV(&w_, Mul_SV(Rho * Delta / Epsilon, &q)));
1042 }
1043 Asgn_VV(&p_, Mul_QV(A, &p));
1044 Epsilon = Mul_VV(&q, &p_);
1045 if (IsZero(Epsilon)) {
1046 LASError(LASBreakdownErr, "QMRIter", "Epsilon", NULL, NULL);
1047 break;
1048 }
1049 Beta = Epsilon / Delta;
1050 if (IsZero(Beta)) {
1051 LASError(LASBreakdownErr, "QMRIter", "Beta", NULL, NULL);
1052 break;
1053 }
1054 Asgn_VV(&v_, Sub_VV(&p_, Mul_SV(Beta, &v)));
1055 RhoNew = l2Norm_V(&v);
1056 Asgn_VV(&w_, Sub_VV(Mul_QV(Transp_Q(A), &q), Mul_SV(Beta, &w)));
1057 Xi = l2Norm_V(&w_);
1058 Theta = RhoNew / (GammaOld * fabs(Beta));
1059 Gamma = 1.0 / sqrt(1.0 + Theta * Theta);
1060 if (IsZero(Gamma)) {
1061 LASError(LASBreakdownErr, "QMRIter", "Gamma", NULL, NULL);
1062 break;
1063 }
1064 Eta = - Eta * Rho * Gamma * Gamma / (Beta * GammaOld * GammaOld);
1065 if (Iter == 1) {
1066 Asgn_VV(&d, Mul_SV(Eta, &p));
1067 Asgn_VV(&s, Mul_SV(Eta, &p_));
1068 } else {
1069 Asgn_VV(&d, Add_VV(Mul_SV(Eta, &p),
1070 Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &d)));
1071 Asgn_VV(&s, Add_VV(Mul_SV(Eta, &p_),
1072 Mul_SV(ThetaOld * ThetaOld * Gamma * Gamma, &s)));
1073 }
1074 AddAsgn_VV(x, &d);
1075 SubAsgn_VV(&r, &s);
1076 GammaOld = Gamma;
1077 ThetaOld = Theta;
1078 }
1079 }
1080 if (Q_KerDefined(A))
1081 OrthoRightKer_VQ(x, A);
1082 }
1083
1084 V_Destr(&r);
1085 V_Destr(&p);
1086 V_Destr(&p_);
1087 V_Destr(&q);
1088 V_Destr(&v);
1089 V_Destr(&v_);
1090 V_Destr(&w);
1091 V_Destr(&w_);
1092 V_Destr(&s);
1093 V_Destr(&d);
1094 if (PrecondProc != NULL || Q_KerDefined(A)) {
1095 V_Destr(&y);
1096 V_Destr(&y_);
1097 V_Destr(&z);
1098 V_Destr(&z_);
1099 }
1100
1101 Q_Unlock(A);
1102 V_Unlock(x);
1103 V_Unlock(b);
1104
1105 return(x);
1106 }
1107
1108 Vector *CGSIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
1109 PrecondProcType PrecondProc, double OmegaPrecond)
1110 {
1111 /*
1112 * for details to the algorithm see
1113 *
1114 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
1115 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
1116 * Templates for the Solution of Linear Systems: Building Blocks
1117 * for Iterative Solvers;
1118 * SIAM, Philadelphia, 1994
1119 *
1120 */
1121
1122 int Iter;
1123 double Alpha, Beta, Rho, RhoOld = 0.0;
1124 double bNorm;
1125 size_t Dim;
1126 Vector r, r_, p, p_, q, u, u_, v_;
1127
1128 Q_Lock(A);
1129 V_Lock(x);
1130 V_Lock(b);
1131
1132 Dim = Q_GetDim(A);
1133 V_Constr(&r, "r", Dim, Normal, True);
1134 V_Constr(&r_, "r_", Dim, Normal, True);
1135 V_Constr(&p, "p", Dim, Normal, True);
1136 V_Constr(&q, "q", Dim, Normal, True);
1137 V_Constr(&u, "u", Dim, Normal, True);
1138 V_Constr(&u_, "u_", Dim, Normal, True);
1139 V_Constr(&v_, "v_", Dim, Normal, True);
1140 if (PrecondProc != NULL || Q_KerDefined(A))
1141 V_Constr(&p_, "p_", Dim, Normal, True);
1142
1143 if (LASResult() == LASOK) {
1144 bNorm = l2Norm_V(b);
1145
1146 Iter = 0;
1147 /* r = b - A * x(i) */
1148 if (!IsZero(l1Norm_V(x) / Dim)) {
1149 if (Q_KerDefined(A))
1150 OrthoRightKer_VQ(x, A);
1151 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
1152 } else {
1153 Asgn_VV(&r, b);
1154 }
1155 if (PrecondProc != NULL || Q_KerDefined(A)) {
1156 /* preconditioned CGS */
1157 Asgn_VV(&r_, &r);
1158 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGSIterId)
1159 && Iter < MaxIter) {
1160 Iter++;
1161 Rho = Mul_VV(&r_, &r);
1162 if (IsZero(Rho)) {
1163 LASError(LASBreakdownErr, "CGSIter", "Rho", NULL, NULL);
1164 break;
1165 }
1166 if (Iter == 1) {
1167 Asgn_VV(&u, &r);
1168 Asgn_VV(&p, &u);
1169 } else {
1170 Beta = Rho / RhoOld;
1171 Asgn_VV(&u, Add_VV(&r, Mul_SV(Beta, &q)));
1172 Asgn_VV(&p, Add_VV(&u, Mul_SV(Beta, Add_VV(&q, Mul_SV(Beta, &p)))));
1173 }
1174 if (PrecondProc != NULL)
1175 (*PrecondProc)(A, &p_, &p, OmegaPrecond);
1176 else
1177 Asgn_VV(&p_, &p);
1178 if (Q_KerDefined(A))
1179 OrthoRightKer_VQ(&p_, A);
1180 Asgn_VV(&v_, Mul_QV(A, &p_));
1181 Alpha = Rho / Mul_VV(&r_, &v_);
1182 Asgn_VV(&q, Sub_VV(&u, Mul_SV(Alpha, &v_)));
1183 if (PrecondProc != NULL)
1184 (*PrecondProc)(A, &u_, Add_VV(&u, &q), OmegaPrecond);
1185 else
1186 Asgn_VV(&u_, Add_VV(&u, &q));
1187 if (Q_KerDefined(A))
1188 OrthoRightKer_VQ(&u_, A);
1189 AddAsgn_VV(x, Mul_SV(Alpha, &u_));
1190 SubAsgn_VV(&r, Mul_SV(Alpha, Mul_QV(A, &u_)));
1191 RhoOld = Rho;
1192 }
1193 } else {
1194 /* plain CGS (p_ = p) */
1195 Asgn_VV(&r_, &r);
1196 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, CGSIterId)
1197 && Iter < MaxIter) {
1198 Iter++;
1199 Rho = Mul_VV(&r_, &r);
1200 if (IsZero(Rho)) {
1201 LASError(LASBreakdownErr, "CGSIter", "Rho", NULL, NULL);
1202 break;
1203 }
1204 if (Iter == 1) {
1205 Asgn_VV(&u, &r);
1206 Asgn_VV(&p, &u);
1207 } else {
1208 Beta = Rho / RhoOld;
1209 Asgn_VV(&u, Add_VV(&r, Mul_SV(Beta, &q)));
1210 Asgn_VV(&p, Add_VV(&u, Mul_SV(Beta, Add_VV(&q, Mul_SV(Beta, &p)))));
1211 }
1212 Asgn_VV(&v_, Mul_QV(A, &p));
1213 Alpha = Rho / Mul_VV(&r_, &v_);
1214 Asgn_VV(&q, Sub_VV(&u, Mul_SV(Alpha, &v_)));
1215 Asgn_VV(&u_, Add_VV(&u, &q));
1216 AddAsgn_VV(x, Mul_SV(Alpha, &u_));
1217 SubAsgn_VV(&r, Mul_SV(Alpha, Mul_QV(A, &u_)));
1218 RhoOld = Rho;
1219 }
1220 }
1221 if (Q_KerDefined(A))
1222 OrthoRightKer_VQ(x, A);
1223 }
1224
1225 V_Destr(&r);
1226 V_Destr(&r_);
1227 V_Destr(&p);
1228 V_Destr(&q);
1229 V_Destr(&u);
1230 V_Destr(&u_);
1231 V_Destr(&v_);
1232 if (PrecondProc != NULL || Q_KerDefined(A))
1233 V_Destr(&p_);
1234
1235 Q_Unlock(A);
1236 V_Unlock(x);
1237 V_Unlock(b);
1238
1239 return(x);
1240 }
1241
1242 Vector *BiCGSTABIter(QMatrix *A, Vector *x, Vector *b, int MaxIter,
1243 PrecondProcType PrecondProc, double OmegaPrecond)
1244 {
1245 /*
1246 * for details to the algorithm see
1247 *
1248 * R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
1249 * V. Eijkhout, R. Pozo, Ch. Romine, H. van der Vorst:
1250 * Templates for the Solution of Linear Systems: Building Blocks
1251 * for Iterative Solvers;
1252 * SIAM, Philadelphia, 1994
1253 *
1254 */
1255
1256 int Iter;
1257 double Alpha = 0.0, Beta, Rho, RhoOld = 0.0, Omega = 0.0;
1258 double bNorm;
1259 size_t Dim;
1260 Vector r, r_, p, p_, v, s, s_, t;
1261
1262 Q_Lock(A);
1263 V_Lock(x);
1264 V_Lock(b);
1265
1266 Dim = Q_GetDim(A);
1267 V_Constr(&r, "r", Dim, Normal, True);
1268 V_Constr(&r_, "r_", Dim, Normal, True);
1269 V_Constr(&p, "p", Dim, Normal, True);
1270 V_Constr(&v, "v", Dim, Normal, True);
1271 V_Constr(&s, "s", Dim, Normal, True);
1272 V_Constr(&t, "t", Dim, Normal, True);
1273 if (PrecondProc != NULL || Q_KerDefined(A)) {
1274 V_Constr(&p_, "p_", Dim, Normal, True);
1275 V_Constr(&s_, "s_", Dim, Normal, True);
1276 }
1277
1278 if (LASResult() == LASOK) {
1279 bNorm = l2Norm_V(b);
1280
1281 Iter = 0;
1282 /* r = b - A * x(i) */
1283 if (!IsZero(l1Norm_V(x) / Dim)) {
1284 if (Q_KerDefined(A))
1285 OrthoRightKer_VQ(x, A);
1286 Asgn_VV(&r, Sub_VV(b, Mul_QV(A, x)));
1287 } else {
1288 Asgn_VV(&r, b);
1289 }
1290 if (PrecondProc != NULL || Q_KerDefined(A)) {
1291 /* preconditioned BiCGStab */
1292 Asgn_VV(&r_, &r);
1293 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGSTABIterId)
1294 && Iter < MaxIter) {
1295 Iter++;
1296 Rho = Mul_VV(&r_, &r);
1297 if (IsZero(Rho)) {
1298 LASError(LASBreakdownErr, "BiCGSTABIter", "Rho", NULL, NULL);
1299 break;
1300 }
1301 if (Iter == 1) {
1302 Asgn_VV(&p, &r);
1303 } else {
1304 Beta = (Rho / RhoOld) * (Alpha / Omega);
1305 Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, Sub_VV(&p, Mul_SV(Omega, &v)))));
1306 }
1307 if (PrecondProc != NULL)
1308 (*PrecondProc)(A, &p_, &p, OmegaPrecond);
1309 else
1310 Asgn_VV(&p_, &p);
1311 if (Q_KerDefined(A))
1312 OrthoRightKer_VQ(&p_, A);
1313 Asgn_VV(&v, Mul_QV(A, &p_));
1314 Alpha = Rho / Mul_VV(&r_, &v);
1315 Asgn_VV(&s, Sub_VV(&r, Mul_SV(Alpha, &v)));
1316 if (PrecondProc != NULL)
1317 (*PrecondProc)(A, &s_, &s, OmegaPrecond);
1318 else
1319 Asgn_VV(&s_, &s);
1320 if (Q_KerDefined(A))
1321 OrthoRightKer_VQ(&s_, A);
1322 Asgn_VV(&t, Mul_QV(A, &s_));
1323 Omega = Mul_VV(&t, &s) / pow(l2Norm_V(&t), 2.0);
1324 AddAsgn_VV(x, Add_VV(Mul_SV(Alpha, &p_), Mul_SV(Omega, &s_)));
1325 Asgn_VV(&r, Sub_VV(&s, Mul_SV(Omega, &t)));
1326 RhoOld = Rho;
1327 if (IsZero(Omega))
1328 break;
1329 }
1330 } else {
1331 /* plain BiCGStab (p_ = p) */
1332 Asgn_VV(&r_, &r);
1333 while (!RTCResult(Iter, l2Norm_V(&r), bNorm, BiCGSTABIterId)
1334 && Iter < MaxIter) {
1335 Iter++;
1336 Rho = Mul_VV(&r_, &r);
1337 if (IsZero(Rho)) {
1338 LASError(LASBreakdownErr, "BiCGSTABIter", "Rho", NULL, NULL);
1339 break;
1340 }
1341 if (Iter == 1) {
1342 Asgn_VV(&p, &r);
1343 } else {
1344 Beta = (Rho / RhoOld) * (Alpha / Omega);
1345 Asgn_VV(&p, Add_VV(&r, Mul_SV(Beta, Sub_VV(&p, Mul_SV(Omega, &v)))));
1346 }
1347 Asgn_VV(&v, Mul_QV(A, &p));
1348 Alpha = Rho / Mul_VV(&r_, &v);
1349 Asgn_VV(&s, Sub_VV(&r, Mul_SV(Alpha, &v)));
1350 Asgn_VV(&t, Mul_QV(A, &s));
1351 Omega = Mul_VV(&t, &s) / pow(l2Norm_V(&t), 2.0);
1352 AddAsgn_VV(x, Add_VV(Mul_SV(Alpha, &p), Mul_SV(Omega, &s)));
1353 Asgn_VV(&r, Sub_VV(&s, Mul_SV(Omega, &t)));
1354 RhoOld = Rho;
1355 if (IsZero(Omega))
1356 break;
1357 }
1358 }
1359 if (Q_KerDefined(A))
1360 OrthoRightKer_VQ(x, A);
1361 }
1362
1363 V_Destr(&r);
1364 V_Destr(&r_);
1365 V_Destr(&p);
1366 V_Destr(&v);
1367 V_Destr(&s);
1368 V_Destr(&t);
1369 if (PrecondProc != NULL || Q_KerDefined(A)) {
1370 V_Destr(&p_);
1371 V_Destr(&s_);
1372 }
1373
1374 Q_Unlock(A);
1375 V_Unlock(x);
1376 V_Unlock(b);
1377
1378 return(x);
1379 }
1380
1381 void SetGMRESRestart(int Steps)
1382 {
1383 GMRESSteps = Steps;
1384 }