Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KrySPvgmres.c
Go to the documentation of this file.
1
27#include <math.h>
28
29#include "fasp.h"
30#include "fasp_functs.h"
31
32/*---------------------------------*/
33/*-- Declare Private Functions --*/
34/*---------------------------------*/
35
36#include "KryUtil.inl"
37
38/*---------------------------------*/
39/*-- Public Functions --*/
40/*---------------------------------*/
41
68INT fasp_solver_dcsr_spvgmres (const dCSRmat *A,
69 const dvector *b,
70 dvector *x,
71 precond *pc,
72 const REAL tol,
73 const INT MaxIt,
74 SHORT restart,
75 const SHORT StopType,
76 const SHORT PrtLvl)
77{
78 const INT n = b->row;
79 const INT MIN_ITER = 0;
80 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
81 const REAL epsmac = SMALLREAL;
82
83 //--------------------------------------------//
84 // Newly added parameters to monitor when //
85 // to change the restart parameter //
86 //--------------------------------------------//
87 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
88 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
89
90 // local variables
91 INT iter = 0;
92 INT restart1 = restart + 1;
93 int i, j, k; // must be signed! -zcs
94
95 REAL r_norm, r_normb, gamma, t;
96 REAL normr0 = BIGREAL, absres = BIGREAL;
97 REAL relres = BIGREAL, normu = BIGREAL;
98
99 REAL cr = 1.0; // convergence rate
100 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
101 INT d = 3; // reduction for restart parameter
102 INT restart_max = restart; // upper bound for restart in each restart cycle
103 INT restart_min = 3; // lower bound for restart (should be small)
104 INT Restart = restart; // real restart in some fixed restarted cycle
105
106 INT iter_best = 0; // initial best known iteration
107 REAL absres_best = BIGREAL; // initial best known residual
108
109 // allocate temp memory (need about (restart+4)*n REAL numbers)
110 REAL *c = NULL, *s = NULL, *rs = NULL;
111 REAL *norms = NULL, *r = NULL, *w = NULL;
112 REAL *work = NULL, *x_best = NULL;
113 REAL **p = NULL, **hh = NULL;
114
115 // Output some info for debuging
116 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe VGMRes solver (CSR) ...\n");
117
118#if DEBUG_MODE > 0
119 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
120 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
121#endif
122
123 /* allocate memory and setup temp work space */
124 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
125
126 /* check whether memory is enough for GMRES */
127 while ( (work == NULL) && (restart > 5 ) ) {
128 restart = restart - 5 ;
129 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
130 printf("### WARNING: vGMRES restart number set to %d!\n", restart);
131 restart1 = restart + 1;
132 }
133
134 if ( work == NULL ) {
135 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
136 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
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;
144 x_best = c + restart; s = x_best + n;
145
146 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
147
148 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
149
150 // r = b-A*x
151 fasp_darray_cp(n, b->val, p[0]);
152 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
153
154 r_norm = fasp_blas_darray_norm2(n, p[0]);
155
156 // compute initial residuals
157 switch (StopType) {
158 case STOP_REL_RES:
159 normr0 = MAX(SMALLREAL,r_norm);
160 relres = r_norm/normr0;
161 break;
162 case STOP_REL_PRECRES:
163 if ( pc == NULL )
164 fasp_darray_cp(n, p[0], r);
165 else
166 pc->fct(p[0], r, pc->data);
167 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
168 normr0 = MAX(SMALLREAL,r_normb);
169 relres = r_normb/normr0;
170 break;
171 case STOP_MOD_REL_RES:
172 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
173 normr0 = r_norm;
174 relres = normr0/normu;
175 break;
176 default:
177 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
178 goto FINISHED;
179 }
180
181 // if initial residual is small, no need to iterate!
182 if ( relres < tol ) goto FINISHED;
183
184 // output iteration information if needed
185 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
186
187 // store initial residual
188 norms[0] = relres;
189
190 /* outer iteration cycle */
191 while ( iter < MaxIt ) {
192
193 rs[0] = r_norm_old = r_norm;
194
195 t = 1.0 / r_norm;
196
197 fasp_blas_darray_ax(n, t, p[0]);
198
199 //-----------------------------------//
200 // adjust the restart parameter //
201 //-----------------------------------//
202 if ( cr > cr_max || iter == 0 ) {
203 Restart = restart_max;
204 }
205 else if ( cr < cr_min ) {
206 // Restart = Restart;
207 }
208 else {
209 if ( Restart - d > restart_min ) {
210 Restart -= d;
211 }
212 else {
213 Restart = restart_max;
214 }
215 }
216
217 /* RESTART CYCLE (right-preconditioning) */
218 i = 0;
219 while ( i < Restart && iter < MaxIt ) {
220
221 i++; iter++;
222
223 /* apply preconditioner */
224 if (pc == NULL)
225 fasp_darray_cp(n, p[i-1], r);
226 else
227 pc->fct(p[i-1], r, pc->data);
228
229 fasp_blas_dcsr_mxv(A, r, p[i]);
230
231 /* modified Gram_Schmidt */
232 for (j = 0; j < i; j ++) {
233 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
234 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
235 }
236 t = fasp_blas_darray_norm2(n, p[i]);
237 hh[i][i-1] = t;
238 if (t != 0.0) {
239 t = 1.0/t;
240 fasp_blas_darray_ax(n, t, p[i]);
241 }
242
243 for (j = 1; j < i; ++j) {
244 t = hh[j-1][i-1];
245 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
246 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
247 }
248 t= hh[i][i-1]*hh[i][i-1];
249 t+= hh[i-1][i-1]*hh[i-1][i-1];
250
251 gamma = sqrt(t);
252 if (gamma == 0.0) gamma = epsmac;
253 c[i-1] = hh[i-1][i-1] / gamma;
254 s[i-1] = hh[i][i-1] / gamma;
255 rs[i] = -s[i-1]*rs[i-1];
256 rs[i-1] = c[i-1]*rs[i-1];
257 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
258
259 absres = r_norm = fabs(rs[i]);
260
261 relres = absres/normr0;
262
263 norms[iter] = relres;
264
265 // output iteration information if needed
266 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
267 norms[iter]/norms[iter-1]);
268
269 // should we exit restart cycle
270 if ( relres <= tol && iter >= MIN_ITER ) break;
271
272 } /* end of restart cycle */
273
274 /* now compute solution, first solve upper triangular system */
275 rs[i-1] = rs[i-1] / hh[i-1][i-1];
276 for (k = i-2; k >= 0; k --) {
277 t = 0.0;
278 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
279
280 t += rs[k];
281 rs[k] = t / hh[k][k];
282 }
283
284 fasp_darray_cp(n, p[i-1], w);
285
286 fasp_blas_darray_ax(n, rs[i-1], w);
287
288 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
289
290 /* apply preconditioner */
291 if ( pc == NULL )
292 fasp_darray_cp(n, w, r);
293 else
294 pc->fct(w, r, pc->data);
295
296 fasp_blas_darray_axpy(n, 1.0, r, x->val);
297
298 // safety net check: save the best-so-far solution
299 if ( fasp_dvec_isnan(x) ) {
300 // If the solution is NAN, restore the best solution
301 absres = BIGREAL;
302 goto RESTORE_BESTSOL;
303 }
304
305 if ( absres < absres_best - maxdiff) {
306 absres_best = absres;
307 iter_best = iter;
308 fasp_darray_cp(n,x->val,x_best);
309 }
310
311 // Check: prevent false convergence
312 if ( relres <= tol && iter >= MIN_ITER ) {
313
314 fasp_darray_cp(n, b->val, r);
315 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
316
317 r_norm = fasp_blas_darray_norm2(n, r);
318
319 switch ( StopType ) {
320 case STOP_REL_RES:
321 absres = r_norm;
322 relres = absres/normr0;
323 break;
324 case STOP_REL_PRECRES:
325 if ( pc == NULL )
326 fasp_darray_cp(n, r, w);
327 else
328 pc->fct(r, w, pc->data);
329 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
330 relres = absres/normr0;
331 break;
332 case STOP_MOD_REL_RES:
333 absres = r_norm;
334 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
335 relres = absres/normu;
336 break;
337 }
338
339 norms[iter] = relres;
340
341 if ( relres <= tol ) {
342 break;
343 }
344 else {
345 // Need to restart
346 fasp_darray_cp(n, r, p[0]); i = 0;
347 }
348
349 } /* end of convergence check */
350
351 /* compute residual vector and continue loop */
352 for ( j = i; j > 0; j-- ) {
353 rs[j-1] = -s[j-1]*rs[j];
354 rs[j] = c[j-1]*rs[j];
355 }
356
357 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
358
359 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
360
361 if ( i ) {
362 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
363 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
364 }
365
366 //-----------------------------------//
367 // compute the convergence rate //
368 //-----------------------------------//
369 cr = r_norm / r_norm_old;
370
371 } /* end of iteration while loop */
372
373RESTORE_BESTSOL: // restore the best-so-far solution if necessary
374 if ( iter != iter_best ) {
375
376 // compute best residual
377 fasp_darray_cp(n,b->val,r);
378 fasp_blas_dcsr_aAxpy(-1.0,A,x_best,r);
379
380 switch ( StopType ) {
381 case STOP_REL_RES:
382 absres_best = fasp_blas_darray_norm2(n,r);
383 break;
384 case STOP_REL_PRECRES:
385 // z = B(r)
386 if ( pc != NULL )
387 pc->fct(r,w,pc->data); /* Apply preconditioner */
388 else
389 fasp_darray_cp(n,r,w); /* No preconditioner */
390 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
391 break;
392 case STOP_MOD_REL_RES:
393 absres_best = fasp_blas_darray_norm2(n,r);
394 break;
395 }
396
397 if ( absres > absres_best + maxdiff || isnan(absres) ) {
398 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
399 fasp_darray_cp(n,x_best,x->val);
400 relres = absres_best / normr0;
401 }
402 }
403
404FINISHED:
405 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
406
407 /*-------------------------------------------
408 * Free some stuff
409 *------------------------------------------*/
410 fasp_mem_free(work); work = NULL;
411 fasp_mem_free(p); p = NULL;
412 fasp_mem_free(hh); hh = NULL;
413 fasp_mem_free(norms); norms = NULL;
414
415#if DEBUG_MODE > 0
416 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
417#endif
418
419 if (iter>=MaxIt)
420 return ERROR_SOLVER_MAXIT;
421 else
422 return iter;
423}
424
449INT fasp_solver_dbsr_spvgmres (const dBSRmat *A,
450 const dvector *b,
451 dvector *x,
452 precond *pc,
453 const REAL tol,
454 const INT MaxIt,
455 SHORT restart,
456 const SHORT StopType,
457 const SHORT PrtLvl)
458{
459 const INT n = b->row;
460 const INT MIN_ITER = 0;
461 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
462 const REAL epsmac = SMALLREAL;
463
464 //--------------------------------------------//
465 // Newly added parameters to monitor when //
466 // to change the restart parameter //
467 //--------------------------------------------//
468 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
469 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
470
471 // local variables
472 INT iter = 0;
473 INT restart1 = restart + 1;
474 int i, j, k; // must be signed! -zcs
475
476 REAL r_norm, r_normb, gamma, t;
477 REAL normr0 = BIGREAL, absres = BIGREAL;
478 REAL relres = BIGREAL, normu = BIGREAL;
479
480 REAL cr = 1.0; // convergence rate
481 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
482 INT d = 3; // reduction for restart parameter
483 INT restart_max = restart; // upper bound for restart in each restart cycle
484 INT restart_min = 3; // lower bound for restart (should be small)
485 INT Restart = restart; // real restart in some fixed restarted cycle
486
487 INT iter_best = 0; // initial best known iteration
488 REAL absres_best = BIGREAL; // initial best known residual
489
490 // allocate temp memory (need about (restart+4)*n REAL numbers)
491 REAL *c = NULL, *s = NULL, *rs = NULL;
492 REAL *norms = NULL, *r = NULL, *w = NULL;
493 REAL *work = NULL, *x_best = NULL;
494 REAL **p = NULL, **hh = NULL;
495
496 // Output some info for debuging
497 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe VGMRes solver (BSR) ...\n");
498
499#if DEBUG_MODE > 0
500 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
501 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
502#endif
503
504 /* allocate memory and setup temp work space */
505 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
506
507 /* check whether memory is enough for GMRES */
508 while ( (work == NULL) && (restart > 5 ) ) {
509 restart = restart - 5 ;
510 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
511 printf("### WARNING: vGMRES restart number set to %d!\n", restart);
512 restart1 = restart + 1;
513 }
514
515 if ( work == NULL ) {
516 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
517 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
518 }
519
520 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
521 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
522 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
523
524 r = work; w = r + n; rs = w + n; c = rs + restart1;
525 x_best = c + restart; s = x_best + n;
526
527 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
528
529 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
530
531 // r = b-A*x
532 fasp_darray_cp(n, b->val, p[0]);
533 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
534
535 r_norm = fasp_blas_darray_norm2(n, p[0]);
536
537 // compute initial residuals
538 switch (StopType) {
539 case STOP_REL_RES:
540 normr0 = MAX(SMALLREAL,r_norm);
541 relres = r_norm/normr0;
542 break;
543 case STOP_REL_PRECRES:
544 if ( pc == NULL )
545 fasp_darray_cp(n, p[0], r);
546 else
547 pc->fct(p[0], r, pc->data);
548 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
549 normr0 = MAX(SMALLREAL,r_normb);
550 relres = r_normb/normr0;
551 break;
552 case STOP_MOD_REL_RES:
553 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
554 normr0 = r_norm;
555 relres = normr0/normu;
556 break;
557 default:
558 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
559 goto FINISHED;
560 }
561
562 // if initial residual is small, no need to iterate!
563 if ( relres < tol ) goto FINISHED;
564
565 // output iteration information if needed
566 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
567
568 // store initial residual
569 norms[0] = relres;
570
571 /* outer iteration cycle */
572 while ( iter < MaxIt ) {
573
574 rs[0] = r_norm_old = r_norm;
575
576 t = 1.0 / r_norm;
577
578 fasp_blas_darray_ax(n, t, p[0]);
579
580 //-----------------------------------//
581 // adjust the restart parameter //
582 //-----------------------------------//
583 if ( cr > cr_max || iter == 0 ) {
584 Restart = restart_max;
585 }
586 else if ( cr < cr_min ) {
587 // Restart = Restart;
588 }
589 else {
590 if ( Restart - d > restart_min ) {
591 Restart -= d;
592 }
593 else {
594 Restart = restart_max;
595 }
596 }
597
598 /* RESTART CYCLE (right-preconditioning) */
599 i = 0;
600 while ( i < Restart && iter < MaxIt ) {
601
602 i++; iter++;
603
604 /* apply preconditioner */
605 if (pc == NULL)
606 fasp_darray_cp(n, p[i-1], r);
607 else
608 pc->fct(p[i-1], r, pc->data);
609
610 fasp_blas_dbsr_mxv(A, r, p[i]);
611
612 /* modified Gram_Schmidt */
613 for (j = 0; j < i; j ++) {
614 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
615 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
616 }
617 t = fasp_blas_darray_norm2(n, p[i]);
618 hh[i][i-1] = t;
619 if (t != 0.0) {
620 t = 1.0/t;
621 fasp_blas_darray_ax(n, t, p[i]);
622 }
623
624 for (j = 1; j < i; ++j) {
625 t = hh[j-1][i-1];
626 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
627 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
628 }
629 t= hh[i][i-1]*hh[i][i-1];
630 t+= hh[i-1][i-1]*hh[i-1][i-1];
631
632 gamma = sqrt(t);
633 if (gamma == 0.0) gamma = epsmac;
634 c[i-1] = hh[i-1][i-1] / gamma;
635 s[i-1] = hh[i][i-1] / gamma;
636 rs[i] = -s[i-1]*rs[i-1];
637 rs[i-1] = c[i-1]*rs[i-1];
638 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
639
640 absres = r_norm = fabs(rs[i]);
641
642 relres = absres/normr0;
643
644 norms[iter] = relres;
645
646 // output iteration information if needed
647 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
648 norms[iter]/norms[iter-1]);
649
650 // should we exit restart cycle
651 if ( relres <= tol && iter >= MIN_ITER ) break;
652
653 } /* end of restart cycle */
654
655 /* now compute solution, first solve upper triangular system */
656 rs[i-1] = rs[i-1] / hh[i-1][i-1];
657 for (k = i-2; k >= 0; k --) {
658 t = 0.0;
659 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
660
661 t += rs[k];
662 rs[k] = t / hh[k][k];
663 }
664
665 fasp_darray_cp(n, p[i-1], w);
666
667 fasp_blas_darray_ax(n, rs[i-1], w);
668
669 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
670
671 /* apply preconditioner */
672 if ( pc == NULL )
673 fasp_darray_cp(n, w, r);
674 else
675 pc->fct(w, r, pc->data);
676
677 fasp_blas_darray_axpy(n, 1.0, r, x->val);
678
679 // safety net check: save the best-so-far solution
680 if ( fasp_dvec_isnan(x) ) {
681 // If the solution is NAN, restore the best solution
682 absres = BIGREAL;
683 goto RESTORE_BESTSOL;
684 }
685
686 if ( absres < absres_best - maxdiff) {
687 absres_best = absres;
688 iter_best = iter;
689 fasp_darray_cp(n,x->val,x_best);
690 }
691
692 // Check: prevent false convergence
693 if ( relres <= tol && iter >= MIN_ITER ) {
694
695 fasp_darray_cp(n, b->val, r);
696 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
697
698 r_norm = fasp_blas_darray_norm2(n, r);
699
700 switch ( StopType ) {
701 case STOP_REL_RES:
702 absres = r_norm;
703 relres = absres/normr0;
704 break;
705 case STOP_REL_PRECRES:
706 if ( pc == NULL )
707 fasp_darray_cp(n, r, w);
708 else
709 pc->fct(r, w, pc->data);
710 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
711 relres = absres/normr0;
712 break;
713 case STOP_MOD_REL_RES:
714 absres = r_norm;
715 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
716 relres = absres/normu;
717 break;
718 }
719
720 norms[iter] = relres;
721
722 if ( relres <= tol ) {
723 break;
724 }
725 else {
726 // Need to restart
727 fasp_darray_cp(n, r, p[0]); i = 0;
728 }
729
730 } /* end of convergence check */
731
732 /* compute residual vector and continue loop */
733 for ( j = i; j > 0; j-- ) {
734 rs[j-1] = -s[j-1]*rs[j];
735 rs[j] = c[j-1]*rs[j];
736 }
737
738 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
739
740 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
741
742 if ( i ) {
743 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
744 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
745 }
746
747 //-----------------------------------//
748 // compute the convergence rate //
749 //-----------------------------------//
750 cr = r_norm / r_norm_old;
751
752 } /* end of iteration while loop */
753
754RESTORE_BESTSOL: // restore the best-so-far solution if necessary
755 if ( iter != iter_best ) {
756
757 // compute best residual
758 fasp_darray_cp(n,b->val,r);
759 fasp_blas_dbsr_aAxpy(-1.0,A,x_best,r);
760
761 switch ( StopType ) {
762 case STOP_REL_RES:
763 absres_best = fasp_blas_darray_norm2(n,r);
764 break;
765 case STOP_REL_PRECRES:
766 // z = B(r)
767 if ( pc != NULL )
768 pc->fct(r,w,pc->data); /* Apply preconditioner */
769 else
770 fasp_darray_cp(n,r,w); /* No preconditioner */
771 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
772 break;
773 case STOP_MOD_REL_RES:
774 absres_best = fasp_blas_darray_norm2(n,r);
775 break;
776 }
777
778 if ( absres > absres_best + maxdiff || isnan(absres) ) {
779 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
780 fasp_darray_cp(n,x_best,x->val);
781 relres = absres_best / normr0;
782 }
783 }
784
785FINISHED:
786 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
787
788 /*-------------------------------------------
789 * Free some stuff
790 *------------------------------------------*/
791 fasp_mem_free(work); work = NULL;
792 fasp_mem_free(p); p = NULL;
793 fasp_mem_free(hh); hh = NULL;
794 fasp_mem_free(norms); norms = NULL;
795
796#if DEBUG_MODE > 0
797 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
798#endif
799
800 if (iter>=MaxIt)
801 return ERROR_SOLVER_MAXIT;
802 else
803 return iter;
804}
805
829INT fasp_solver_dblc_spvgmres (const dBLCmat *A,
830 const dvector *b,
831 dvector *x,
832 precond *pc,
833 const REAL tol,
834 const INT MaxIt,
835 SHORT restart,
836 const SHORT StopType,
837 const SHORT PrtLvl)
838{
839 const INT n = b->row;
840 const INT MIN_ITER = 0;
841 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
842 const REAL epsmac = SMALLREAL;
843
844 //--------------------------------------------//
845 // Newly added parameters to monitor when //
846 // to change the restart parameter //
847 //--------------------------------------------//
848 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
849 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
850
851 // local variables
852 INT iter = 0;
853 INT restart1 = restart + 1;
854 int i, j, k; // must be signed! -zcs
855
856 REAL r_norm, r_normb, gamma, t;
857 REAL normr0 = BIGREAL, absres = BIGREAL;
858 REAL relres = BIGREAL, normu = BIGREAL;
859
860 REAL cr = 1.0; // convergence rate
861 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
862 INT d = 3; // reduction for restart parameter
863 INT restart_max = restart; // upper bound for restart in each restart cycle
864 INT restart_min = 3; // lower bound for restart (should be small)
865 INT Restart = restart; // real restart in some fixed restarted cycle
866
867 INT iter_best = 0; // initial best known iteration
868 REAL absres_best = BIGREAL; // initial best known residual
869
870 // allocate temp memory (need about (restart+4)*n REAL numbers)
871 REAL *c = NULL, *s = NULL, *rs = NULL;
872 REAL *norms = NULL, *r = NULL, *w = NULL;
873 REAL *work = NULL, *x_best = NULL;
874 REAL **p = NULL, **hh = NULL;
875
876 // Output some info for debuging
877 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe VGMRes solver (BLC) ...\n");
878
879#if DEBUG_MODE > 0
880 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
881 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
882#endif
883
884 /* allocate memory and setup temp work space */
885 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
886
887 /* check whether memory is enough for GMRES */
888 while ( (work == NULL) && (restart > 5 ) ) {
889 restart = restart - 5 ;
890 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
891 printf("### WARNING: vGMRES restart number set to %d!\n", restart);
892 restart1 = restart + 1;
893 }
894
895 if ( work == NULL ) {
896 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
897 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
898 }
899
900 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
901 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
902 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
903
904 r = work; w = r + n; rs = w + n; c = rs + restart1;
905 x_best = c + restart; s = x_best + n;
906
907 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
908
909 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
910
911 // r = b-A*x
912 fasp_darray_cp(n, b->val, p[0]);
913 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
914
915 r_norm = fasp_blas_darray_norm2(n, p[0]);
916
917 // compute initial residuals
918 switch (StopType) {
919 case STOP_REL_RES:
920 normr0 = MAX(SMALLREAL,r_norm);
921 relres = r_norm/normr0;
922 break;
923 case STOP_REL_PRECRES:
924 if ( pc == NULL )
925 fasp_darray_cp(n, p[0], r);
926 else
927 pc->fct(p[0], r, pc->data);
928 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
929 normr0 = MAX(SMALLREAL,r_normb);
930 relres = r_normb/normr0;
931 break;
932 case STOP_MOD_REL_RES:
933 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
934 normr0 = r_norm;
935 relres = normr0/normu;
936 break;
937 default:
938 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
939 goto FINISHED;
940 }
941
942 // if initial residual is small, no need to iterate!
943 if ( relres < tol ) goto FINISHED;
944
945 // output iteration information if needed
946 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
947
948 // store initial residual
949 norms[0] = relres;
950
951 /* outer iteration cycle */
952 while ( iter < MaxIt ) {
953
954 rs[0] = r_norm_old = r_norm;
955
956 t = 1.0 / r_norm;
957
958 fasp_blas_darray_ax(n, t, p[0]);
959
960 //-----------------------------------//
961 // adjust the restart parameter //
962 //-----------------------------------//
963 if ( cr > cr_max || iter == 0 ) {
964 Restart = restart_max;
965 }
966 else if ( cr < cr_min ) {
967 // Restart = Restart;
968 }
969 else {
970 if ( Restart - d > restart_min ) {
971 Restart -= d;
972 }
973 else {
974 Restart = restart_max;
975 }
976 }
977
978 /* RESTART CYCLE (right-preconditioning) */
979 i = 0;
980 while ( i < Restart && iter < MaxIt ) {
981
982 i++; iter++;
983
984 /* apply preconditioner */
985 if (pc == NULL)
986 fasp_darray_cp(n, p[i-1], r);
987 else
988 pc->fct(p[i-1], r, pc->data);
989
990 fasp_blas_dblc_mxv(A, r, p[i]);
991
992 /* modified Gram_Schmidt */
993 for (j = 0; j < i; j ++) {
994 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
995 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
996 }
997 t = fasp_blas_darray_norm2(n, p[i]);
998 hh[i][i-1] = t;
999 if (t != 0.0) {
1000 t = 1.0/t;
1001 fasp_blas_darray_ax(n, t, p[i]);
1002 }
1003
1004 for (j = 1; j < i; ++j) {
1005 t = hh[j-1][i-1];
1006 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1007 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1008 }
1009 t= hh[i][i-1]*hh[i][i-1];
1010 t+= hh[i-1][i-1]*hh[i-1][i-1];
1011
1012 gamma = sqrt(t);
1013 if (gamma == 0.0) gamma = epsmac;
1014 c[i-1] = hh[i-1][i-1] / gamma;
1015 s[i-1] = hh[i][i-1] / gamma;
1016 rs[i] = -s[i-1]*rs[i-1];
1017 rs[i-1] = c[i-1]*rs[i-1];
1018 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1019
1020 absres = r_norm = fabs(rs[i]);
1021
1022 relres = absres/normr0;
1023
1024 norms[iter] = relres;
1025
1026 // output iteration information if needed
1027 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1028 norms[iter]/norms[iter-1]);
1029
1030 // should we exit restart cycle
1031 if ( relres <= tol && iter >= MIN_ITER ) break;
1032
1033 } /* end of restart cycle */
1034
1035 /* now compute solution, first solve upper triangular system */
1036 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1037 for (k = i-2; k >= 0; k --) {
1038 t = 0.0;
1039 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
1040
1041 t += rs[k];
1042 rs[k] = t / hh[k][k];
1043 }
1044
1045 fasp_darray_cp(n, p[i-1], w);
1046
1047 fasp_blas_darray_ax(n, rs[i-1], w);
1048
1049 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1050
1051 /* apply preconditioner */
1052 if ( pc == NULL )
1053 fasp_darray_cp(n, w, r);
1054 else
1055 pc->fct(w, r, pc->data);
1056
1057 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1058
1059 // safety net check: save the best-so-far solution
1060 if ( fasp_dvec_isnan(x) ) {
1061 // If the solution is NAN, restore the best solution
1062 absres = BIGREAL;
1063 goto RESTORE_BESTSOL;
1064 }
1065
1066 if ( absres < absres_best - maxdiff) {
1067 absres_best = absres;
1068 iter_best = iter;
1069 fasp_darray_cp(n,x->val,x_best);
1070 }
1071
1072 // Check: prevent false convergence
1073 if ( relres <= tol && iter >= MIN_ITER ) {
1074
1075 fasp_darray_cp(n, b->val, r);
1076 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
1077
1078 r_norm = fasp_blas_darray_norm2(n, r);
1079
1080 switch ( StopType ) {
1081 case STOP_REL_RES:
1082 absres = r_norm;
1083 relres = absres/normr0;
1084 break;
1085 case STOP_REL_PRECRES:
1086 if ( pc == NULL )
1087 fasp_darray_cp(n, r, w);
1088 else
1089 pc->fct(r, w, pc->data);
1090 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
1091 relres = absres/normr0;
1092 break;
1093 case STOP_MOD_REL_RES:
1094 absres = r_norm;
1095 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1096 relres = absres/normu;
1097 break;
1098 }
1099
1100 norms[iter] = relres;
1101
1102 if ( relres <= tol ) {
1103 break;
1104 }
1105 else {
1106 // Need to restart
1107 fasp_darray_cp(n, r, p[0]); i = 0;
1108 }
1109
1110 } /* end of convergence check */
1111
1112 /* compute residual vector and continue loop */
1113 for ( j = i; j > 0; j-- ) {
1114 rs[j-1] = -s[j-1]*rs[j];
1115 rs[j] = c[j-1]*rs[j];
1116 }
1117
1118 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1119
1120 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1121
1122 if ( i ) {
1123 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1124 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1125 }
1126
1127 //-----------------------------------//
1128 // compute the convergence rate //
1129 //-----------------------------------//
1130 cr = r_norm / r_norm_old;
1131
1132 } /* end of iteration while loop */
1133
1134RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1135 if ( iter != iter_best ) {
1136
1137 // compute best residual
1138 fasp_darray_cp(n,b->val,r);
1139 fasp_blas_dblc_aAxpy(-1.0,A,x_best,r);
1140
1141 switch ( StopType ) {
1142 case STOP_REL_RES:
1143 absres_best = fasp_blas_darray_norm2(n,r);
1144 break;
1145 case STOP_REL_PRECRES:
1146 // z = B(r)
1147 if ( pc != NULL )
1148 pc->fct(r,w,pc->data); /* Apply preconditioner */
1149 else
1150 fasp_darray_cp(n,r,w); /* No preconditioner */
1151 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
1152 break;
1153 case STOP_MOD_REL_RES:
1154 absres_best = fasp_blas_darray_norm2(n,r);
1155 break;
1156 }
1157
1158 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1159 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1160 fasp_darray_cp(n,x_best,x->val);
1161 relres = absres_best / normr0;
1162 }
1163 }
1164
1165FINISHED:
1166 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1167
1168 /*-------------------------------------------
1169 * Free some stuff
1170 *------------------------------------------*/
1171 fasp_mem_free(work); work = NULL;
1172 fasp_mem_free(p); p = NULL;
1173 fasp_mem_free(hh); hh = NULL;
1174 fasp_mem_free(norms); norms = NULL;
1175
1176#if DEBUG_MODE > 0
1177 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1178#endif
1179
1180 if (iter>=MaxIt)
1181 return ERROR_SOLVER_MAXIT;
1182 else
1183 return iter;
1184}
1185
1210INT fasp_solver_dstr_spvgmres (const dSTRmat *A,
1211 const dvector *b,
1212 dvector *x,
1213 precond *pc,
1214 const REAL tol,
1215 const INT MaxIt,
1216 SHORT restart,
1217 const SHORT StopType,
1218 const SHORT PrtLvl)
1219{
1220 const INT n = b->row;
1221 const INT MIN_ITER = 0;
1222 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1223 const REAL epsmac = SMALLREAL;
1224
1225 //--------------------------------------------//
1226 // Newly added parameters to monitor when //
1227 // to change the restart parameter //
1228 //--------------------------------------------//
1229 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
1230 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
1231
1232 // local variables
1233 INT iter = 0;
1234 INT restart1 = restart + 1;
1235 int i, j, k; // must be signed! -zcs
1236
1237 REAL r_norm, r_normb, gamma, t;
1238 REAL normr0 = BIGREAL, absres = BIGREAL;
1239 REAL relres = BIGREAL, normu = BIGREAL;
1240
1241 REAL cr = 1.0; // convergence rate
1242 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
1243 INT d = 3; // reduction for restart parameter
1244 INT restart_max = restart; // upper bound for restart in each restart cycle
1245 INT restart_min = 3; // lower bound for restart (should be small)
1246 INT Restart = restart; // real restart in some fixed restarted cycle
1247
1248 INT iter_best = 0; // initial best known iteration
1249 REAL absres_best = BIGREAL; // initial best known residual
1250
1251 // allocate temp memory (need about (restart+4)*n REAL numbers)
1252 REAL *c = NULL, *s = NULL, *rs = NULL;
1253 REAL *norms = NULL, *r = NULL, *w = NULL;
1254 REAL *work = NULL, *x_best = NULL;
1255 REAL **p = NULL, **hh = NULL;
1256
1257 // Output some info for debuging
1258 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe VGMRes solver (STR) ...\n");
1259
1260#if DEBUG_MODE > 0
1261 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1262 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1263#endif
1264
1265 /* allocate memory and setup temp work space */
1266 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1267
1268 /* check whether memory is enough for GMRES */
1269 while ( (work == NULL) && (restart > 5 ) ) {
1270 restart = restart - 5 ;
1271 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1272 printf("### WARNING: vGMRES restart number set to %d!\n", restart);
1273 restart1 = restart + 1;
1274 }
1275
1276 if ( work == NULL ) {
1277 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1278 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1279 }
1280
1281 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1282 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1283 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1284
1285 r = work; w = r + n; rs = w + n; c = rs + restart1;
1286 x_best = c + restart; s = x_best + n;
1287
1288 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
1289
1290 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
1291
1292 // r = b-A*x
1293 fasp_darray_cp(n, b->val, p[0]);
1294 fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1295
1296 r_norm = fasp_blas_darray_norm2(n, p[0]);
1297
1298 // compute initial residuals
1299 switch (StopType) {
1300 case STOP_REL_RES:
1301 normr0 = MAX(SMALLREAL,r_norm);
1302 relres = r_norm/normr0;
1303 break;
1304 case STOP_REL_PRECRES:
1305 if ( pc == NULL )
1306 fasp_darray_cp(n, p[0], r);
1307 else
1308 pc->fct(p[0], r, pc->data);
1309 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
1310 normr0 = MAX(SMALLREAL,r_normb);
1311 relres = r_normb/normr0;
1312 break;
1313 case STOP_MOD_REL_RES:
1314 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1315 normr0 = r_norm;
1316 relres = normr0/normu;
1317 break;
1318 default:
1319 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1320 goto FINISHED;
1321 }
1322
1323 // if initial residual is small, no need to iterate!
1324 if ( relres < tol ) goto FINISHED;
1325
1326 // output iteration information if needed
1327 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
1328
1329 // store initial residual
1330 norms[0] = relres;
1331
1332 /* outer iteration cycle */
1333 while ( iter < MaxIt ) {
1334
1335 rs[0] = r_norm_old = r_norm;
1336
1337 t = 1.0 / r_norm;
1338
1339 fasp_blas_darray_ax(n, t, p[0]);
1340
1341 //-----------------------------------//
1342 // adjust the restart parameter //
1343 //-----------------------------------//
1344 if ( cr > cr_max || iter == 0 ) {
1345 Restart = restart_max;
1346 }
1347 else if ( cr < cr_min ) {
1348 // Restart = Restart;
1349 }
1350 else {
1351 if ( Restart - d > restart_min ) {
1352 Restart -= d;
1353 }
1354 else {
1355 Restart = restart_max;
1356 }
1357 }
1358
1359 /* RESTART CYCLE (right-preconditioning) */
1360 i = 0;
1361 while ( i < Restart && iter < MaxIt ) {
1362
1363 i++; iter++;
1364
1365 /* apply preconditioner */
1366 if (pc == NULL)
1367 fasp_darray_cp(n, p[i-1], r);
1368 else
1369 pc->fct(p[i-1], r, pc->data);
1370
1371 fasp_blas_dstr_mxv(A, r, p[i]);
1372
1373 /* modified Gram_Schmidt */
1374 for (j = 0; j < i; j ++) {
1375 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1376 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
1377 }
1378 t = fasp_blas_darray_norm2(n, p[i]);
1379 hh[i][i-1] = t;
1380 if (t != 0.0) {
1381 t = 1.0/t;
1382 fasp_blas_darray_ax(n, t, p[i]);
1383 }
1384
1385 for (j = 1; j < i; ++j) {
1386 t = hh[j-1][i-1];
1387 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1388 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1389 }
1390 t= hh[i][i-1]*hh[i][i-1];
1391 t+= hh[i-1][i-1]*hh[i-1][i-1];
1392
1393 gamma = sqrt(t);
1394 if (gamma == 0.0) gamma = epsmac;
1395 c[i-1] = hh[i-1][i-1] / gamma;
1396 s[i-1] = hh[i][i-1] / gamma;
1397 rs[i] = -s[i-1]*rs[i-1];
1398 rs[i-1] = c[i-1]*rs[i-1];
1399 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1400
1401 absres = r_norm = fabs(rs[i]);
1402
1403 relres = absres/normr0;
1404
1405 norms[iter] = relres;
1406
1407 // output iteration information if needed
1408 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1409 norms[iter]/norms[iter-1]);
1410
1411 // should we exit restart cycle
1412 if ( relres <= tol && iter >= MIN_ITER ) break;
1413
1414 } /* end of restart cycle */
1415
1416 /* now compute solution, first solve upper triangular system */
1417 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1418 for (k = i-2; k >= 0; k --) {
1419 t = 0.0;
1420 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
1421
1422 t += rs[k];
1423 rs[k] = t / hh[k][k];
1424 }
1425
1426 fasp_darray_cp(n, p[i-1], w);
1427
1428 fasp_blas_darray_ax(n, rs[i-1], w);
1429
1430 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1431
1432 /* apply preconditioner */
1433 if ( pc == NULL )
1434 fasp_darray_cp(n, w, r);
1435 else
1436 pc->fct(w, r, pc->data);
1437
1438 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1439
1440 // safety net check: save the best-so-far solution
1441 if ( fasp_dvec_isnan(x) ) {
1442 // If the solution is NAN, restore the best solution
1443 absres = BIGREAL;
1444 goto RESTORE_BESTSOL;
1445 }
1446
1447 if ( absres < absres_best - maxdiff) {
1448 absres_best = absres;
1449 iter_best = iter;
1450 fasp_darray_cp(n,x->val,x_best);
1451 }
1452
1453 // Check: prevent false convergence
1454 if ( relres <= tol && iter >= MIN_ITER ) {
1455
1456 fasp_darray_cp(n, b->val, r);
1457 fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1458
1459 r_norm = fasp_blas_darray_norm2(n, r);
1460
1461 switch ( StopType ) {
1462 case STOP_REL_RES:
1463 absres = r_norm;
1464 relres = absres/normr0;
1465 break;
1466 case STOP_REL_PRECRES:
1467 if ( pc == NULL )
1468 fasp_darray_cp(n, r, w);
1469 else
1470 pc->fct(r, w, pc->data);
1471 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
1472 relres = absres/normr0;
1473 break;
1474 case STOP_MOD_REL_RES:
1475 absres = r_norm;
1476 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1477 relres = absres/normu;
1478 break;
1479 }
1480
1481 norms[iter] = relres;
1482
1483 if ( relres <= tol ) {
1484 break;
1485 }
1486 else {
1487 // Need to restart
1488 fasp_darray_cp(n, r, p[0]); i = 0;
1489 }
1490
1491 } /* end of convergence check */
1492
1493 /* compute residual vector and continue loop */
1494 for ( j = i; j > 0; j-- ) {
1495 rs[j-1] = -s[j-1]*rs[j];
1496 rs[j] = c[j-1]*rs[j];
1497 }
1498
1499 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1500
1501 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1502
1503 if ( i ) {
1504 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1505 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1506 }
1507
1508 //-----------------------------------//
1509 // compute the convergence rate //
1510 //-----------------------------------//
1511 cr = r_norm / r_norm_old;
1512
1513 } /* end of iteration while loop */
1514
1515RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1516 if ( iter != iter_best ) {
1517
1518 // compute best residual
1519 fasp_darray_cp(n,b->val,r);
1520 fasp_blas_dstr_aAxpy(-1.0,A,x_best,r);
1521
1522 switch ( StopType ) {
1523 case STOP_REL_RES:
1524 absres_best = fasp_blas_darray_norm2(n,r);
1525 break;
1526 case STOP_REL_PRECRES:
1527 // z = B(r)
1528 if ( pc != NULL )
1529 pc->fct(r,w,pc->data); /* Apply preconditioner */
1530 else
1531 fasp_darray_cp(n,r,w); /* No preconditioner */
1532 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
1533 break;
1534 case STOP_MOD_REL_RES:
1535 absres_best = fasp_blas_darray_norm2(n,r);
1536 break;
1537 }
1538
1539 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1540 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1541 fasp_darray_cp(n,x_best,x->val);
1542 relres = absres_best / normr0;
1543 }
1544 }
1545
1546FINISHED:
1547 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1548
1549 /*-------------------------------------------
1550 * Free some stuff
1551 *------------------------------------------*/
1552 fasp_mem_free(work); work = NULL;
1553 fasp_mem_free(p); p = NULL;
1554 fasp_mem_free(hh); hh = NULL;
1555 fasp_mem_free(norms); norms = NULL;
1556
1557#if DEBUG_MODE > 0
1558 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1559#endif
1560
1561 if (iter>=MaxIt)
1562 return ERROR_SOLVER_MAXIT;
1563 else
1564 return iter;
1565}
1566
1567/*---------------------------------*/
1568/*-- End of File --*/
1569/*---------------------------------*/
fasp_darray_cp
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:210
fasp_mem_free
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:152
fasp_mem_calloc
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
fasp_itinfo
void fasp_itinfo(const INT ptrlvl, const INT stop_type, const INT iter, const REAL relres, const REAL absres, const REAL factor)
Print out iteration information for iterative solvers.
Definition: AuxMessage.c:41
fasp_chkerr
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
fasp_dvec_isnan
SHORT fasp_dvec_isnan(const dvector *u)
Check a dvector whether there is NAN.
Definition: AuxVector.c:39
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_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_dcsr_spvgmres
INT fasp_solver_dcsr_spvgmres(const dCSRmat *A, const 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: KrySPvgmres.c:68
fasp_solver_dblc_spvgmres
INT fasp_solver_dblc_spvgmres(const dBLCmat *A, const dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KrySPvgmres.c:829
fasp_solver_dbsr_spvgmres
INT fasp_solver_dbsr_spvgmres(const dBSRmat *A, const 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: KrySPvgmres.c:449
fasp_solver_dstr_spvgmres
INT fasp_solver_dstr_spvgmres(const dSTRmat *A, const 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: KrySPvgmres.c:1210
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
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
STAG_RATIO
#define STAG_RATIO
Definition: fasp_const.h:265
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
ERROR_ALLOC_MEM
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
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
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
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