Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
spgmres.c
Go to the documentation of this file.
1 
13 #include <math.h>
14 
15 #include "fasp.h"
16 #include "fasp_functs.h"
17 #include "itsolver_util.inl"
18 
19 /*---------------------------------*/
20 /*-- Public Functions --*/
21 /*---------------------------------*/
22 
47  dvector *b,
48  dvector *x,
49  precond *pc,
50  const REAL tol,
51  const INT MaxIt,
52  SHORT restart,
53  const SHORT stop_type,
54  const SHORT prtlvl)
55 {
56  const INT n = b->row;
57  const INT MIN_ITER = 0;
58  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
59  const REAL epsmac = SMALLREAL;
60 
61  // local variables
62  INT iter = 0;
63  INT restart1 = restart + 1;
64  INT i, j, k;
65 
66  REAL r_norm, r_normb, gamma, t;
67  REAL normr0 = BIGREAL, absres = BIGREAL;
68  REAL relres = BIGREAL, normu = BIGREAL;
69 
70  INT iter_best = 0; // initial best known iteration
71  REAL absres_best = BIGREAL; // initial best known residual
72 
73  // allocate temp memory (need about (restart+4)*n REAL numbers)
74  REAL *c = NULL, *s = NULL, *rs = NULL;
75  REAL *norms = NULL, *r = NULL, *w = NULL;
76  REAL *work = NULL, *x_best = NULL;
77  REAL **p = NULL, **hh = NULL;
78 
79 #if DEBUG_MODE > 0
80  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
81  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
82 #endif
83 
84  /* allocate memory and setup temp work space */
85  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
86 
87  /* check whether memory is enough for GMRES */
88  while ( (work == NULL) && (restart > 5 ) ) {
89  restart = restart - 5 ;
90  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
91  printf("### WARNING: GMRES restart number set to %d!\n", restart );
92  restart1 = restart + 1;
93  }
94 
95  if ( work == NULL ) {
96  printf("### ERROR: No enough memory for GMRES %s : %s: %d !\n",
97  __FILE__, __FUNCTION__, __LINE__ );
98  exit(ERROR_ALLOC_MEM);
99  }
100 
101  p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
102  hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
103  norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
104 
105  r = work; w = r + n; rs = w + n; c = rs + restart1;
106  x_best = c + restart; s = x_best + n;
107 
108  for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
109 
110  for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
111 
112  // r = b-A*x
113  fasp_array_cp(n, b->val, p[0]);
114  fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
115 
116  r_norm = fasp_blas_array_norm2(n,p[0]);
117 
118  // compute initial residuals
119  switch (stop_type) {
120  case STOP_REL_RES:
121  normr0 = MAX(SMALLREAL,r_norm);
122  relres = r_norm/normr0;
123  break;
124  case STOP_REL_PRECRES:
125  if ( pc == NULL )
126  fasp_array_cp(n, p[0], r);
127  else
128  pc->fct(p[0], r, pc->data);
129  r_normb = sqrt(fasp_blas_array_dotprod(n,p[0],r));
130  normr0 = MAX(SMALLREAL,r_normb);
131  relres = r_normb/normr0;
132  break;
133  case STOP_MOD_REL_RES:
134  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
135  normr0 = r_norm;
136  relres = normr0/normu;
137  break;
138  default:
139  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
140  goto FINISHED;
141  }
142 
143  // if initial residual is small, no need to iterate!
144  if ( relres < tol ) goto FINISHED;
145 
146  // output iteration information if needed
147  print_itinfo(prtlvl,stop_type,0,relres,normr0,0.0);
148 
149  // store initial residual
150  norms[0] = relres;
151 
152  /* outer iteration cycle */
153  while ( iter < MaxIt ) {
154 
155  rs[0] = r_norm;
156 
157  t = 1.0 / r_norm;
158 
159  fasp_blas_array_ax(n, t, p[0]);
160 
161  /* RESTART CYCLE (right-preconditioning) */
162  i = 0;
163  while ( i < restart && iter < MaxIt ) {
164 
165  i++; iter++;
166 
167  /* apply the preconditioner */
168  if ( pc == NULL )
169  fasp_array_cp(n, p[i-1], r);
170  else
171  pc->fct(p[i-1], r, pc->data);
172 
173  fasp_blas_dcsr_mxv(A, r, p[i]);
174 
175  /* modified Gram_Schmidt */
176  for ( j = 0; j < i; j++ ) {
177  hh[j][i-1] = fasp_blas_array_dotprod(n, p[j], p[i]);
178  fasp_blas_array_axpy(n, -hh[j][i-1], p[j], p[i]);
179  }
180  t = fasp_blas_array_norm2(n, p[i]);
181  hh[i][i-1] = t;
182  if (t != 0.0) {
183  t = 1.0/t;
184  fasp_blas_array_ax(n, t, p[i]);
185  }
186 
187  for (j = 1; j < i; ++j) {
188  t = hh[j-1][i-1];
189  hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
190  hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
191  }
192  t= hh[i][i-1]*hh[i][i-1];
193  t+= hh[i-1][i-1]*hh[i-1][i-1];
194 
195  gamma = sqrt(t);
196  if (gamma == 0.0) gamma = epsmac;
197  c[i-1] = hh[i-1][i-1] / gamma;
198  s[i-1] = hh[i][i-1] / gamma;
199  rs[i] = -s[i-1]*rs[i-1];
200  rs[i-1] = c[i-1]*rs[i-1];
201  hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
202 
203  absres = r_norm = fabs(rs[i]);
204 
205  relres = absres/normr0;
206 
207  norms[iter] = relres;
208 
209  // output iteration information if needed
210  print_itinfo(prtlvl, stop_type, iter, relres, absres,
211  norms[iter]/norms[iter-1]);
212 
213  // should we exit the restart cycle
214  if ( relres <= tol && iter >= MIN_ITER ) break;
215 
216  } /* end of restart cycle */
217 
218  /* compute solution, first solve upper triangular system */
219  rs[i-1] = rs[i-1] / hh[i-1][i-1];
220  for ( k = i-2; k >= 0; k-- ) {
221  t = 0.0;
222  for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
223 
224  t += rs[k];
225  rs[k] = t / hh[k][k];
226  }
227 
228  fasp_array_cp(n, p[i-1], w);
229 
230  fasp_blas_array_ax(n, rs[i-1], w);
231 
232  for ( j = i-2; j >= 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], w);
233 
234  /* apply the preconditioner */
235  if ( pc == NULL )
236  fasp_array_cp(n, w, r);
237  else
238  pc->fct(w, r, pc->data);
239 
240  fasp_blas_array_axpy(n, 1.0, r, x->val);
241 
242  // safety net check: save the best-so-far solution
243  if ( fasp_dvec_isnan(x) ) {
244  // If the solution is NAN, restrore the best solution
245  absres = BIGREAL;
246  goto RESTORE_BESTSOL;
247  }
248 
249  if ( absres < absres_best - maxdiff) {
250  absres_best = absres;
251  iter_best = iter;
252  fasp_array_cp(n,x->val,x_best);
253  }
254 
255  // Check: prevent false convergence
256  if ( relres <= tol && iter >= MIN_ITER ) {
257 
258  fasp_array_cp(n, b->val, r);
259  fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
260 
261  r_norm = fasp_blas_array_norm2(n, r);
262 
263  switch ( stop_type ) {
264  case STOP_REL_RES:
265  absres = r_norm;
266  relres = absres/normr0;
267  break;
268  case STOP_REL_PRECRES:
269  if ( pc == NULL )
270  fasp_array_cp(n, r, w);
271  else
272  pc->fct(r, w, pc->data);
273  absres = sqrt(fasp_blas_array_dotprod(n,w,r));
274  relres = absres/normr0;
275  break;
276  case STOP_MOD_REL_RES:
277  absres = r_norm;
278  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
279  relres = absres/normu;
280  break;
281  }
282 
283  norms[iter] = relres;
284 
285  if ( relres <= tol ) {
286  break;
287  }
288  else {
289  // Need to restart
290  fasp_array_cp(n, r, p[0]); i = 0;
291  }
292 
293  } /* end of convergence check */
294 
295  /* compute residual vector and continue loop */
296  for (j = i; j > 0; j--) {
297  rs[j-1] = -s[j-1]*rs[j];
298  rs[j] = c[j-1]*rs[j];
299  }
300 
301  if ( i ) fasp_blas_array_axpy(n, rs[i]-1.0, p[i], p[i]);
302 
303  for ( j = i-1 ; j > 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], p[i]);
304 
305  if ( i ) {
306  fasp_blas_array_axpy(n, rs[0]-1.0, p[0], p[0]);
307  fasp_blas_array_axpy(n, 1.0, p[i], p[0]);
308  }
309 
310  } /* end of main while loop */
311 
312 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
313  if ( iter != iter_best ) {
314 
315  // compute best residual
316  fasp_array_cp(n,b->val,r);
317  fasp_blas_dcsr_aAxpy(-1.0,A,x_best,r);
318 
319  switch ( stop_type ) {
320  case STOP_REL_RES:
321  absres_best = fasp_blas_array_norm2(n,r);
322  break;
323  case STOP_REL_PRECRES:
324  // z = B(r)
325  if ( pc != NULL )
326  pc->fct(r,w,pc->data); /* Apply preconditioner */
327  else
328  fasp_array_cp(n,r,w); /* No preconditioner */
329  absres_best = sqrt(ABS(fasp_blas_array_dotprod(n,w,r)));
330  break;
331  case STOP_MOD_REL_RES:
332  absres_best = fasp_blas_array_norm2(n,r);
333  break;
334  }
335 
336  if ( absres > absres_best + maxdiff ) {
337  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
338  fasp_array_cp(n,x_best,x->val);
339  relres = absres_best / normr0;
340  }
341  }
342 
343 FINISHED:
344  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
345 
346  /*-------------------------------------------
347  * Clean up workspace
348  *------------------------------------------*/
349  fasp_mem_free(work);
350  fasp_mem_free(p);
351  fasp_mem_free(hh);
352  fasp_mem_free(norms);
353 
354 #if DEBUG_MODE > 0
355  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
356 #endif
357 
358  if ( iter >= MaxIt )
359  return ERROR_SOLVER_MAXIT;
360  else
361  return iter;
362 }
363 
387  dvector *b,
388  dvector *x,
389  precond *pc,
390  const REAL tol,
391  const INT MaxIt,
392  SHORT restart,
393  const SHORT stop_type,
394  const SHORT prtlvl)
395 {
396  const INT n = b->row;
397  const INT MIN_ITER = 0;
398  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
399  const REAL epsmac = SMALLREAL;
400 
401  // local variables
402  INT iter = 0;
403  INT restart1 = restart + 1;
404  INT i, j, k;
405 
406  REAL r_norm, r_normb, gamma, t;
407  REAL normr0 = BIGREAL, absres = BIGREAL;
408  REAL relres = BIGREAL, normu = BIGREAL;
409 
410  INT iter_best = 0; // initial best known iteration
411  REAL absres_best = BIGREAL; // initial best known residual
412 
413  // allocate temp memory (need about (restart+4)*n REAL numbers)
414  REAL *c = NULL, *s = NULL, *rs = NULL;
415  REAL *norms = NULL, *r = NULL, *w = NULL;
416  REAL *work = NULL, *x_best = NULL;
417  REAL **p = NULL, **hh = NULL;
418 
419 #if DEBUG_MODE > 0
420  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
421  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
422 #endif
423 
424  /* allocate memory and setup temp work space */
425  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
426 
427  /* check whether memory is enough for GMRES */
428  while ( (work == NULL) && (restart > 5 ) ) {
429  restart = restart - 5 ;
430  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
431  printf("### WARNING: GMRES restart number set to %d!\n", restart );
432  restart1 = restart + 1;
433  }
434 
435  if ( work == NULL ) {
436  printf("### ERROR: No enough memory for GMRES %s : %s: %d !\n",
437  __FILE__, __FUNCTION__, __LINE__ );
438  exit(ERROR_ALLOC_MEM);
439  }
440 
441  p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
442  hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
443  norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
444 
445  r = work; w = r + n; rs = w + n; c = rs + restart1;
446  x_best = c + restart; s = x_best + n;
447 
448  for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
449 
450  for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
451 
452  // r = b-A*x
453  fasp_array_cp(n, b->val, p[0]);
454  fasp_blas_bdcsr_aAxpy(-1.0, A, x->val, p[0]);
455 
456  r_norm = fasp_blas_array_norm2(n,p[0]);
457 
458  // compute initial residuals
459  switch (stop_type) {
460  case STOP_REL_RES:
461  normr0 = MAX(SMALLREAL,r_norm);
462  relres = r_norm/normr0;
463  break;
464  case STOP_REL_PRECRES:
465  if ( pc == NULL )
466  fasp_array_cp(n, p[0], r);
467  else
468  pc->fct(p[0], r, pc->data);
469  r_normb = sqrt(fasp_blas_array_dotprod(n,p[0],r));
470  normr0 = MAX(SMALLREAL,r_normb);
471  relres = r_normb/normr0;
472  break;
473  case STOP_MOD_REL_RES:
474  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
475  normr0 = r_norm;
476  relres = normr0/normu;
477  break;
478  default:
479  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
480  goto FINISHED;
481  }
482 
483  // if initial residual is small, no need to iterate!
484  if ( relres < tol ) goto FINISHED;
485 
486  // output iteration information if needed
487  print_itinfo(prtlvl,stop_type,0,relres,normr0,0.0);
488 
489  // store initial residual
490  norms[0] = relres;
491 
492  /* outer iteration cycle */
493  while ( iter < MaxIt ) {
494 
495  rs[0] = r_norm;
496 
497  t = 1.0 / r_norm;
498 
499  fasp_blas_array_ax(n, t, p[0]);
500 
501  /* RESTART CYCLE (right-preconditioning) */
502  i = 0;
503  while ( i < restart && iter < MaxIt ) {
504 
505  i++; iter++;
506 
507  /* apply the preconditioner */
508  if ( pc == NULL )
509  fasp_array_cp(n, p[i-1], r);
510  else
511  pc->fct(p[i-1], r, pc->data);
512 
513  fasp_blas_bdcsr_mxv(A, r, p[i]);
514 
515  /* modified Gram_Schmidt */
516  for ( j = 0; j < i; j++ ) {
517  hh[j][i-1] = fasp_blas_array_dotprod(n, p[j], p[i]);
518  fasp_blas_array_axpy(n, -hh[j][i-1], p[j], p[i]);
519  }
520  t = fasp_blas_array_norm2(n, p[i]);
521  hh[i][i-1] = t;
522  if (t != 0.0) {
523  t = 1.0/t;
524  fasp_blas_array_ax(n, t, p[i]);
525  }
526 
527  for (j = 1; j < i; ++j) {
528  t = hh[j-1][i-1];
529  hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
530  hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
531  }
532  t= hh[i][i-1]*hh[i][i-1];
533  t+= hh[i-1][i-1]*hh[i-1][i-1];
534 
535  gamma = sqrt(t);
536  if (gamma == 0.0) gamma = epsmac;
537  c[i-1] = hh[i-1][i-1] / gamma;
538  s[i-1] = hh[i][i-1] / gamma;
539  rs[i] = -s[i-1]*rs[i-1];
540  rs[i-1] = c[i-1]*rs[i-1];
541  hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
542 
543  absres = r_norm = fabs(rs[i]);
544 
545  relres = absres/normr0;
546 
547  norms[iter] = relres;
548 
549  // output iteration information if needed
550  print_itinfo(prtlvl, stop_type, iter, relres, absres,
551  norms[iter]/norms[iter-1]);
552 
553  // should we exit the restart cycle
554  if ( relres <= tol && iter >= MIN_ITER ) break;
555 
556  } /* end of restart cycle */
557 
558  /* compute solution, first solve upper triangular system */
559  rs[i-1] = rs[i-1] / hh[i-1][i-1];
560  for ( k = i-2; k >= 0; k-- ) {
561  t = 0.0;
562  for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
563 
564  t += rs[k];
565  rs[k] = t / hh[k][k];
566  }
567 
568  fasp_array_cp(n, p[i-1], w);
569 
570  fasp_blas_array_ax(n, rs[i-1], w);
571 
572  for ( j = i-2; j >= 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], w);
573 
574  /* apply the preconditioner */
575  if ( pc == NULL )
576  fasp_array_cp(n, w, r);
577  else
578  pc->fct(w, r, pc->data);
579 
580  fasp_blas_array_axpy(n, 1.0, r, x->val);
581 
582  // safety net check: save the best-so-far solution
583  if ( fasp_dvec_isnan(x) ) {
584  // If the solution is NAN, restrore the best solution
585  absres = BIGREAL;
586  goto RESTORE_BESTSOL;
587  }
588 
589  if ( absres < absres_best - maxdiff) {
590  absres_best = absres;
591  iter_best = iter;
592  fasp_array_cp(n,x->val,x_best);
593  }
594 
595  // Check: prevent false convergence
596  if ( relres <= tol && iter >= MIN_ITER ) {
597 
598  fasp_array_cp(n, b->val, r);
599  fasp_blas_bdcsr_aAxpy(-1.0, A, x->val, r);
600 
601  r_norm = fasp_blas_array_norm2(n, r);
602 
603  switch ( stop_type ) {
604  case STOP_REL_RES:
605  absres = r_norm;
606  relres = absres/normr0;
607  break;
608  case STOP_REL_PRECRES:
609  if ( pc == NULL )
610  fasp_array_cp(n, r, w);
611  else
612  pc->fct(r, w, pc->data);
613  absres = sqrt(fasp_blas_array_dotprod(n,w,r));
614  relres = absres/normr0;
615  break;
616  case STOP_MOD_REL_RES:
617  absres = r_norm;
618  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
619  relres = absres/normu;
620  break;
621  }
622 
623  norms[iter] = relres;
624 
625  if ( relres <= tol ) {
626  break;
627  }
628  else {
629  // Need to restart
630  fasp_array_cp(n, r, p[0]); i = 0;
631  }
632 
633  } /* end of convergence check */
634 
635  /* compute residual vector and continue loop */
636  for (j = i; j > 0; j--) {
637  rs[j-1] = -s[j-1]*rs[j];
638  rs[j] = c[j-1]*rs[j];
639  }
640 
641  if ( i ) fasp_blas_array_axpy(n, rs[i]-1.0, p[i], p[i]);
642 
643  for ( j = i-1 ; j > 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], p[i]);
644 
645  if ( i ) {
646  fasp_blas_array_axpy(n, rs[0]-1.0, p[0], p[0]);
647  fasp_blas_array_axpy(n, 1.0, p[i], p[0]);
648  }
649 
650  } /* end of main while loop */
651 
652 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
653  if ( iter != iter_best ) {
654 
655  // compute best residual
656  fasp_array_cp(n,b->val,r);
657  fasp_blas_bdcsr_aAxpy(-1.0,A,x_best,r);
658 
659  switch ( stop_type ) {
660  case STOP_REL_RES:
661  absres_best = fasp_blas_array_norm2(n,r);
662  break;
663  case STOP_REL_PRECRES:
664  // z = B(r)
665  if ( pc != NULL )
666  pc->fct(r,w,pc->data); /* Apply preconditioner */
667  else
668  fasp_array_cp(n,r,w); /* No preconditioner */
669  absres_best = sqrt(ABS(fasp_blas_array_dotprod(n,w,r)));
670  break;
671  case STOP_MOD_REL_RES:
672  absres_best = fasp_blas_array_norm2(n,r);
673  break;
674  }
675 
676  if ( absres > absres_best + maxdiff ) {
677  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
678  fasp_array_cp(n,x_best,x->val);
679  relres = absres_best / normr0;
680  }
681  }
682 
683 FINISHED:
684  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
685 
686  /*-------------------------------------------
687  * Clean up workspace
688  *------------------------------------------*/
689  fasp_mem_free(work);
690  fasp_mem_free(p);
691  fasp_mem_free(hh);
692  fasp_mem_free(norms);
693 
694 #if DEBUG_MODE > 0
695  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
696 #endif
697 
698  if ( iter >= MaxIt )
699  return ERROR_SOLVER_MAXIT;
700  else
701  return iter;
702 }
703 
727  dvector *b,
728  dvector *x,
729  precond *pc,
730  const REAL tol,
731  const INT MaxIt,
732  SHORT restart,
733  const SHORT stop_type,
734  const SHORT prtlvl)
735 {
736  const INT n = b->row;
737  const INT MIN_ITER = 0;
738  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
739  const REAL epsmac = SMALLREAL;
740 
741  // local variables
742  INT iter = 0;
743  INT restart1 = restart + 1;
744  INT i, j, k;
745 
746  REAL r_norm, r_normb, gamma, t;
747  REAL normr0 = BIGREAL, absres = BIGREAL;
748  REAL relres = BIGREAL, normu = BIGREAL;
749 
750  INT iter_best = 0; // initial best known iteration
751  REAL absres_best = BIGREAL; // initial best known residual
752 
753  // allocate temp memory (need about (restart+4)*n REAL numbers)
754  REAL *c = NULL, *s = NULL, *rs = NULL;
755  REAL *norms = NULL, *r = NULL, *w = NULL;
756  REAL *work = NULL, *x_best = NULL;
757  REAL **p = NULL, **hh = NULL;
758 
759 #if DEBUG_MODE > 0
760  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
761  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
762 #endif
763 
764  /* allocate memory and setup temp work space */
765  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
766 
767  /* check whether memory is enough for GMRES */
768  while ( (work == NULL) && (restart > 5 ) ) {
769  restart = restart - 5 ;
770  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
771  printf("### WARNING: GMRES restart number set to %d!\n", restart );
772  restart1 = restart + 1;
773  }
774 
775  if ( work == NULL ) {
776  printf("### ERROR: No enough memory for GMRES %s : %s: %d !\n",
777  __FILE__, __FUNCTION__, __LINE__ );
778  exit(ERROR_ALLOC_MEM);
779  }
780 
781  p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
782  hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
783  norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
784 
785  r = work; w = r + n; rs = w + n; c = rs + restart1;
786  x_best = c + restart; s = x_best + n;
787 
788  for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
789 
790  for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
791 
792  // r = b-A*x
793  fasp_array_cp(n, b->val, p[0]);
794  fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
795 
796  r_norm = fasp_blas_array_norm2(n,p[0]);
797 
798  // compute initial residuals
799  switch (stop_type) {
800  case STOP_REL_RES:
801  normr0 = MAX(SMALLREAL,r_norm);
802  relres = r_norm/normr0;
803  break;
804  case STOP_REL_PRECRES:
805  if ( pc == NULL )
806  fasp_array_cp(n, p[0], r);
807  else
808  pc->fct(p[0], r, pc->data);
809  r_normb = sqrt(fasp_blas_array_dotprod(n,p[0],r));
810  normr0 = MAX(SMALLREAL,r_normb);
811  relres = r_normb/normr0;
812  break;
813  case STOP_MOD_REL_RES:
814  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
815  normr0 = r_norm;
816  relres = normr0/normu;
817  break;
818  default:
819  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
820  goto FINISHED;
821  }
822 
823  // if initial residual is small, no need to iterate!
824  if ( relres < tol ) goto FINISHED;
825 
826  // output iteration information if needed
827  print_itinfo(prtlvl,stop_type,0,relres,normr0,0.0);
828 
829  // store initial residual
830  norms[0] = relres;
831 
832  /* outer iteration cycle */
833  while ( iter < MaxIt ) {
834 
835  rs[0] = r_norm;
836 
837  t = 1.0 / r_norm;
838 
839  fasp_blas_array_ax(n, t, p[0]);
840 
841  /* RESTART CYCLE (right-preconditioning) */
842  i = 0;
843  while ( i < restart && iter < MaxIt ) {
844 
845  i++; iter++;
846 
847  /* apply the preconditioner */
848  if ( pc == NULL )
849  fasp_array_cp(n, p[i-1], r);
850  else
851  pc->fct(p[i-1], r, pc->data);
852 
853  fasp_blas_dbsr_mxv(A, r, p[i]);
854 
855  /* modified Gram_Schmidt */
856  for ( j = 0; j < i; j++ ) {
857  hh[j][i-1] = fasp_blas_array_dotprod(n, p[j], p[i]);
858  fasp_blas_array_axpy(n, -hh[j][i-1], p[j], p[i]);
859  }
860  t = fasp_blas_array_norm2(n, p[i]);
861  hh[i][i-1] = t;
862  if (t != 0.0) {
863  t = 1.0/t;
864  fasp_blas_array_ax(n, t, p[i]);
865  }
866 
867  for (j = 1; j < i; ++j) {
868  t = hh[j-1][i-1];
869  hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
870  hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
871  }
872  t= hh[i][i-1]*hh[i][i-1];
873  t+= hh[i-1][i-1]*hh[i-1][i-1];
874 
875  gamma = sqrt(t);
876  if (gamma == 0.0) gamma = epsmac;
877  c[i-1] = hh[i-1][i-1] / gamma;
878  s[i-1] = hh[i][i-1] / gamma;
879  rs[i] = -s[i-1]*rs[i-1];
880  rs[i-1] = c[i-1]*rs[i-1];
881  hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
882 
883  absres = r_norm = fabs(rs[i]);
884 
885  relres = absres/normr0;
886 
887  norms[iter] = relres;
888 
889  // output iteration information if needed
890  print_itinfo(prtlvl, stop_type, iter, relres, absres,
891  norms[iter]/norms[iter-1]);
892 
893  // should we exit the restart cycle
894  if ( relres <= tol && iter >= MIN_ITER ) break;
895 
896  } /* end of restart cycle */
897 
898  /* compute solution, first solve upper triangular system */
899  rs[i-1] = rs[i-1] / hh[i-1][i-1];
900  for ( k = i-2; k >= 0; k-- ) {
901  t = 0.0;
902  for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
903 
904  t += rs[k];
905  rs[k] = t / hh[k][k];
906  }
907 
908  fasp_array_cp(n, p[i-1], w);
909 
910  fasp_blas_array_ax(n, rs[i-1], w);
911 
912  for ( j = i-2; j >= 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], w);
913 
914  /* apply the preconditioner */
915  if ( pc == NULL )
916  fasp_array_cp(n, w, r);
917  else
918  pc->fct(w, r, pc->data);
919 
920  fasp_blas_array_axpy(n, 1.0, r, x->val);
921 
922  // safety net check: save the best-so-far solution
923  if ( fasp_dvec_isnan(x) ) {
924  // If the solution is NAN, restrore the best solution
925  absres = BIGREAL;
926  goto RESTORE_BESTSOL;
927  }
928 
929  if ( absres < absres_best - maxdiff) {
930  absres_best = absres;
931  iter_best = iter;
932  fasp_array_cp(n,x->val,x_best);
933  }
934 
935  // Check: prevent false convergence
936  if ( relres <= tol && iter >= MIN_ITER ) {
937 
938  fasp_array_cp(n, b->val, r);
939  fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
940 
941  r_norm = fasp_blas_array_norm2(n, r);
942 
943  switch ( stop_type ) {
944  case STOP_REL_RES:
945  absres = r_norm;
946  relres = absres/normr0;
947  break;
948  case STOP_REL_PRECRES:
949  if ( pc == NULL )
950  fasp_array_cp(n, r, w);
951  else
952  pc->fct(r, w, pc->data);
953  absres = sqrt(fasp_blas_array_dotprod(n,w,r));
954  relres = absres/normr0;
955  break;
956  case STOP_MOD_REL_RES:
957  absres = r_norm;
958  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
959  relres = absres/normu;
960  break;
961  }
962 
963  norms[iter] = relres;
964 
965  if ( relres <= tol ) {
966  break;
967  }
968  else {
969  // Need to restart
970  fasp_array_cp(n, r, p[0]); i = 0;
971  }
972 
973  } /* end of convergence check */
974 
975  /* compute residual vector and continue loop */
976  for (j = i; j > 0; j--) {
977  rs[j-1] = -s[j-1]*rs[j];
978  rs[j] = c[j-1]*rs[j];
979  }
980 
981  if ( i ) fasp_blas_array_axpy(n, rs[i]-1.0, p[i], p[i]);
982 
983  for ( j = i-1 ; j > 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], p[i]);
984 
985  if ( i ) {
986  fasp_blas_array_axpy(n, rs[0]-1.0, p[0], p[0]);
987  fasp_blas_array_axpy(n, 1.0, p[i], p[0]);
988  }
989 
990  } /* end of main while loop */
991 
992 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
993  if ( iter != iter_best ) {
994 
995  // compute best residual
996  fasp_array_cp(n,b->val,r);
997  fasp_blas_dbsr_aAxpy(-1.0,A,x_best,r);
998 
999  switch ( stop_type ) {
1000  case STOP_REL_RES:
1001  absres_best = fasp_blas_array_norm2(n,r);
1002  break;
1003  case STOP_REL_PRECRES:
1004  // z = B(r)
1005  if ( pc != NULL )
1006  pc->fct(r,w,pc->data); /* Apply preconditioner */
1007  else
1008  fasp_array_cp(n,r,w); /* No preconditioner */
1009  absres_best = sqrt(ABS(fasp_blas_array_dotprod(n,w,r)));
1010  break;
1011  case STOP_MOD_REL_RES:
1012  absres_best = fasp_blas_array_norm2(n,r);
1013  break;
1014  }
1015 
1016  if ( absres > absres_best + maxdiff ) {
1017  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1018  fasp_array_cp(n,x_best,x->val);
1019  relres = absres_best / normr0;
1020  }
1021  }
1022 
1023 FINISHED:
1024  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1025 
1026  /*-------------------------------------------
1027  * Clean up workspace
1028  *------------------------------------------*/
1029  fasp_mem_free(work);
1030  fasp_mem_free(p);
1031  fasp_mem_free(hh);
1032  fasp_mem_free(norms);
1033 
1034 #if DEBUG_MODE > 0
1035  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1036 #endif
1037 
1038  if ( iter >= MaxIt )
1039  return ERROR_SOLVER_MAXIT;
1040  else
1041  return iter;
1042 }
1043 
1067  dvector *b,
1068  dvector *x,
1069  precond *pc,
1070  const REAL tol,
1071  const INT MaxIt,
1072  SHORT restart,
1073  const SHORT stop_type,
1074  const SHORT prtlvl)
1075 {
1076  const INT n = b->row;
1077  const INT MIN_ITER = 0;
1078  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1079  const REAL epsmac = SMALLREAL;
1080 
1081  // local variables
1082  INT iter = 0;
1083  INT restart1 = restart + 1;
1084  INT i, j, k;
1085 
1086  REAL r_norm, r_normb, gamma, t;
1087  REAL normr0 = BIGREAL, absres = BIGREAL;
1088  REAL relres = BIGREAL, normu = BIGREAL;
1089 
1090  INT iter_best = 0; // initial best known iteration
1091  REAL absres_best = BIGREAL; // initial best known residual
1092 
1093  // allocate temp memory (need about (restart+4)*n REAL numbers)
1094  REAL *c = NULL, *s = NULL, *rs = NULL;
1095  REAL *norms = NULL, *r = NULL, *w = NULL;
1096  REAL *work = NULL, *x_best = NULL;
1097  REAL **p = NULL, **hh = NULL;
1098 
1099 #if DEBUG_MODE > 0
1100  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
1101  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1102 #endif
1103 
1104  /* allocate memory and setup temp work space */
1105  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1106 
1107  /* check whether memory is enough for GMRES */
1108  while ( (work == NULL) && (restart > 5 ) ) {
1109  restart = restart - 5 ;
1110  work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1111  printf("### WARNING: GMRES restart number set to %d!\n", restart );
1112  restart1 = restart + 1;
1113  }
1114 
1115  if ( work == NULL ) {
1116  printf("### ERROR: No enough memory for GMRES %s : %s: %d !\n",
1117  __FILE__, __FUNCTION__, __LINE__ );
1118  exit(ERROR_ALLOC_MEM);
1119  }
1120 
1121  p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1122  hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1123  norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1124 
1125  r = work; w = r + n; rs = w + n; c = rs + restart1;
1126  x_best = c + restart; s = x_best + n;
1127 
1128  for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
1129 
1130  for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
1131 
1132  // r = b-A*x
1133  fasp_array_cp(n, b->val, p[0]);
1134  fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1135 
1136  r_norm = fasp_blas_array_norm2(n,p[0]);
1137 
1138  // compute initial residuals
1139  switch (stop_type) {
1140  case STOP_REL_RES:
1141  normr0 = MAX(SMALLREAL,r_norm);
1142  relres = r_norm/normr0;
1143  break;
1144  case STOP_REL_PRECRES:
1145  if ( pc == NULL )
1146  fasp_array_cp(n, p[0], r);
1147  else
1148  pc->fct(p[0], r, pc->data);
1149  r_normb = sqrt(fasp_blas_array_dotprod(n,p[0],r));
1150  normr0 = MAX(SMALLREAL,r_normb);
1151  relres = r_normb/normr0;
1152  break;
1153  case STOP_MOD_REL_RES:
1154  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
1155  normr0 = r_norm;
1156  relres = normr0/normu;
1157  break;
1158  default:
1159  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
1160  goto FINISHED;
1161  }
1162 
1163  // if initial residual is small, no need to iterate!
1164  if ( relres < tol ) goto FINISHED;
1165 
1166  // output iteration information if needed
1167  print_itinfo(prtlvl,stop_type,0,relres,normr0,0.0);
1168 
1169  // store initial residual
1170  norms[0] = relres;
1171 
1172  /* outer iteration cycle */
1173  while ( iter < MaxIt ) {
1174 
1175  rs[0] = r_norm;
1176 
1177  t = 1.0 / r_norm;
1178 
1179  fasp_blas_array_ax(n, t, p[0]);
1180 
1181  /* RESTART CYCLE (right-preconditioning) */
1182  i = 0;
1183  while ( i < restart && iter < MaxIt ) {
1184 
1185  i++; iter++;
1186 
1187  /* apply the preconditioner */
1188  if ( pc == NULL )
1189  fasp_array_cp(n, p[i-1], r);
1190  else
1191  pc->fct(p[i-1], r, pc->data);
1192 
1193  fasp_blas_dstr_mxv(A, r, p[i]);
1194 
1195  /* modified Gram_Schmidt */
1196  for ( j = 0; j < i; j++ ) {
1197  hh[j][i-1] = fasp_blas_array_dotprod(n, p[j], p[i]);
1198  fasp_blas_array_axpy(n, -hh[j][i-1], p[j], p[i]);
1199  }
1200  t = fasp_blas_array_norm2(n, p[i]);
1201  hh[i][i-1] = t;
1202  if (t != 0.0) {
1203  t = 1.0/t;
1204  fasp_blas_array_ax(n, t, p[i]);
1205  }
1206 
1207  for (j = 1; j < i; ++j) {
1208  t = hh[j-1][i-1];
1209  hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1210  hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1211  }
1212  t= hh[i][i-1]*hh[i][i-1];
1213  t+= hh[i-1][i-1]*hh[i-1][i-1];
1214 
1215  gamma = sqrt(t);
1216  if (gamma == 0.0) gamma = epsmac;
1217  c[i-1] = hh[i-1][i-1] / gamma;
1218  s[i-1] = hh[i][i-1] / gamma;
1219  rs[i] = -s[i-1]*rs[i-1];
1220  rs[i-1] = c[i-1]*rs[i-1];
1221  hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1222 
1223  absres = r_norm = fabs(rs[i]);
1224 
1225  relres = absres/normr0;
1226 
1227  norms[iter] = relres;
1228 
1229  // output iteration information if needed
1230  print_itinfo(prtlvl, stop_type, iter, relres, absres,
1231  norms[iter]/norms[iter-1]);
1232 
1233  // should we exit the restart cycle
1234  if ( relres <= tol && iter >= MIN_ITER ) break;
1235 
1236  } /* end of restart cycle */
1237 
1238  /* compute solution, first solve upper triangular system */
1239  rs[i-1] = rs[i-1] / hh[i-1][i-1];
1240  for ( k = i-2; k >= 0; k-- ) {
1241  t = 0.0;
1242  for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
1243 
1244  t += rs[k];
1245  rs[k] = t / hh[k][k];
1246  }
1247 
1248  fasp_array_cp(n, p[i-1], w);
1249 
1250  fasp_blas_array_ax(n, rs[i-1], w);
1251 
1252  for ( j = i-2; j >= 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], w);
1253 
1254  /* apply the preconditioner */
1255  if ( pc == NULL )
1256  fasp_array_cp(n, w, r);
1257  else
1258  pc->fct(w, r, pc->data);
1259 
1260  fasp_blas_array_axpy(n, 1.0, r, x->val);
1261 
1262  // safety net check: save the best-so-far solution
1263  if ( fasp_dvec_isnan(x) ) {
1264  // If the solution is NAN, restrore the best solution
1265  absres = BIGREAL;
1266  goto RESTORE_BESTSOL;
1267  }
1268 
1269  if ( absres < absres_best - maxdiff) {
1270  absres_best = absres;
1271  iter_best = iter;
1272  fasp_array_cp(n,x->val,x_best);
1273  }
1274 
1275  // Check: prevent false convergence
1276  if ( relres <= tol && iter >= MIN_ITER ) {
1277 
1278  fasp_array_cp(n, b->val, r);
1279  fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1280 
1281  r_norm = fasp_blas_array_norm2(n, r);
1282 
1283  switch ( stop_type ) {
1284  case STOP_REL_RES:
1285  absres = r_norm;
1286  relres = absres/normr0;
1287  break;
1288  case STOP_REL_PRECRES:
1289  if ( pc == NULL )
1290  fasp_array_cp(n, r, w);
1291  else
1292  pc->fct(r, w, pc->data);
1293  absres = sqrt(fasp_blas_array_dotprod(n,w,r));
1294  relres = absres/normr0;
1295  break;
1296  case STOP_MOD_REL_RES:
1297  absres = r_norm;
1298  normu = MAX(SMALLREAL,fasp_blas_array_norm2(n,x->val));
1299  relres = absres/normu;
1300  break;
1301  }
1302 
1303  norms[iter] = relres;
1304 
1305  if ( relres <= tol ) {
1306  break;
1307  }
1308  else {
1309  // Need to restart
1310  fasp_array_cp(n, r, p[0]); i = 0;
1311  }
1312 
1313  } /* end of convergence check */
1314 
1315  /* compute residual vector and continue loop */
1316  for (j = i; j > 0; j--) {
1317  rs[j-1] = -s[j-1]*rs[j];
1318  rs[j] = c[j-1]*rs[j];
1319  }
1320 
1321  if ( i ) fasp_blas_array_axpy(n, rs[i]-1.0, p[i], p[i]);
1322 
1323  for ( j = i-1 ; j > 0; j-- ) fasp_blas_array_axpy(n, rs[j], p[j], p[i]);
1324 
1325  if ( i ) {
1326  fasp_blas_array_axpy(n, rs[0]-1.0, p[0], p[0]);
1327  fasp_blas_array_axpy(n, 1.0, p[i], p[0]);
1328  }
1329 
1330  } /* end of main while loop */
1331 
1332 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1333  if ( iter != iter_best ) {
1334 
1335  // compute best residual
1336  fasp_array_cp(n,b->val,r);
1337  fasp_blas_dstr_aAxpy(-1.0,A,x_best,r);
1338 
1339  switch ( stop_type ) {
1340  case STOP_REL_RES:
1341  absres_best = fasp_blas_array_norm2(n,r);
1342  break;
1343  case STOP_REL_PRECRES:
1344  // z = B(r)
1345  if ( pc != NULL )
1346  pc->fct(r,w,pc->data); /* Apply preconditioner */
1347  else
1348  fasp_array_cp(n,r,w); /* No preconditioner */
1349  absres_best = sqrt(ABS(fasp_blas_array_dotprod(n,w,r)));
1350  break;
1351  case STOP_MOD_REL_RES:
1352  absres_best = fasp_blas_array_norm2(n,r);
1353  break;
1354  }
1355 
1356  if ( absres > absres_best + maxdiff ) {
1357  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1358  fasp_array_cp(n,x_best,x->val);
1359  relres = absres_best / normr0;
1360  }
1361  }
1362 
1363 FINISHED:
1364  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1365 
1366  /*-------------------------------------------
1367  * Clean up workspace
1368  *------------------------------------------*/
1369  fasp_mem_free(work);
1370  fasp_mem_free(p);
1371  fasp_mem_free(hh);
1372  fasp_mem_free(norms);
1373 
1374 #if DEBUG_MODE > 0
1375  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1376 #endif
1377 
1378  if ( iter >= MaxIt )
1379  return ERROR_SOLVER_MAXIT;
1380  else
1381  return iter;
1382 }
1383 
1384 /*---------------------------------*/
1385 /*-- End of File --*/
1386 /*---------------------------------*/
INT fasp_solver_dstr_spgmres(dSTRmat *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 with safe-guard.
Definition: spgmres.c:1066
#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
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
INT fasp_solver_bdcsr_spgmres(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 with safe-guard.
Definition: spgmres.c:386
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
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
#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
INT fasp_solver_dbsr_spgmres(dBSRmat *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 with safe-guard.
Definition: spgmres.c:726
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_dcsr_spgmres(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 with safe-guard.
Definition: spgmres.c:46
#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