Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KryPgmres.c
Go to the documentation of this file.
1
24#include "fasp.h"
25#include "fasp_functs.h"
26#include <math.h>
27
28/*---------------------------------*/
29/*-- Declare Private Functions --*/
30/*---------------------------------*/
31
32#include "KryUtil.inl"
33
34/*---------------------------------*/
35/*-- Public Functions --*/
36/*---------------------------------*/
37
66INT fasp_solver_dcsr_pgmres(dCSRmat* A, dvector* b, dvector* x, precond* pc,
67 const REAL tol, const REAL abstol, const INT MaxIt,
68 const SHORT restart, const SHORT StopType,
69 const SHORT PrtLvl)
70{
71 const INT n = b->row;
72 const INT MIN_ITER = 0;
73
74 // local variables
75 INT iter = 0;
76 int i, j, k; // must be signed! -zcs
77
78 REAL r_norm, r_normb, gamma, t;
79 REAL absres0 = BIGREAL, absres = BIGREAL;
80 REAL relres = BIGREAL, normu = BIGREAL;
81
82 // allocate temp memory (need about (restart+4)*n REAL numbers)
83 REAL * c = NULL, *s = NULL, *rs = NULL;
84 REAL * norms = NULL, *r = NULL, *w = NULL;
85 REAL* work = NULL;
86 REAL **p = NULL, **hh = NULL;
87
88 INT Restart = MIN(restart, MaxIt);
89 INT Restart1 = Restart + 1;
90 LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
91
92 /* allocate memory and setup temp work space */
93 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
94
95 // Output some info for debugging
96 if (PrtLvl > PRINT_NONE) printf("\nCalling GMRes solver (CSR) ...\n");
97
98#if DEBUG_MODE > 0
99 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
100 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
101#endif
102
103 /* check whether memory is enough for GMRES */
104 while ((work == NULL) && (Restart > 5)) {
105 Restart = Restart - 5;
106 Restart1 = Restart + 1;
107 worksize = (Restart + 4) * (Restart + n) + 1 - n;
108 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
109 }
110
111 if (work == NULL) {
112 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
113 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
114 }
115
116 if (PrtLvl > PRINT_MIN && Restart < restart) {
117 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
118 }
119
120 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
121 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
122 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
123
124 r = work;
125 w = r + n;
126 rs = w + n;
127 c = rs + Restart1;
128 s = c + Restart;
129
130 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
131
132 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
133
134 // compute initial residual: r = b-A*x
135 fasp_darray_cp(n, b->val, p[0]);
136 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
137 r_norm = fasp_blas_darray_norm2(n, p[0]);
138
139 // compute stopping criteria
140 switch (StopType) {
141 case STOP_REL_RES:
142 absres0 = MAX(SMALLREAL, r_norm);
143 relres = r_norm / absres0;
144 break;
145 case STOP_REL_PRECRES:
146 if (pc == NULL)
147 fasp_darray_cp(n, p[0], r);
148 else
149 pc->fct(p[0], r, pc->data);
150 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
151 absres0 = MAX(SMALLREAL, r_normb);
152 relres = r_normb / absres0;
153 break;
154 case STOP_MOD_REL_RES:
155 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
156 absres0 = r_norm;
157 relres = absres0 / normu;
158 break;
159 default:
160 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
161 goto FINISHED;
162 }
163
164 // if initial residual is small, no need to iterate!
165 if (relres < tol || absres0 < abstol) goto FINISHED;
166
167 // output iteration information if needed
168 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0.0);
169
170 // store initial residual
171 norms[0] = relres;
172
173 /* GMRES(M) outer iteration */
174 while (iter < MaxIt && relres > tol) {
175
176 rs[0] = r_norm;
177
178 t = 1.0 / r_norm;
179
180 fasp_blas_darray_ax(n, t, p[0]);
181
182 /* RESTART CYCLE (right-preconditioning) */
183 i = 0;
184 while (i < Restart && iter < MaxIt) {
185
186 i++;
187 iter++;
188
189 /* apply preconditioner */
190 if (pc == NULL)
191 fasp_darray_cp(n, p[i - 1], r);
192 else
193 pc->fct(p[i - 1], r, pc->data);
194
195 fasp_blas_dcsr_mxv(A, r, p[i]);
196
197 /* Modified Gram_Schmidt orthogonalization */
198 for (j = 0; j < i; j++) {
199 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
200 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
201 }
202 t = fasp_blas_darray_norm2(n, p[i]);
203 hh[i][i - 1] = t;
204
205 if (ABS(t) > SMALLREAL) { // If t=0, we get solution subspace
206 t = 1.0 / t;
207 fasp_blas_darray_ax(n, t, p[i]);
208 }
209
210 for (j = 1; j < i; ++j) {
211 t = hh[j - 1][i - 1];
212 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
213 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
214 }
215 t = hh[i][i - 1] * hh[i][i - 1];
216 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
217
218 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
219 c[i - 1] = hh[i - 1][i - 1] / gamma;
220 s[i - 1] = hh[i][i - 1] / gamma;
221 rs[i] = -s[i - 1] * rs[i - 1];
222 rs[i - 1] = c[i - 1] * rs[i - 1];
223 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
224
225 absres = r_norm = fabs(rs[i]);
226
227 relres = absres / absres0;
228
229 norms[iter] = relres;
230
231 // output iteration information if needed
232 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
233 norms[iter] / norms[iter - 1]);
234
235 // exit restart cycle if reaches tolerance
236 if (relres < tol && iter >= MIN_ITER) break;
237
238 } /* end of restart cycle */
239
240 /* compute solution, first solve upper triangular system */
241 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
242 for (k = i - 2; k >= 0; k--) {
243 t = 0.0;
244 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
245 t += rs[k];
246 rs[k] = t / hh[k][k];
247 }
248
249 fasp_darray_cp(n, p[i - 1], w);
250
251 fasp_blas_darray_ax(n, rs[i - 1], w);
252
253 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
254
255 /* apply preconditioner */
256 if (pc == NULL)
257 fasp_darray_cp(n, w, r);
258 else
259 pc->fct(w, r, pc->data);
260
261 fasp_blas_darray_axpy(n, 1.0, r, x->val);
262
263 // Check: prevent false convergence
264 if (relres < tol && iter >= MIN_ITER) {
265
266 REAL computed_relres = relres;
267
268 // compute residual
269 fasp_darray_cp(n, b->val, r);
270 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
271 r_norm = fasp_blas_darray_norm2(n, r);
272
273 switch (StopType) {
274 case STOP_REL_RES:
275 absres = r_norm;
276 relres = absres / absres0;
277 break;
278 case STOP_REL_PRECRES:
279 if (pc == NULL)
280 fasp_darray_cp(n, r, w);
281 else
282 pc->fct(r, w, pc->data);
283 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
284 relres = absres / absres0;
285 break;
286 case STOP_MOD_REL_RES:
287 absres = r_norm;
288 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
289 relres = absres / normu;
290 break;
291 }
292
293 norms[iter] = relres;
294
295 if (relres < tol) {
296 break;
297 } else { // Need to restart
298 fasp_darray_cp(n, r, p[0]);
299 i = 0;
300 }
301
302 if (PrtLvl >= PRINT_MORE) {
303 ITS_COMPRES(computed_relres);
304 ITS_REALRES(relres);
305 }
306
307 } /* end of convergence check */
308
309 /* compute residual vector and continue loop */
310 for (j = i; j > 0; j--) {
311 rs[j - 1] = -s[j - 1] * rs[j];
312 rs[j] = c[j - 1] * rs[j];
313 }
314
315 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
316
317 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
318
319 if (i) {
320 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
321 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
322 }
323
324 } /* end of main while loop */
325FINISHED:
326 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
327
328 /*-------------------------------------------
329 * Clean up workspace
330 *------------------------------------------*/
331 fasp_mem_free(work);
332 work = NULL;
333 fasp_mem_free(p);
334 p = NULL;
335 fasp_mem_free(hh);
336 hh = NULL;
337 fasp_mem_free(norms);
338 norms = NULL;
339
340#if DEBUG_MODE > 0
341 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
342#endif
343
344 if (iter >= MaxIt)
345 return ERROR_SOLVER_MAXIT;
346 else
347 return iter;
348}
349
376INT fasp_solver_dbsr_pgmres(dBSRmat* A, dvector* b, dvector* x, precond* pc,
377 const REAL tol, const REAL abstol, const INT MaxIt,
378 const SHORT restart, const SHORT StopType,
379 const SHORT PrtLvl)
380{
381 const INT n = b->row;
382 const INT MIN_ITER = 0;
383
384 // local variables
385 INT iter = 0;
386 int i, j, k; // must be signed! -zcs
387
388 REAL r_norm, r_normb, gamma, t;
389 REAL absres0 = BIGREAL, absres = BIGREAL;
390 REAL relres = BIGREAL, normu = BIGREAL;
391
392 // allocate temp memory (need about (restart+4)*n REAL numbers)
393 REAL * c = NULL, *s = NULL, *rs = NULL;
394 REAL * norms = NULL, *r = NULL, *w = NULL;
395 REAL* work = NULL;
396 REAL **p = NULL, **hh = NULL;
397
398 INT Restart = MIN(restart, MaxIt);
399 INT Restart1 = Restart + 1;
400 LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
401
402 /* allocate memory and setup temp work space */
403 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
404
405 // Output some info for debugging
406 if (PrtLvl > PRINT_NONE) printf("\nCalling GMRes solver (BSR) ...\n");
407
408#if DEBUG_MODE > 0
409 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
410 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
411#endif
412
413 /* check whether memory is enough for GMRES */
414 while ((work == NULL) && (Restart > 5)) {
415 Restart = Restart - 5;
416 Restart1 = Restart + 1;
417 worksize = (Restart + 4) * (Restart + n) + 1 - n;
418 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
419 }
420
421 if (work == NULL) {
422 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
423 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
424 }
425
426 if (PrtLvl > PRINT_MIN && Restart < restart) {
427 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
428 }
429
430 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
431 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
432 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
433
434 r = work;
435 w = r + n;
436 rs = w + n;
437 c = rs + Restart1;
438 s = c + Restart;
439
440 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
441
442 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
443
444 // compute initial residual: r = b-A*x
445 fasp_darray_cp(n, b->val, p[0]);
446 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
447 r_norm = fasp_blas_darray_norm2(n, p[0]);
448
449 // compute stopping criteria
450 switch (StopType) {
451 case STOP_REL_RES:
452 absres0 = MAX(SMALLREAL, r_norm);
453 relres = r_norm / absres0;
454 break;
455 case STOP_REL_PRECRES:
456 if (pc == NULL)
457 fasp_darray_cp(n, p[0], r);
458 else
459 pc->fct(p[0], r, pc->data);
460 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
461 absres0 = MAX(SMALLREAL, r_normb);
462 relres = r_normb / absres0;
463 break;
464 case STOP_MOD_REL_RES:
465 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
466 absres0 = r_norm;
467 relres = absres0 / normu;
468 break;
469 default:
470 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
471 goto FINISHED;
472 }
473
474 // if initial residual is small, no need to iterate!
475 if (relres < tol || absres0 < abstol) goto FINISHED;
476
477 // output iteration information if needed
478 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0.0);
479
480 // store initial residual
481 norms[0] = relres;
482
483 /* GMRES(M) outer iteration */
484 while (iter < MaxIt && relres > tol) {
485
486 rs[0] = r_norm;
487
488 t = 1.0 / r_norm;
489
490 fasp_blas_darray_ax(n, t, p[0]);
491
492 /* RESTART CYCLE (right-preconditioning) */
493 i = 0;
494 while (i < Restart && iter < MaxIt) {
495
496 i++;
497 iter++;
498
499 /* apply preconditioner */
500 if (pc == NULL)
501 fasp_darray_cp(n, p[i - 1], r);
502 else
503 pc->fct(p[i - 1], r, pc->data);
504
505 fasp_blas_dbsr_mxv(A, r, p[i]);
506
507 /* Modified Gram_Schmidt orthogonalization */
508 for (j = 0; j < i; j++) {
509 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
510 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
511 }
512 t = fasp_blas_darray_norm2(n, p[i]);
513 hh[i][i - 1] = t;
514
515 if (ABS(t) > SMALLREAL) { // If t=0, we get solution subspace
516 t = 1.0 / t;
517 fasp_blas_darray_ax(n, t, p[i]);
518 }
519
520 for (j = 1; j < i; ++j) {
521 t = hh[j - 1][i - 1];
522 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
523 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
524 }
525 t = hh[i][i - 1] * hh[i][i - 1];
526 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
527
528 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
529 c[i - 1] = hh[i - 1][i - 1] / gamma;
530 s[i - 1] = hh[i][i - 1] / gamma;
531 rs[i] = -s[i - 1] * rs[i - 1];
532 rs[i - 1] = c[i - 1] * rs[i - 1];
533 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
534
535 absres = r_norm = fabs(rs[i]);
536
537 relres = absres / absres0;
538
539 norms[iter] = relres;
540
541 // output iteration information if needed
542 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
543 norms[iter] / norms[iter - 1]);
544
545 // exit restart cycle if reaches tolerance
546 if (relres < tol && iter >= MIN_ITER) break;
547
548 } /* end of restart cycle */
549
550 /* compute solution, first solve upper triangular system */
551 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
552 for (k = i - 2; k >= 0; k--) {
553 t = 0.0;
554 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
555 t += rs[k];
556 rs[k] = t / hh[k][k];
557 }
558
559 fasp_darray_cp(n, p[i - 1], w);
560
561 fasp_blas_darray_ax(n, rs[i - 1], w);
562
563 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
564
565 /* apply preconditioner */
566 if (pc == NULL)
567 fasp_darray_cp(n, w, r);
568 else
569 pc->fct(w, r, pc->data);
570
571 fasp_blas_darray_axpy(n, 1.0, r, x->val);
572
573 // Check: prevent false convergence
574 if (relres < tol && iter >= MIN_ITER) {
575
576 REAL computed_relres = relres;
577
578 // compute residual
579 fasp_darray_cp(n, b->val, r);
580 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
581 r_norm = fasp_blas_darray_norm2(n, r);
582
583 switch (StopType) {
584 case STOP_REL_RES:
585 absres = r_norm;
586 relres = absres / absres0;
587 break;
588 case STOP_REL_PRECRES:
589 if (pc == NULL)
590 fasp_darray_cp(n, r, w);
591 else
592 pc->fct(r, w, pc->data);
593 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
594 relres = absres / absres0;
595 break;
596 case STOP_MOD_REL_RES:
597 absres = r_norm;
598 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
599 relres = absres / normu;
600 break;
601 }
602
603 norms[iter] = relres;
604
605 if (relres < tol) {
606 break;
607 } else { // Need to restart
608 fasp_darray_cp(n, r, p[0]);
609 i = 0;
610 }
611
612 if (PrtLvl >= PRINT_MORE) {
613 ITS_COMPRES(computed_relres);
614 ITS_REALRES(relres);
615 }
616
617 } /* end of convergence check */
618
619 /* compute residual vector and continue loop */
620 for (j = i; j > 0; j--) {
621 rs[j - 1] = -s[j - 1] * rs[j];
622 rs[j] = c[j - 1] * rs[j];
623 }
624
625 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
626
627 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
628
629 if (i) {
630 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
631 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
632 }
633
634 } /* end of main while loop */
635
636FINISHED:
637 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
638
639 /*-------------------------------------------
640 * Clean up workspace
641 *------------------------------------------*/
642 fasp_mem_free(work);
643 work = NULL;
644 fasp_mem_free(p);
645 p = NULL;
646 fasp_mem_free(hh);
647 hh = NULL;
648 fasp_mem_free(norms);
649 norms = NULL;
650
651#if DEBUG_MODE > 0
652 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
653#endif
654
655 if (iter >= MaxIt)
656 return ERROR_SOLVER_MAXIT;
657 else
658 return iter;
659}
660
687INT fasp_solver_dblc_pgmres(dBLCmat* A, dvector* b, dvector* x, precond* pc,
688 const REAL tol, const REAL abstol, const INT MaxIt,
689 const SHORT restart, const SHORT StopType,
690 const SHORT PrtLvl)
691{
692 const INT n = b->row;
693 const INT MIN_ITER = 0;
694
695 // local variables
696 INT iter = 0;
697 int i, j, k; // must be signed! -zcs
698
699 REAL r_norm, r_normb, gamma, t;
700 REAL absres0 = BIGREAL, absres = BIGREAL;
701 REAL relres = BIGREAL, normu = BIGREAL;
702
703 // allocate temp memory (need about (restart+4)*n REAL numbers)
704 REAL * c = NULL, *s = NULL, *rs = NULL;
705 REAL * norms = NULL, *r = NULL, *w = NULL;
706 REAL* work = NULL;
707 REAL **p = NULL, **hh = NULL;
708
709 INT Restart = MIN(restart, MaxIt);
710 INT Restart1 = Restart + 1;
711 LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
712
713 /* allocate memory and setup temp work space */
714 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
715
716 // Output some info for debugging
717 if (PrtLvl > PRINT_NONE) printf("\nCalling GMRes solver (BLC) ...\n");
718
719#if DEBUG_MODE > 0
720 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
721 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
722#endif
723
724 /* check whether memory is enough for GMRES */
725 while ((work == NULL) && (Restart > 5)) {
726 Restart = Restart - 5;
727 Restart1 = Restart + 1;
728 worksize = (Restart + 4) * (Restart + n) + 1 - n;
729 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
730 }
731
732 if (work == NULL) {
733 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
734 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
735 }
736
737 if (PrtLvl > PRINT_MIN && Restart < restart) {
738 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
739 }
740
741 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
742 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
743 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
744
745 r = work;
746 w = r + n;
747 rs = w + n;
748 c = rs + Restart1;
749 s = c + Restart;
750
751 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
752
753 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
754
755 // compute initial residual: r = b-A*x
756 fasp_darray_cp(n, b->val, p[0]);
757 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
758 r_norm = fasp_blas_darray_norm2(n, p[0]);
759
760 // compute stopping criteria
761 switch (StopType) {
762 case STOP_REL_RES:
763 absres0 = MAX(SMALLREAL, r_norm);
764 relres = r_norm / absres0;
765 break;
766 case STOP_REL_PRECRES:
767 if (pc == NULL)
768 fasp_darray_cp(n, p[0], r);
769 else
770 pc->fct(p[0], r, pc->data);
771 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
772 absres0 = MAX(SMALLREAL, r_normb);
773 relres = r_normb / absres0;
774 break;
775 case STOP_MOD_REL_RES:
776 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
777 absres0 = r_norm;
778 relres = absres0 / normu;
779 break;
780 default:
781 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
782 goto FINISHED;
783 }
784
785 // if initial residual is small, no need to iterate!
786 if (relres < tol || absres0 < abstol) goto FINISHED;
787
788 // output iteration information if needed
789 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0.0);
790
791 // store initial residual
792 norms[0] = relres;
793
794 /* GMRES(M) outer iteration */
795 while (iter < MaxIt && relres > tol) {
796
797 rs[0] = r_norm;
798
799 t = 1.0 / r_norm;
800
801 fasp_blas_darray_ax(n, t, p[0]);
802
803 /* RESTART CYCLE (right-preconditioning) */
804 i = 0;
805 while (i < Restart && iter < MaxIt) {
806
807 i++;
808 iter++;
809
810 /* apply preconditioner */
811 if (pc == NULL)
812 fasp_darray_cp(n, p[i - 1], r);
813 else
814 pc->fct(p[i - 1], r, pc->data);
815
816 fasp_blas_dblc_mxv(A, r, p[i]);
817
818 /* Modified Gram_Schmidt orthogonalization */
819 for (j = 0; j < i; j++) {
820 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
821 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
822 }
823 t = fasp_blas_darray_norm2(n, p[i]);
824 hh[i][i - 1] = t;
825
826 if (ABS(t) > SMALLREAL) { // If t=0, we get solution subspace
827 t = 1.0 / t;
828 fasp_blas_darray_ax(n, t, p[i]);
829 }
830
831 for (j = 1; j < i; ++j) {
832 t = hh[j - 1][i - 1];
833 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
834 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
835 }
836 t = hh[i][i - 1] * hh[i][i - 1];
837 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
838
839 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
840 c[i - 1] = hh[i - 1][i - 1] / gamma;
841 s[i - 1] = hh[i][i - 1] / gamma;
842 rs[i] = -s[i - 1] * rs[i - 1];
843 rs[i - 1] = c[i - 1] * rs[i - 1];
844 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
845
846 absres = r_norm = fabs(rs[i]);
847
848 relres = absres / absres0;
849
850 norms[iter] = relres;
851
852 // output iteration information if needed
853 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
854 norms[iter] / norms[iter - 1]);
855
856 // exit restart cycle if reaches tolerance
857 if (relres < tol && iter >= MIN_ITER) break;
858
859 } /* end of restart cycle */
860
861 /* compute solution, first solve upper triangular system */
862 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
863 for (k = i - 2; k >= 0; k--) {
864 t = 0.0;
865 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
866 t += rs[k];
867 rs[k] = t / hh[k][k];
868 }
869
870 fasp_darray_cp(n, p[i - 1], w);
871
872 fasp_blas_darray_ax(n, rs[i - 1], w);
873
874 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
875
876 /* apply preconditioner */
877 if (pc == NULL)
878 fasp_darray_cp(n, w, r);
879 else
880 pc->fct(w, r, pc->data);
881
882 fasp_blas_darray_axpy(n, 1.0, r, x->val);
883
884 // Check: prevent false convergence
885 if (relres < tol && iter >= MIN_ITER) {
886
887 REAL computed_relres = relres;
888
889 // compute residual
890 fasp_darray_cp(n, b->val, r);
891 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
892 r_norm = fasp_blas_darray_norm2(n, r);
893
894 switch (StopType) {
895 case STOP_REL_RES:
896 absres = r_norm;
897 relres = absres / absres0;
898 break;
899 case STOP_REL_PRECRES:
900 if (pc == NULL)
901 fasp_darray_cp(n, r, w);
902 else
903 pc->fct(r, w, pc->data);
904 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
905 relres = absres / absres0;
906 break;
907 case STOP_MOD_REL_RES:
908 absres = r_norm;
909 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
910 relres = absres / normu;
911 break;
912 }
913
914 norms[iter] = relres;
915
916 if (relres < tol) {
917 break;
918 } else { // Need to restart
919 fasp_darray_cp(n, r, p[0]);
920 i = 0;
921 }
922
923 if (PrtLvl >= PRINT_MORE) {
924 ITS_COMPRES(computed_relres);
925 ITS_REALRES(relres);
926 }
927
928 } /* end of convergence check */
929
930 /* compute residual vector and continue loop */
931 for (j = i; j > 0; j--) {
932 rs[j - 1] = -s[j - 1] * rs[j];
933 rs[j] = c[j - 1] * rs[j];
934 }
935
936 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
937
938 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
939
940 if (i) {
941 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
942 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
943 }
944
945 } /* end of main while loop */
946
947FINISHED:
948 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
949
950 /*-------------------------------------------
951 * Clean up workspace
952 *------------------------------------------*/
953 fasp_mem_free(work);
954 work = NULL;
955 fasp_mem_free(p);
956 p = NULL;
957 fasp_mem_free(hh);
958 hh = NULL;
959 fasp_mem_free(norms);
960 norms = NULL;
961
962#if DEBUG_MODE > 0
963 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
964#endif
965
966 if (iter >= MaxIt)
967 return ERROR_SOLVER_MAXIT;
968 else
969 return iter;
970}
971
998INT fasp_solver_dstr_pgmres(dSTRmat* A, dvector* b, dvector* x, precond* pc,
999 const REAL tol, const REAL abstol, const INT MaxIt,
1000 const SHORT restart, const SHORT StopType,
1001 const SHORT PrtLvl)
1002{
1003 const INT n = b->row;
1004 const INT MIN_ITER = 0;
1005
1006 // local variables
1007 INT iter = 0;
1008 int i, j, k; // must be signed! -zcs
1009
1010 REAL r_norm, r_normb, gamma, t;
1011 REAL absres0 = BIGREAL, absres = BIGREAL;
1012 REAL relres = BIGREAL, normu = BIGREAL;
1013
1014 // allocate temp memory (need about (restart+4)*n REAL numbers)
1015 REAL * c = NULL, *s = NULL, *rs = NULL;
1016 REAL * norms = NULL, *r = NULL, *w = NULL;
1017 REAL* work = NULL;
1018 REAL **p = NULL, **hh = NULL;
1019
1020 INT Restart = MIN(restart, MaxIt);
1021 INT Restart1 = Restart + 1;
1022 LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
1023
1024 /* allocate memory and setup temp work space */
1025 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1026
1027 // Output some info for debugging
1028 if (PrtLvl > PRINT_NONE) printf("\nCalling GMRes solver (STR) ...\n");
1029
1030#if DEBUG_MODE > 0
1031 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1032 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1033#endif
1034
1035 /* check whether memory is enough for GMRES */
1036 while ((work == NULL) && (Restart > 5)) {
1037 Restart = Restart - 5;
1038 Restart1 = Restart + 1;
1039 worksize = (Restart + 4) * (Restart + n) + 1 - n;
1040 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1041 }
1042
1043 if (work == NULL) {
1044 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1045 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1046 }
1047
1048 if (PrtLvl > PRINT_MIN && Restart < restart) {
1049 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
1050 }
1051
1052 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1053 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1054 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
1055
1056 r = work;
1057 w = r + n;
1058 rs = w + n;
1059 c = rs + Restart1;
1060 s = c + Restart;
1061
1062 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
1063
1064 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
1065
1066 // compute initial residual: r = b-A*x
1067 fasp_darray_cp(n, b->val, p[0]);
1068 fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1069 r_norm = fasp_blas_darray_norm2(n, p[0]);
1070
1071 // compute stopping criteria
1072 switch (StopType) {
1073 case STOP_REL_RES:
1074 absres0 = MAX(SMALLREAL, r_norm);
1075 relres = r_norm / absres0;
1076 break;
1077 case STOP_REL_PRECRES:
1078 if (pc == NULL)
1079 fasp_darray_cp(n, p[0], r);
1080 else
1081 pc->fct(p[0], r, pc->data);
1082 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
1083 absres0 = MAX(SMALLREAL, r_normb);
1084 relres = r_normb / absres0;
1085 break;
1086 case STOP_MOD_REL_RES:
1087 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
1088 absres0 = r_norm;
1089 relres = absres0 / normu;
1090 break;
1091 default:
1092 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1093 goto FINISHED;
1094 }
1095
1096 // if initial residual is small, no need to iterate!
1097 if (relres < tol || absres0 < abstol) goto FINISHED;
1098
1099 // output iteration information if needed
1100 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0.0);
1101
1102 // store initial residual
1103 norms[0] = relres;
1104
1105 /* GMRES(M) outer iteration */
1106 while (iter < MaxIt && relres > tol) {
1107
1108 rs[0] = r_norm;
1109
1110 t = 1.0 / r_norm;
1111
1112 fasp_blas_darray_ax(n, t, p[0]);
1113
1114 /* RESTART CYCLE (right-preconditioning) */
1115 i = 0;
1116 while (i < Restart && iter < MaxIt) {
1117
1118 i++;
1119 iter++;
1120
1121 /* apply preconditioner */
1122 if (pc == NULL)
1123 fasp_darray_cp(n, p[i - 1], r);
1124 else
1125 pc->fct(p[i - 1], r, pc->data);
1126
1127 fasp_blas_dstr_mxv(A, r, p[i]);
1128
1129 /* Modified Gram_Schmidt orthogonalization */
1130 for (j = 0; j < i; j++) {
1131 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1132 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
1133 }
1134 t = fasp_blas_darray_norm2(n, p[i]);
1135 hh[i][i - 1] = t;
1136
1137 if (ABS(t) > SMALLREAL) { // If t=0, we get solution subspace
1138 t = 1.0 / t;
1139 fasp_blas_darray_ax(n, t, p[i]);
1140 }
1141
1142 for (j = 1; j < i; ++j) {
1143 t = hh[j - 1][i - 1];
1144 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
1145 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
1146 }
1147 t = hh[i][i - 1] * hh[i][i - 1];
1148 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
1149
1150 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
1151 c[i - 1] = hh[i - 1][i - 1] / gamma;
1152 s[i - 1] = hh[i][i - 1] / gamma;
1153 rs[i] = -s[i - 1] * rs[i - 1];
1154 rs[i - 1] = c[i - 1] * rs[i - 1];
1155 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
1156
1157 absres = r_norm = fabs(rs[i]);
1158
1159 relres = absres / absres0;
1160
1161 norms[iter] = relres;
1162
1163 // output iteration information if needed
1164 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1165 norms[iter] / norms[iter - 1]);
1166
1167 // exit restart cycle if reaches tolerance
1168 if (relres < tol && iter >= MIN_ITER) break;
1169
1170 } /* end of restart cycle */
1171
1172 /* compute solution, first solve upper triangular system */
1173 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
1174 for (k = i - 2; k >= 0; k--) {
1175 t = 0.0;
1176 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
1177 t += rs[k];
1178 rs[k] = t / hh[k][k];
1179 }
1180
1181 fasp_darray_cp(n, p[i - 1], w);
1182
1183 fasp_blas_darray_ax(n, rs[i - 1], w);
1184
1185 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1186
1187 /* apply preconditioner */
1188 if (pc == NULL)
1189 fasp_darray_cp(n, w, r);
1190 else
1191 pc->fct(w, r, pc->data);
1192
1193 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1194
1195 // Check: prevent false convergence
1196 if (relres < tol && iter >= MIN_ITER) {
1197
1198 REAL computed_relres = relres;
1199
1200 // compute residual
1201 fasp_darray_cp(n, b->val, r);
1202 fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1203 r_norm = fasp_blas_darray_norm2(n, r);
1204
1205 switch (StopType) {
1206 case STOP_REL_RES:
1207 absres = r_norm;
1208 relres = absres / absres0;
1209 break;
1210 case STOP_REL_PRECRES:
1211 if (pc == NULL)
1212 fasp_darray_cp(n, r, w);
1213 else
1214 pc->fct(r, w, pc->data);
1215 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
1216 relres = absres / absres0;
1217 break;
1218 case STOP_MOD_REL_RES:
1219 absres = r_norm;
1220 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
1221 relres = absres / normu;
1222 break;
1223 }
1224
1225 norms[iter] = relres;
1226
1227 if (relres < tol) {
1228 break;
1229 } else { // Need to restart
1230 fasp_darray_cp(n, r, p[0]);
1231 i = 0;
1232 }
1233
1234 if (PrtLvl >= PRINT_MORE) {
1235 ITS_COMPRES(computed_relres);
1236 ITS_REALRES(relres);
1237 }
1238
1239 } /* end of convergence check */
1240
1241 /* compute residual vector and continue loop */
1242 for (j = i; j > 0; j--) {
1243 rs[j - 1] = -s[j - 1] * rs[j];
1244 rs[j] = c[j - 1] * rs[j];
1245 }
1246
1247 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
1248
1249 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1250
1251 if (i) {
1252 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
1253 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1254 }
1255
1256 } /* end of main while loop */
1257
1258FINISHED:
1259 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
1260
1261 /*-------------------------------------------
1262 * Clean up workspace
1263 *------------------------------------------*/
1264 fasp_mem_free(work);
1265 work = NULL;
1266 fasp_mem_free(p);
1267 p = NULL;
1268 fasp_mem_free(hh);
1269 hh = NULL;
1270 fasp_mem_free(norms);
1271 norms = NULL;
1272
1273#if DEBUG_MODE > 0
1274 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1275#endif
1276
1277 if (iter >= MaxIt)
1278 return ERROR_SOLVER_MAXIT;
1279 else
1280 return iter;
1281}
1282
1309INT fasp_solver_pgmres(mxv_matfree* mf, dvector* b, dvector* x, precond* pc,
1310 const REAL tol, const REAL abstol, const INT MaxIt,
1311 const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
1312{
1313 const INT n = b->row;
1314 const INT min_iter = 0;
1315
1316 // local variables
1317 INT iter = 0;
1318 int i, j, k; // must be signed! -zcs
1319
1320 REAL epsmac = SMALLREAL;
1321 REAL r_norm, b_norm, den_norm;
1322 REAL epsilon, gamma, t;
1323
1324 // allocate temp memory (need about (restart+4)*n REAL numbers)
1325 REAL * c = NULL, *s = NULL, *rs = NULL;
1326 REAL * norms = NULL, *r = NULL, *w = NULL;
1327 REAL* work = NULL;
1328 REAL **p = NULL, **hh = NULL;
1329
1330 INT Restart = restart;
1331 INT Restart1 = Restart + 1;
1332 LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
1333
1334 // Output some info for debugging
1335 if (PrtLvl > PRINT_NONE) printf("\nCalling GMRes solver (MatFree) ...\n");
1336
1337#if DEBUG_MODE > 0
1338 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1339 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1340#endif
1341
1342 /* allocate memory and setup temp work space */
1343 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1344
1345 /* check whether memory is enough for GMRES */
1346 while ((work == NULL) && (Restart > 5)) {
1347 Restart = Restart - 5;
1348 worksize = (Restart + 4) * (Restart + n) + 1 - n;
1349 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1350 Restart1 = Restart + 1;
1351 }
1352
1353 if (work == NULL) {
1354 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1355 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1356 }
1357
1358 if (PrtLvl > PRINT_MIN && Restart < restart) {
1359 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
1360 }
1361
1362 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1363 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1364 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
1365
1366 r = work;
1367 w = r + n;
1368 rs = w + n;
1369 c = rs + Restart1;
1370 s = c + Restart;
1371
1372 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
1373
1374 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
1375
1376 /* initialization */
1377 mf->fct(mf->data, x->val, p[0]);
1378 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, p[0]);
1379
1380 b_norm = fasp_blas_darray_norm2(n, b->val);
1381 r_norm = fasp_blas_darray_norm2(n, p[0]);
1382
1383 if (PrtLvl > PRINT_NONE) {
1384 norms[0] = r_norm;
1385 if (PrtLvl >= PRINT_SOME) {
1386 ITS_PUTNORM("right-hand side", b_norm);
1387 ITS_PUTNORM("residual", r_norm);
1388 }
1389 }
1390
1391 if (b_norm > 0.0)
1392 den_norm = b_norm;
1393 else
1394 den_norm = r_norm;
1395
1396 epsilon = tol * den_norm;
1397
1398 /* outer iteration cycle */
1399 while (iter < MaxIt) {
1400
1401 rs[0] = r_norm;
1402 if (r_norm == 0.0) {
1403 fasp_mem_free(work);
1404 work = NULL;
1405 fasp_mem_free(p);
1406 p = NULL;
1407 fasp_mem_free(hh);
1408 hh = NULL;
1409 fasp_mem_free(norms);
1410 norms = NULL;
1411 return iter;
1412 }
1413
1414 if (r_norm <= epsilon && iter >= min_iter) {
1415 mf->fct(mf->data, x->val, r);
1416 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1417 r_norm = fasp_blas_darray_norm2(n, r);
1418
1419 if (r_norm <= epsilon) {
1420 break;
1421 } else {
1422 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1423 }
1424 }
1425
1426 t = 1.0 / r_norm;
1427 // for (j = 0; j < n; j ++) p[0][j] *= t;
1428 fasp_blas_darray_ax(n, t, p[0]);
1429
1430 /* RESTART CYCLE (right-preconditioning) */
1431 i = 0;
1432 while (i < Restart && iter < MaxIt) {
1433
1434 i++;
1435 iter++;
1436
1437 /* apply preconditioner */
1438 if (pc == NULL)
1439 fasp_darray_cp(n, p[i - 1], r);
1440 else
1441 pc->fct(p[i - 1], r, pc->data);
1442
1443 mf->fct(mf->data, r, p[i]);
1444
1445 /* modified Gram_Schmidt */
1446 for (j = 0; j < i; j++) {
1447 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1448 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
1449 }
1450 t = fasp_blas_darray_norm2(n, p[i]);
1451 hh[i][i - 1] = t;
1452 if (t != 0.0) {
1453 t = 1.0 / t;
1454 // for (j = 0; j < n; j ++) p[i][j] *= t;
1455 fasp_blas_darray_ax(n, t, p[i]);
1456 }
1457
1458 for (j = 1; j < i; ++j) {
1459 t = hh[j - 1][i - 1];
1460 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
1461 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
1462 }
1463 t = hh[i][i - 1] * hh[i][i - 1];
1464 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
1465 gamma = sqrt(t);
1466 if (gamma == 0.0) gamma = epsmac;
1467 c[i - 1] = hh[i - 1][i - 1] / gamma;
1468 s[i - 1] = hh[i][i - 1] / gamma;
1469 rs[i] = -s[i - 1] * rs[i - 1];
1470 rs[i - 1] = c[i - 1] * rs[i - 1];
1471 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
1472 r_norm = fabs(rs[i]);
1473
1474 norms[iter] = r_norm;
1475
1476 if (b_norm > 0) {
1477 fasp_itinfo(PrtLvl, StopType, iter, norms[iter] / b_norm, norms[iter],
1478 norms[iter] / norms[iter - 1]);
1479 } else {
1480 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
1481 norms[iter] / norms[iter - 1]);
1482 }
1483
1484 /* should we exit restart cycle? */
1485 if (r_norm <= epsilon && iter >= min_iter) {
1486 break;
1487 }
1488 } /* end of restart cycle */
1489
1490 /* now compute solution, first solve upper triangular system */
1491 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
1492 for (k = i - 2; k >= 0; k--) {
1493 t = 0.0;
1494 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
1495
1496 t += rs[k];
1497 rs[k] = t / hh[k][k];
1498 }
1499 fasp_darray_cp(n, p[i - 1], w);
1500 // for (j = 0; j < n; j ++) w[j] *= rs[i-1];
1501 fasp_blas_darray_ax(n, rs[i - 1], w);
1502 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1503
1504 /* apply preconditioner */
1505 if (pc == NULL)
1506 fasp_darray_cp(n, w, r);
1507 else
1508 pc->fct(w, r, pc->data);
1509
1510 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1511
1512 if (r_norm <= epsilon && iter >= min_iter) {
1513 mf->fct(mf->data, x->val, r);
1514 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1515 r_norm = fasp_blas_darray_norm2(n, r);
1516
1517 if (r_norm <= epsilon) {
1518 break;
1519 } else {
1520 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1521 fasp_darray_cp(n, r, p[0]);
1522 i = 0;
1523 }
1524 } /* end of convergence check */
1525
1526 /* compute residual vector and continue loop */
1527 for (j = i; j > 0; j--) {
1528 rs[j - 1] = -s[j - 1] * rs[j];
1529 rs[j] = c[j - 1] * rs[j];
1530 }
1531
1532 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
1533
1534 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1535
1536 if (i) {
1537 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
1538 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1539 }
1540 } /* end of iteration while loop */
1541
1542 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, r_norm);
1543
1544 /*-------------------------------------------
1545 * Clean up workspace
1546 *------------------------------------------*/
1547 fasp_mem_free(work);
1548 work = NULL;
1549 fasp_mem_free(p);
1550 p = NULL;
1551 fasp_mem_free(hh);
1552 hh = NULL;
1553 fasp_mem_free(norms);
1554 norms = NULL;
1555
1556#if DEBUG_MODE > 0
1557 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1558#endif
1559
1560 if (iter >= MaxIt)
1561 return ERROR_SOLVER_MAXIT;
1562 else
1563 return iter;
1564}
1565
1566#if 0
1567static double estimate_spectral_radius (const double **A, int n, size_t k = 20)
1568{
1569 double *x = (double *)malloc(n* sizeof(double));
1570 double *y = (double *)malloc(n* sizeof(double));
1571 double *z = (double *)malloc(n* sizeof(double));
1572 double t;
1573 int i1,j1;
1574
1575 // initialize x to random values in [0,1)
1576 // cusp::copy(cusp::detail::random_reals<ValueType>(N), x);
1577 dvector px;
1578 px.row = n;
1579 px.val = x;
1580
1581 fasp_dvec_rand(n, &px);
1582
1583 for(size_t i = 0; i < k; i++)
1584 {
1585 //cusp::blas::scal(x, ValueType(1.0) / cusp::blas::nrmmax(x));
1586 t= 1.0/ fasp_blas_darray_norminf(n, px);
1587 for(i1= 0; i1 <n; i1++) x[i1] *= t;
1588
1589 //cusp::multiply(A, x, y);
1590
1591 for(i1= 0; i1 <n; i1++) {
1592 t= 0.0
1593 for(j1= 0; j1 <n; j1++) t += A[i1][j1] * x[j1];
1594 y[i1] = t; }
1595 // x.swap(y);
1596 for(i1= 0; i1 <n; i1++) z[i1] = x[i1];
1597 for(i1= 0; i1 <n; i1++) x[i1] = y[i1];
1598 for(i1= 0; i1 <n; i1++) y[i1] = z[i1];
1599 }
1600
1601 free(x);
1602 free(y);
1603 free(z);
1604
1605 if (k == 0)
1606 return 0;
1607 else
1608 //return cusp::blas::nrm2(x) / cusp::blas::nrm2(y);
1609 return fasp_blas_darray_norm2(n,x) / fasp_blas_darray_norm2(n,y) ;
1610}
1611
1612static double spectral_radius (dCSRmat *A,
1613 const SHORT restart)
1614{
1615 const INT n = A->row;
1616 const INT MIN_ITER = 0;
1617
1618 // local variables
1619 INT iter = 0;
1620 INT Restart1 = restart + 1;
1621 INT i, j, k;
1622
1623 REAL r_norm, den_norm;
1624 REAL epsilon, gamma, t;
1625
1626 REAL *c = NULL, *s = NULL, *rs = NULL;
1627 REAL *norms = NULL, *r = NULL, *w = NULL;
1628 REAL **p = NULL, **hh = NULL;
1629 REAL *work = NULL;
1630
1631 /* allocate memory */
1632 work = (REAL *)fasp_mem_calloc((restart+4)*(restart+n)+1-n, sizeof(REAL));
1633 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1634 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1635
1636 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1637
1638 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + restart;
1639
1640 for (i = 0; i < Restart1; i ++) p[i] = s + restart + i*n;
1641 for (i = 0; i < Restart1; i ++) hh[i] = p[restart] + n + i*restart;
1642
1643 /* initialization */
1644 dvector p0;
1645 p0.row = n;
1646 p0.val = p[0];
1647 fasp_dvec_rand(n, &p0);
1648
1649 r_norm = fasp_blas_darray_norm2(n, p[0]);
1650 t = 1.0 / r_norm;
1651 for (j = 0; j < n; j ++) p[0][j] *= t;
1652
1653 int maxiter = MIN(n, restart) ;
1654 for ( j = 0; j < maxiter; j++ ) {
1655 fasp_blas_bdbsr_mxv(A, p[j], p[j+1]);
1656
1657 for( i = 0; i <= j; i++ ) {
1658 hh[i][j] = fasp_blas_darray_dotprod(n, p[i], p[j+1]);
1659 fasp_blas_darray_axpy(n, -hh[i][j], p[i], p[ j+1 ]);
1660 }
1661
1662 hh[j+1][j] = fasp_blas_darray_norm2 (n, p[j+1]);
1663 if ( hh[j+1][j] < 1e-10) break;
1664 t = 1.0/hh[j+1][j];
1665 for (k = 0; k < n; k ++) p[j+1][k] *= t;
1666 }
1667
1668 H = (REAL **)fasp_mem_calloc(j, sizeof(REAL *));
1669 H[0] = (REAL *)fasp_mem_calloc(j*j, sizeof(REAL));
1670 for (i = 1; i < j; i ++) H[i] = H[i-1] + j;
1671
1672
1673 for( size_t row = 0; row < j; row++ )
1674 for( size_t col = 0; col < j; col++ )
1675 H[row][col] = hh[row][col];
1676
1677 double spectral_radius = estimate_spectral_radius( H, j, 20);
1678
1679 /*-------------------------------------------
1680 * Clean up workspace
1681 *------------------------------------------*/
1682 fasp_mem_free(work); work = NULL;
1683 fasp_mem_free(p); p = NULL;
1684 fasp_mem_free(hh); hh = NULL;
1685 fasp_mem_free(norms); norms = NULL;
1686 fasp_mem_free(H[0]); H[0] = NULL;
1687 fasp_mem_free(H); H = NULL;
1688
1689 return spectral_radius;
1690}
1691#endif
1692
1693/*---------------------------------*/
1694/*-- End of File --*/
1695/*---------------------------------*/
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_itinfo
void fasp_itinfo(const INT ptrlvl, const INT stop_type, const INT iter, const REAL relres, const REAL absres, const REAL factor)
Print out iteration information for iterative solvers.
Definition: AuxMessage.c:41
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_dvec_rand
void fasp_dvec_rand(const INT n, dvector *x)
Generate fake random REAL vector in the range from 0 to 1.
Definition: AuxVector.c:192
fasp_blas_darray_dotprod
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
Definition: BlaArray.c:771
fasp_blas_darray_norminf
REAL fasp_blas_darray_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: BlaArray.c:719
fasp_blas_darray_axpby
void fasp_blas_darray_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: BlaArray.c:620
fasp_blas_darray_norm2
REAL fasp_blas_darray_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: BlaArray.c:691
fasp_blas_darray_ax
void fasp_blas_darray_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: BlaArray.c:43
fasp_blas_darray_axpy
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:90
fasp_blas_dblc_mxv
void fasp_blas_dblc_mxv(const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvBLC.c:164
fasp_blas_dblc_aAxpy
void fasp_blas_dblc_aAxpy(const REAL alpha, const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvBLC.c:38
fasp_blas_dbsr_aAxpy
void fasp_blas_dbsr_aAxpy(const REAL alpha, const dBSRmat *A, const REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: BlaSpmvBSR.c:514
fasp_blas_dbsr_mxv
void fasp_blas_dbsr_mxv(const dBSRmat *A, const REAL *x, REAL *y)
Compute y := A*x.
Definition: BlaSpmvBSR.c:1055
fasp_blas_dcsr_mxv
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:242
fasp_blas_dcsr_aAxpy
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:494
fasp_blas_dstr_aAxpy
void fasp_blas_dstr_aAxpy(const REAL alpha, const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvSTR.c:61
fasp_blas_dstr_mxv
void fasp_blas_dstr_mxv(const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvSTR.c:131
fasp_solver_dbsr_pgmres
INT fasp_solver_dbsr_pgmres(dBSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:376
fasp_solver_dcsr_pgmres
INT fasp_solver_dcsr_pgmres(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Right preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:66
fasp_solver_dstr_pgmres
INT fasp_solver_dstr_pgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:998
fasp_solver_dblc_pgmres
INT fasp_solver_dblc_pgmres(dBLCmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:687
fasp_solver_pgmres
INT fasp_solver_pgmres(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Solve "Ax=b" using PGMRES (right preconditioned) iterative method.
Definition: KryPgmres.c:1309
fasp.h
Main header file for the FASP project.
MIN
#define MIN(a, b)
Definition: fasp.h:83
REAL
#define REAL
Definition: fasp.h:75
SHORT
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:71
LONG
#define LONG
Definition: fasp.h:73
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
ERROR_SOLVER_MAXIT
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
STOP_REL_RES
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:132
STOP_MOD_REL_RES
#define STOP_MOD_REL_RES
Definition: fasp_const.h:134
STOP_REL_PRECRES
#define STOP_REL_PRECRES
Definition: fasp_const.h:133
BIGREAL
#define BIGREAL
Some global constants.
Definition: fasp_const.h:255
PRINT_SOME
#define PRINT_SOME
Definition: fasp_const.h:75
ERROR_ALLOC_MEM
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
PRINT_MORE
#define PRINT_MORE
Definition: fasp_const.h:76
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
SMALLREAL
#define SMALLREAL
Definition: fasp_const.h:256
PRINT_MIN
#define PRINT_MIN
Definition: fasp_const.h:74
dBLCmat
Block REAL CSR matrix format.
Definition: fasp_block.h:74
dBSRmat
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
dCSRmat
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:151
dCSRmat::row
INT row
row number of matrix A, m
Definition: fasp.h:154
dSTRmat
Structure matrix of REAL type.
Definition: fasp.h:316
dvector
Vector with n entries of REAL type.
Definition: fasp.h:354
dvector::val
REAL * val
actual vector entries
Definition: fasp.h:360
dvector::row
INT row
number of rows
Definition: fasp.h:357
mxv_matfree
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1109
mxv_matfree::data
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1112
mxv_matfree::fct
void(* fct)(const void *, const REAL *, REAL *)
action for MxV, void function pointer
Definition: fasp.h:1115
precond
Preconditioner data and action.
Definition: fasp.h:1095
precond::data
void * data
data for preconditioner, void pointer
Definition: fasp.h:1098
precond::fct
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1101