Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KryPvgmres.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_pvgmres(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 const REAL epsmac = SMALLREAL;
74
75 //--------------------------------------------//
76 // Newly added parameters to monitor when //
77 // to change the restart parameter //
78 //--------------------------------------------//
79 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
80 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
81
82 // local variables
83 INT iter = 0;
84 int i, j, k; // must be signed! -zcs
85
86 REAL r_norm, r_normb, gamma, t;
87 REAL absres0 = BIGREAL, absres = BIGREAL;
88 REAL relres = BIGREAL, normu = BIGREAL;
89
90 REAL cr = 1.0; // convergence rate
91 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
92 INT d = 3; // reduction for restart parameter
93 INT restart_max = restart; // upper bound for restart in each restart cycle
94 INT restart_min = 3; // lower bound for restart in each restart cycle
95
96 INT Restart = restart; // real restart in some fixed restarted cycle
97 INT Restart1 = Restart + 1;
98 unsigned LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
99
100 // allocate temp memory (need about (restart+4)*n REAL numbers)
101 REAL * c = NULL, *s = NULL, *rs = NULL;
102 REAL * norms = NULL, *r = NULL, *w = NULL;
103 REAL* work = NULL;
104 REAL **p = NULL, **hh = NULL;
105
106 // Output some info for debuging
107 if (PrtLvl > PRINT_NONE) printf("\nCalling VGMRes solver (CSR) ...\n");
108
109#if DEBUG_MODE > 0
110 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
111 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
112#endif
113
114 /* allocate memory and setup temp work space */
115 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
116
117 /* check whether memory is enough for GMRES */
118 while ((work == NULL) && (Restart > 5)) {
119 Restart = Restart - 5;
120 worksize = (Restart + 4) * (Restart + n) + 1 - n;
121 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
122 Restart1 = Restart + 1;
123 }
124
125 if (work == NULL) {
126 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
127 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
128 }
129
130 if (PrtLvl > PRINT_MIN && Restart < restart) {
131 printf("### WARNING: vGMRES restart number set to %d!\n", Restart);
132 }
133
134 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
135 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
136 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
137
138 r = work;
139 w = r + n;
140 rs = w + n;
141 c = rs + Restart1;
142 s = c + Restart;
143
144 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
145
146 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
147
148 // r = b-A*x
149 fasp_darray_cp(n, b->val, p[0]);
150 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
151
152 r_norm = fasp_blas_darray_norm2(n, p[0]);
153
154 // compute initial residuals
155 switch (StopType) {
156 case STOP_REL_RES:
157 absres0 = MAX(SMALLREAL, r_norm);
158 relres = r_norm / absres0;
159 break;
160 case STOP_REL_PRECRES:
161 if (pc == NULL)
162 fasp_darray_cp(n, p[0], r);
163 else
164 pc->fct(p[0], r, pc->data);
165 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
166 absres0 = MAX(SMALLREAL, r_normb);
167 relres = r_normb / absres0;
168 break;
169 case STOP_MOD_REL_RES:
170 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
171 absres0 = r_norm;
172 relres = absres0 / normu;
173 break;
174 default:
175 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
176 goto FINISHED;
177 }
178
179 // if initial residual is small, no need to iterate!
180 if (relres < tol || absres0 < abstol) goto FINISHED;
181
182 // output iteration information if needed
183 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0);
184
185 // store initial residual
186 norms[0] = relres;
187
188 /* outer iteration cycle */
189 while (iter < MaxIt) {
190
191 rs[0] = r_norm_old = r_norm;
192
193 t = 1.0 / r_norm;
194
195 fasp_blas_darray_ax(n, t, p[0]);
196
197 //-----------------------------------//
198 // adjust the restart parameter //
199 //-----------------------------------//
200 if (cr > cr_max || iter == 0) {
201 Restart = restart_max;
202 } else if (cr < cr_min) {
203 // Restart = Restart;
204 } else {
205 if (Restart - d > restart_min) {
206 Restart -= d;
207 } else {
208 Restart = restart_max;
209 }
210 }
211
212 /* RESTART CYCLE (right-preconditioning) */
213 i = 0;
214 while (i < Restart && iter < MaxIt) {
215
216 i++;
217 iter++;
218
219 /* apply preconditioner */
220 if (pc == NULL)
221 fasp_darray_cp(n, p[i - 1], r);
222 else
223 pc->fct(p[i - 1], r, pc->data);
224
225 fasp_blas_dcsr_mxv(A, r, p[i]);
226
227 /* modified Gram_Schmidt */
228 for (j = 0; j < i; j++) {
229 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
230 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
231 }
232 t = fasp_blas_darray_norm2(n, p[i]);
233 hh[i][i - 1] = t;
234 if (t != 0.0) {
235 t = 1.0 / t;
236 fasp_blas_darray_ax(n, t, p[i]);
237 }
238
239 for (j = 1; j < i; ++j) {
240 t = hh[j - 1][i - 1];
241 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
242 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
243 }
244 t = hh[i][i - 1] * hh[i][i - 1];
245 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
246
247 gamma = sqrt(t);
248 if (gamma == 0.0) gamma = epsmac;
249 c[i - 1] = hh[i - 1][i - 1] / gamma;
250 s[i - 1] = hh[i][i - 1] / gamma;
251 rs[i] = -s[i - 1] * rs[i - 1];
252 rs[i - 1] = c[i - 1] * rs[i - 1];
253 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
254
255 absres = r_norm = fabs(rs[i]);
256
257 relres = absres / absres0;
258
259 norms[iter] = relres;
260
261 // output iteration information if needed
262 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
263 norms[iter] / norms[iter - 1]);
264
265 // should we exit restart cycle
266 if (relres < tol && iter >= MIN_ITER) break;
267
268 } /* end of restart cycle */
269
270 /* now compute solution, first solve upper triangular system */
271 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
272 for (k = i - 2; k >= 0; k--) {
273 t = 0.0;
274 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
275
276 t += rs[k];
277 rs[k] = t / hh[k][k];
278 }
279
280 fasp_darray_cp(n, p[i - 1], w);
281
282 fasp_blas_darray_ax(n, rs[i - 1], w);
283
284 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
285
286 /* apply preconditioner */
287 if (pc == NULL)
288 fasp_darray_cp(n, w, r);
289 else
290 pc->fct(w, r, pc->data);
291
292 fasp_blas_darray_axpy(n, 1.0, r, x->val);
293
294 // Check: prevent false convergence
295 if (relres < tol && iter >= MIN_ITER) {
296
297 REAL computed_relres = relres;
298
299 // compute current residual
300 fasp_darray_cp(n, b->val, r);
301 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
302
303 r_norm = fasp_blas_darray_norm2(n, r);
304
305 switch (StopType) {
306 case STOP_REL_RES:
307 absres = r_norm;
308 relres = absres / absres0;
309 break;
310 case STOP_REL_PRECRES:
311 if (pc == NULL)
312 fasp_darray_cp(n, r, w);
313 else
314 pc->fct(r, w, pc->data);
315 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
316 relres = absres / absres0;
317 break;
318 case STOP_MOD_REL_RES:
319 absres = r_norm;
320 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
321 relres = absres / normu;
322 break;
323 }
324
325 norms[iter] = relres;
326
327 if (relres < tol) {
328 break;
329 } else {
330 // Need to restart
331 fasp_darray_cp(n, r, p[0]);
332 i = 0;
333 }
334
335 if (PrtLvl >= PRINT_MORE) {
336 ITS_COMPRES(computed_relres);
337 ITS_REALRES(relres);
338 }
339
340 } /* end of convergence check */
341
342 /* compute residual vector and continue loop */
343 for (j = i; j > 0; j--) {
344 rs[j - 1] = -s[j - 1] * rs[j];
345 rs[j] = c[j - 1] * rs[j];
346 }
347
348 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
349
350 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
351
352 if (i) {
353 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
354 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
355 }
356
357 //-----------------------------------//
358 // compute the convergence rate //
359 //-----------------------------------//
360 cr = r_norm / r_norm_old;
361
362 } /* end of iteration while loop */
363
364FINISHED:
365 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
366
367 /*-------------------------------------------
368 * Free some stuff
369 *------------------------------------------*/
370 fasp_mem_free(work);
371 work = NULL;
372 fasp_mem_free(p);
373 p = NULL;
374 fasp_mem_free(hh);
375 hh = NULL;
376 fasp_mem_free(norms);
377 norms = NULL;
378
379#if DEBUG_MODE > 0
380 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
381#endif
382
383 if (iter >= MaxIt)
384 return ERROR_SOLVER_MAXIT;
385 else
386 return iter;
387}
388
416INT fasp_solver_dbsr_pvgmres(dBSRmat* A, dvector* b, dvector* x, precond* pc,
417 const REAL tol, const REAL abstol, const INT MaxIt,
418 const SHORT restart, const SHORT StopType,
419 const SHORT PrtLvl)
420{
421 const INT n = b->row;
422 const INT MIN_ITER = 0;
423 const REAL epsmac = SMALLREAL;
424
425 //--------------------------------------------//
426 // Newly added parameters to monitor when //
427 // to change the restart parameter //
428 //--------------------------------------------//
429 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
430 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
431
432 // local variables
433 INT iter = 0;
434 int i, j, k; // must be signed! -zcs
435
436 REAL r_norm, r_normb, gamma, t;
437 REAL absres0 = BIGREAL, absres = BIGREAL;
438 REAL relres = BIGREAL, normu = BIGREAL;
439
440 REAL cr = 1.0; // convergence rate
441 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
442 INT d = 3; // reduction for restart parameter
443 INT restart_max = restart; // upper bound for restart in each restart cycle
444 INT restart_min =
445 3; // lower bound for restart in each restart cycle (should be small)
446
447 INT Restart = restart; // real restart in some fixed restarted cycle
448 INT Restart1 = Restart + 1;
449 unsigned LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
450
451 // allocate temp memory (need about (restart+4)*n REAL numbers)
452 REAL * c = NULL, *s = NULL, *rs = NULL;
453 REAL * norms = NULL, *r = NULL, *w = NULL;
454 REAL* work = NULL;
455 REAL **p = NULL, **hh = NULL;
456
457 // Output some info for debuging
458 if (PrtLvl > PRINT_NONE) printf("\nCalling VGMRes solver (BSR) ...\n");
459
460#if DEBUG_MODE > 0
461 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
462 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
463#endif
464
465 /* allocate memory and setup temp work space */
466 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
467
468 /* check whether memory is enough for GMRES */
469 while ((work == NULL) && (Restart > 5)) {
470 Restart = Restart - 5;
471 worksize = (Restart + 4) * (Restart + n) + 1 - n;
472 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
473 Restart1 = Restart + 1;
474 }
475
476 if (work == NULL) {
477 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
478 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
479 }
480
481 if (PrtLvl > PRINT_MIN && Restart < restart) {
482 printf("### WARNING: vGMRES restart number set to %d!\n", Restart);
483 }
484
485 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
486 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
487 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
488
489 r = work;
490 w = r + n;
491 rs = w + n;
492 c = rs + Restart1;
493 s = c + Restart;
494
495 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
496
497 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
498
499 // r = b-A*x
500 fasp_darray_cp(n, b->val, p[0]);
501 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
502
503 r_norm = fasp_blas_darray_norm2(n, p[0]);
504
505 // compute initial residuals
506 switch (StopType) {
507 case STOP_REL_RES:
508 absres0 = MAX(SMALLREAL, r_norm);
509 relres = r_norm / absres0;
510 break;
511 case STOP_REL_PRECRES:
512 if (pc == NULL)
513 fasp_darray_cp(n, p[0], r);
514 else
515 pc->fct(p[0], r, pc->data);
516 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
517 absres0 = MAX(SMALLREAL, r_normb);
518 relres = r_normb / absres0;
519 break;
520 case STOP_MOD_REL_RES:
521 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
522 absres0 = r_norm;
523 relres = absres0 / normu;
524 break;
525 default:
526 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
527 goto FINISHED;
528 }
529
530 // if initial residual is small, no need to iterate!
531 if (relres < tol || absres0 < abstol) goto FINISHED;
532
533 // output iteration information if needed
534 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0);
535
536 // store initial residual
537 norms[0] = relres;
538
539 /* outer iteration cycle */
540 while (iter < MaxIt) {
541
542 rs[0] = r_norm_old = r_norm;
543
544 t = 1.0 / r_norm;
545
546 fasp_blas_darray_ax(n, t, p[0]);
547
548 //-----------------------------------//
549 // adjust the restart parameter //
550 //-----------------------------------//
551 if (cr > cr_max || iter == 0) {
552 Restart = restart_max;
553 } else if (cr < cr_min) {
554 // Restart = Restart;
555 } else {
556 if (Restart - d > restart_min) {
557 Restart -= d;
558 } else {
559 Restart = restart_max;
560 }
561 }
562
563 /* RESTART CYCLE (right-preconditioning) */
564 i = 0;
565 while (i < Restart && iter < MaxIt) {
566
567 i++;
568 iter++;
569
570 /* apply preconditioner */
571 if (pc == NULL)
572 fasp_darray_cp(n, p[i - 1], r);
573 else
574 pc->fct(p[i - 1], r, pc->data);
575
576 fasp_blas_dbsr_mxv(A, r, p[i]);
577
578 /* modified Gram_Schmidt */
579 for (j = 0; j < i; j++) {
580 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
581 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
582 }
583 t = fasp_blas_darray_norm2(n, p[i]);
584 hh[i][i - 1] = t;
585 if (t != 0.0) {
586 t = 1.0 / t;
587 fasp_blas_darray_ax(n, t, p[i]);
588 }
589
590 for (j = 1; j < i; ++j) {
591 t = hh[j - 1][i - 1];
592 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
593 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
594 }
595 t = hh[i][i - 1] * hh[i][i - 1];
596 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
597
598 gamma = sqrt(t);
599 if (gamma == 0.0) gamma = epsmac;
600 c[i - 1] = hh[i - 1][i - 1] / gamma;
601 s[i - 1] = hh[i][i - 1] / gamma;
602 rs[i] = -s[i - 1] * rs[i - 1];
603 rs[i - 1] = c[i - 1] * rs[i - 1];
604 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
605
606 absres = r_norm = fabs(rs[i]);
607
608 relres = absres / absres0;
609
610 norms[iter] = relres;
611
612 // output iteration information if needed
613 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
614 norms[iter] / norms[iter - 1]);
615
616 // should we exit restart cycle
617 if (relres < tol && iter >= MIN_ITER) break;
618
619 } /* end of restart cycle */
620
621 /* now compute solution, first solve upper triangular system */
622 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
623 for (k = i - 2; k >= 0; k--) {
624 t = 0.0;
625 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
626
627 t += rs[k];
628 rs[k] = t / hh[k][k];
629 }
630
631 fasp_darray_cp(n, p[i - 1], w);
632
633 fasp_blas_darray_ax(n, rs[i - 1], w);
634
635 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
636
637 /* apply preconditioner */
638 if (pc == NULL)
639 fasp_darray_cp(n, w, r);
640 else
641 pc->fct(w, r, pc->data);
642
643 fasp_blas_darray_axpy(n, 1.0, r, x->val);
644
645 // Check: prevent false convergence
646 if (relres < tol && iter >= MIN_ITER) {
647
648 REAL computed_relres = relres;
649
650 // compute current residual
651 fasp_darray_cp(n, b->val, r);
652 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
653
654 r_norm = fasp_blas_darray_norm2(n, r);
655
656 switch (StopType) {
657 case STOP_REL_RES:
658 absres = r_norm;
659 relres = absres / absres0;
660 break;
661 case STOP_REL_PRECRES:
662 if (pc == NULL)
663 fasp_darray_cp(n, r, w);
664 else
665 pc->fct(r, w, pc->data);
666 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
667 relres = absres / absres0;
668 break;
669 case STOP_MOD_REL_RES:
670 absres = r_norm;
671 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
672 relres = absres / normu;
673 break;
674 }
675
676 norms[iter] = relres;
677
678 if (relres < tol) {
679 break;
680 } else {
681 // Need to restart
682 fasp_darray_cp(n, r, p[0]);
683 i = 0;
684 }
685
686 if (PrtLvl >= PRINT_MORE) {
687 ITS_COMPRES(computed_relres);
688 ITS_REALRES(relres);
689 }
690
691 } /* end of convergence check */
692
693 /* compute residual vector and continue loop */
694 for (j = i; j > 0; j--) {
695 rs[j - 1] = -s[j - 1] * rs[j];
696 rs[j] = c[j - 1] * rs[j];
697 }
698
699 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
700
701 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
702
703 if (i) {
704 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
705 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
706 }
707
708 //-----------------------------------//
709 // compute the convergence rate //
710 //-----------------------------------//
711 cr = r_norm / r_norm_old;
712
713 } /* end of iteration while loop */
714
715FINISHED:
716 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
717
718 /*-------------------------------------------
719 * Free some stuff
720 *------------------------------------------*/
721 fasp_mem_free(work);
722 work = NULL;
723 fasp_mem_free(p);
724 p = NULL;
725 fasp_mem_free(hh);
726 hh = NULL;
727 fasp_mem_free(norms);
728 norms = NULL;
729
730#if DEBUG_MODE > 0
731 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
732#endif
733
734 if (iter >= MaxIt)
735 return ERROR_SOLVER_MAXIT;
736 else
737 return iter;
738}
739
765INT fasp_solver_dblc_pvgmres(dBLCmat* A, dvector* b, dvector* x, precond* pc,
766 const REAL tol, const REAL abstol, const INT MaxIt,
767 const SHORT restart, const SHORT StopType,
768 const SHORT PrtLvl)
769{
770 const INT n = b->row;
771 const INT MIN_ITER = 0;
772 const REAL epsmac = SMALLREAL;
773
774 //--------------------------------------------//
775 // Newly added parameters to monitor when //
776 // to change the restart parameter //
777 //--------------------------------------------//
778 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
779 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
780
781 // local variables
782 INT iter = 0;
783 int i, j, k; // must be signed! -zcs
784
785 REAL r_norm, r_normb, gamma, t;
786 REAL absres0 = BIGREAL, absres = BIGREAL;
787 REAL relres = BIGREAL, normu = BIGREAL;
788
789 REAL cr = 1.0; // convergence rate
790 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
791 INT d = 3; // reduction for restart parameter
792 INT restart_max = restart; // upper bound for restart in each restart cycle
793 INT restart_min =
794 3; // lower bound for restart in each restart cycle (should be small)
795
796 INT Restart = restart; // real restart in some fixed restarted cycle
797 INT Restart1 = Restart + 1;
798 unsigned LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
799
800 // allocate temp memory (need about (restart+4)*n REAL numbers)
801 REAL * c = NULL, *s = NULL, *rs = NULL;
802 REAL * norms = NULL, *r = NULL, *w = NULL;
803 REAL* work = NULL;
804 REAL **p = NULL, **hh = NULL;
805
806 // Output some info for debuging
807 if (PrtLvl > PRINT_NONE) printf("\nCalling VGMRes solver (BLC) ...\n");
808
809#if DEBUG_MODE > 0
810 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
811 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
812#endif
813
814 /* allocate memory and setup temp work space */
815 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
816
817 /* check whether memory is enough for GMRES */
818 while ((work == NULL) && (Restart > 5)) {
819 Restart = Restart - 5;
820 worksize = (Restart + 4) * (Restart + n) + 1 - n;
821 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
822 Restart1 = Restart + 1;
823 }
824
825 if (work == NULL) {
826 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
827 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
828 }
829
830 if (PrtLvl > PRINT_MIN && Restart < restart) {
831 printf("### WARNING: vGMRES restart number set to %d!\n", Restart);
832 }
833
834 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
835 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
836 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
837
838 r = work;
839 w = r + n;
840 rs = w + n;
841 c = rs + Restart1;
842 s = c + Restart;
843
844 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
845
846 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
847
848 // r = b-A*x
849 fasp_darray_cp(n, b->val, p[0]);
850 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
851
852 r_norm = fasp_blas_darray_norm2(n, p[0]);
853
854 // compute initial residuals
855 switch (StopType) {
856 case STOP_REL_RES:
857 absres0 = MAX(SMALLREAL, r_norm);
858 relres = r_norm / absres0;
859 break;
860 case STOP_REL_PRECRES:
861 if (pc == NULL)
862 fasp_darray_cp(n, p[0], r);
863 else
864 pc->fct(p[0], r, pc->data);
865 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
866 absres0 = MAX(SMALLREAL, r_normb);
867 relres = r_normb / absres0;
868 break;
869 case STOP_MOD_REL_RES:
870 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
871 absres0 = r_norm;
872 relres = absres0 / normu;
873 break;
874 default:
875 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
876 goto FINISHED;
877 }
878
879 // if initial residual is small, no need to iterate!
880 if (relres < tol || absres0 < abstol) goto FINISHED;
881
882 // output iteration information if needed
883 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0);
884
885 // store initial residual
886 norms[0] = relres;
887
888 /* outer iteration cycle */
889 while (iter < MaxIt) {
890
891 rs[0] = r_norm_old = r_norm;
892
893 t = 1.0 / r_norm;
894
895 fasp_blas_darray_ax(n, t, p[0]);
896
897 //-----------------------------------//
898 // adjust the restart parameter //
899 //-----------------------------------//
900 if (cr > cr_max || iter == 0) {
901 Restart = restart_max;
902 } else if (cr < cr_min) {
903 // Restart = Restart;
904 } else {
905 if (Restart - d > restart_min) {
906 Restart -= d;
907 } else {
908 Restart = restart_max;
909 }
910 }
911
912 /* RESTART CYCLE (right-preconditioning) */
913 i = 0;
914 while (i < Restart && iter < MaxIt) {
915
916 i++;
917 iter++;
918
919 /* apply preconditioner */
920 if (pc == NULL)
921 fasp_darray_cp(n, p[i - 1], r);
922 else
923 pc->fct(p[i - 1], r, pc->data);
924
925 fasp_blas_dblc_mxv(A, r, p[i]);
926
927 /* modified Gram_Schmidt */
928 for (j = 0; j < i; j++) {
929 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
930 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
931 }
932 t = fasp_blas_darray_norm2(n, p[i]);
933 hh[i][i - 1] = t;
934 if (t != 0.0) {
935 t = 1.0 / t;
936 fasp_blas_darray_ax(n, t, p[i]);
937 }
938
939 for (j = 1; j < i; ++j) {
940 t = hh[j - 1][i - 1];
941 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
942 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
943 }
944 t = hh[i][i - 1] * hh[i][i - 1];
945 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
946
947 gamma = sqrt(t);
948 if (gamma == 0.0) gamma = epsmac;
949 c[i - 1] = hh[i - 1][i - 1] / gamma;
950 s[i - 1] = hh[i][i - 1] / gamma;
951 rs[i] = -s[i - 1] * rs[i - 1];
952 rs[i - 1] = c[i - 1] * rs[i - 1];
953 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
954
955 absres = r_norm = fabs(rs[i]);
956
957 relres = absres / absres0;
958
959 norms[iter] = relres;
960
961 // output iteration information if needed
962 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
963 norms[iter] / norms[iter - 1]);
964
965 // should we exit restart cycle
966 if (relres < tol && iter >= MIN_ITER) break;
967
968 } /* end of restart cycle */
969
970 /* now compute solution, first solve upper triangular system */
971 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
972 for (k = i - 2; k >= 0; k--) {
973 t = 0.0;
974 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
975
976 t += rs[k];
977 rs[k] = t / hh[k][k];
978 }
979
980 fasp_darray_cp(n, p[i - 1], w);
981
982 fasp_blas_darray_ax(n, rs[i - 1], w);
983
984 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
985
986 /* apply preconditioner */
987 if (pc == NULL)
988 fasp_darray_cp(n, w, r);
989 else
990 pc->fct(w, r, pc->data);
991
992 fasp_blas_darray_axpy(n, 1.0, r, x->val);
993
994 // Check: prevent false convergence
995 if (relres < tol && iter >= MIN_ITER) {
996
997 REAL computed_relres = relres;
998
999 // compute current residual
1000 fasp_darray_cp(n, b->val, r);
1001 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
1002
1003 r_norm = fasp_blas_darray_norm2(n, r);
1004
1005 switch (StopType) {
1006 case STOP_REL_RES:
1007 absres = r_norm;
1008 relres = absres / absres0;
1009 break;
1010 case STOP_REL_PRECRES:
1011 if (pc == NULL)
1012 fasp_darray_cp(n, r, w);
1013 else
1014 pc->fct(r, w, pc->data);
1015 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
1016 relres = absres / absres0;
1017 break;
1018 case STOP_MOD_REL_RES:
1019 absres = r_norm;
1020 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
1021 relres = absres / normu;
1022 break;
1023 }
1024
1025 norms[iter] = relres;
1026
1027 if (relres < tol) {
1028 break;
1029 } else {
1030 // Need to restart
1031 fasp_darray_cp(n, r, p[0]);
1032 i = 0;
1033 }
1034
1035 if (PrtLvl >= PRINT_MORE) {
1036 ITS_COMPRES(computed_relres);
1037 ITS_REALRES(relres);
1038 }
1039
1040 } /* end of convergence check */
1041
1042 /* compute residual vector and continue loop */
1043 for (j = i; j > 0; j--) {
1044 rs[j - 1] = -s[j - 1] * rs[j];
1045 rs[j] = c[j - 1] * rs[j];
1046 }
1047
1048 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
1049
1050 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1051
1052 if (i) {
1053 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
1054 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1055 }
1056
1057 //-----------------------------------//
1058 // compute the convergence rate //
1059 //-----------------------------------//
1060 cr = r_norm / r_norm_old;
1061
1062 } /* end of iteration while loop */
1063
1064FINISHED:
1065 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
1066
1067 /*-------------------------------------------
1068 * Free some stuff
1069 *------------------------------------------*/
1070 fasp_mem_free(work);
1071 work = NULL;
1072 fasp_mem_free(p);
1073 p = NULL;
1074 fasp_mem_free(hh);
1075 hh = NULL;
1076 fasp_mem_free(norms);
1077 norms = NULL;
1078
1079#if DEBUG_MODE > 0
1080 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1081#endif
1082
1083 if (iter >= MaxIt)
1084 return ERROR_SOLVER_MAXIT;
1085 else
1086 return iter;
1087}
1088
1116INT fasp_solver_dstr_pvgmres(dSTRmat* A, dvector* b, dvector* x, precond* pc,
1117 const REAL tol, const REAL abstol, const INT MaxIt,
1118 const SHORT restart, const SHORT StopType,
1119 const SHORT PrtLvl)
1120{
1121 const INT n = b->row;
1122 const INT MIN_ITER = 0;
1123 const REAL epsmac = SMALLREAL;
1124
1125 //--------------------------------------------//
1126 // Newly added parameters to monitor when //
1127 // to change the restart parameter //
1128 //--------------------------------------------//
1129 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
1130 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
1131
1132 // local variables
1133 INT iter = 0;
1134 int i, j, k; // must be signed! -zcs
1135
1136 REAL r_norm, r_normb, gamma, t;
1137 REAL absres0 = BIGREAL, absres = BIGREAL;
1138 REAL relres = BIGREAL, normu = BIGREAL;
1139
1140 REAL cr = 1.0; // convergence rate
1141 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
1142 INT d = 3; // reduction for restart parameter
1143 INT restart_max = restart; // upper bound for restart in each restart cycle
1144 INT restart_min =
1145 3; // lower bound for restart in each restart cycle (should be small)
1146
1147 INT Restart = restart; // real restart in some fixed restarted cycle
1148 INT Restart1 = Restart + 1;
1149 unsigned LONG worksize = (Restart + 4) * (Restart + n) + 1 - n;
1150
1151 // allocate temp memory (need about (restart+4)*n REAL numbers)
1152 REAL * c = NULL, *s = NULL, *rs = NULL;
1153 REAL * norms = NULL, *r = NULL, *w = NULL;
1154 REAL* work = NULL;
1155 REAL **p = NULL, **hh = NULL;
1156
1157 // Output some info for debuging
1158 if (PrtLvl > PRINT_NONE) printf("\nCalling VGMRes solver (STR) ...\n");
1159
1160#if DEBUG_MODE > 0
1161 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1162 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1163#endif
1164
1165 /* allocate memory and setup temp work space */
1166 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1167
1168 /* check whether memory is enough for GMRES */
1169 while ((work == NULL) && (Restart > 5)) {
1170 Restart = Restart - 5;
1171 worksize = (Restart + 4) * (Restart + n) + 1 - n;
1172 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1173 Restart1 = Restart + 1;
1174 }
1175
1176 if (work == NULL) {
1177 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1178 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1179 }
1180
1181 if (PrtLvl > PRINT_MIN && Restart < restart) {
1182 printf("### WARNING: vGMRES restart number set to %d!\n", Restart);
1183 }
1184
1185 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1186 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1187 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
1188
1189 r = work;
1190 w = r + n;
1191 rs = w + n;
1192 c = rs + Restart1;
1193 s = c + Restart;
1194
1195 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
1196
1197 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
1198
1199 // r = b-A*x
1200 fasp_darray_cp(n, b->val, p[0]);
1201 fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1202
1203 r_norm = fasp_blas_darray_norm2(n, p[0]);
1204
1205 // compute initial residuals
1206 switch (StopType) {
1207 case STOP_REL_RES:
1208 absres0 = MAX(SMALLREAL, r_norm);
1209 relres = r_norm / absres0;
1210 break;
1211 case STOP_REL_PRECRES:
1212 if (pc == NULL)
1213 fasp_darray_cp(n, p[0], r);
1214 else
1215 pc->fct(p[0], r, pc->data);
1216 r_normb = sqrt(fasp_blas_darray_dotprod(n, p[0], r));
1217 absres0 = MAX(SMALLREAL, r_normb);
1218 relres = r_normb / absres0;
1219 break;
1220 case STOP_MOD_REL_RES:
1221 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
1222 absres0 = r_norm;
1223 relres = absres0 / normu;
1224 break;
1225 default:
1226 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1227 goto FINISHED;
1228 }
1229
1230 // if initial residual is small, no need to iterate!
1231 if (relres < tol || absres0 < abstol) goto FINISHED;
1232
1233 // output iteration information if needed
1234 fasp_itinfo(PrtLvl, StopType, 0, relres, absres0, 0);
1235
1236 // store initial residual
1237 norms[0] = relres;
1238
1239 /* outer iteration cycle */
1240 while (iter < MaxIt) {
1241
1242 rs[0] = r_norm_old = r_norm;
1243
1244 t = 1.0 / r_norm;
1245
1246 fasp_blas_darray_ax(n, t, p[0]);
1247
1248 //-----------------------------------//
1249 // adjust the restart parameter //
1250 //-----------------------------------//
1251 if (cr > cr_max || iter == 0) {
1252 Restart = restart_max;
1253 } else if (cr < cr_min) {
1254 // Restart = Restart;
1255 } else {
1256 if (Restart - d > restart_min) {
1257 Restart -= d;
1258 } else {
1259 Restart = restart_max;
1260 }
1261 }
1262
1263 /* RESTART CYCLE (right-preconditioning) */
1264 i = 0;
1265 while (i < Restart && iter < MaxIt) {
1266
1267 i++;
1268 iter++;
1269
1270 /* apply preconditioner */
1271 if (pc == NULL)
1272 fasp_darray_cp(n, p[i - 1], r);
1273 else
1274 pc->fct(p[i - 1], r, pc->data);
1275
1276 fasp_blas_dstr_mxv(A, r, p[i]);
1277
1278 /* modified Gram_Schmidt */
1279 for (j = 0; j < i; j++) {
1280 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1281 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
1282 }
1283 t = fasp_blas_darray_norm2(n, p[i]);
1284 hh[i][i - 1] = t;
1285 if (t != 0.0) {
1286 t = 1.0 / t;
1287 fasp_blas_darray_ax(n, t, p[i]);
1288 }
1289
1290 for (j = 1; j < i; ++j) {
1291 t = hh[j - 1][i - 1];
1292 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
1293 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
1294 }
1295 t = hh[i][i - 1] * hh[i][i - 1];
1296 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
1297
1298 gamma = sqrt(t);
1299 if (gamma == 0.0) gamma = epsmac;
1300 c[i - 1] = hh[i - 1][i - 1] / gamma;
1301 s[i - 1] = hh[i][i - 1] / gamma;
1302 rs[i] = -s[i - 1] * rs[i - 1];
1303 rs[i - 1] = c[i - 1] * rs[i - 1];
1304 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
1305
1306 absres = r_norm = fabs(rs[i]);
1307
1308 relres = absres / absres0;
1309
1310 norms[iter] = relres;
1311
1312 // output iteration information if needed
1313 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1314 norms[iter] / norms[iter - 1]);
1315
1316 // should we exit restart cycle
1317 if (relres < tol && iter >= MIN_ITER) break;
1318
1319 } /* end of restart cycle */
1320
1321 /* now compute solution, first solve upper triangular system */
1322 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
1323 for (k = i - 2; k >= 0; k--) {
1324 t = 0.0;
1325 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
1326
1327 t += rs[k];
1328 rs[k] = t / hh[k][k];
1329 }
1330
1331 fasp_darray_cp(n, p[i - 1], w);
1332
1333 fasp_blas_darray_ax(n, rs[i - 1], w);
1334
1335 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1336
1337 /* apply preconditioner */
1338 if (pc == NULL)
1339 fasp_darray_cp(n, w, r);
1340 else
1341 pc->fct(w, r, pc->data);
1342
1343 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1344
1345 // Check: prevent false convergence
1346 if (relres < tol && iter >= MIN_ITER) {
1347
1348 REAL computed_relres = relres;
1349
1350 // compute current residual
1351 fasp_darray_cp(n, b->val, r);
1352 fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1353
1354 r_norm = fasp_blas_darray_norm2(n, r);
1355
1356 switch (StopType) {
1357 case STOP_REL_RES:
1358 absres = r_norm;
1359 relres = absres / absres0;
1360 break;
1361 case STOP_REL_PRECRES:
1362 if (pc == NULL)
1363 fasp_darray_cp(n, r, w);
1364 else
1365 pc->fct(r, w, pc->data);
1366 absres = sqrt(fasp_blas_darray_dotprod(n, w, r));
1367 relres = absres / absres0;
1368 break;
1369 case STOP_MOD_REL_RES:
1370 absres = r_norm;
1371 normu = MAX(SMALLREAL, fasp_blas_darray_norm2(n, x->val));
1372 relres = absres / normu;
1373 break;
1374 }
1375
1376 norms[iter] = relres;
1377
1378 if (relres < tol) {
1379 break;
1380 } else {
1381 // Need to restart
1382 fasp_darray_cp(n, r, p[0]);
1383 i = 0;
1384 }
1385
1386 if (PrtLvl >= PRINT_MORE) {
1387 ITS_COMPRES(computed_relres);
1388 ITS_REALRES(relres);
1389 }
1390
1391 } /* end of convergence check */
1392
1393 /* compute residual vector and continue loop */
1394 for (j = i; j > 0; j--) {
1395 rs[j - 1] = -s[j - 1] * rs[j];
1396 rs[j] = c[j - 1] * rs[j];
1397 }
1398
1399 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
1400
1401 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1402
1403 if (i) {
1404 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
1405 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1406 }
1407
1408 //-----------------------------------//
1409 // compute the convergence rate //
1410 //-----------------------------------//
1411 cr = r_norm / r_norm_old;
1412
1413 } /* end of iteration while loop */
1414
1415FINISHED:
1416 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, relres);
1417
1418 /*-------------------------------------------
1419 * Free some stuff
1420 *------------------------------------------*/
1421 fasp_mem_free(work);
1422 work = NULL;
1423 fasp_mem_free(p);
1424 p = NULL;
1425 fasp_mem_free(hh);
1426 hh = NULL;
1427 fasp_mem_free(norms);
1428 norms = NULL;
1429
1430#if DEBUG_MODE > 0
1431 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1432#endif
1433
1434 if (iter >= MaxIt)
1435 return ERROR_SOLVER_MAXIT;
1436 else
1437 return iter;
1438}
1439
1468INT fasp_solver_pvgmres(mxv_matfree* mf, dvector* b, dvector* x, precond* pc,
1469 const REAL tol, const REAL abstol, const INT MaxIt,
1470 SHORT restart, const SHORT StopType, const SHORT PrtLvl)
1471{
1472 const INT n = b->row;
1473 const INT min_iter = 0;
1474
1475 //--------------------------------------------//
1476 // Newly added parameters to monitor when //
1477 // to change the restart parameter //
1478 //--------------------------------------------//
1479 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
1480 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
1481
1482 // local variables
1483 INT iter = 0;
1484 int i, j, k; // must be signed! -zcs
1485
1486 REAL epsmac = SMALLREAL;
1487 REAL r_norm, b_norm, den_norm;
1488 REAL epsilon, gamma, t;
1489
1490 REAL * c = NULL, *s = NULL, *rs = NULL;
1491 REAL * norms = NULL, *r = NULL, *w = NULL;
1492 REAL **p = NULL, **hh = NULL;
1493 REAL* work = NULL;
1494
1495 REAL cr = 1.0; // convergence rate
1496 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
1497 INT d = 3; // reduction for restart parameter
1498 INT restart_max = restart; // upper bound for restart in each restart cycle
1499 INT restart_min = 3; // lower bound for restart in each restart cycle
1500
1501 INT Restart = restart; // real restart in some fixed restarted cycle
1502 INT Restart1 = Restart + 1;
1503 unsigned LONG worksize = (restart + 4) * (restart + n) + 1 - n;
1504
1505 // Output some info for debuging
1506 if (PrtLvl > PRINT_NONE) printf("\nCalling VGMRes solver (MatFree) ...\n");
1507
1508#if DEBUG_MODE > 0
1509 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1510 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1511#endif
1512
1513 /* allocate memory and setup temp work space */
1514 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1515
1516 /* check whether memory is enough for GMRES */
1517 while ((work == NULL) && (Restart > 5)) {
1518 Restart = Restart - 5;
1519 worksize = (Restart + 4) * (Restart + n) + 1 - n;
1520 work = (REAL*)fasp_mem_calloc(worksize, sizeof(REAL));
1521 Restart1 = Restart + 1;
1522 }
1523
1524 if (work == NULL) {
1525 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1526 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1527 }
1528
1529 if (PrtLvl > PRINT_MIN && Restart < restart) {
1530 printf("### WARNING: vGMRES restart number set to %d!\n", Restart);
1531 }
1532
1533 p = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1534 hh = (REAL**)fasp_mem_calloc(Restart1, sizeof(REAL*));
1535 norms = (REAL*)fasp_mem_calloc(MaxIt + 1, sizeof(REAL));
1536
1537 r = work;
1538 w = r + n;
1539 rs = w + n;
1540 c = rs + Restart1;
1541 s = c + Restart;
1542 for (i = 0; i < Restart1; i++) p[i] = s + Restart + i * n;
1543 for (i = 0; i < Restart1; i++) hh[i] = p[Restart] + n + i * Restart;
1544
1545 /* initialization */
1546 mf->fct(mf->data, x->val, p[0]);
1547 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, p[0]);
1548
1549 b_norm = fasp_blas_darray_norm2(n, b->val);
1550 r_norm = fasp_blas_darray_norm2(n, p[0]);
1551 norms[0] = r_norm;
1552
1553 if (PrtLvl >= PRINT_SOME) {
1554 ITS_PUTNORM("right-hand side", b_norm);
1555 ITS_PUTNORM("residual", r_norm);
1556 }
1557
1558 if (b_norm > 0.0)
1559 den_norm = b_norm;
1560 else
1561 den_norm = r_norm;
1562
1563 epsilon = tol * den_norm;
1564
1565 /* outer iteration cycle */
1566 while (iter < MaxIt) {
1567 rs[0] = r_norm;
1568 r_norm_old = r_norm;
1569 if (r_norm == 0.0) {
1570 fasp_mem_free(work);
1571 work = NULL;
1572 fasp_mem_free(p);
1573 p = NULL;
1574 fasp_mem_free(hh);
1575 hh = NULL;
1576 fasp_mem_free(norms);
1577 norms = NULL;
1578 return iter;
1579 }
1580
1581 //-----------------------------------//
1582 // adjust the restart parameter //
1583 //-----------------------------------//
1584
1585 if (cr > cr_max || iter == 0) {
1586 Restart = restart_max;
1587 } else if (cr < cr_min) {
1588 // Restart = Restart;
1589 } else {
1590 if (Restart - d > restart_min) {
1591 Restart -= d;
1592 } else {
1593 Restart = restart_max;
1594 }
1595 }
1596
1597 if (r_norm <= epsilon && iter >= min_iter) {
1598 mf->fct(mf->data, x->val, r);
1599 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1600 r_norm = fasp_blas_darray_norm2(n, r);
1601
1602 if (r_norm <= epsilon) {
1603 break;
1604 } else {
1605 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1606 }
1607 }
1608
1609 t = 1.0 / r_norm;
1610
1611 // for (j = 0; j < n; j ++) p[0][j] *= t;
1612 fasp_blas_darray_ax(n, t, p[0]);
1613
1614 /* RESTART CYCLE (right-preconditioning) */
1615 i = 0;
1616 while (i < Restart && iter < MaxIt) {
1617
1618 i++;
1619 iter++;
1620
1621 /* apply preconditioner */
1622 if (pc == NULL)
1623 fasp_darray_cp(n, p[i - 1], r);
1624 else
1625 pc->fct(p[i - 1], r, pc->data);
1626
1627 mf->fct(mf->data, r, p[i]);
1628
1629 /* modified Gram_Schmidt */
1630 for (j = 0; j < i; j++) {
1631 hh[j][i - 1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1632 fasp_blas_darray_axpy(n, -hh[j][i - 1], p[j], p[i]);
1633 }
1634 t = fasp_blas_darray_norm2(n, p[i]);
1635 hh[i][i - 1] = t;
1636 if (t != 0.0) {
1637 t = 1.0 / t;
1638 // for (j = 0; j < n; j ++) p[i][j] *= t;
1639 fasp_blas_darray_ax(n, t, p[i]);
1640 }
1641
1642 for (j = 1; j < i; ++j) {
1643 t = hh[j - 1][i - 1];
1644 hh[j - 1][i - 1] = s[j - 1] * hh[j][i - 1] + c[j - 1] * t;
1645 hh[j][i - 1] = -s[j - 1] * t + c[j - 1] * hh[j][i - 1];
1646 }
1647 t = hh[i][i - 1] * hh[i][i - 1];
1648 t += hh[i - 1][i - 1] * hh[i - 1][i - 1];
1649 gamma = sqrt(t);
1650 if (gamma == 0.0) gamma = epsmac;
1651 c[i - 1] = hh[i - 1][i - 1] / gamma;
1652 s[i - 1] = hh[i][i - 1] / gamma;
1653 rs[i] = -s[i - 1] * rs[i - 1];
1654 rs[i - 1] = c[i - 1] * rs[i - 1];
1655 hh[i - 1][i - 1] = s[i - 1] * hh[i][i - 1] + c[i - 1] * hh[i - 1][i - 1];
1656 r_norm = fabs(rs[i]);
1657
1658 norms[iter] = r_norm;
1659
1660 if (b_norm > 0) {
1661 fasp_itinfo(PrtLvl, StopType, iter, norms[iter] / b_norm, norms[iter],
1662 norms[iter] / norms[iter - 1]);
1663 } else {
1664 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
1665 norms[iter] / norms[iter - 1]);
1666 }
1667
1668 /* should we exit restart cycle? */
1669 if (r_norm <= epsilon && iter >= min_iter) break;
1670
1671 } /* end of restart cycle */
1672
1673 /* now compute solution, first solve upper triangular system */
1674
1675 rs[i - 1] = rs[i - 1] / hh[i - 1][i - 1];
1676 for (k = i - 2; k >= 0; k--) {
1677 t = 0.0;
1678 for (j = k + 1; j < i; j++) t -= hh[k][j] * rs[j];
1679
1680 t += rs[k];
1681 rs[k] = t / hh[k][k];
1682 }
1683 fasp_darray_cp(n, p[i - 1], w);
1684 // for (j = 0; j < n; j ++) w[j] *= rs[i-1];
1685 fasp_blas_darray_ax(n, rs[i - 1], w);
1686 for (j = i - 2; j >= 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1687
1688 /* apply preconditioner */
1689 if (pc == NULL)
1690 fasp_darray_cp(n, w, r);
1691 else
1692 pc->fct(w, r, pc->data);
1693
1694 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1695
1696 if (r_norm <= epsilon && iter >= min_iter) {
1697 mf->fct(mf->data, x->val, r);
1698 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1699 r_norm = fasp_blas_darray_norm2(n, r);
1700
1701 if (r_norm <= epsilon) {
1702 break;
1703 } else {
1704 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1705 fasp_darray_cp(n, r, p[0]);
1706 i = 0;
1707 }
1708 } /* end of convergence check */
1709
1710 /* compute residual vector and continue loop */
1711 for (j = i; j > 0; j--) {
1712 rs[j - 1] = -s[j - 1] * rs[j];
1713 rs[j] = c[j - 1] * rs[j];
1714 }
1715
1716 if (i) fasp_blas_darray_axpy(n, rs[i] - 1.0, p[i], p[i]);
1717
1718 for (j = i - 1; j > 0; j--) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1719
1720 if (i) {
1721 fasp_blas_darray_axpy(n, rs[0] - 1.0, p[0], p[0]);
1722 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1723 }
1724
1725 //-----------------------------------//
1726 // compute the convergence rate //
1727 //-----------------------------------//
1728 cr = r_norm / r_norm_old;
1729
1730 } /* end of iteration while loop */
1731
1732 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter, MaxIt, r_norm);
1733
1734 /*-------------------------------------------
1735 * Free some stuff
1736 *------------------------------------------*/
1737 fasp_mem_free(work);
1738 work = NULL;
1739 fasp_mem_free(p);
1740 p = NULL;
1741 fasp_mem_free(hh);
1742 hh = NULL;
1743 fasp_mem_free(norms);
1744 norms = NULL;
1745
1746#if DEBUG_MODE > 0
1747 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1748#endif
1749
1750 if (iter >= MaxIt)
1751 return ERROR_SOLVER_MAXIT;
1752 else
1753 return iter;
1754}
1755
1756/*---------------------------------*/
1757/*-- End of File --*/
1758/*---------------------------------*/
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_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_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_dblc_pvgmres
INT fasp_solver_dblc_pvgmres(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)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:765
fasp_solver_dcsr_pvgmres
INT fasp_solver_dcsr_pvgmres(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 in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:66
fasp_solver_pvgmres
INT fasp_solver_pvgmres(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, const REAL tol, const REAL abstol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Solve "Ax=b" using PGMRES(right preconditioned) iterative method in which the restart parameter can b...
Definition: KryPvgmres.c:1468
fasp_solver_dbsr_pvgmres
INT fasp_solver_dbsr_pvgmres(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)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:416
fasp_solver_dstr_pvgmres
INT fasp_solver_dstr_pvgmres(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)
Right preconditioned GMRES method in which the restart parameter can be adaptively modified during it...
Definition: KryPvgmres.c:1116
fasp.h
Main header file for the FASP project.
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
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
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