Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
spcg.c
Go to the documentation of this file.
1 
57 #include <math.h>
58 
59 #include "fasp.h"
60 #include "fasp_functs.h"
61 #include "itsolver_util.inl"
62 
63 /*---------------------------------*/
64 /*-- Public Functions --*/
65 /*---------------------------------*/
66 
89  dvector *b,
90  dvector *u,
91  precond *pc,
92  const REAL tol,
93  const INT MaxIt,
94  const SHORT stop_type,
95  const SHORT prtlvl)
96 {
97  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
98  const INT m = b->row;
99  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
100  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
101 
102  // local variables
103  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
104  REAL absres0 = BIGREAL, absres = BIGREAL;
105  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
106  REAL reldiff, factor, infnormu;
107  REAL alpha, beta, temp1, temp2;
108 
109  INT iter_best = 0; // initial best known iteration
110  REAL absres_best = BIGREAL; // initial best known residual
111 
112  // allocate temp memory (need 5*m REAL numbers)
113  REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
114  REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
115 
116 #if DEBUG_MODE > 0
117  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
118  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
119 #endif
120 
121  // r = b-A*u
122  fasp_array_cp(m,b->val,r);
123  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
124 
125  if (pc != NULL)
126  pc->fct(r,z,pc->data); /* Apply preconditioner */
127  else
128  fasp_array_cp(m,r,z); /* No preconditioner */
129 
130  // compute initial residuals
131  switch (stop_type) {
132  case STOP_REL_RES:
133  absres0 = fasp_blas_array_norm2(m,r);
134  normr0 = MAX(SMALLREAL,absres0);
135  relres = absres0/normr0;
136  break;
137  case STOP_REL_PRECRES:
138  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
139  normr0 = MAX(SMALLREAL,absres0);
140  relres = absres0/normr0;
141  break;
142  case STOP_MOD_REL_RES:
143  absres0 = fasp_blas_array_norm2(m,r);
144  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
145  relres = absres0/normu;
146  break;
147  default:
148  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
149  goto FINISHED;
150  }
151 
152  // if initial residual is small, no need to iterate!
153  if ( relres < tol ) goto FINISHED;
154 
155  // output iteration information if needed
156  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
157 
158  fasp_array_cp(m,z,p);
159  temp1 = fasp_blas_array_dotprod(m,z,r);
160 
161  // main PCG loop
162  while ( iter++ < MaxIt ) {
163 
164  // t=A*p
165  fasp_blas_dcsr_mxv(A,p,t);
166 
167  // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
168  temp2 = fasp_blas_array_dotprod(m,t,p);
169  if ( ABS(temp2) > SMALLREAL2 ) {
170  alpha = temp1/temp2;
171  }
172  else { // Possible breakdown
173  goto RESTORE_BESTSOL;
174  }
175 
176  // u_k=u_{k-1} + alpha_k*p_{k-1}
177  fasp_blas_array_axpy(m,alpha,p,u->val);
178 
179  // r_k=r_{k-1} - alpha_k*A*p_{k-1}
180  fasp_blas_array_axpy(m,-alpha,t,r);
181 
182  // compute residuals
183  switch ( stop_type ) {
184  case STOP_REL_RES:
185  absres = fasp_blas_array_norm2(m,r);
186  relres = absres/normr0;
187  break;
188  case STOP_REL_PRECRES:
189  // z = B(r)
190  if ( pc != NULL )
191  pc->fct(r,z,pc->data); /* Apply preconditioner */
192  else
193  fasp_array_cp(m,r,z); /* No preconditioner */
194  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
195  relres = absres/normr0;
196  break;
197  case STOP_MOD_REL_RES:
198  absres = fasp_blas_array_norm2(m,r);
199  relres = absres/normu;
200  break;
201  }
202 
203  // compute reducation factor of residual ||r||
204  factor = absres/absres0;
205 
206  // output iteration information if needed
207  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
208 
209  // If the solution is NAN, restrore the best solution
210  if ( fasp_dvec_isnan(u) ) {
211  absres = BIGREAL;
212  goto RESTORE_BESTSOL;
213  }
214 
215  // safety net check: save the best-so-far solution
216  if ( absres < absres_best - maxdiff) {
217  absres_best = absres;
218  iter_best = iter;
219  fasp_array_cp(m,u->val,u_best);
220  }
221 
222  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
223  infnormu = fasp_blas_array_norminf(m, u->val);
224  if ( infnormu <= sol_inf_tol ) {
225  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
226  iter = ERROR_SOLVER_SOLSTAG;
227  break;
228  }
229 
230  // Check II: if staggenated, try to restart
231  normu = fasp_blas_dvec_norm2(u);
232 
233  // compute relative difference
234  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/normu;
235  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
236 
237  if ( prtlvl >= PRINT_MORE ) {
238  ITS_DIFFRES(reldiff,relres);
239  ITS_RESTART;
240  }
241 
242  fasp_array_cp(m,b->val,r);
243  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
244 
245  // compute residuals
246  switch ( stop_type ) {
247  case STOP_REL_RES:
248  absres = fasp_blas_array_norm2(m,r);
249  relres = absres/normr0;
250  break;
251  case STOP_REL_PRECRES:
252  // z = B(r)
253  if ( pc != NULL )
254  pc->fct(r,z,pc->data); /* Apply preconditioner */
255  else
256  fasp_array_cp(m,r,z); /* No preconditioner */
257  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
258  relres = absres/normr0;
259  break;
260  case STOP_MOD_REL_RES:
261  absres = fasp_blas_array_norm2(m,r);
262  relres = absres/normu;
263  break;
264  }
265 
266  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
267 
268  if ( relres < tol )
269  break;
270  else {
271  if ( stag >= MaxStag ) {
272  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
273  iter = ERROR_SOLVER_STAG;
274  break;
275  }
276  fasp_array_set(m,p,0.0);
277  ++stag;
278  ++restart_step;
279  }
280  } // end of staggnation check!
281 
282  // Check III: prevent false convergence
283  if ( relres < tol ) {
284 
285  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
286 
287  // compute residual r = b - Ax again
288  fasp_array_cp(m,b->val,r);
289  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
290 
291  // compute residuals
292  switch ( stop_type ) {
293  case STOP_REL_RES:
294  absres = fasp_blas_array_norm2(m,r);
295  relres = absres/normr0;
296  break;
297  case STOP_REL_PRECRES:
298  // z = B(r)
299  if ( pc != NULL )
300  pc->fct(r,z,pc->data); /* Apply preconditioner */
301  else
302  fasp_array_cp(m,r,z); /* No preconditioner */
303  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
304  relres = absres/normr0;
305  break;
306  case STOP_MOD_REL_RES:
307  absres = fasp_blas_array_norm2(m,r);
308  relres = absres/normu;
309  break;
310  }
311 
312  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
313 
314  // check convergence
315  if ( relres < tol ) break;
316 
317  if ( more_step >= MaxRestartStep ) {
318  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
319  iter = ERROR_SOLVER_TOLSMALL;
320  break;
321  }
322 
323  // prepare for restarting the method
324  fasp_array_set(m,p,0.0);
325  ++more_step;
326  ++restart_step;
327 
328  } // end of safe-guard check!
329 
330  // save residual for next iteration
331  absres0 = absres;
332 
333  // compute z_k = B(r_k)
334  if ( stop_type != STOP_REL_PRECRES ) {
335  if ( pc != NULL )
336  pc->fct(r,z,pc->data); /* Apply preconditioner */
337  else
338  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
339  }
340 
341  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
342  temp2 = fasp_blas_array_dotprod(m,z,r);
343  beta = temp2/temp1;
344  temp1 = temp2;
345 
346  // compute p_k = z_k + beta_k*p_{k-1}
347  fasp_blas_array_axpby(m,1.0,z,beta,p);
348 
349  } // end of main PCG loop.
350 
351 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
352  if ( iter != iter_best ) {
353 
354  // compute best residual
355  fasp_array_cp(m,b->val,r);
356  fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
357 
358  switch ( stop_type ) {
359  case STOP_REL_RES:
360  absres_best = fasp_blas_array_norm2(m,r);
361  break;
362  case STOP_REL_PRECRES:
363  // z = B(r)
364  if ( pc != NULL )
365  pc->fct(r,z,pc->data); /* Apply preconditioner */
366  else
367  fasp_array_cp(m,r,z); /* No preconditioner */
368  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
369  break;
370  case STOP_MOD_REL_RES:
371  absres_best = fasp_blas_array_norm2(m,r);
372  break;
373  }
374 
375  if ( absres > absres_best + maxdiff ) {
376  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
377  fasp_array_cp(m,u_best,u->val);
378  relres = absres_best / normr0;
379  }
380  }
381 
382 FINISHED: // finish the iterative method
383  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
384 
385  // clean up temp memory
386  fasp_mem_free(work);
387 
388 #if DEBUG_MODE > 0
389  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
390 #endif
391 
392  if ( iter > MaxIt )
393  return ERROR_SOLVER_MAXIT;
394  else
395  return iter;
396 }
397 
398 
421  dvector *b,
422  dvector *u,
423  precond *pc,
424  const REAL tol,
425  const INT MaxIt,
426  const SHORT stop_type,
427  const SHORT prtlvl)
428 {
429  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
430  const INT m = b->row;
431  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
432  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
433 
434  // local variables
435  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
436  REAL absres0 = BIGREAL, absres = BIGREAL;
437  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
438  REAL reldiff, factor, infnormu;
439  REAL alpha, beta, temp1, temp2;
440 
441  INT iter_best = 0; // initial best known iteration
442  REAL absres_best = BIGREAL; // initial best known residual
443 
444  // allocate temp memory (need 5*m REAL numbers)
445  REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
446  REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
447 
448 #if DEBUG_MODE > 0
449  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
450  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
451 #endif
452 
453  // r = b-A*u
454  fasp_array_cp(m,b->val,r);
455  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
456 
457  if (pc != NULL)
458  pc->fct(r,z,pc->data); /* Apply preconditioner */
459  else
460  fasp_array_cp(m,r,z); /* No preconditioner */
461 
462  // compute initial residuals
463  switch (stop_type) {
464  case STOP_REL_RES:
465  absres0 = fasp_blas_array_norm2(m,r);
466  normr0 = MAX(SMALLREAL,absres0);
467  relres = absres0/normr0;
468  break;
469  case STOP_REL_PRECRES:
470  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
471  normr0 = MAX(SMALLREAL,absres0);
472  relres = absres0/normr0;
473  break;
474  case STOP_MOD_REL_RES:
475  absres0 = fasp_blas_array_norm2(m,r);
476  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
477  relres = absres0/normu;
478  break;
479  default:
480  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
481  goto FINISHED;
482  }
483 
484  // if initial residual is small, no need to iterate!
485  if ( relres < tol ) goto FINISHED;
486 
487  // output iteration information if needed
488  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
489 
490  fasp_array_cp(m,z,p);
491  temp1 = fasp_blas_array_dotprod(m,z,r);
492 
493  // main PCG loop
494  while ( iter++ < MaxIt ) {
495 
496  // t=A*p
497  fasp_blas_bdcsr_mxv(A,p,t);
498 
499  // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
500  temp2 = fasp_blas_array_dotprod(m,t,p);
501  if ( ABS(temp2) > SMALLREAL2 ) {
502  alpha = temp1/temp2;
503  }
504  else { // Possible breakdown
505  goto RESTORE_BESTSOL;
506  }
507 
508  // u_k=u_{k-1} + alpha_k*p_{k-1}
509  fasp_blas_array_axpy(m,alpha,p,u->val);
510 
511  // r_k=r_{k-1} - alpha_k*A*p_{k-1}
512  fasp_blas_array_axpy(m,-alpha,t,r);
513 
514  // compute residuals
515  switch ( stop_type ) {
516  case STOP_REL_RES:
517  absres = fasp_blas_array_norm2(m,r);
518  relres = absres/normr0;
519  break;
520  case STOP_REL_PRECRES:
521  // z = B(r)
522  if ( pc != NULL )
523  pc->fct(r,z,pc->data); /* Apply preconditioner */
524  else
525  fasp_array_cp(m,r,z); /* No preconditioner */
526  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
527  relres = absres/normr0;
528  break;
529  case STOP_MOD_REL_RES:
530  absres = fasp_blas_array_norm2(m,r);
531  relres = absres/normu;
532  break;
533  }
534 
535  // compute reducation factor of residual ||r||
536  factor = absres/absres0;
537 
538  // output iteration information if needed
539  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
540 
541  // If the solution is NAN, restrore the best solution
542  if ( fasp_dvec_isnan(u) ) {
543  absres = BIGREAL;
544  goto RESTORE_BESTSOL;
545  }
546 
547  // safety net check: save the best-so-far solution
548  if ( absres < absres_best - maxdiff) {
549  absres_best = absres;
550  iter_best = iter;
551  fasp_array_cp(m,u->val,u_best);
552  }
553 
554  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
555  infnormu = fasp_blas_array_norminf(m, u->val);
556  if ( infnormu <= sol_inf_tol ) {
557  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
558  iter = ERROR_SOLVER_SOLSTAG;
559  break;
560  }
561 
562  // Check II: if staggenated, try to restart
563  normu = fasp_blas_dvec_norm2(u);
564 
565  // compute relative difference
566  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/normu;
567  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
568 
569  if ( prtlvl >= PRINT_MORE ) {
570  ITS_DIFFRES(reldiff,relres);
571  ITS_RESTART;
572  }
573 
574  fasp_array_cp(m,b->val,r);
575  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
576 
577  // compute residuals
578  switch ( stop_type ) {
579  case STOP_REL_RES:
580  absres = fasp_blas_array_norm2(m,r);
581  relres = absres/normr0;
582  break;
583  case STOP_REL_PRECRES:
584  // z = B(r)
585  if ( pc != NULL )
586  pc->fct(r,z,pc->data); /* Apply preconditioner */
587  else
588  fasp_array_cp(m,r,z); /* No preconditioner */
589  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
590  relres = absres/normr0;
591  break;
592  case STOP_MOD_REL_RES:
593  absres = fasp_blas_array_norm2(m,r);
594  relres = absres/normu;
595  break;
596  }
597 
598  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
599 
600  if ( relres < tol )
601  break;
602  else {
603  if ( stag >= MaxStag ) {
604  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
605  iter = ERROR_SOLVER_STAG;
606  break;
607  }
608  fasp_array_set(m,p,0.0);
609  ++stag;
610  ++restart_step;
611  }
612  } // end of staggnation check!
613 
614  // Check III: prevent false convergence
615  if ( relres < tol ) {
616 
617  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
618 
619  // compute residual r = b - Ax again
620  fasp_array_cp(m,b->val,r);
621  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
622 
623  // compute residuals
624  switch ( stop_type ) {
625  case STOP_REL_RES:
626  absres = fasp_blas_array_norm2(m,r);
627  relres = absres/normr0;
628  break;
629  case STOP_REL_PRECRES:
630  // z = B(r)
631  if ( pc != NULL )
632  pc->fct(r,z,pc->data); /* Apply preconditioner */
633  else
634  fasp_array_cp(m,r,z); /* No preconditioner */
635  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
636  relres = absres/normr0;
637  break;
638  case STOP_MOD_REL_RES:
639  absres = fasp_blas_array_norm2(m,r);
640  relres = absres/normu;
641  break;
642  }
643 
644  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
645 
646  // check convergence
647  if ( relres < tol ) break;
648 
649  if ( more_step >= MaxRestartStep ) {
650  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
651  iter = ERROR_SOLVER_TOLSMALL;
652  break;
653  }
654 
655  // prepare for restarting the method
656  fasp_array_set(m,p,0.0);
657  ++more_step;
658  ++restart_step;
659 
660  } // end of safe-guard check!
661 
662  // save residual for next iteration
663  absres0 = absres;
664 
665  // compute z_k = B(r_k)
666  if ( stop_type != STOP_REL_PRECRES ) {
667  if ( pc != NULL )
668  pc->fct(r,z,pc->data); /* Apply preconditioner */
669  else
670  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
671  }
672 
673  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
674  temp2 = fasp_blas_array_dotprod(m,z,r);
675  beta = temp2/temp1;
676  temp1 = temp2;
677 
678  // compute p_k = z_k + beta_k*p_{k-1}
679  fasp_blas_array_axpby(m,1.0,z,beta,p);
680 
681  } // end of main PCG loop.
682 
683 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
684  if ( iter != iter_best ) {
685 
686  // compute best residual
687  fasp_array_cp(m,b->val,r);
688  fasp_blas_bdcsr_aAxpy(-1.0,A,u_best,r);
689 
690  switch ( stop_type ) {
691  case STOP_REL_RES:
692  absres_best = fasp_blas_array_norm2(m,r);
693  break;
694  case STOP_REL_PRECRES:
695  // z = B(r)
696  if ( pc != NULL )
697  pc->fct(r,z,pc->data); /* Apply preconditioner */
698  else
699  fasp_array_cp(m,r,z); /* No preconditioner */
700  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
701  break;
702  case STOP_MOD_REL_RES:
703  absres_best = fasp_blas_array_norm2(m,r);
704  break;
705  }
706 
707  if ( absres > absres_best + maxdiff ) {
708  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
709  fasp_array_cp(m,u_best,u->val);
710  relres = absres_best / normr0;
711  }
712  }
713 
714 FINISHED: // finish the iterative method
715  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
716 
717  // clean up temp memory
718  fasp_mem_free(work);
719 
720 #if DEBUG_MODE > 0
721  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
722 #endif
723 
724  if ( iter > MaxIt )
725  return ERROR_SOLVER_MAXIT;
726  else
727  return iter;
728 }
729 
752  dvector *b,
753  dvector *u,
754  precond *pc,
755  const REAL tol,
756  const INT MaxIt,
757  const SHORT stop_type,
758  const SHORT prtlvl)
759 {
760  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
761  const INT m = b->row;
762  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
763  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
764 
765  // local variables
766  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
767  REAL absres0 = BIGREAL, absres = BIGREAL;
768  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
769  REAL reldiff, factor, infnormu;
770  REAL alpha, beta, temp1, temp2;
771 
772  INT iter_best = 0; // initial best known iteration
773  REAL absres_best = BIGREAL; // initial best known residual
774 
775  // allocate temp memory (need 5*m REAL numbers)
776  REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
777  REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
778 
779 #if DEBUG_MODE > 0
780  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
781  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
782 #endif
783 
784  // r = b-A*u
785  fasp_array_cp(m,b->val,r);
786  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
787 
788  if (pc != NULL)
789  pc->fct(r,z,pc->data); /* Apply preconditioner */
790  else
791  fasp_array_cp(m,r,z); /* No preconditioner */
792 
793  // compute initial residuals
794  switch (stop_type) {
795  case STOP_REL_RES:
796  absres0 = fasp_blas_array_norm2(m,r);
797  normr0 = MAX(SMALLREAL,absres0);
798  relres = absres0/normr0;
799  break;
800  case STOP_REL_PRECRES:
801  absres0 = sqrt(fasp_blas_array_dotprod(m,r,z));
802  normr0 = MAX(SMALLREAL,absres0);
803  relres = absres0/normr0;
804  break;
805  case STOP_MOD_REL_RES:
806  absres0 = fasp_blas_array_norm2(m,r);
807  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
808  relres = absres0/normu;
809  break;
810  default:
811  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
812  goto FINISHED;
813  }
814 
815  // if initial residual is small, no need to iterate!
816  if ( relres < tol ) goto FINISHED;
817 
818  // output iteration information if needed
819  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
820 
821  fasp_array_cp(m,z,p);
822  temp1 = fasp_blas_array_dotprod(m,z,r);
823 
824  // main PCG loop
825  while ( iter++ < MaxIt ) {
826 
827  // t=A*p
828  fasp_blas_dstr_mxv(A,p,t);
829 
830  // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
831  temp2 = fasp_blas_array_dotprod(m,t,p);
832  if ( ABS(temp2) > SMALLREAL2 ) {
833  alpha = temp1/temp2;
834  }
835  else { // Possible breakdown
836  goto RESTORE_BESTSOL;
837  }
838 
839  // u_k=u_{k-1} + alpha_k*p_{k-1}
840  fasp_blas_array_axpy(m,alpha,p,u->val);
841 
842  // r_k=r_{k-1} - alpha_k*A*p_{k-1}
843  fasp_blas_array_axpy(m,-alpha,t,r);
844 
845  // compute residuals
846  switch ( stop_type ) {
847  case STOP_REL_RES:
848  absres = fasp_blas_array_norm2(m,r);
849  relres = absres/normr0;
850  break;
851  case STOP_REL_PRECRES:
852  // z = B(r)
853  if ( pc != NULL )
854  pc->fct(r,z,pc->data); /* Apply preconditioner */
855  else
856  fasp_array_cp(m,r,z); /* No preconditioner */
857  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
858  relres = absres/normr0;
859  break;
860  case STOP_MOD_REL_RES:
861  absres = fasp_blas_array_norm2(m,r);
862  relres = absres/normu;
863  break;
864  }
865 
866  // compute reducation factor of residual ||r||
867  factor = absres/absres0;
868 
869  // output iteration information if needed
870  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
871 
872  // If the solution is NAN, restrore the best solution
873  if ( fasp_dvec_isnan(u) ) {
874  absres = BIGREAL;
875  goto RESTORE_BESTSOL;
876  }
877 
878  // safety net check: save the best-so-far solution
879  if ( absres < absres_best - maxdiff) {
880  absres_best = absres;
881  iter_best = iter;
882  fasp_array_cp(m,u->val,u_best);
883  }
884 
885  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
886  infnormu = fasp_blas_array_norminf(m, u->val);
887  if ( infnormu <= sol_inf_tol ) {
888  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
889  iter = ERROR_SOLVER_SOLSTAG;
890  break;
891  }
892 
893  // Check II: if staggenated, try to restart
894  normu = fasp_blas_dvec_norm2(u);
895 
896  // compute relative difference
897  reldiff = ABS(alpha)*fasp_blas_array_norm2(m,p)/normu;
898  if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
899 
900  if ( prtlvl >= PRINT_MORE ) {
901  ITS_DIFFRES(reldiff,relres);
902  ITS_RESTART;
903  }
904 
905  fasp_array_cp(m,b->val,r);
906  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
907 
908  // compute residuals
909  switch ( stop_type ) {
910  case STOP_REL_RES:
911  absres = fasp_blas_array_norm2(m,r);
912  relres = absres/normr0;
913  break;
914  case STOP_REL_PRECRES:
915  // z = B(r)
916  if ( pc != NULL )
917  pc->fct(r,z,pc->data); /* Apply preconditioner */
918  else
919  fasp_array_cp(m,r,z); /* No preconditioner */
920  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
921  relres = absres/normr0;
922  break;
923  case STOP_MOD_REL_RES:
924  absres = fasp_blas_array_norm2(m,r);
925  relres = absres/normu;
926  break;
927  }
928 
929  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
930 
931  if ( relres < tol )
932  break;
933  else {
934  if ( stag >= MaxStag ) {
935  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
936  iter = ERROR_SOLVER_STAG;
937  break;
938  }
939  fasp_array_set(m,p,0.0);
940  ++stag;
941  ++restart_step;
942  }
943  } // end of staggnation check!
944 
945  // Check III: prevent false convergence
946  if ( relres < tol ) {
947 
948  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
949 
950  // compute residual r = b - Ax again
951  fasp_array_cp(m,b->val,r);
952  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
953 
954  // compute residuals
955  switch ( stop_type ) {
956  case STOP_REL_RES:
957  absres = fasp_blas_array_norm2(m,r);
958  relres = absres/normr0;
959  break;
960  case STOP_REL_PRECRES:
961  // z = B(r)
962  if ( pc != NULL )
963  pc->fct(r,z,pc->data); /* Apply preconditioner */
964  else
965  fasp_array_cp(m,r,z); /* No preconditioner */
966  absres = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
967  relres = absres/normr0;
968  break;
969  case STOP_MOD_REL_RES:
970  absres = fasp_blas_array_norm2(m,r);
971  relres = absres/normu;
972  break;
973  }
974 
975  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
976 
977  // check convergence
978  if ( relres < tol ) break;
979 
980  if ( more_step >= MaxRestartStep ) {
981  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
982  iter = ERROR_SOLVER_TOLSMALL;
983  break;
984  }
985 
986  // prepare for restarting the method
987  fasp_array_set(m,p,0.0);
988  ++more_step;
989  ++restart_step;
990 
991  } // end of safe-guard check!
992 
993  // save residual for next iteration
994  absres0 = absres;
995 
996  // compute z_k = B(r_k)
997  if ( stop_type != STOP_REL_PRECRES ) {
998  if ( pc != NULL )
999  pc->fct(r,z,pc->data); /* Apply preconditioner */
1000  else
1001  fasp_array_cp(m,r,z); /* No preconditioner, B=I */
1002  }
1003 
1004  // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
1005  temp2 = fasp_blas_array_dotprod(m,z,r);
1006  beta = temp2/temp1;
1007  temp1 = temp2;
1008 
1009  // compute p_k = z_k + beta_k*p_{k-1}
1010  fasp_blas_array_axpby(m,1.0,z,beta,p);
1011 
1012  } // end of main PCG loop.
1013 
1014 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1015  if ( iter != iter_best ) {
1016 
1017  // compute best residual
1018  fasp_array_cp(m,b->val,r);
1019  fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
1020 
1021  switch ( stop_type ) {
1022  case STOP_REL_RES:
1023  absres_best = fasp_blas_array_norm2(m,r);
1024  break;
1025  case STOP_REL_PRECRES:
1026  // z = B(r)
1027  if ( pc != NULL )
1028  pc->fct(r,z,pc->data); /* Apply preconditioner */
1029  else
1030  fasp_array_cp(m,r,z); /* No preconditioner */
1031  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1032  break;
1033  case STOP_MOD_REL_RES:
1034  absres_best = fasp_blas_array_norm2(m,r);
1035  break;
1036  }
1037 
1038  if ( absres > absres_best + maxdiff ) {
1039  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1040  fasp_array_cp(m,u_best,u->val);
1041  relres = absres_best / normr0;
1042  }
1043  }
1044 
1045 FINISHED: // finish the iterative method
1046  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1047 
1048  // clean up temp memory
1049  fasp_mem_free(work);
1050 
1051 #if DEBUG_MODE > 0
1052  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1053 #endif
1054 
1055  if ( iter > MaxIt )
1056  return ERROR_SOLVER_MAXIT;
1057  else
1058  return iter;
1059 }
1060 
1061 /*---------------------------------*/
1062 /*-- End of File --*/
1063 /*---------------------------------*/
INT fasp_solver_bdcsr_spcg(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 with safety net.
Definition: spcg.c:420
#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_dstr_spcg(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 with safety net.
Definition: spcg.c:751
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
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1005
Vector with n entries of REAL type.
Definition: fasp.h:342
INT fasp_dvec_isnan(dvector *u)
Check a dvector whether there is NAN.
Definition: vec.c:33
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:54
#define INT
Definition: fasp.h:64
void print_itinfo(const INT ptrlvl, const INT stop_type, const INT iter, const REAL relres, const REAL absres, const REAL factor)
Print out iteration information for iterative solvers.
Definition: message.c:36
Preconditioner data and action.
Definition: fasp.h:999
void fasp_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
#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
#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
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 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
INT fasp_solver_dcsr_spcg(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 with safety net.
Definition: spcg.c:88
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
#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