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