Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
spbcgs.c
Go to the documentation of this file.
1 
59 #include <math.h>
60 
61 #include "fasp.h"
62 #include "fasp_functs.h"
63 #include "itsolver_util.inl"
64 
65 /*---------------------------------*/
66 /*-- Public Functions --*/
67 /*---------------------------------*/
68 
91  dvector *b,
92  dvector *u,
93  precond *pc,
94  const REAL tol,
95  const INT MaxIt,
96  const SHORT stop_type,
97  const SHORT prtlvl)
98 {
99  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
100  const INT m = b->row;
101  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
102  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
103  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
104 
105  // local variables
106  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
107  REAL alpha, beta, omega, temp1, temp2;
108  REAL absres0 = BIGREAL, absres = BIGREAL;
109  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
110  REAL reldiff, factor, normd, tempr, infnormu;
111 
112  REAL *uval = u->val, *bval = b->val;
113  INT iter_best = 0; // initial best known iteration
114  REAL absres_best = BIGREAL; // initial best known residual
115 
116  // allocate temp memory (need 8*m REAL)
117  REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
118  REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
119  REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
120 
121 #if DEBUG_MODE > 0
122  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
123  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
124 #endif
125 
126  // r = b-A*u
127  fasp_array_cp(m,bval,r);
128  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
129  absres0 = fasp_blas_array_norm2(m,r);
130 
131  // compute initial relative residual
132  switch (stop_type) {
133  case STOP_REL_RES:
134  normr0 = MAX(SMALLREAL,absres0);
135  relres = absres0/normr0;
136  break;
137  case STOP_REL_PRECRES:
138  normr0 = MAX(SMALLREAL,absres0);
139  relres = absres0/normr0;
140  break;
141  case STOP_MOD_REL_RES:
142  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
143  relres = absres0/normu;
144  break;
145  default:
146  printf("### ERROR: Unrecognized 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) goto FINISHED;
152 
153  // output iteration information if needed
154  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
155 
156  // rho = r* := r
157  fasp_array_cp(m,r,rho);
158  temp1 = fasp_blas_array_dotprod(m,r,rho);
159 
160  // p = r
161  fasp_array_cp(m,r,p);
162 
163  // main BiCGstab loop
164  while ( iter++ < MaxIt ) {
165 
166  // pp = precond(p)
167  if ( pc != NULL )
168  pc->fct(p,pp,pc->data); /* Apply preconditioner */
169  else
170  fasp_array_cp(m,p,pp); /* No preconditioner */
171 
172  // z = A*pp
173  fasp_blas_dcsr_mxv(A,pp,z);
174 
175  // alpha = (r,rho)/(A*p,rho)
176  temp2 = fasp_blas_array_dotprod(m,z,rho);
177  if ( ABS(temp2) > SMALLREAL ) {
178  alpha = temp1/temp2;
179  }
180  else {
181  ITS_DIVZERO; goto FINISHED;
182  }
183 
184  // s = r - alpha z
185  fasp_array_cp(m,r,s);
186  fasp_blas_array_axpy(m,-alpha,z,s);
187 
188  // sp = precond(s)
189  if ( pc != NULL )
190  pc->fct(s,sp,pc->data); /* Apply preconditioner */
191  else
192  fasp_array_cp(m,s,sp); /* No preconditioner */
193 
194  // t = A*sp;
195  fasp_blas_dcsr_mxv(A,sp,t);
196 
197  // omega = (t,s)/(t,t)
198  tempr = fasp_blas_array_dotprod(m,t,t);
199 
200  if ( ABS(tempr) > SMALLREAL ) {
201  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
202  }
203  else {
204  omega = 0.0;
205  if ( prtlvl >= PRINT_SOME ) ITS_DIVZERO;
206  }
207 
208  // delu = alpha pp + omega sp
209  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
210 
211  // u = u + delu
212  fasp_blas_array_axpy(m,1.0,sp,uval);
213 
214  // r = s - omega t
215  fasp_blas_array_axpy(m,-omega,t,s);
216  fasp_array_cp(m,s,r);
217 
218  // beta = (r,rho)/(rp,rho)
219  temp2 = temp1;
220  temp1 = fasp_blas_array_dotprod(m,r,rho);
221 
222  if ( ABS(temp2) > SMALLREAL ) {
223  beta = (temp1*alpha)/(temp2*omega);
224  }
225  else {
226  ITS_DIVZERO; goto RESTORE_BESTSOL;
227  }
228 
229  // p = p - omega z
230  fasp_blas_array_axpy(m,-omega,z,p);
231 
232  // p = r + beta p
233  fasp_blas_array_axpby(m,1.0,r,beta,p);
234 
235  // compute difference
236  normd = fasp_blas_array_norm2(m,sp);
237  normu = fasp_blas_array_norm2(m,uval);
238  reldiff = normd/normu;
239 
240  if ( normd < TOL_s ) {
241  ITS_SMALLSP; goto FINISHED;
242  }
243 
244  // compute residuals
245  switch (stop_type) {
246  case STOP_REL_RES:
247  absres = fasp_blas_array_norm2(m,r);
248  relres = absres/normr0;
249  break;
250  case STOP_REL_PRECRES:
251  if ( pc == NULL )
252  fasp_array_cp(m,r,z);
253  else
254  pc->fct(r,z,pc->data);
255  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
256  relres = absres/normr0;
257  break;
258  case STOP_MOD_REL_RES:
259  absres = fasp_blas_array_norm2(m,r);
260  relres = absres/normu;
261  break;
262  }
263 
264  // safety net check: save the best-so-far solution
265  if ( fasp_dvec_isnan(u) ) {
266  // If the solution is NAN, restrore the best solution
267  absres = BIGREAL;
268  goto RESTORE_BESTSOL;
269  }
270 
271  if ( absres < absres_best - maxdiff) {
272  absres_best = absres;
273  iter_best = iter;
274  fasp_array_cp(m,uval,u_best);
275  }
276 
277  // compute reducation factor of residual ||r||
278  factor = absres/absres0;
279 
280  // output iteration information if needed
281  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
282 
283  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
284  infnormu = fasp_blas_array_norminf(m, uval);
285  if ( infnormu <= sol_inf_tol ) {
286  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
287  iter = ERROR_SOLVER_SOLSTAG;
288  goto FINISHED;
289  }
290 
291  // Check II: if staggenated, try to restart
292  if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
293 
294  if ( prtlvl >= PRINT_MORE ) {
295  ITS_DIFFRES(reldiff,relres);
296  ITS_RESTART;
297  }
298 
299  // re-init iteration param
300  fasp_array_cp(m,bval,r);
301  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
302 
303  // pp = precond(p)
304  fasp_array_cp(m,r,p);
305  if ( pc != NULL )
306  pc->fct(p,pp,pc->data); /* Apply preconditioner */
307  else
308  fasp_array_cp(m,p,pp); /* No preconditioner */
309 
310  // rho = r* := r
311  fasp_array_cp(m,r,rho);
312  temp1 = fasp_blas_array_dotprod(m,r,rho);
313 
314  // compute residuals
315  switch (stop_type) {
316  case STOP_REL_RES:
317  absres = fasp_blas_array_norm2(m,r);
318  relres = absres/normr0;
319  break;
320  case STOP_REL_PRECRES:
321  if ( pc != NULL )
322  pc->fct(r,z,pc->data);
323  else
324  fasp_array_cp(m,r,z);
325  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
326  relres = absres/normr0;
327  break;
328  case STOP_MOD_REL_RES:
329  absres = fasp_blas_array_norm2(m,r);
330  relres = absres/normu;
331  break;
332  }
333 
334  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
335 
336  if ( relres < tol )
337  break;
338  else {
339  if ( stag >= MaxStag ) {
340  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
341  iter = ERROR_SOLVER_STAG;
342  goto FINISHED;
343  }
344  ++stag;
345  ++restart_step;
346  }
347 
348  } // end of stagnation check!
349 
350  // Check III: prevent false convergence
351  if ( relres < tol ) {
352  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
353 
354  // re-init iteration param
355  fasp_array_cp(m,bval,r);
356  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
357 
358  // pp = precond(p)
359  fasp_array_cp(m,r,p);
360  if ( pc != NULL )
361  pc->fct(p,pp,pc->data); /* Apply preconditioner */
362  else
363  fasp_array_cp(m,p,pp); /* No preconditioner */
364 
365  // rho = r* := r
366  fasp_array_cp(m,r,rho);
367  temp1 = fasp_blas_array_dotprod(m,r,rho);
368 
369  // compute residuals
370  switch (stop_type) {
371  case STOP_REL_RES:
372  absres = fasp_blas_array_norm2(m,r);
373  relres = absres/normr0;
374  break;
375  case STOP_REL_PRECRES:
376  if ( pc != NULL )
377  pc->fct(r,z,pc->data);
378  else
379  fasp_array_cp(m,r,z);
380  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
381  relres = tempr/normr0;
382  break;
383  case STOP_MOD_REL_RES:
384  absres = fasp_blas_array_norm2(m,r);
385  relres = absres/normu;
386  break;
387  }
388 
389  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
390 
391  // check convergence
392  if ( relres < tol ) break;
393 
394  if ( more_step >= MaxRestartStep ) {
395  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
396  iter = ERROR_SOLVER_TOLSMALL;
397  goto FINISHED;
398  }
399  else {
400  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
401  }
402 
403  ++more_step;
404  ++restart_step;
405  } // end if safe guard
406 
407  absres0 = absres;
408 
409  } // end of main BiCGstab loop
410 
411 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
412  if ( iter != iter_best ) {
413 
414  // compute best residual
415  fasp_array_cp(m,b->val,r);
416  fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
417 
418  switch ( stop_type ) {
419  case STOP_REL_RES:
420  absres_best = fasp_blas_array_norm2(m,r);
421  break;
422  case STOP_REL_PRECRES:
423  // z = B(r)
424  if ( pc != NULL )
425  pc->fct(r,z,pc->data); /* Apply preconditioner */
426  else
427  fasp_array_cp(m,r,z); /* No preconditioner */
428  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
429  break;
430  case STOP_MOD_REL_RES:
431  absres_best = fasp_blas_array_norm2(m,r);
432  break;
433  }
434 
435  if ( absres > absres_best + maxdiff ) {
436  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
437  fasp_array_cp(m,u_best,u->val);
438  relres = absres_best / normr0;
439  }
440  }
441 
442 FINISHED: // finish the iterative method
443  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
444 
445  // clean up temp memory
446  fasp_mem_free(work);
447 
448 #if DEBUG_MODE > 0
449  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
450 #endif
451 
452  if ( iter > MaxIt )
453  return ERROR_SOLVER_MAXIT;
454  else
455  return iter;
456 }
457 
480  dvector *b,
481  dvector *u,
482  precond *pc,
483  const REAL tol,
484  const INT MaxIt,
485  const SHORT stop_type,
486  const SHORT prtlvl)
487 {
488  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
489  const INT m = b->row;
490  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
491  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
492  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
493 
494  // local variables
495  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
496  REAL alpha, beta, omega, temp1, temp2;
497  REAL absres0 = BIGREAL, absres = BIGREAL;
498  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
499  REAL reldiff, factor, normd, tempr, infnormu;
500 
501  REAL *uval = u->val, *bval = b->val;
502  INT iter_best = 0; // initial best known iteration
503  REAL absres_best = BIGREAL; // initial best known residual
504 
505  // allocate temp memory (need 8*m REAL)
506  REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
507  REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
508  REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
509 
510 #if DEBUG_MODE > 0
511  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
512  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
513 #endif
514 
515  // r = b-A*u
516  fasp_array_cp(m,bval,r);
517  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
518  absres0 = fasp_blas_array_norm2(m,r);
519 
520  // compute initial relative residual
521  switch (stop_type) {
522  case STOP_REL_RES:
523  normr0 = MAX(SMALLREAL,absres0);
524  relres = absres0/normr0;
525  break;
526  case STOP_REL_PRECRES:
527  normr0 = MAX(SMALLREAL,absres0);
528  relres = absres0/normr0;
529  break;
530  case STOP_MOD_REL_RES:
531  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
532  relres = absres0/normu;
533  break;
534  default:
535  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
536  goto FINISHED;
537  }
538 
539  // if initial residual is small, no need to iterate!
540  if (relres<tol) goto FINISHED;
541 
542  // output iteration information if needed
543  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
544 
545  // rho = r* := r
546  fasp_array_cp(m,r,rho);
547  temp1 = fasp_blas_array_dotprod(m,r,rho);
548 
549  // p = r
550  fasp_array_cp(m,r,p);
551 
552  // main BiCGstab loop
553  while ( iter++ < MaxIt ) {
554 
555  // pp = precond(p)
556  if ( pc != NULL )
557  pc->fct(p,pp,pc->data); /* Apply preconditioner */
558  else
559  fasp_array_cp(m,p,pp); /* No preconditioner */
560 
561  // z = A*pp
562  fasp_blas_dbsr_mxv(A,pp,z);
563 
564  // alpha = (r,rho)/(A*p,rho)
565  temp2 = fasp_blas_array_dotprod(m,z,rho);
566  if ( ABS(temp2) > SMALLREAL ) {
567  alpha = temp1/temp2;
568  }
569  else {
570  ITS_DIVZERO; goto FINISHED;
571  }
572 
573  // s = r - alpha z
574  fasp_array_cp(m,r,s);
575  fasp_blas_array_axpy(m,-alpha,z,s);
576 
577  // sp = precond(s)
578  if ( pc != NULL )
579  pc->fct(s,sp,pc->data); /* Apply preconditioner */
580  else
581  fasp_array_cp(m,s,sp); /* No preconditioner */
582 
583  // t = A*sp;
584  fasp_blas_dbsr_mxv(A,sp,t);
585 
586  // omega = (t,s)/(t,t)
587  tempr = fasp_blas_array_dotprod(m,t,t);
588 
589  if ( ABS(tempr) > SMALLREAL ) {
590  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
591  }
592  else {
593  omega = 0.0;
594  if ( prtlvl >= PRINT_SOME ) ITS_DIVZERO;
595  }
596 
597  // delu = alpha pp + omega sp
598  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
599 
600  // u = u + delu
601  fasp_blas_array_axpy(m,1.0,sp,uval);
602 
603  // r = s - omega t
604  fasp_blas_array_axpy(m,-omega,t,s);
605  fasp_array_cp(m,s,r);
606 
607  // beta = (r,rho)/(rp,rho)
608  temp2 = temp1;
609  temp1 = fasp_blas_array_dotprod(m,r,rho);
610 
611  if ( ABS(temp2) > SMALLREAL ) {
612  beta = (temp1*alpha)/(temp2*omega);
613  }
614  else {
615  ITS_DIVZERO; goto RESTORE_BESTSOL;
616  }
617 
618  // p = p - omega z
619  fasp_blas_array_axpy(m,-omega,z,p);
620 
621  // p = r + beta p
622  fasp_blas_array_axpby(m,1.0,r,beta,p);
623 
624  // compute difference
625  normd = fasp_blas_array_norm2(m,sp);
626  normu = fasp_blas_array_norm2(m,uval);
627  reldiff = normd/normu;
628 
629  if ( normd < TOL_s ) {
630  ITS_SMALLSP; goto FINISHED;
631  }
632 
633  // compute residuals
634  switch (stop_type) {
635  case STOP_REL_RES:
636  absres = fasp_blas_array_norm2(m,r);
637  relres = absres/normr0;
638  break;
639  case STOP_REL_PRECRES:
640  if ( pc == NULL )
641  fasp_array_cp(m,r,z);
642  else
643  pc->fct(r,z,pc->data);
644  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
645  relres = absres/normr0;
646  break;
647  case STOP_MOD_REL_RES:
648  absres = fasp_blas_array_norm2(m,r);
649  relres = absres/normu;
650  break;
651  }
652 
653  // safety net check: save the best-so-far solution
654  if ( fasp_dvec_isnan(u) ) {
655  // If the solution is NAN, restrore the best solution
656  absres = BIGREAL;
657  goto RESTORE_BESTSOL;
658  }
659 
660  if ( absres < absres_best - maxdiff) {
661  absres_best = absres;
662  iter_best = iter;
663  fasp_array_cp(m,uval,u_best);
664  }
665 
666  // compute reducation factor of residual ||r||
667  factor = absres/absres0;
668 
669  // output iteration information if needed
670  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
671 
672  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
673  infnormu = fasp_blas_array_norminf(m, uval);
674  if ( infnormu <= sol_inf_tol ) {
675  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
676  iter = ERROR_SOLVER_SOLSTAG;
677  goto FINISHED;
678  }
679 
680  // Check II: if staggenated, try to restart
681  if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
682 
683  if ( prtlvl >= PRINT_MORE ) {
684  ITS_DIFFRES(reldiff,relres);
685  ITS_RESTART;
686  }
687 
688  // re-init iteration param
689  fasp_array_cp(m,bval,r);
690  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
691 
692  // pp = precond(p)
693  fasp_array_cp(m,r,p);
694  if ( pc != NULL )
695  pc->fct(p,pp,pc->data); /* Apply preconditioner */
696  else
697  fasp_array_cp(m,p,pp); /* No preconditioner */
698 
699  // rho = r* := r
700  fasp_array_cp(m,r,rho);
701  temp1 = fasp_blas_array_dotprod(m,r,rho);
702 
703  // compute residuals
704  switch (stop_type) {
705  case STOP_REL_RES:
706  absres = fasp_blas_array_norm2(m,r);
707  relres = absres/normr0;
708  break;
709  case STOP_REL_PRECRES:
710  if ( pc != NULL )
711  pc->fct(r,z,pc->data);
712  else
713  fasp_array_cp(m,r,z);
714  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
715  relres = absres/normr0;
716  break;
717  case STOP_MOD_REL_RES:
718  absres = fasp_blas_array_norm2(m,r);
719  relres = absres/normu;
720  break;
721  }
722 
723  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
724 
725  if ( relres < tol )
726  break;
727  else {
728  if ( stag >= MaxStag ) {
729  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
730  iter = ERROR_SOLVER_STAG;
731  goto FINISHED;
732  }
733  ++stag;
734  ++restart_step;
735  }
736 
737  } // end of stagnation check!
738 
739  // Check III: prevent false convergence
740  if ( relres < tol ) {
741  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
742 
743  // re-init iteration param
744  fasp_array_cp(m,bval,r);
745  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
746 
747  // pp = precond(p)
748  fasp_array_cp(m,r,p);
749  if ( pc != NULL )
750  pc->fct(p,pp,pc->data); /* Apply preconditioner */
751  else
752  fasp_array_cp(m,p,pp); /* No preconditioner */
753 
754  // rho = r* := r
755  fasp_array_cp(m,r,rho);
756  temp1 = fasp_blas_array_dotprod(m,r,rho);
757 
758  // compute residuals
759  switch (stop_type) {
760  case STOP_REL_RES:
761  absres = fasp_blas_array_norm2(m,r);
762  relres = absres/normr0;
763  break;
764  case STOP_REL_PRECRES:
765  if ( pc != NULL )
766  pc->fct(r,z,pc->data);
767  else
768  fasp_array_cp(m,r,z);
769  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
770  relres = tempr/normr0;
771  break;
772  case STOP_MOD_REL_RES:
773  absres = fasp_blas_array_norm2(m,r);
774  relres = absres/normu;
775  break;
776  }
777 
778  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
779 
780  // check convergence
781  if ( relres < tol ) break;
782 
783  if ( more_step >= MaxRestartStep ) {
784  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
785  iter = ERROR_SOLVER_TOLSMALL;
786  goto FINISHED;
787  }
788  else {
789  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
790  }
791 
792  ++more_step;
793  ++restart_step;
794  } // end if safe guard
795 
796  absres0 = absres;
797 
798  } // end of main BiCGstab loop
799 
800 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
801  if ( iter != iter_best ) {
802 
803  // compute best residual
804  fasp_array_cp(m,b->val,r);
805  fasp_blas_dbsr_aAxpy(-1.0,A,u_best,r);
806 
807  switch ( stop_type ) {
808  case STOP_REL_RES:
809  absres_best = fasp_blas_array_norm2(m,r);
810  break;
811  case STOP_REL_PRECRES:
812  // z = B(r)
813  if ( pc != NULL )
814  pc->fct(r,z,pc->data); /* Apply preconditioner */
815  else
816  fasp_array_cp(m,r,z); /* No preconditioner */
817  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
818  break;
819  case STOP_MOD_REL_RES:
820  absres_best = fasp_blas_array_norm2(m,r);
821  break;
822  }
823 
824  if ( absres > absres_best + maxdiff ) {
825  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
826  fasp_array_cp(m,u_best,u->val);
827  relres = absres_best / normr0;
828  }
829  }
830 
831 FINISHED: // finish the iterative method
832  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
833 
834  // clean up temp memory
835  fasp_mem_free(work);
836 
837 #if DEBUG_MODE > 0
838  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
839 #endif
840 
841  if ( iter > MaxIt )
842  return ERROR_SOLVER_MAXIT;
843  else
844  return iter;
845 }
846 
869  dvector *b,
870  dvector *u,
871  precond *pc,
872  const REAL tol,
873  const INT MaxIt,
874  const SHORT stop_type,
875  const SHORT prtlvl)
876 {
877  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
878  const INT m = b->row;
879  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
880  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
881  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
882 
883  // local variables
884  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
885  REAL alpha, beta, omega, temp1, temp2;
886  REAL absres0 = BIGREAL, absres = BIGREAL;
887  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
888  REAL reldiff, factor, normd, tempr, infnormu;
889 
890  REAL *uval = u->val, *bval = b->val;
891  INT iter_best = 0; // initial best known iteration
892  REAL absres_best = BIGREAL; // initial best known residual
893 
894  // allocate temp memory (need 8*m REAL)
895  REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
896  REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
897  REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
898 
899 #if DEBUG_MODE > 0
900  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
901  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
902 #endif
903 
904  // r = b-A*u
905  fasp_array_cp(m,bval,r);
906  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
907  absres0 = fasp_blas_array_norm2(m,r);
908 
909  // compute initial relative residual
910  switch (stop_type) {
911  case STOP_REL_RES:
912  normr0 = MAX(SMALLREAL,absres0);
913  relres = absres0/normr0;
914  break;
915  case STOP_REL_PRECRES:
916  normr0 = MAX(SMALLREAL,absres0);
917  relres = absres0/normr0;
918  break;
919  case STOP_MOD_REL_RES:
920  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
921  relres = absres0/normu;
922  break;
923  default:
924  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
925  goto FINISHED;
926  }
927 
928  // if initial residual is small, no need to iterate!
929  if (relres<tol) goto FINISHED;
930 
931  // output iteration information if needed
932  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
933 
934  // rho = r* := r
935  fasp_array_cp(m,r,rho);
936  temp1 = fasp_blas_array_dotprod(m,r,rho);
937 
938  // p = r
939  fasp_array_cp(m,r,p);
940 
941  // main BiCGstab loop
942  while ( iter++ < MaxIt ) {
943 
944  // pp = precond(p)
945  if ( pc != NULL )
946  pc->fct(p,pp,pc->data); /* Apply preconditioner */
947  else
948  fasp_array_cp(m,p,pp); /* No preconditioner */
949 
950  // z = A*pp
951  fasp_blas_bdcsr_mxv(A,pp,z);
952 
953  // alpha = (r,rho)/(A*p,rho)
954  temp2 = fasp_blas_array_dotprod(m,z,rho);
955  if ( ABS(temp2) > SMALLREAL ) {
956  alpha = temp1/temp2;
957  }
958  else {
959  ITS_DIVZERO; goto FINISHED;
960  }
961 
962  // s = r - alpha z
963  fasp_array_cp(m,r,s);
964  fasp_blas_array_axpy(m,-alpha,z,s);
965 
966  // sp = precond(s)
967  if ( pc != NULL )
968  pc->fct(s,sp,pc->data); /* Apply preconditioner */
969  else
970  fasp_array_cp(m,s,sp); /* No preconditioner */
971 
972  // t = A*sp;
973  fasp_blas_bdcsr_mxv(A,sp,t);
974 
975  // omega = (t,s)/(t,t)
976  tempr = fasp_blas_array_dotprod(m,t,t);
977 
978  if ( ABS(tempr) > SMALLREAL ) {
979  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
980  }
981  else {
982  omega = 0.0;
983  if ( prtlvl >= PRINT_SOME ) ITS_DIVZERO;
984  }
985 
986  // delu = alpha pp + omega sp
987  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
988 
989  // u = u + delu
990  fasp_blas_array_axpy(m,1.0,sp,uval);
991 
992  // r = s - omega t
993  fasp_blas_array_axpy(m,-omega,t,s);
994  fasp_array_cp(m,s,r);
995 
996  // beta = (r,rho)/(rp,rho)
997  temp2 = temp1;
998  temp1 = fasp_blas_array_dotprod(m,r,rho);
999 
1000  if ( ABS(temp2) > SMALLREAL ) {
1001  beta = (temp1*alpha)/(temp2*omega);
1002  }
1003  else {
1004  ITS_DIVZERO; goto RESTORE_BESTSOL;
1005  }
1006 
1007  // p = p - omega z
1008  fasp_blas_array_axpy(m,-omega,z,p);
1009 
1010  // p = r + beta p
1011  fasp_blas_array_axpby(m,1.0,r,beta,p);
1012 
1013  // compute difference
1014  normd = fasp_blas_array_norm2(m,sp);
1015  normu = fasp_blas_array_norm2(m,uval);
1016  reldiff = normd/normu;
1017 
1018  if ( normd < TOL_s ) {
1019  ITS_SMALLSP; goto FINISHED;
1020  }
1021 
1022  // compute residuals
1023  switch (stop_type) {
1024  case STOP_REL_RES:
1025  absres = fasp_blas_array_norm2(m,r);
1026  relres = absres/normr0;
1027  break;
1028  case STOP_REL_PRECRES:
1029  if ( pc == NULL )
1030  fasp_array_cp(m,r,z);
1031  else
1032  pc->fct(r,z,pc->data);
1033  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1034  relres = absres/normr0;
1035  break;
1036  case STOP_MOD_REL_RES:
1037  absres = fasp_blas_array_norm2(m,r);
1038  relres = absres/normu;
1039  break;
1040  }
1041 
1042  // safety net check: save the best-so-far solution
1043  if ( fasp_dvec_isnan(u) ) {
1044  // If the solution is NAN, restrore the best solution
1045  absres = BIGREAL;
1046  goto RESTORE_BESTSOL;
1047  }
1048 
1049  if ( absres < absres_best - maxdiff) {
1050  absres_best = absres;
1051  iter_best = iter;
1052  fasp_array_cp(m,uval,u_best);
1053  }
1054 
1055  // compute reducation factor of residual ||r||
1056  factor = absres/absres0;
1057 
1058  // output iteration information if needed
1059  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1060 
1061  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1062  infnormu = fasp_blas_array_norminf(m, uval);
1063  if ( infnormu <= sol_inf_tol ) {
1064  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1065  iter = ERROR_SOLVER_SOLSTAG;
1066  goto FINISHED;
1067  }
1068 
1069  // Check II: if staggenated, try to restart
1070  if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
1071 
1072  if ( prtlvl >= PRINT_MORE ) {
1073  ITS_DIFFRES(reldiff,relres);
1074  ITS_RESTART;
1075  }
1076 
1077  // re-init iteration param
1078  fasp_array_cp(m,bval,r);
1079  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
1080 
1081  // pp = precond(p)
1082  fasp_array_cp(m,r,p);
1083  if ( pc != NULL )
1084  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1085  else
1086  fasp_array_cp(m,p,pp); /* No preconditioner */
1087 
1088  // rho = r* := r
1089  fasp_array_cp(m,r,rho);
1090  temp1 = fasp_blas_array_dotprod(m,r,rho);
1091 
1092  // compute residuals
1093  switch (stop_type) {
1094  case STOP_REL_RES:
1095  absres = fasp_blas_array_norm2(m,r);
1096  relres = absres/normr0;
1097  break;
1098  case STOP_REL_PRECRES:
1099  if ( pc != NULL )
1100  pc->fct(r,z,pc->data);
1101  else
1102  fasp_array_cp(m,r,z);
1103  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1104  relres = absres/normr0;
1105  break;
1106  case STOP_MOD_REL_RES:
1107  absres = fasp_blas_array_norm2(m,r);
1108  relres = absres/normu;
1109  break;
1110  }
1111 
1112  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1113 
1114  if ( relres < tol )
1115  break;
1116  else {
1117  if ( stag >= MaxStag ) {
1118  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1119  iter = ERROR_SOLVER_STAG;
1120  goto FINISHED;
1121  }
1122  ++stag;
1123  ++restart_step;
1124  }
1125 
1126  } // end of stagnation check!
1127 
1128  // Check III: prevent false convergence
1129  if ( relres < tol ) {
1130  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
1131 
1132  // re-init iteration param
1133  fasp_array_cp(m,bval,r);
1134  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
1135 
1136  // pp = precond(p)
1137  fasp_array_cp(m,r,p);
1138  if ( pc != NULL )
1139  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1140  else
1141  fasp_array_cp(m,p,pp); /* No preconditioner */
1142 
1143  // rho = r* := r
1144  fasp_array_cp(m,r,rho);
1145  temp1 = fasp_blas_array_dotprod(m,r,rho);
1146 
1147  // compute residuals
1148  switch (stop_type) {
1149  case STOP_REL_RES:
1150  absres = fasp_blas_array_norm2(m,r);
1151  relres = absres/normr0;
1152  break;
1153  case STOP_REL_PRECRES:
1154  if ( pc != NULL )
1155  pc->fct(r,z,pc->data);
1156  else
1157  fasp_array_cp(m,r,z);
1158  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1159  relres = tempr/normr0;
1160  break;
1161  case STOP_MOD_REL_RES:
1162  absres = fasp_blas_array_norm2(m,r);
1163  relres = absres/normu;
1164  break;
1165  }
1166 
1167  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1168 
1169  // check convergence
1170  if ( relres < tol ) break;
1171 
1172  if ( more_step >= MaxRestartStep ) {
1173  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1174  iter = ERROR_SOLVER_TOLSMALL;
1175  goto FINISHED;
1176  }
1177  else {
1178  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
1179  }
1180 
1181  ++more_step;
1182  ++restart_step;
1183  } // end if safe guard
1184 
1185  absres0 = absres;
1186 
1187  } // end of main BiCGstab loop
1188 
1189 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1190  if ( iter != iter_best ) {
1191 
1192  // compute best residual
1193  fasp_array_cp(m,b->val,r);
1194  fasp_blas_bdcsr_aAxpy(-1.0,A,u_best,r);
1195 
1196  switch ( stop_type ) {
1197  case STOP_REL_RES:
1198  absres_best = fasp_blas_array_norm2(m,r);
1199  break;
1200  case STOP_REL_PRECRES:
1201  // z = B(r)
1202  if ( pc != NULL )
1203  pc->fct(r,z,pc->data); /* Apply preconditioner */
1204  else
1205  fasp_array_cp(m,r,z); /* No preconditioner */
1206  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1207  break;
1208  case STOP_MOD_REL_RES:
1209  absres_best = fasp_blas_array_norm2(m,r);
1210  break;
1211  }
1212 
1213  if ( absres > absres_best + maxdiff ) {
1214  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1215  fasp_array_cp(m,u_best,u->val);
1216  relres = absres_best / normr0;
1217  }
1218  }
1219 
1220 FINISHED: // finish the iterative method
1221  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1222 
1223  // clean up temp memory
1224  fasp_mem_free(work);
1225 
1226 #if DEBUG_MODE > 0
1227  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1228 #endif
1229 
1230  if ( iter > MaxIt )
1231  return ERROR_SOLVER_MAXIT;
1232  else
1233  return iter;
1234 }
1235 
1258  dvector *b,
1259  dvector *u,
1260  precond *pc,
1261  const REAL tol,
1262  const INT MaxIt,
1263  const SHORT stop_type,
1264  const SHORT prtlvl)
1265 {
1266  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
1267  const INT m = b->row;
1268  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1269  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
1270  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
1271 
1272  // local variables
1273  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
1274  REAL alpha, beta, omega, temp1, temp2;
1275  REAL absres0 = BIGREAL, absres = BIGREAL;
1276  REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
1277  REAL reldiff, factor, normd, tempr, infnormu;
1278 
1279  REAL *uval = u->val, *bval = b->val;
1280  INT iter_best = 0; // initial best known iteration
1281  REAL absres_best = BIGREAL; // initial best known residual
1282 
1283  // allocate temp memory (need 8*m REAL)
1284  REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
1285  REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
1286  REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
1287 
1288 #if DEBUG_MODE > 0
1289  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
1290  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1291 #endif
1292 
1293  // r = b-A*u
1294  fasp_array_cp(m,bval,r);
1295  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1296  absres0 = fasp_blas_array_norm2(m,r);
1297 
1298  // compute initial relative residual
1299  switch (stop_type) {
1300  case STOP_REL_RES:
1301  normr0 = MAX(SMALLREAL,absres0);
1302  relres = absres0/normr0;
1303  break;
1304  case STOP_REL_PRECRES:
1305  normr0 = MAX(SMALLREAL,absres0);
1306  relres = absres0/normr0;
1307  break;
1308  case STOP_MOD_REL_RES:
1309  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
1310  relres = absres0/normu;
1311  break;
1312  default:
1313  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
1314  goto FINISHED;
1315  }
1316 
1317  // if initial residual is small, no need to iterate!
1318  if (relres<tol) goto FINISHED;
1319 
1320  // output iteration information if needed
1321  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
1322 
1323  // rho = r* := r
1324  fasp_array_cp(m,r,rho);
1325  temp1 = fasp_blas_array_dotprod(m,r,rho);
1326 
1327  // p = r
1328  fasp_array_cp(m,r,p);
1329 
1330  // main BiCGstab loop
1331  while ( iter++ < MaxIt ) {
1332 
1333  // pp = precond(p)
1334  if ( pc != NULL )
1335  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1336  else
1337  fasp_array_cp(m,p,pp); /* No preconditioner */
1338 
1339  // z = A*pp
1340  fasp_blas_dstr_mxv(A,pp,z);
1341 
1342  // alpha = (r,rho)/(A*p,rho)
1343  temp2 = fasp_blas_array_dotprod(m,z,rho);
1344  if ( ABS(temp2) > SMALLREAL ) {
1345  alpha = temp1/temp2;
1346  }
1347  else {
1348  ITS_DIVZERO; goto FINISHED;
1349  }
1350 
1351  // s = r - alpha z
1352  fasp_array_cp(m,r,s);
1353  fasp_blas_array_axpy(m,-alpha,z,s);
1354 
1355  // sp = precond(s)
1356  if ( pc != NULL )
1357  pc->fct(s,sp,pc->data); /* Apply preconditioner */
1358  else
1359  fasp_array_cp(m,s,sp); /* No preconditioner */
1360 
1361  // t = A*sp;
1362  fasp_blas_dstr_mxv(A,sp,t);
1363 
1364  // omega = (t,s)/(t,t)
1365  tempr = fasp_blas_array_dotprod(m,t,t);
1366 
1367  if ( ABS(tempr) > SMALLREAL ) {
1368  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
1369  }
1370  else {
1371  omega = 0.0;
1372  if ( prtlvl >= PRINT_SOME ) ITS_DIVZERO;
1373  }
1374 
1375  // delu = alpha pp + omega sp
1376  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
1377 
1378  // u = u + delu
1379  fasp_blas_array_axpy(m,1.0,sp,uval);
1380 
1381  // r = s - omega t
1382  fasp_blas_array_axpy(m,-omega,t,s);
1383  fasp_array_cp(m,s,r);
1384 
1385  // beta = (r,rho)/(rp,rho)
1386  temp2 = temp1;
1387  temp1 = fasp_blas_array_dotprod(m,r,rho);
1388 
1389  if ( ABS(temp2) > SMALLREAL ) {
1390  beta = (temp1*alpha)/(temp2*omega);
1391  }
1392  else {
1393  ITS_DIVZERO; goto RESTORE_BESTSOL;
1394  }
1395 
1396  // p = p - omega z
1397  fasp_blas_array_axpy(m,-omega,z,p);
1398 
1399  // p = r + beta p
1400  fasp_blas_array_axpby(m,1.0,r,beta,p);
1401 
1402  // compute difference
1403  normd = fasp_blas_array_norm2(m,sp);
1404  normu = fasp_blas_array_norm2(m,uval);
1405  reldiff = normd/normu;
1406 
1407  if ( normd < TOL_s ) {
1408  ITS_SMALLSP; goto FINISHED;
1409  }
1410 
1411  // compute residuals
1412  switch (stop_type) {
1413  case STOP_REL_RES:
1414  absres = fasp_blas_array_norm2(m,r);
1415  relres = absres/normr0;
1416  break;
1417  case STOP_REL_PRECRES:
1418  if ( pc == NULL )
1419  fasp_array_cp(m,r,z);
1420  else
1421  pc->fct(r,z,pc->data);
1422  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1423  relres = absres/normr0;
1424  break;
1425  case STOP_MOD_REL_RES:
1426  absres = fasp_blas_array_norm2(m,r);
1427  relres = absres/normu;
1428  break;
1429  }
1430 
1431  // safety net check: save the best-so-far solution
1432  if ( fasp_dvec_isnan(u) ) {
1433  // If the solution is NAN, restrore the best solution
1434  absres = BIGREAL;
1435  goto RESTORE_BESTSOL;
1436  }
1437 
1438  if ( absres < absres_best - maxdiff) {
1439  absres_best = absres;
1440  iter_best = iter;
1441  fasp_array_cp(m,uval,u_best);
1442  }
1443 
1444  // compute reducation factor of residual ||r||
1445  factor = absres/absres0;
1446 
1447  // output iteration information if needed
1448  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1449 
1450  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1451  infnormu = fasp_blas_array_norminf(m, uval);
1452  if ( infnormu <= sol_inf_tol ) {
1453  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1454  iter = ERROR_SOLVER_SOLSTAG;
1455  goto FINISHED;
1456  }
1457 
1458  // Check II: if staggenated, try to restart
1459  if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
1460 
1461  if ( prtlvl >= PRINT_MORE ) {
1462  ITS_DIFFRES(reldiff,relres);
1463  ITS_RESTART;
1464  }
1465 
1466  // re-init iteration param
1467  fasp_array_cp(m,bval,r);
1468  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1469 
1470  // pp = precond(p)
1471  fasp_array_cp(m,r,p);
1472  if ( pc != NULL )
1473  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1474  else
1475  fasp_array_cp(m,p,pp); /* No preconditioner */
1476 
1477  // rho = r* := r
1478  fasp_array_cp(m,r,rho);
1479  temp1 = fasp_blas_array_dotprod(m,r,rho);
1480 
1481  // compute residuals
1482  switch (stop_type) {
1483  case STOP_REL_RES:
1484  absres = fasp_blas_array_norm2(m,r);
1485  relres = absres/normr0;
1486  break;
1487  case STOP_REL_PRECRES:
1488  if ( pc != NULL )
1489  pc->fct(r,z,pc->data);
1490  else
1491  fasp_array_cp(m,r,z);
1492  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1493  relres = absres/normr0;
1494  break;
1495  case STOP_MOD_REL_RES:
1496  absres = fasp_blas_array_norm2(m,r);
1497  relres = absres/normu;
1498  break;
1499  }
1500 
1501  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1502 
1503  if ( relres < tol )
1504  break;
1505  else {
1506  if ( stag >= MaxStag ) {
1507  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1508  iter = ERROR_SOLVER_STAG;
1509  goto FINISHED;
1510  }
1511  ++stag;
1512  ++restart_step;
1513  }
1514 
1515  } // end of stagnation check!
1516 
1517  // Check III: prevent false convergence
1518  if ( relres < tol ) {
1519  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
1520 
1521  // re-init iteration param
1522  fasp_array_cp(m,bval,r);
1523  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1524 
1525  // pp = precond(p)
1526  fasp_array_cp(m,r,p);
1527  if ( pc != NULL )
1528  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1529  else
1530  fasp_array_cp(m,p,pp); /* No preconditioner */
1531 
1532  // rho = r* := r
1533  fasp_array_cp(m,r,rho);
1534  temp1 = fasp_blas_array_dotprod(m,r,rho);
1535 
1536  // compute residuals
1537  switch (stop_type) {
1538  case STOP_REL_RES:
1539  absres = fasp_blas_array_norm2(m,r);
1540  relres = absres/normr0;
1541  break;
1542  case STOP_REL_PRECRES:
1543  if ( pc != NULL )
1544  pc->fct(r,z,pc->data);
1545  else
1546  fasp_array_cp(m,r,z);
1547  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1548  relres = tempr/normr0;
1549  break;
1550  case STOP_MOD_REL_RES:
1551  absres = fasp_blas_array_norm2(m,r);
1552  relres = absres/normu;
1553  break;
1554  }
1555 
1556  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1557 
1558  // check convergence
1559  if ( relres < tol ) break;
1560 
1561  if ( more_step >= MaxRestartStep ) {
1562  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1563  iter = ERROR_SOLVER_TOLSMALL;
1564  goto FINISHED;
1565  }
1566  else {
1567  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
1568  }
1569 
1570  ++more_step;
1571  ++restart_step;
1572  } // end if safe guard
1573 
1574  absres0 = absres;
1575 
1576  } // end of main BiCGstab loop
1577 
1578 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1579  if ( iter != iter_best ) {
1580 
1581  // compute best residual
1582  fasp_array_cp(m,b->val,r);
1583  fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
1584 
1585  switch ( stop_type ) {
1586  case STOP_REL_RES:
1587  absres_best = fasp_blas_array_norm2(m,r);
1588  break;
1589  case STOP_REL_PRECRES:
1590  // z = B(r)
1591  if ( pc != NULL )
1592  pc->fct(r,z,pc->data); /* Apply preconditioner */
1593  else
1594  fasp_array_cp(m,r,z); /* No preconditioner */
1595  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,z,r)));
1596  break;
1597  case STOP_MOD_REL_RES:
1598  absres_best = fasp_blas_array_norm2(m,r);
1599  break;
1600  }
1601 
1602  if ( absres > absres_best + maxdiff ) {
1603  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1604  fasp_array_cp(m,u_best,u->val);
1605  relres = absres_best / normr0;
1606  }
1607  }
1608 
1609 FINISHED: // finish the iterative method
1610  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1611 
1612  // clean up temp memory
1613  fasp_mem_free(work);
1614 
1615 #if DEBUG_MODE > 0
1616  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1617 #endif
1618 
1619  if ( iter > MaxIt )
1620  return ERROR_SOLVER_MAXIT;
1621  else
1622  return iter;
1623 }
1624 
1625 /*---------------------------------*/
1626 /*-- End of File --*/
1627 /*---------------------------------*/
INT fasp_solver_dbsr_spbcgs(dBSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: spbcgs.c:479
#define REAL
Definition: fasp.h:67
void fasp_blas_dstr_aAxpy(const REAL alpha, dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_str.c:47
void fasp_blas_bdcsr_mxv(block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:155
void fasp_blas_bdcsr_aAxpy(const REAL alpha, block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_bcsr.c:30
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
Structure matrix of REAL type.
Definition: fasp.h:304
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:44
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1005
Vector with n entries of REAL type.
Definition: fasp.h:342
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 PRINT_SOME
Definition: fasp_const.h:81
#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
INT fasp_solver_dstr_spbcgs(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: spbcgs.c:1257
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
void fasp_blas_dstr_mxv(dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_str.c:117
Main header file for FASP.
#define 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
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
INT fasp_solver_bdcsr_spbcgs(block_dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: spbcgs.c:868
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
void fasp_blas_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_spbcgs(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: spbcgs.c:90
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