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