Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pcg.c
Go to the documentation of this file.
1 
50 #include <math.h>
51 
52 #include "fasp.h"
53 #include "fasp_functs.h"
54 #include "itsolver_util.inl"
55 
56 /*---------------------------------*/
57 /*-- Public Functions --*/
58 /*---------------------------------*/
59 
85  dvector *b,
86  dvector *u,
87  precond *pc,
88  const REAL tol,
89  const INT MaxIt,
90  const SHORT stop_type,
91  const SHORT prtlvl)
92 {
93  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
94  const INT m = b->row;
95  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
96  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
97 
98  // local variables
99  INT iter = 0, stag = 1, more_step = 1;
100  REAL absres0 = BIGREAL, absres = BIGREAL;
101  REAL relres = BIGREAL, unorm2 = BIGREAL, normr0 = BIGREAL;
102  REAL reldiff, factor, unorminf;
103  REAL alpha, beta, temp1, temp2;
104 
105  // allocate temp memory (need 4*m REAL numbers)
106  REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
107  REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
108 
109 #if DEBUG_MODE > 0
110  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
111  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
112 #endif
113 
114  // r = b-A*u
115  fasp_array_cp(m,b->val,r);
116  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
117 
118  if ( pc != NULL )
119  pc->fct(r,z,pc->data); /* Apply preconditioner */
120  else
121  fasp_array_cp(m,r,z); /* No preconditioner */
122 
123  // compute initial residuals
124  switch ( stop_type ) {
125  case STOP_REL_RES:
126  absres0 = fasp_blas_array_norm2(m,r);
127  normr0 = MAX(SMALLREAL,absres0);
128  relres = absres0/normr0;
129  break;
130  case STOP_REL_PRECRES:
131  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
132  normr0 = MAX(SMALLREAL,absres0);
133  relres = absres0/normr0;
134  break;
135  case STOP_MOD_REL_RES:
136  absres0 = fasp_blas_array_norm2(m,r);
137  unorm2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
138  relres = absres0/unorm2;
139  break;
140  default:
141  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
142  goto FINISHED;
143  }
144 
145  // if initial residual is small, no need to iterate!
146  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
147 
148  // output iteration information if needed
149  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
150 
151  fasp_array_cp(m,z,p);
152  temp1 = fasp_blas_array_dotprod(m,z,r);
153 
154  // main PCG loop
155  while ( iter++ < MaxIt ) {
156 
157  // t = A*p
158  fasp_blas_dcsr_mxv(A,p,t);
159 
160  // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
161  temp2 = fasp_blas_array_dotprod(m,t,p);
162  if ( ABS(temp2) > SMALLREAL2 ) {
163  alpha = temp1/temp2;
164  }
165  else { // Possible breakdown
166  ITS_DIVZERO; goto FINISHED;
167  }
168 
169  // u_k = u_{k-1} + alpha_k*p_{k-1}
170  fasp_blas_array_axpy(m,alpha,p,u->val);
171 
172  // r_k = r_{k-1} - alpha_k*A*p_{k-1}
173  fasp_blas_array_axpy(m,-alpha,t,r);
174 
175  // compute norm of residual
176  switch ( stop_type ) {
177  case STOP_REL_RES:
178  absres = fasp_blas_array_norm2(m,r);
179  relres = absres/normr0;
180  break;
181  case STOP_REL_PRECRES:
182  // z = B(r)
183  if ( pc != NULL )
184  pc->fct(r,z,pc->data); /* Apply preconditioner */
185  else
186  fasp_array_cp(m,r,z); /* No preconditioner */
187  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
188  relres = absres/normr0;
189  break;
190  case STOP_MOD_REL_RES:
191  absres = fasp_blas_array_norm2(m,r);
192  relres = absres/unorm2;
193  break;
194  }
195 
196  // compute reduction factor of residual ||r||
197  factor = absres/absres0;
198 
199  // output iteration information if needed
200  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
201 
202  if ( factor > 0.9 ) {
203 
204  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
205  unorminf = fasp_blas_array_norminf(m, u->val);
206  if ( unorminf <= sol_inf_tol ) {
207  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
208  iter = ERROR_SOLVER_SOLSTAG;
209  break;
210  }
211 
212  // Check II: if stagnated, try to restart
213  unorm2 = fasp_blas_dvec_norm2(u);
214 
215  // compute relative difference
216  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/unorm2;
217  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
218 
219  if ( prtlvl >= PRINT_MORE ) {
220  ITS_DIFFRES(reldiff,relres);
221  ITS_RESTART;
222  }
223 
224  fasp_array_cp(m,b->val,r);
225  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
226 
227  // compute residual norms
228  switch ( stop_type ) {
229  case STOP_REL_RES:
230  absres = fasp_blas_array_norm2(m,r);
231  relres = absres/normr0;
232  break;
233  case STOP_REL_PRECRES:
234  // z = B(r)
235  if ( pc != NULL )
236  pc->fct(r,z,pc->data); /* Apply preconditioner */
237  else
238  fasp_array_cp(m,r,z); /* No preconditioner */
239  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
240  relres = absres/normr0;
241  break;
242  case STOP_MOD_REL_RES:
243  absres = fasp_blas_array_norm2(m,r);
244  relres = absres/unorm2;
245  break;
246  }
247 
248  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
249 
250  if ( relres < tol )
251  break;
252  else {
253  if ( stag >= MaxStag ) {
254  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
255  iter = ERROR_SOLVER_STAG;
256  break;
257  }
258  fasp_array_set(m,p,0.0);
259  ++stag;
260  }
261 
262  } // end of stagnation check!
263 
264  } // end of check I and II
265 
266  // Check III: prevent false convergence
267  if ( relres < tol ) {
268 
269  REAL updated_relres = relres;
270 
271  // compute true residual r = b - Ax and update residual
272  fasp_array_cp(m,b->val,r);
273  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
274 
275  // compute residual norms
276  switch ( stop_type ) {
277  case STOP_REL_RES:
278  absres = fasp_blas_array_norm2(m,r);
279  relres = absres/normr0;
280  break;
281  case STOP_REL_PRECRES:
282  // z = B(r)
283  if ( pc != NULL )
284  pc->fct(r,z,pc->data); /* Apply preconditioner */
285  else
286  fasp_array_cp(m,r,z); /* No preconditioner */
287  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
288  relres = absres/normr0;
289  break;
290  case STOP_MOD_REL_RES:
291  absres = fasp_blas_array_norm2(m,r);
292  relres = absres/unorm2;
293  break;
294  }
295 
296  // check convergence
297  if ( relres < tol ) break;
298 
299  if ( prtlvl >= PRINT_MORE ) {
300  ITS_COMPRES(updated_relres); ITS_REALRES(relres);
301  }
302 
303  if ( more_step >= MaxRestartStep ) {
304  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
305  iter = ERROR_SOLVER_TOLSMALL;
306  break;
307  }
308 
309  // prepare for restarting the method
310  fasp_array_set(m,p,0.0);
311  ++more_step;
312 
313  } // end of safe-guard check!
314 
315  // save residual for next iteration
316  absres0 = absres;
317 
318  // compute z_k = B(r_k)
319  if ( stop_type != STOP_REL_PRECRES ) {
320  if ( pc != NULL )
321  pc->fct(r,z,pc->data); /* Apply preconditioner */
322  else
323  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
324  }
325 
326  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
327  temp2 = fasp_blas_array_dotprod(m,z,r);
328  beta = temp2/temp1;
329  temp1 = temp2;
330 
331  // compute p_k = z_k + beta_k*p_{k-1}
332  fasp_blas_array_axpby(m,1.0,z,beta,p);
333 
334  } // end of main PCG loop.
335 
336 FINISHED: // finish the iterative method
337  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
338 
339  // clean up temp memory
340  fasp_mem_free(work);
341 
342 #if DEBUG_MODE > 0
343  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
344 #endif
345 
346  if ( iter > MaxIt )
347  return ERROR_SOLVER_MAXIT;
348  else
349  return iter;
350 }
351 
374  dvector *b,
375  dvector *u,
376  precond *pc,
377  const REAL tol,
378  const INT MaxIt,
379  const SHORT stop_type,
380  const SHORT prtlvl)
381 {
382  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
383  const INT m = b->row;
384  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
385  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
386 
387  // local variables
388  INT iter = 0, stag = 1, more_step = 1;
389  REAL absres0 = BIGREAL, absres = BIGREAL;
390  REAL relres = BIGREAL, unorm2 = BIGREAL, normr0 = BIGREAL;
391  REAL reldiff, factor, unorminf;
392  REAL alpha, beta, temp1, temp2;
393 
394  // allocate temp memory (need 4*m REAL numbers)
395  REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
396  REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
397 
398 #if DEBUG_MODE > 0
399  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
400  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
401 #endif
402 
403  // r = b-A*u
404  fasp_array_cp(m,b->val,r);
405  fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
406 
407  if ( pc != NULL )
408  pc->fct(r,z,pc->data); /* Apply preconditioner */
409  else
410  fasp_array_cp(m,r,z); /* No preconditioner */
411 
412  // compute initial residuals
413  switch ( stop_type ) {
414  case STOP_REL_RES:
415  absres0 = fasp_blas_array_norm2(m,r);
416  normr0 = MAX(SMALLREAL,absres0);
417  relres = absres0/normr0;
418  break;
419  case STOP_REL_PRECRES:
420  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
421  normr0 = MAX(SMALLREAL,absres0);
422  relres = absres0/normr0;
423  break;
424  case STOP_MOD_REL_RES:
425  absres0 = fasp_blas_array_norm2(m,r);
426  unorm2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
427  relres = absres0/unorm2;
428  break;
429  default:
430  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
431  goto FINISHED;
432  }
433 
434  // if initial residual is small, no need to iterate!
435  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
436 
437  // output iteration information if needed
438  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
439 
440  fasp_array_cp(m,z,p);
441  temp1 = fasp_blas_array_dotprod(m,z,r);
442 
443  // main PCG loop
444  while ( iter++ < MaxIt ) {
445 
446  // t = A*p
447  fasp_blas_dbsr_mxv(A,p,t);
448 
449  // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
450  temp2 = fasp_blas_array_dotprod(m,t,p);
451  if ( ABS(temp2) > SMALLREAL2 ) {
452  alpha = temp1/temp2;
453  }
454  else { // Possible breakdown
455  ITS_DIVZERO; goto FINISHED;
456  }
457 
458  // u_k = u_{k-1} + alpha_k*p_{k-1}
459  fasp_blas_array_axpy(m,alpha,p,u->val);
460 
461  // r_k = r_{k-1} - alpha_k*A*p_{k-1}
462  fasp_blas_array_axpy(m,-alpha,t,r);
463 
464  // compute norm of residual
465  switch ( stop_type ) {
466  case STOP_REL_RES:
467  absres = fasp_blas_array_norm2(m,r);
468  relres = absres/normr0;
469  break;
470  case STOP_REL_PRECRES:
471  // z = B(r)
472  if ( pc != NULL )
473  pc->fct(r,z,pc->data); /* Apply preconditioner */
474  else
475  fasp_array_cp(m,r,z); /* No preconditioner */
476  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
477  relres = absres/normr0;
478  break;
479  case STOP_MOD_REL_RES:
480  absres = fasp_blas_array_norm2(m,r);
481  relres = absres/unorm2;
482  break;
483  }
484 
485  // compute reduction factor of residual ||r||
486  factor = absres/absres0;
487 
488  // output iteration information if needed
489  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
490 
491  if ( factor > 0.9 ) {
492 
493  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
494  unorminf = fasp_blas_array_norminf(m, u->val);
495  if ( unorminf <= sol_inf_tol ) {
496  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
497  iter = ERROR_SOLVER_SOLSTAG;
498  break;
499  }
500 
501  // Check II: if stagnated, try to restart
502  unorm2 = fasp_blas_dvec_norm2(u);
503 
504  // compute relative difference
505  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/unorm2;
506  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
507 
508  if ( prtlvl >= PRINT_MORE ) {
509  ITS_DIFFRES(reldiff,relres);
510  ITS_RESTART;
511  }
512 
513  fasp_array_cp(m,b->val,r);
514  fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
515 
516  // compute residual norms
517  switch ( stop_type ) {
518  case STOP_REL_RES:
519  absres = fasp_blas_array_norm2(m,r);
520  relres = absres/normr0;
521  break;
522  case STOP_REL_PRECRES:
523  // z = B(r)
524  if ( pc != NULL )
525  pc->fct(r,z,pc->data); /* Apply preconditioner */
526  else
527  fasp_array_cp(m,r,z); /* No preconditioner */
528  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
529  relres = absres/normr0;
530  break;
531  case STOP_MOD_REL_RES:
532  absres = fasp_blas_array_norm2(m,r);
533  relres = absres/unorm2;
534  break;
535  }
536 
537  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
538 
539  if ( relres < tol )
540  break;
541  else {
542  if ( stag >= MaxStag ) {
543  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
544  iter = ERROR_SOLVER_STAG;
545  break;
546  }
547  fasp_array_set(m,p,0.0);
548  ++stag;
549  }
550 
551  } // end of stagnation check!
552 
553  } // end of check I and II
554 
555  // Check III: prevent false convergence
556  if ( relres < tol ) {
557 
558  REAL updated_relres = relres;
559 
560  // compute true residual r = b - Ax and update residual
561  fasp_array_cp(m,b->val,r);
562  fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
563 
564  // compute residual norms
565  switch ( stop_type ) {
566  case STOP_REL_RES:
567  absres = fasp_blas_array_norm2(m,r);
568  relres = absres/normr0;
569  break;
570  case STOP_REL_PRECRES:
571  // z = B(r)
572  if ( pc != NULL )
573  pc->fct(r,z,pc->data); /* Apply preconditioner */
574  else
575  fasp_array_cp(m,r,z); /* No preconditioner */
576  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
577  relres = absres/normr0;
578  break;
579  case STOP_MOD_REL_RES:
580  absres = fasp_blas_array_norm2(m,r);
581  relres = absres/unorm2;
582  break;
583  }
584 
585  // check convergence
586  if ( relres < tol ) break;
587 
588  if ( prtlvl >= PRINT_MORE ) {
589  ITS_COMPRES(updated_relres); ITS_REALRES(relres);
590  }
591 
592  if ( more_step >= MaxRestartStep ) {
593  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
594  iter = ERROR_SOLVER_TOLSMALL;
595  break;
596  }
597 
598  // prepare for restarting the method
599  fasp_array_set(m,p,0.0);
600  ++more_step;
601 
602  } // end of safe-guard check!
603 
604  // save residual for next iteration
605  absres0 = absres;
606 
607  // compute z_k = B(r_k)
608  if ( stop_type != STOP_REL_PRECRES ) {
609  if ( pc != NULL )
610  pc->fct(r,z,pc->data); /* Apply preconditioner */
611  else
612  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
613  }
614 
615  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
616  temp2 = fasp_blas_array_dotprod(m,z,r);
617  beta = temp2/temp1;
618  temp1 = temp2;
619 
620  // compute p_k = z_k + beta_k*p_{k-1}
621  fasp_blas_array_axpby(m,1.0,z,beta,p);
622 
623  } // end of main PCG loop.
624 
625 FINISHED: // finish the iterative method
626  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
627 
628  // clean up temp memory
629  fasp_mem_free(work);
630 
631 #if DEBUG_MODE > 0
632  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
633 #endif
634 
635  if ( iter > MaxIt )
636  return ERROR_SOLVER_MAXIT;
637  else
638  return iter;
639 }
640 
666  dvector *b,
667  dvector *u,
668  precond *pc,
669  const REAL tol,
670  const INT MaxIt,
671  const SHORT stop_type,
672  const SHORT prtlvl)
673 {
674  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
675  const INT m = b->row;
676  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
677  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
678 
679  // local variables
680  INT iter = 0, stag = 1, more_step = 1;
681  REAL absres0 = BIGREAL, absres = BIGREAL;
682  REAL relres = BIGREAL, unorm2 = BIGREAL, normr0 = BIGREAL;
683  REAL reldiff, factor, unorminf;
684  REAL alpha, beta, temp1, temp2;
685 
686  // allocate temp memory (need 4*m REAL numbers)
687  REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
688  REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
689 
690 #if DEBUG_MODE > 0
691  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
692  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
693 #endif
694 
695  // r = b-A*u
696  fasp_array_cp(m,b->val,r);
697  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
698 
699  if ( pc != NULL )
700  pc->fct(r,z,pc->data); /* Apply preconditioner */
701  else
702  fasp_array_cp(m,r,z); /* No preconditioner */
703 
704  // compute initial residuals
705  switch ( stop_type ) {
706  case STOP_REL_RES:
707  absres0 = fasp_blas_array_norm2(m,r);
708  normr0 = MAX(SMALLREAL,absres0);
709  relres = absres0/normr0;
710  break;
711  case STOP_REL_PRECRES:
712  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
713  normr0 = MAX(SMALLREAL,absres0);
714  relres = absres0/normr0;
715  break;
716  case STOP_MOD_REL_RES:
717  absres0 = fasp_blas_array_norm2(m,r);
718  unorm2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
719  relres = absres0/unorm2;
720  break;
721  default:
722  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
723  goto FINISHED;
724  }
725 
726  // if initial residual is small, no need to iterate!
727  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
728 
729  // output iteration information if needed
730  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
731 
732  fasp_array_cp(m,z,p);
733  temp1 = fasp_blas_array_dotprod(m,z,r);
734 
735  // main PCG loop
736  while ( iter++ < MaxIt ) {
737 
738  // t = A*p
739  fasp_blas_bdcsr_mxv(A,p,t);
740 
741  // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
742  temp2 = fasp_blas_array_dotprod(m,t,p);
743  if ( ABS(temp2) > SMALLREAL2 ) {
744  alpha = temp1/temp2;
745  }
746  else { // Possible breakdown
747  ITS_DIVZERO; goto FINISHED;
748  }
749 
750  // u_k = u_{k-1} + alpha_k*p_{k-1}
751  fasp_blas_array_axpy(m,alpha,p,u->val);
752 
753  // r_k = r_{k-1} - alpha_k*A*p_{k-1}
754  fasp_blas_array_axpy(m,-alpha,t,r);
755 
756  // compute norm of residual
757  switch ( stop_type ) {
758  case STOP_REL_RES:
759  absres = fasp_blas_array_norm2(m,r);
760  relres = absres/normr0;
761  break;
762  case STOP_REL_PRECRES:
763  // z = B(r)
764  if ( pc != NULL )
765  pc->fct(r,z,pc->data); /* Apply preconditioner */
766  else
767  fasp_array_cp(m,r,z); /* No preconditioner */
768  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
769  relres = absres/normr0;
770  break;
771  case STOP_MOD_REL_RES:
772  absres = fasp_blas_array_norm2(m,r);
773  relres = absres/unorm2;
774  break;
775  }
776 
777  // compute reduction factor of residual ||r||
778  factor = absres/absres0;
779 
780  // output iteration information if needed
781  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
782 
783  if ( factor > 0.9 ) {
784 
785  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
786  unorminf = fasp_blas_array_norminf(m, u->val);
787  if ( unorminf <= sol_inf_tol ) {
788  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
789  iter = ERROR_SOLVER_SOLSTAG;
790  break;
791  }
792 
793  // Check II: if stagnated, try to restart
794  unorm2 = fasp_blas_dvec_norm2(u);
795 
796  // compute relative difference
797  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/unorm2;
798  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
799 
800  if ( prtlvl >= PRINT_MORE ) {
801  ITS_DIFFRES(reldiff,relres);
802  ITS_RESTART;
803  }
804 
805  fasp_array_cp(m,b->val,r);
806  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
807 
808  // compute residual norms
809  switch ( stop_type ) {
810  case STOP_REL_RES:
811  absres = fasp_blas_array_norm2(m,r);
812  relres = absres/normr0;
813  break;
814  case STOP_REL_PRECRES:
815  // z = B(r)
816  if ( pc != NULL )
817  pc->fct(r,z,pc->data); /* Apply preconditioner */
818  else
819  fasp_array_cp(m,r,z); /* No preconditioner */
820  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
821  relres = absres/normr0;
822  break;
823  case STOP_MOD_REL_RES:
824  absres = fasp_blas_array_norm2(m,r);
825  relres = absres/unorm2;
826  break;
827  }
828 
829  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
830 
831  if ( relres < tol )
832  break;
833  else {
834  if ( stag >= MaxStag ) {
835  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
836  iter = ERROR_SOLVER_STAG;
837  break;
838  }
839  fasp_array_set(m,p,0.0);
840  ++stag;
841  }
842 
843  } // end of stagnation check!
844 
845  } // end of check I and II
846 
847  // Check III: prevent false convergence
848  if ( relres < tol ) {
849 
850  REAL updated_relres = relres;
851 
852  // compute true residual r = b - Ax and update residual
853  fasp_array_cp(m,b->val,r);
854  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
855 
856  // compute residual norms
857  switch ( stop_type ) {
858  case STOP_REL_RES:
859  absres = fasp_blas_array_norm2(m,r);
860  relres = absres/normr0;
861  break;
862  case STOP_REL_PRECRES:
863  // z = B(r)
864  if ( pc != NULL )
865  pc->fct(r,z,pc->data); /* Apply preconditioner */
866  else
867  fasp_array_cp(m,r,z); /* No preconditioner */
868  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
869  relres = absres/normr0;
870  break;
871  case STOP_MOD_REL_RES:
872  absres = fasp_blas_array_norm2(m,r);
873  relres = absres/unorm2;
874  break;
875  }
876 
877  // check convergence
878  if ( relres < tol ) break;
879 
880  if ( prtlvl >= PRINT_MORE ) {
881  ITS_COMPRES(updated_relres); ITS_REALRES(relres);
882  }
883 
884  if ( more_step >= MaxRestartStep ) {
885  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
886  iter = ERROR_SOLVER_TOLSMALL;
887  break;
888  }
889 
890  // prepare for restarting the method
891  fasp_array_set(m,p,0.0);
892  ++more_step;
893 
894  } // end of safe-guard check!
895 
896  // save residual for next iteration
897  absres0 = absres;
898 
899  // compute z_k = B(r_k)
900  if ( stop_type != STOP_REL_PRECRES ) {
901  if ( pc != NULL )
902  pc->fct(r,z,pc->data); /* Apply preconditioner */
903  else
904  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
905  }
906 
907  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
908  temp2 = fasp_blas_array_dotprod(m,z,r);
909  beta = temp2/temp1;
910  temp1 = temp2;
911 
912  // compute p_k = z_k + beta_k*p_{k-1}
913  fasp_blas_array_axpby(m,1.0,z,beta,p);
914 
915  } // end of main PCG loop.
916 
917 FINISHED: // finish the iterative method
918  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
919 
920  // clean up temp memory
921  fasp_mem_free(work);
922 
923 #if DEBUG_MODE > 0
924  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
925 #endif
926 
927  if ( iter > MaxIt )
928  return ERROR_SOLVER_MAXIT;
929  else
930  return iter;
931 }
932 
958  dvector *b,
959  dvector *u,
960  precond *pc,
961  const REAL tol,
962  const INT MaxIt,
963  const SHORT stop_type,
964  const SHORT prtlvl)
965 {
966  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
967  const INT m = b->row;
968  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
969  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
970 
971  // local variables
972  INT iter = 0, stag = 1, more_step = 1;
973  REAL absres0 = BIGREAL, absres = BIGREAL;
974  REAL relres = BIGREAL, unorm2 = BIGREAL, normr0 = BIGREAL;
975  REAL reldiff, factor, unorminf;
976  REAL alpha, beta, temp1, temp2;
977 
978  // allocate temp memory (need 4*m REAL numbers)
979  REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
980  REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
981 
982 #if DEBUG_MODE > 0
983  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
984  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
985 #endif
986 
987  // r = b-A*u
988  fasp_array_cp(m,b->val,r);
989  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
990 
991  if ( pc != NULL )
992  pc->fct(r,z,pc->data); /* Apply preconditioner */
993  else
994  fasp_array_cp(m,r,z); /* No preconditioner */
995 
996  // compute initial residuals
997  switch ( stop_type ) {
998  case STOP_REL_RES:
999  absres0 = fasp_blas_array_norm2(m,r);
1000  normr0 = MAX(SMALLREAL,absres0);
1001  relres = absres0/normr0;
1002  break;
1003  case STOP_REL_PRECRES:
1004  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
1005  normr0 = MAX(SMALLREAL,absres0);
1006  relres = absres0/normr0;
1007  break;
1008  case STOP_MOD_REL_RES:
1009  absres0 = fasp_blas_array_norm2(m,r);
1010  unorm2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
1011  relres = absres0/unorm2;
1012  break;
1013  default:
1014  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
1015  goto FINISHED;
1016  }
1017 
1018  // if initial residual is small, no need to iterate!
1019  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
1020 
1021  // output iteration information if needed
1022  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
1023 
1024  fasp_array_cp(m,z,p);
1025  temp1 = fasp_blas_array_dotprod(m,z,r);
1026 
1027  // main PCG loop
1028  while ( iter++ < MaxIt ) {
1029 
1030  // t = A*p
1031  fasp_blas_dstr_mxv(A,p,t);
1032 
1033  // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
1034  temp2 = fasp_blas_array_dotprod(m,t,p);
1035  if ( ABS(temp2) > SMALLREAL2 ) {
1036  alpha = temp1/temp2;
1037  }
1038  else { // Possible breakdown
1039  ITS_DIVZERO; goto FINISHED;
1040  }
1041 
1042  // u_k = u_{k-1} + alpha_k*p_{k-1}
1043  fasp_blas_array_axpy(m,alpha,p,u->val);
1044 
1045  // r_k = r_{k-1} - alpha_k*A*p_{k-1}
1046  fasp_blas_array_axpy(m,-alpha,t,r);
1047 
1048  // compute norm of residual
1049  switch ( stop_type ) {
1050  case STOP_REL_RES:
1051  absres = fasp_blas_array_norm2(m,r);
1052  relres = absres/normr0;
1053  break;
1054  case STOP_REL_PRECRES:
1055  // z = B(r)
1056  if ( pc != NULL )
1057  pc->fct(r,z,pc->data); /* Apply preconditioner */
1058  else
1059  fasp_array_cp(m,r,z); /* No preconditioner */
1060  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1061  relres = absres/normr0;
1062  break;
1063  case STOP_MOD_REL_RES:
1064  absres = fasp_blas_array_norm2(m,r);
1065  relres = absres/unorm2;
1066  break;
1067  }
1068 
1069  // compute reduction factor of residual ||r||
1070  factor = absres/absres0;
1071 
1072  // output iteration information if needed
1073  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1074 
1075  if ( factor > 0.9 ) {
1076 
1077  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
1078  unorminf = fasp_blas_array_norminf(m, u->val);
1079  if ( unorminf <= sol_inf_tol ) {
1080  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1081  iter = ERROR_SOLVER_SOLSTAG;
1082  break;
1083  }
1084 
1085  // Check II: if stagnated, try to restart
1086  unorm2 = fasp_blas_dvec_norm2(u);
1087 
1088  // compute relative difference
1089  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/unorm2;
1090  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
1091 
1092  if ( prtlvl >= PRINT_MORE ) {
1093  ITS_DIFFRES(reldiff,relres);
1094  ITS_RESTART;
1095  }
1096 
1097  fasp_array_cp(m,b->val,r);
1098  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1099 
1100  // compute residual norms
1101  switch ( stop_type ) {
1102  case STOP_REL_RES:
1103  absres = fasp_blas_array_norm2(m,r);
1104  relres = absres/normr0;
1105  break;
1106  case STOP_REL_PRECRES:
1107  // z = B(r)
1108  if ( pc != NULL )
1109  pc->fct(r,z,pc->data); /* Apply preconditioner */
1110  else
1111  fasp_array_cp(m,r,z); /* No preconditioner */
1112  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1113  relres = absres/normr0;
1114  break;
1115  case STOP_MOD_REL_RES:
1116  absres = fasp_blas_array_norm2(m,r);
1117  relres = absres/unorm2;
1118  break;
1119  }
1120 
1121  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1122 
1123  if ( relres < tol )
1124  break;
1125  else {
1126  if ( stag >= MaxStag ) {
1127  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1128  iter = ERROR_SOLVER_STAG;
1129  break;
1130  }
1131  fasp_array_set(m,p,0.0);
1132  ++stag;
1133  }
1134 
1135  } // end of stagnation check!
1136 
1137  } // end of check I and II
1138 
1139  // Check III: prevent false convergence
1140  if ( relres < tol ) {
1141 
1142  REAL updated_relres = relres;
1143 
1144  // compute true residual r = b - Ax and update residual
1145  fasp_array_cp(m,b->val,r);
1146  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1147 
1148  // compute residual norms
1149  switch ( stop_type ) {
1150  case STOP_REL_RES:
1151  absres = fasp_blas_array_norm2(m,r);
1152  relres = absres/normr0;
1153  break;
1154  case STOP_REL_PRECRES:
1155  // z = B(r)
1156  if ( pc != NULL )
1157  pc->fct(r,z,pc->data); /* Apply preconditioner */
1158  else
1159  fasp_array_cp(m,r,z); /* No preconditioner */
1160  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1161  relres = absres/normr0;
1162  break;
1163  case STOP_MOD_REL_RES:
1164  absres = fasp_blas_array_norm2(m,r);
1165  relres = absres/unorm2;
1166  break;
1167  }
1168 
1169  // check convergence
1170  if ( relres < tol ) break;
1171 
1172  if ( prtlvl >= PRINT_MORE ) {
1173  ITS_COMPRES(updated_relres); ITS_REALRES(relres);
1174  }
1175 
1176  if ( more_step >= MaxRestartStep ) {
1177  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1178  iter = ERROR_SOLVER_TOLSMALL;
1179  break;
1180  }
1181 
1182  // prepare for restarting the method
1183  fasp_array_set(m,p,0.0);
1184  ++more_step;
1185 
1186  } // end of safe-guard check!
1187 
1188  // save residual for next iteration
1189  absres0 = absres;
1190 
1191  // compute z_k = B(r_k)
1192  if ( stop_type != STOP_REL_PRECRES ) {
1193  if ( pc != NULL )
1194  pc->fct(r,z,pc->data); /* Apply preconditioner */
1195  else
1196  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
1197  }
1198 
1199  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
1200  temp2 = fasp_blas_array_dotprod(m,z,r);
1201  beta = temp2/temp1;
1202  temp1 = temp2;
1203 
1204  // compute p_k = z_k + beta_k*p_{k-1}
1205  fasp_blas_array_axpby(m,1.0,z,beta,p);
1206 
1207  } // end of main PCG loop.
1208 
1209 FINISHED: // finish the iterative method
1210  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1211 
1212  // clean up temp memory
1213  fasp_mem_free(work);
1214 
1215 #if DEBUG_MODE > 0
1216  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1217 #endif
1218 
1219  if ( iter > MaxIt )
1220  return ERROR_SOLVER_MAXIT;
1221  else
1222  return iter;
1223 }
1224 
1225 /*---------------------------------*/
1226 /*-- End of File --*/
1227 /*---------------------------------*/
#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_dcsr_pcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: pcg.c:84
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_blas_array_axpby(const INT n, const REAL a, REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: blas_array.c:218
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
REAL fasp_blas_dvec_norm2(dvector *x)
L2 norm of dvector x.
Definition: blas_vec.c:265
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 PRINT_MORE
Definition: fasp_const.h:82
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
#define ERROR_SOLVER_TOLSMALL
Definition: fasp_const.h:51
INT fasp_solver_bdcsr_pcg(block_dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: pcg.c:665
#define MAX_RESTART
Definition: fasp_const.h:245
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
INT fasp_solver_dstr_pcg(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: pcg.c:957
#define ERROR_SOLVER_SOLSTAG
Definition: fasp_const.h:50
#define ERROR_SOLVER_STAG
Definition: fasp_const.h:49
#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
Block REAL CSR matrix format.
Definition: fasp_block.h:84
#define ABS(a)
Definition: fasp.h:74
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
void fasp_blas_dbsr_mxv(dBSRmat *A, REAL *x, REAL *y)
Compute y := A*x.
Definition: blas_bsr.c:895
REAL fasp_blas_array_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies (x,y)
Definition: blas_array.c:267
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
INT fasp_solver_dbsr_pcg(dBSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: pcg.c:373
#define SMALLREAL2
Definition: fasp_const.h:239
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
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_STAG
Definition: fasp_const.h:246
#define STAG_RATIO
Definition: fasp_const.h:247
#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