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