Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pbcgs.c
Go to the documentation of this file.
1 
53 #include <math.h>
54 
55 #include "fasp.h"
56 #include "fasp_functs.h"
57 #include "itsolver_util.inl"
58 
59 /*---------------------------------*/
60 /*-- Public Functions --*/
61 /*---------------------------------*/
62 
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; // stagnation tolerance
100  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
101  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
102 
103  // local variables
104  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
105  REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
106  REAL normd = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
107  REAL reldiff, factor, infnormu;
108  REAL alpha, beta, omega, temp1, temp2, tempr;
109 
110  REAL *uval=u->val, *bval=b->val;
111 
112  // allocate temp memory (need 8*m REAL)
113  REAL *work=(REAL *)fasp_mem_calloc(8*m,sizeof(REAL));
114  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
115  REAL *rho=t+m, *pp=rho+m, *s=pp+m, *sp=s+m;
116 
117 #if DEBUG_MODE > 0
118  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
119  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
120 #endif
121 
122  // r = b-A*u
123  fasp_array_cp(m,bval,r);
124  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
125  absres0 = fasp_blas_array_norm2(m,r);
126 
127  // compute initial relative residual
128  switch (stop_type) {
129  case STOP_REL_RES:
130  normr0 = MAX(SMALLREAL,absres0);
131  relres = absres0/normr0;
132  break;
133  case STOP_REL_PRECRES:
134  normr0 = MAX(SMALLREAL,absres0);
135  relres = absres0/normr0;
136  break;
137  case STOP_MOD_REL_RES:
138  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
139  relres = absres0/normu;
140  break;
141  default:
142  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
143  goto FINISHED;
144  }
145 
146  // if initial residual is small, no need to iterate!
147  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
148 
149  // output iteration information if needed
150  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
151 
152  // shadow residual rho = r* := r
153  fasp_array_cp(m,r,rho);
154  temp1 = fasp_blas_array_dotprod(m,r,rho);
155 
156  // p = r
157  fasp_array_cp(m,r,p);
158 
159  // main BiCGstab loop
160  while ( iter++ < MaxIt ) {
161 
162  // pp = precond(p)
163  if ( pc != NULL )
164  pc->fct(p,pp,pc->data); /* Apply preconditioner */
165  else
166  fasp_array_cp(m,p,pp); /* No preconditioner */
167 
168  // z = A*pp
169  fasp_blas_dcsr_mxv(A,pp,z);
170 
171  // alpha = (r,rho)/(A*p,rho)
172  temp2 = fasp_blas_array_dotprod(m,z,rho);
173  if ( ABS(temp2) > SMALLREAL2 ) {
174  alpha = temp1/temp2;
175  }
176  else { // Possible breakdown
177  ITS_DIVZERO; goto FINISHED;
178  }
179 
180  // s = r - alpha z
181  fasp_array_cp(m,r,s);
182  fasp_blas_array_axpy(m,-alpha,z,s);
183 
184  // sp = precond(s)
185  if ( pc != NULL )
186  pc->fct(s,sp,pc->data); /* Apply preconditioner */
187  else
188  fasp_array_cp(m,s,sp); /* No preconditioner */
189 
190  // t = A*sp;
191  fasp_blas_dcsr_mxv(A,sp,t);
192 
193  // omega = (t,s)/(t,t)
194  tempr = fasp_blas_array_dotprod(m,t,t);
195  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
196 
197  // diffu = alpha pp + omega sp
198  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
199 
200  // u = u + diffu
201  fasp_blas_array_axpy(m,1.0,sp,uval);
202 
203  // r = s - omega t
204  fasp_blas_array_axpy(m,-omega,t,s);
205  fasp_array_cp(m,s,r);
206 
207  // beta = (r,rho)/(rp,rho)
208  temp2 = temp1;
209  temp1 = fasp_blas_array_dotprod(m,r,rho);
210 
211  if ( ABS(temp2) > SMALLREAL2 || ABS(omega) > SMALLREAL2 ) {
212  beta = (temp1*alpha)/(temp2*omega);
213  }
214  else { // Possible breakdown
215  ITS_DIVZERO; goto FINISHED;
216  }
217 
218  // p = p - omega z
219  fasp_blas_array_axpy(m,-omega,z,p);
220 
221  // p = r + beta p
222  fasp_blas_array_axpby(m,1.0,r,beta,p);
223 
224  // compute difference
225  normd = fasp_blas_array_norm2(m,sp);
226  if ( normd < TOL_s ) { // Possible breakdown?
227  ITS_SMALLSP; goto FINISHED;
228  }
229 
230  normu = fasp_blas_array_norm2(m,uval);
231  reldiff = normd/normu;
232 
233  // compute residuals
234  switch (stop_type) {
235  case STOP_REL_RES:
236  absres = fasp_blas_array_norm2(m,r);
237  relres = absres/normr0;
238  break;
239  case STOP_REL_PRECRES:
240  if ( pc == NULL )
241  fasp_array_cp(m,r,z);
242  else
243  pc->fct(r,z,pc->data);
244  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
245  relres = absres/normr0;
246  break;
247  case STOP_MOD_REL_RES:
248  absres = fasp_blas_array_norm2(m,r);
249  relres = absres/normu;
250  break;
251  }
252 
253  // compute reduction factor of residual ||r||
254  factor = absres/absres0;
255 
256  // output iteration information if needed
257  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
258 
259  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
260  infnormu = fasp_blas_array_norminf(m, uval);
261  if ( infnormu <= sol_inf_tol ) {
262  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
263  iter = ERROR_SOLVER_SOLSTAG;
264  goto FINISHED;
265  }
266 
267  // Check II: if stagnated, try to restart
268  if ( (stag<=MaxStag) && (reldiff<maxdiff) ) {
269 
270  if ( prtlvl >= PRINT_MORE ) {
271  ITS_DIFFRES(reldiff,relres);
272  ITS_RESTART;
273  }
274 
275  // re-init iteration param
276  fasp_array_cp(m,bval,r);
277  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
278 
279  // pp = precond(p)
280  fasp_array_cp(m,r,p);
281  if ( pc != NULL )
282  pc->fct(p,pp,pc->data); /* Apply preconditioner */
283  else
284  fasp_array_cp(m,p,pp); /* No preconditioner */
285 
286  // rho = r* := r
287  fasp_array_cp(m,r,rho);
288  temp1 = fasp_blas_array_dotprod(m,r,rho);
289 
290  // compute residuals
291  switch (stop_type) {
292  case STOP_REL_RES:
293  absres = fasp_blas_array_norm2(m,r);
294  relres = absres/normr0;
295  break;
296  case STOP_REL_PRECRES:
297  if ( pc != NULL )
298  pc->fct(r,z,pc->data);
299  else
300  fasp_array_cp(m,r,z);
301  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
302  relres = absres/normr0;
303  break;
304  case STOP_MOD_REL_RES:
305  absres = fasp_blas_array_norm2(m,r);
306  relres = absres/normu;
307  break;
308  }
309 
310  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
311 
312  if ( relres < tol )
313  break;
314  else {
315  if ( stag >= MaxStag ) {
316  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
317  iter = ERROR_SOLVER_STAG;
318  goto FINISHED;
319  }
320  ++stag;
321  ++restart_step;
322  }
323 
324  } // end of stagnation check!
325 
326  // Check III: prevent false convergence
327  if ( relres < tol ) {
328 
329  REAL computed_relres = relres;
330 
331  // compute current residual
332  fasp_array_cp(m,bval,r);
333  fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
334 
335  // pp = precond(p)
336  fasp_array_cp(m,r,p);
337  if ( pc != NULL )
338  pc->fct(p,pp,pc->data); /* Apply preconditioner */
339  else
340  fasp_array_cp(m,p,pp); /* No preconditioner */
341 
342  // rho = r* := r
343  fasp_array_cp(m,r,rho);
344  temp1 = fasp_blas_array_dotprod(m,r,rho);
345 
346  // compute residuals
347  switch (stop_type) {
348  case STOP_REL_RES:
349  absres = fasp_blas_array_norm2(m,r);
350  relres = absres/normr0;
351  break;
352  case STOP_REL_PRECRES:
353  if ( pc != NULL )
354  pc->fct(r,z,pc->data);
355  else
356  fasp_array_cp(m,r,z);
357  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
358  relres = tempr/normr0;
359  break;
360  case STOP_MOD_REL_RES:
361  absres = fasp_blas_array_norm2(m,r);
362  relres = absres/normu;
363  break;
364  }
365 
366  // check convergence
367  if ( relres < tol ) break;
368 
369  if ( prtlvl >= PRINT_MORE ) {
370  ITS_COMPRES(computed_relres); ITS_REALRES(relres);
371  }
372 
373  if ( more_step >= MaxRestartStep ) {
374  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
375  iter = ERROR_SOLVER_TOLSMALL;
376  goto FINISHED;
377  }
378  else {
379  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
380  }
381 
382  ++more_step;
383  ++restart_step;
384  } // end of false convergence check
385 
386  absres0 = absres;
387 
388  } // end of main BiCGstab loop
389 
390 FINISHED: // finish the iterative method
391  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
392 
393  // clean up temp memory
394  fasp_mem_free(work);
395 
396 #if DEBUG_MODE > 0
397  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
398 #endif
399 
400  if ( iter > MaxIt )
401  return ERROR_SOLVER_MAXIT;
402  else
403  return iter;
404 }
405 
432  dvector *b,
433  dvector *u,
434  precond *pc,
435  const REAL tol,
436  const INT MaxIt,
437  const SHORT stop_type,
438  const SHORT prtlvl)
439 {
440  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
441  const INT m = b->row;
442  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
443  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
444  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
445 
446  // local variables
447  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
448  REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
449  REAL normd = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
450  REAL reldiff, factor, infnormu;
451  REAL alpha, beta, omega, temp1, temp2, tempr;
452 
453  REAL *uval=u->val, *bval=b->val;
454 
455  // allocate temp memory (need 8*m REAL)
456  REAL *work=(REAL *)fasp_mem_calloc(8*m,sizeof(REAL));
457  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
458  REAL *rho=t+m, *pp=rho+m, *s=pp+m, *sp=s+m;
459 
460 #if DEBUG_MODE > 0
461  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
462  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
463 #endif
464 
465  // r = b-A*u
466  fasp_array_cp(m,bval,r);
467  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
468  absres0 = fasp_blas_array_norm2(m,r);
469 
470  // compute initial relative residual
471  switch (stop_type) {
472  case STOP_REL_RES:
473  normr0 = MAX(SMALLREAL,absres0);
474  relres = absres0/normr0;
475  break;
476  case STOP_REL_PRECRES:
477  normr0 = MAX(SMALLREAL,absres0);
478  relres = absres0/normr0;
479  break;
480  case STOP_MOD_REL_RES:
481  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
482  relres = absres0/normu;
483  break;
484  default:
485  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
486  goto FINISHED;
487  }
488 
489  // if initial residual is small, no need to iterate!
490  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
491 
492  // output iteration information if needed
493  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
494 
495  // shadow residual rho = r* := r
496  fasp_array_cp(m,r,rho);
497  temp1 = fasp_blas_array_dotprod(m,r,rho);
498 
499  // p = r
500  fasp_array_cp(m,r,p);
501 
502  // main BiCGstab loop
503  while ( iter++ < MaxIt ) {
504 
505  // pp = precond(p)
506  if ( pc != NULL )
507  pc->fct(p,pp,pc->data); /* Apply preconditioner */
508  else
509  fasp_array_cp(m,p,pp); /* No preconditioner */
510 
511  // z = A*pp
512  fasp_blas_dbsr_mxv(A,pp,z);
513 
514  // alpha = (r,rho)/(A*p,rho)
515  temp2 = fasp_blas_array_dotprod(m,z,rho);
516  if ( ABS(temp2) > SMALLREAL2 ) {
517  alpha = temp1/temp2;
518  }
519  else { // Possible breakdown
520  ITS_DIVZERO; goto FINISHED;
521  }
522 
523  // s = r - alpha z
524  fasp_array_cp(m,r,s);
525  fasp_blas_array_axpy(m,-alpha,z,s);
526 
527  // sp = precond(s)
528  if ( pc != NULL )
529  pc->fct(s,sp,pc->data); /* Apply preconditioner */
530  else
531  fasp_array_cp(m,s,sp); /* No preconditioner */
532 
533  // t = A*sp;
534  fasp_blas_dbsr_mxv(A,sp,t);
535 
536  // omega = (t,s)/(t,t)
537  tempr = fasp_blas_array_dotprod(m,t,t);
538  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
539 
540  // diffu = alpha pp + omega sp
541  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
542 
543  // u = u + diffu
544  fasp_blas_array_axpy(m,1.0,sp,uval);
545 
546  // r = s - omega t
547  fasp_blas_array_axpy(m,-omega,t,s);
548  fasp_array_cp(m,s,r);
549 
550  // beta = (r,rho)/(rp,rho)
551  temp2 = temp1;
552  temp1 = fasp_blas_array_dotprod(m,r,rho);
553 
554  if ( ABS(temp2) > SMALLREAL2 || ABS(omega) > SMALLREAL2 ) {
555  beta = (temp1*alpha)/(temp2*omega);
556  }
557  else { // Possible breakdown
558  ITS_DIVZERO; goto FINISHED;
559  }
560 
561  // p = p - omega z
562  fasp_blas_array_axpy(m,-omega,z,p);
563 
564  // p = r + beta p
565  fasp_blas_array_axpby(m,1.0,r,beta,p);
566 
567  // compute difference
568  normd = fasp_blas_array_norm2(m,sp);
569  if ( normd < TOL_s ) { // Possible breakdown?
570  ITS_SMALLSP; goto FINISHED;
571  }
572 
573  normu = fasp_blas_array_norm2(m,uval);
574  reldiff = normd/normu;
575 
576  // compute residuals
577  switch (stop_type) {
578  case STOP_REL_RES:
579  absres = fasp_blas_array_norm2(m,r);
580  relres = absres/normr0;
581  break;
582  case STOP_REL_PRECRES:
583  if ( pc == NULL )
584  fasp_array_cp(m,r,z);
585  else
586  pc->fct(r,z,pc->data);
587  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
588  relres = absres/normr0;
589  break;
590  case STOP_MOD_REL_RES:
591  absres = fasp_blas_array_norm2(m,r);
592  relres = absres/normu;
593  break;
594  }
595 
596  // compute reduction factor of residual ||r||
597  factor = absres/absres0;
598 
599  // output iteration information if needed
600  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
601 
602  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
603  infnormu = fasp_blas_array_norminf(m, uval);
604  if ( infnormu <= sol_inf_tol ) {
605  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
606  iter = ERROR_SOLVER_SOLSTAG;
607  goto FINISHED;
608  }
609 
610  // Check II: if stagnated, try to restart
611  if ( (stag<=MaxStag) && (reldiff<maxdiff) ) {
612 
613  if ( prtlvl >= PRINT_MORE ) {
614  ITS_DIFFRES(reldiff,relres);
615  ITS_RESTART;
616  }
617 
618  // re-init iteration param
619  fasp_array_cp(m,bval,r);
620  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
621 
622  // pp = precond(p)
623  fasp_array_cp(m,r,p);
624  if ( pc != NULL )
625  pc->fct(p,pp,pc->data); /* Apply preconditioner */
626  else
627  fasp_array_cp(m,p,pp); /* No preconditioner */
628 
629  // rho = r* := r
630  fasp_array_cp(m,r,rho);
631  temp1 = fasp_blas_array_dotprod(m,r,rho);
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  pc->fct(r,z,pc->data);
642  else
643  fasp_array_cp(m,r,z);
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  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
654 
655  if ( relres < tol )
656  break;
657  else {
658  if ( stag >= MaxStag ) {
659  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
660  iter = ERROR_SOLVER_STAG;
661  goto FINISHED;
662  }
663  ++stag;
664  ++restart_step;
665  }
666 
667  } // end of stagnation check!
668 
669  // Check III: prevent false convergence
670  if ( relres < tol ) {
671 
672  REAL computed_relres = relres;
673 
674  // compute current residual
675  fasp_array_cp(m,bval,r);
676  fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
677 
678  // pp = precond(p)
679  fasp_array_cp(m,r,p);
680  if ( pc != NULL )
681  pc->fct(p,pp,pc->data); /* Apply preconditioner */
682  else
683  fasp_array_cp(m,p,pp); /* No preconditioner */
684 
685  // rho = r* := r
686  fasp_array_cp(m,r,rho);
687  temp1 = fasp_blas_array_dotprod(m,r,rho);
688 
689  // compute residuals
690  switch (stop_type) {
691  case STOP_REL_RES:
692  absres = fasp_blas_array_norm2(m,r);
693  relres = absres/normr0;
694  break;
695  case STOP_REL_PRECRES:
696  if ( pc != NULL )
697  pc->fct(r,z,pc->data);
698  else
699  fasp_array_cp(m,r,z);
700  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
701  relres = tempr/normr0;
702  break;
703  case STOP_MOD_REL_RES:
704  absres = fasp_blas_array_norm2(m,r);
705  relres = absres/normu;
706  break;
707  }
708 
709  // check convergence
710  if ( relres < tol ) break;
711 
712  if ( prtlvl >= PRINT_MORE ) {
713  ITS_COMPRES(computed_relres); ITS_REALRES(relres);
714  }
715 
716  if ( more_step >= MaxRestartStep ) {
717  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
718  iter = ERROR_SOLVER_TOLSMALL;
719  goto FINISHED;
720  }
721  else {
722  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
723  }
724 
725  ++more_step;
726  ++restart_step;
727  } // end of false convergence check
728 
729  absres0 = absres;
730 
731  } // end of main BiCGstab loop
732 
733 FINISHED: // finish the iterative method
734  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
735 
736  // clean up temp memory
737  fasp_mem_free(work);
738 
739 #if DEBUG_MODE > 0
740  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
741 #endif
742 
743  if ( iter > MaxIt )
744  return ERROR_SOLVER_MAXIT;
745  else
746  return iter;
747 }
748 
775  dvector *b,
776  dvector *u,
777  precond *pc,
778  const REAL tol,
779  const INT MaxIt,
780  const SHORT stop_type,
781  const SHORT prtlvl)
782 {
783  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
784  const INT m = b->row;
785  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
786  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
787  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
788 
789  // local variables
790  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
791  REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
792  REAL normd = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
793  REAL reldiff, factor, infnormu;
794  REAL alpha, beta, omega, temp1, temp2, tempr;
795 
796  REAL *uval=u->val, *bval=b->val;
797 
798  // allocate temp memory (need 8*m REAL)
799  REAL *work=(REAL *)fasp_mem_calloc(8*m,sizeof(REAL));
800  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
801  REAL *rho=t+m, *pp=rho+m, *s=pp+m, *sp=s+m;
802 
803 #if DEBUG_MODE > 0
804  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
805  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
806 #endif
807 
808  // r = b-A*u
809  fasp_array_cp(m,bval,r);
810  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
811  absres0 = fasp_blas_array_norm2(m,r);
812 
813  // compute initial relative residual
814  switch (stop_type) {
815  case STOP_REL_RES:
816  normr0 = MAX(SMALLREAL,absres0);
817  relres = absres0/normr0;
818  break;
819  case STOP_REL_PRECRES:
820  normr0 = MAX(SMALLREAL,absres0);
821  relres = absres0/normr0;
822  break;
823  case STOP_MOD_REL_RES:
824  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
825  relres = absres0/normu;
826  break;
827  default:
828  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
829  goto FINISHED;
830  }
831 
832  // if initial residual is small, no need to iterate!
833  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
834 
835  // output iteration information if needed
836  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
837 
838  // shadow residual rho = r* := r
839  fasp_array_cp(m,r,rho);
840  temp1 = fasp_blas_array_dotprod(m,r,rho);
841 
842  // p = r
843  fasp_array_cp(m,r,p);
844 
845  // main BiCGstab loop
846  while ( iter++ < MaxIt ) {
847 
848  // pp = precond(p)
849  if ( pc != NULL )
850  pc->fct(p,pp,pc->data); /* Apply preconditioner */
851  else
852  fasp_array_cp(m,p,pp); /* No preconditioner */
853 
854  // z = A*pp
855  fasp_blas_bdcsr_mxv(A,pp,z);
856 
857  // alpha = (r,rho)/(A*p,rho)
858  temp2 = fasp_blas_array_dotprod(m,z,rho);
859  if ( ABS(temp2) > SMALLREAL2 ) {
860  alpha = temp1/temp2;
861  }
862  else { // Possible breakdown
863  ITS_DIVZERO; goto FINISHED;
864  }
865 
866  // s = r - alpha z
867  fasp_array_cp(m,r,s);
868  fasp_blas_array_axpy(m,-alpha,z,s);
869 
870  // sp = precond(s)
871  if ( pc != NULL )
872  pc->fct(s,sp,pc->data); /* Apply preconditioner */
873  else
874  fasp_array_cp(m,s,sp); /* No preconditioner */
875 
876  // t = A*sp;
877  fasp_blas_bdcsr_mxv(A,sp,t);
878 
879  // omega = (t,s)/(t,t)
880  tempr = fasp_blas_array_dotprod(m,t,t);
881  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
882 
883  // diffu = alpha pp + omega sp
884  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
885 
886  // u = u + diffu
887  fasp_blas_array_axpy(m,1.0,sp,uval);
888 
889  // r = s - omega t
890  fasp_blas_array_axpy(m,-omega,t,s);
891  fasp_array_cp(m,s,r);
892 
893  // beta = (r,rho)/(rp,rho)
894  temp2 = temp1;
895  temp1 = fasp_blas_array_dotprod(m,r,rho);
896 
897  if ( ABS(temp2) > SMALLREAL2 || ABS(omega) > SMALLREAL2 ) {
898  beta = (temp1*alpha)/(temp2*omega);
899  }
900  else { // Possible breakdown
901  ITS_DIVZERO; goto FINISHED;
902  }
903 
904  // p = p - omega z
905  fasp_blas_array_axpy(m,-omega,z,p);
906 
907  // p = r + beta p
908  fasp_blas_array_axpby(m,1.0,r,beta,p);
909 
910  // compute difference
911  normd = fasp_blas_array_norm2(m,sp);
912  if ( normd < TOL_s ) { // Possible breakdown?
913  ITS_SMALLSP; goto FINISHED;
914  }
915 
916  normu = fasp_blas_array_norm2(m,uval);
917  reldiff = normd/normu;
918 
919  // compute residuals
920  switch (stop_type) {
921  case STOP_REL_RES:
922  absres = fasp_blas_array_norm2(m,r);
923  relres = absres/normr0;
924  break;
925  case STOP_REL_PRECRES:
926  if ( pc == NULL )
927  fasp_array_cp(m,r,z);
928  else
929  pc->fct(r,z,pc->data);
930  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
931  relres = absres/normr0;
932  break;
933  case STOP_MOD_REL_RES:
934  absres = fasp_blas_array_norm2(m,r);
935  relres = absres/normu;
936  break;
937  }
938 
939  // compute reduction factor of residual ||r||
940  factor = absres/absres0;
941 
942  // output iteration information if needed
943  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
944 
945  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
946  infnormu = fasp_blas_array_norminf(m, uval);
947  if ( infnormu <= sol_inf_tol ) {
948  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
949  iter = ERROR_SOLVER_SOLSTAG;
950  goto FINISHED;
951  }
952 
953  // Check II: if stagnated, try to restart
954  if ( (stag<=MaxStag) && (reldiff<maxdiff) ) {
955 
956  if ( prtlvl >= PRINT_MORE ) {
957  ITS_DIFFRES(reldiff,relres);
958  ITS_RESTART;
959  }
960 
961  // re-init iteration param
962  fasp_array_cp(m,bval,r);
963  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
964 
965  // pp = precond(p)
966  fasp_array_cp(m,r,p);
967  if ( pc != NULL )
968  pc->fct(p,pp,pc->data); /* Apply preconditioner */
969  else
970  fasp_array_cp(m,p,pp); /* No preconditioner */
971 
972  // rho = r* := r
973  fasp_array_cp(m,r,rho);
974  temp1 = fasp_blas_array_dotprod(m,r,rho);
975 
976  // compute residuals
977  switch (stop_type) {
978  case STOP_REL_RES:
979  absres = fasp_blas_array_norm2(m,r);
980  relres = absres/normr0;
981  break;
982  case STOP_REL_PRECRES:
983  if ( pc != NULL )
984  pc->fct(r,z,pc->data);
985  else
986  fasp_array_cp(m,r,z);
987  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
988  relres = absres/normr0;
989  break;
990  case STOP_MOD_REL_RES:
991  absres = fasp_blas_array_norm2(m,r);
992  relres = absres/normu;
993  break;
994  }
995 
996  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
997 
998  if ( relres < tol )
999  break;
1000  else {
1001  if ( stag >= MaxStag ) {
1002  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1003  iter = ERROR_SOLVER_STAG;
1004  goto FINISHED;
1005  }
1006  ++stag;
1007  ++restart_step;
1008  }
1009 
1010  } // end of stagnation check!
1011 
1012  // Check III: prevent false convergence
1013  if ( relres < tol ) {
1014 
1015  REAL computed_relres = relres;
1016 
1017  // compute current residual
1018  fasp_array_cp(m,bval,r);
1019  fasp_blas_bdcsr_aAxpy(-1.0,A,uval,r);
1020 
1021  // pp = precond(p)
1022  fasp_array_cp(m,r,p);
1023  if ( pc != NULL )
1024  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1025  else
1026  fasp_array_cp(m,p,pp); /* No preconditioner */
1027 
1028  // rho = r* := r
1029  fasp_array_cp(m,r,rho);
1030  temp1 = fasp_blas_array_dotprod(m,r,rho);
1031 
1032  // compute residuals
1033  switch (stop_type) {
1034  case STOP_REL_RES:
1035  absres = fasp_blas_array_norm2(m,r);
1036  relres = absres/normr0;
1037  break;
1038  case STOP_REL_PRECRES:
1039  if ( pc != NULL )
1040  pc->fct(r,z,pc->data);
1041  else
1042  fasp_array_cp(m,r,z);
1043  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1044  relres = tempr/normr0;
1045  break;
1046  case STOP_MOD_REL_RES:
1047  absres = fasp_blas_array_norm2(m,r);
1048  relres = absres/normu;
1049  break;
1050  }
1051 
1052  // check convergence
1053  if ( relres < tol ) break;
1054 
1055  if ( prtlvl >= PRINT_MORE ) {
1056  ITS_COMPRES(computed_relres); ITS_REALRES(relres);
1057  }
1058 
1059  if ( more_step >= MaxRestartStep ) {
1060  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1061  iter = ERROR_SOLVER_TOLSMALL;
1062  goto FINISHED;
1063  }
1064  else {
1065  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
1066  }
1067 
1068  ++more_step;
1069  ++restart_step;
1070  } // end of false convergence check
1071 
1072  absres0 = absres;
1073 
1074  } // end of main BiCGstab loop
1075 
1076 FINISHED: // finish the iterative method
1077  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1078 
1079  // clean up temp memory
1080  fasp_mem_free(work);
1081 
1082 #if DEBUG_MODE > 0
1083  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1084 #endif
1085 
1086  if ( iter > MaxIt )
1087  return ERROR_SOLVER_MAXIT;
1088  else
1089  return iter;
1090 }
1091 
1118  dvector *b,
1119  dvector *u,
1120  precond *pc,
1121  const REAL tol,
1122  const INT MaxIt,
1123  const SHORT stop_type,
1124  const SHORT prtlvl)
1125 {
1126  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
1127  const INT m = b->row;
1128  const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
1129  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
1130  const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
1131 
1132  // local variables
1133  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
1134  REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
1135  REAL normd = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
1136  REAL reldiff, factor, infnormu;
1137  REAL alpha, beta, omega, temp1, temp2, tempr;
1138 
1139  REAL *uval=u->val, *bval=b->val;
1140 
1141  // allocate temp memory (need 8*m REAL)
1142  REAL *work=(REAL *)fasp_mem_calloc(8*m,sizeof(REAL));
1143  REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
1144  REAL *rho=t+m, *pp=rho+m, *s=pp+m, *sp=s+m;
1145 
1146 #if DEBUG_MODE > 0
1147  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
1148  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1149 #endif
1150 
1151  // r = b-A*u
1152  fasp_array_cp(m,bval,r);
1153  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1154  absres0 = fasp_blas_array_norm2(m,r);
1155 
1156  // compute initial relative residual
1157  switch (stop_type) {
1158  case STOP_REL_RES:
1159  normr0 = MAX(SMALLREAL,absres0);
1160  relres = absres0/normr0;
1161  break;
1162  case STOP_REL_PRECRES:
1163  normr0 = MAX(SMALLREAL,absres0);
1164  relres = absres0/normr0;
1165  break;
1166  case STOP_MOD_REL_RES:
1167  normu = MAX(SMALLREAL,fasp_blas_array_norm2(m,uval));
1168  relres = absres0/normu;
1169  break;
1170  default:
1171  printf("### ERROR: Unrecognised stopping type for %s!\n", __FUNCTION__);
1172  goto FINISHED;
1173  }
1174 
1175  // if initial residual is small, no need to iterate!
1176  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
1177 
1178  // output iteration information if needed
1179  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
1180 
1181  // shadow residual rho = r* := r
1182  fasp_array_cp(m,r,rho);
1183  temp1 = fasp_blas_array_dotprod(m,r,rho);
1184 
1185  // p = r
1186  fasp_array_cp(m,r,p);
1187 
1188  // main BiCGstab loop
1189  while ( iter++ < MaxIt ) {
1190 
1191  // pp = precond(p)
1192  if ( pc != NULL )
1193  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1194  else
1195  fasp_array_cp(m,p,pp); /* No preconditioner */
1196 
1197  // z = A*pp
1198  fasp_blas_dstr_mxv(A,pp,z);
1199 
1200  // alpha = (r,rho)/(A*p,rho)
1201  temp2 = fasp_blas_array_dotprod(m,z,rho);
1202  if ( ABS(temp2) > SMALLREAL2 ) {
1203  alpha = temp1/temp2;
1204  }
1205  else { // Possible breakdown
1206  ITS_DIVZERO; goto FINISHED;
1207  }
1208 
1209  // s = r - alpha z
1210  fasp_array_cp(m,r,s);
1211  fasp_blas_array_axpy(m,-alpha,z,s);
1212 
1213  // sp = precond(s)
1214  if ( pc != NULL )
1215  pc->fct(s,sp,pc->data); /* Apply preconditioner */
1216  else
1217  fasp_array_cp(m,s,sp); /* No preconditioner */
1218 
1219  // t = A*sp;
1220  fasp_blas_dstr_mxv(A,sp,t);
1221 
1222  // omega = (t,s)/(t,t)
1223  tempr = fasp_blas_array_dotprod(m,t,t);
1224  omega = fasp_blas_array_dotprod(m,s,t)/tempr;
1225 
1226  // diffu = alpha pp + omega sp
1227  fasp_blas_array_axpby(m,alpha,pp,omega,sp);
1228 
1229  // u = u + diffu
1230  fasp_blas_array_axpy(m,1.0,sp,uval);
1231 
1232  // r = s - omega t
1233  fasp_blas_array_axpy(m,-omega,t,s);
1234  fasp_array_cp(m,s,r);
1235 
1236  // beta = (r,rho)/(rp,rho)
1237  temp2 = temp1;
1238  temp1 = fasp_blas_array_dotprod(m,r,rho);
1239 
1240  if ( ABS(temp2) > SMALLREAL2 || ABS(omega) > SMALLREAL2 ) {
1241  beta = (temp1*alpha)/(temp2*omega);
1242  }
1243  else { // Possible breakdown
1244  ITS_DIVZERO; goto FINISHED;
1245  }
1246 
1247  // p = p - omega z
1248  fasp_blas_array_axpy(m,-omega,z,p);
1249 
1250  // p = r + beta p
1251  fasp_blas_array_axpby(m,1.0,r,beta,p);
1252 
1253  // compute difference
1254  normd = fasp_blas_array_norm2(m,sp);
1255  if ( normd < TOL_s ) { // Possible breakdown?
1256  ITS_SMALLSP; goto FINISHED;
1257  }
1258 
1259  normu = fasp_blas_array_norm2(m,uval);
1260  reldiff = normd/normu;
1261 
1262  // compute residuals
1263  switch (stop_type) {
1264  case STOP_REL_RES:
1265  absres = fasp_blas_array_norm2(m,r);
1266  relres = absres/normr0;
1267  break;
1268  case STOP_REL_PRECRES:
1269  if ( pc == NULL )
1270  fasp_array_cp(m,r,z);
1271  else
1272  pc->fct(r,z,pc->data);
1273  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1274  relres = absres/normr0;
1275  break;
1276  case STOP_MOD_REL_RES:
1277  absres = fasp_blas_array_norm2(m,r);
1278  relres = absres/normu;
1279  break;
1280  }
1281 
1282  // compute reduction factor of residual ||r||
1283  factor = absres/absres0;
1284 
1285  // output iteration information if needed
1286  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1287 
1288  // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
1289  infnormu = fasp_blas_array_norminf(m, uval);
1290  if ( infnormu <= sol_inf_tol ) {
1291  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1292  iter = ERROR_SOLVER_SOLSTAG;
1293  goto FINISHED;
1294  }
1295 
1296  // Check II: if stagnated, try to restart
1297  if ( (stag<=MaxStag) && (reldiff<maxdiff) ) {
1298 
1299  if ( prtlvl >= PRINT_MORE ) {
1300  ITS_DIFFRES(reldiff,relres);
1301  ITS_RESTART;
1302  }
1303 
1304  // re-init iteration param
1305  fasp_array_cp(m,bval,r);
1306  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1307 
1308  // pp = precond(p)
1309  fasp_array_cp(m,r,p);
1310  if ( pc != NULL )
1311  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1312  else
1313  fasp_array_cp(m,p,pp); /* No preconditioner */
1314 
1315  // rho = r* := r
1316  fasp_array_cp(m,r,rho);
1317  temp1 = fasp_blas_array_dotprod(m,r,rho);
1318 
1319  // compute residuals
1320  switch (stop_type) {
1321  case STOP_REL_RES:
1322  absres = fasp_blas_array_norm2(m,r);
1323  relres = absres/normr0;
1324  break;
1325  case STOP_REL_PRECRES:
1326  if ( pc != NULL )
1327  pc->fct(r,z,pc->data);
1328  else
1329  fasp_array_cp(m,r,z);
1330  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1331  relres = absres/normr0;
1332  break;
1333  case STOP_MOD_REL_RES:
1334  absres = fasp_blas_array_norm2(m,r);
1335  relres = absres/normu;
1336  break;
1337  }
1338 
1339  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1340 
1341  if ( relres < tol )
1342  break;
1343  else {
1344  if ( stag >= MaxStag ) {
1345  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1346  iter = ERROR_SOLVER_STAG;
1347  goto FINISHED;
1348  }
1349  ++stag;
1350  ++restart_step;
1351  }
1352 
1353  } // end of stagnation check!
1354 
1355  // Check III: prevent false convergence
1356  if ( relres < tol ) {
1357 
1358  REAL computed_relres = relres;
1359 
1360  // compute current residual
1361  fasp_array_cp(m,bval,r);
1362  fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1363 
1364  // pp = precond(p)
1365  fasp_array_cp(m,r,p);
1366  if ( pc != NULL )
1367  pc->fct(p,pp,pc->data); /* Apply preconditioner */
1368  else
1369  fasp_array_cp(m,p,pp); /* No preconditioner */
1370 
1371  // rho = r* := r
1372  fasp_array_cp(m,r,rho);
1373  temp1 = fasp_blas_array_dotprod(m,r,rho);
1374 
1375  // compute residuals
1376  switch (stop_type) {
1377  case STOP_REL_RES:
1378  absres = fasp_blas_array_norm2(m,r);
1379  relres = absres/normr0;
1380  break;
1381  case STOP_REL_PRECRES:
1382  if ( pc != NULL )
1383  pc->fct(r,z,pc->data);
1384  else
1385  fasp_array_cp(m,r,z);
1386  absres = sqrt(ABS(fasp_blas_array_dotprod(m,r,z)));
1387  relres = tempr/normr0;
1388  break;
1389  case STOP_MOD_REL_RES:
1390  absres = fasp_blas_array_norm2(m,r);
1391  relres = absres/normu;
1392  break;
1393  }
1394 
1395  // check convergence
1396  if ( relres < tol ) break;
1397 
1398  if ( prtlvl >= PRINT_MORE ) {
1399  ITS_COMPRES(computed_relres); ITS_REALRES(relres);
1400  }
1401 
1402  if ( more_step >= MaxRestartStep ) {
1403  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1404  iter = ERROR_SOLVER_TOLSMALL;
1405  goto FINISHED;
1406  }
1407  else {
1408  if ( prtlvl > PRINT_NONE ) ITS_RESTART;
1409  }
1410 
1411  ++more_step;
1412  ++restart_step;
1413  } // end of false convergence check
1414 
1415  absres0 = absres;
1416 
1417  } // end of main BiCGstab loop
1418 
1419 FINISHED: // finish the iterative method
1420  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1421 
1422  // clean up temp memory
1423  fasp_mem_free(work);
1424 
1425 #if DEBUG_MODE > 0
1426  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1427 #endif
1428 
1429  if ( iter > MaxIt )
1430  return ERROR_SOLVER_MAXIT;
1431  else
1432  return iter;
1433 }
1434 
1435 /*---------------------------------*/
1436 /*-- End of File --*/
1437 /*---------------------------------*/
#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_pbcgs(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.
Definition: pbcgs.c:88
INT fasp_solver_dbsr_pbcgs(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.
Definition: pbcgs.c:431
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
INT fasp_solver_dstr_pbcgs(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.
Definition: pbcgs.c:1117
#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 * 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
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_bdcsr_pbcgs(block_dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
A preconditioned BiCGstab method for solving Au=b.
Definition: pbcgs.c:774
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