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