29#include "fasp_functs.h"
75 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
76 INT flag, maxstagsteps, half_step=0;
78 REAL alpha,beta,omega,rho,rho1,rtv,tt;
79 REAL normr,normr_act,normph,normx,imin;
80 REAL norm_sh,norm_xhalf,normrmin,factor;
85 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
86 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
87 REAL *t = sh+m, *xmin = t+m;
90 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling BiCGstab solver (CSR) ...\n");
93 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
94 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
115 if ( normr <= tolb ) {
138 for (iter=1;iter <= MaxIt;iter++) {
143 if ((rho ==0.0 )|| (
ABS(rho) >= DBL_MAX )) {
152 beta = (rho/rho1)*(alpha/omega);
154 if ((beta == 0)||(
ABS(beta) > DBL_MAX )) {
175 if (( rtv==0.0 )||(
ABS(rtv) > DBL_MAX )){
182 if (
ABS(alpha) > DBL_MAX ){
190 if (
ABS(alpha)*normph < DBL_EPSILON*normx )
204 factor = absres/absres0;
205 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
208 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
214 if (normr_act <= tolb) {
221 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
222 flag,stag,imin,half_step);
226 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
228 moresteps = moresteps + 1;
229 if (moresteps >= maxmsteps) {
238 if ( stag >= maxstagsteps ) {
243 if ( normr_act < normrmin )
245 normrmin = normr_act;
250 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
251 flag,stag,imin,half_step);
265 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
272 if (
ABS(omega) > DBL_MAX ) {
280 if (
ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
291 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
296 if ( normr_act <= tolb ) {
301 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
303 moresteps = moresteps + 1;
304 if ( moresteps >= maxmsteps ) {
312 if ( normr_act < normrmin ) {
313 normrmin = normr_act;
318 if ( stag >= maxstagsteps ) {
323 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
331 relres = normr_act / n2b;
337 if ( normr <= normr_act) {
343 relres = normr_act/n2b;
347 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
350 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
351 flag,stag,imin,half_step);
357 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
393 const SHORT StopType,
400 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
401 INT flag, maxstagsteps, half_step=0;
403 REAL alpha,beta,omega,rho,rho1,rtv,tt;
404 REAL normr,normr_act,normph,normx,imin;
405 REAL norm_sh,norm_xhalf,normrmin,factor;
410 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
411 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
412 REAL *t = sh+m, *xmin = t+m;
415 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling BiCGstab solver (BSR) ...\n");
418 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
419 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
440 if ( normr <= tolb ) {
463 for (iter=1;iter <= MaxIt;iter++) {
468 if ((rho ==0.0 )|| (
ABS(rho) >= DBL_MAX )) {
477 beta = (rho/rho1)*(alpha/omega);
479 if ((beta == 0)||(
ABS(beta) > DBL_MAX )) {
500 if (( rtv==0.0 )||(
ABS(rtv) > DBL_MAX )) {
507 if (
ABS(alpha) > DBL_MAX ) {
515 if (
ABS(alpha)*normph < DBL_EPSILON*normx )
529 factor = absres/absres0;
530 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
533 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
539 if (normr_act <= tolb) {
546 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
547 flag,stag,imin,half_step);
551 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
553 moresteps = moresteps + 1;
554 if (moresteps >= maxmsteps){
563 if ( stag >= maxstagsteps ) {
568 if ( normr_act < normrmin )
570 normrmin = normr_act;
575 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
576 flag,stag,imin,half_step);
590 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
597 if (
ABS(omega) > DBL_MAX ) {
605 if (
ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
616 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
621 if ( normr_act <= tolb ) {
626 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
628 moresteps = moresteps + 1;
629 if ( moresteps >= maxmsteps ) {
637 if ( normr_act < normrmin ) {
638 normrmin = normr_act;
643 if ( stag >= maxstagsteps )
649 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
657 relres = normr_act / n2b;
663 if ( normr <= normr_act) {
669 relres = normr_act/n2b;
673 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
676 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
677 flag,stag,imin,half_step);
683 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
719 const SHORT StopType,
726 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
727 INT flag, maxstagsteps, half_step=0;
729 REAL alpha,beta,omega,rho,rho1,rtv,tt;
730 REAL normr,normr_act,normph,normx,imin;
731 REAL norm_sh,norm_xhalf,normrmin,factor;
736 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
737 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
738 REAL *t = sh+m, *xmin = t+m;
741 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling BiCGstab solver (BLC) ...\n");
744 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
745 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
766 if ( normr <= tolb ) {
789 for (iter=1;iter <= MaxIt;iter++) {
794 if ((rho ==0.0 )|| (
ABS(rho) >= DBL_MAX )) {
803 beta = (rho/rho1)*(alpha/omega);
805 if ((beta == 0)||(
ABS(beta) > DBL_MAX )) {
826 if (( rtv==0.0 )||(
ABS(rtv) > DBL_MAX )) {
833 if (
ABS(alpha) > DBL_MAX ) {
841 if (
ABS(alpha)*normph < DBL_EPSILON*normx )
855 factor = absres/absres0;
856 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
859 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
865 if (normr_act <= tolb) {
872 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
873 flag,stag,imin,half_step);
877 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
879 moresteps = moresteps + 1;
880 if (moresteps >= maxmsteps) {
889 if ( stag >= maxstagsteps ) {
894 if ( normr_act < normrmin )
896 normrmin = normr_act;
901 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
902 flag,stag,imin,half_step);
916 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
923 if (
ABS(omega) > DBL_MAX ) {
931 if (
ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
942 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
947 if ( normr_act <= tolb ) {
952 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
954 moresteps = moresteps + 1;
955 if ( moresteps >= maxmsteps ) {
963 if ( normr_act < normrmin ) {
964 normrmin = normr_act;
969 if ( stag >= maxstagsteps )
975 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
983 relres = normr_act / n2b;
989 if ( normr <= normr_act) {
995 relres = normr_act/n2b;
999 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1002 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1003 flag,stag,imin,half_step);
1009 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1045 const SHORT StopType,
1052 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
1053 INT flag, maxstagsteps, half_step=0;
1055 REAL alpha,beta,omega,rho,rho1,rtv,tt;
1056 REAL normr,normr_act,normph,normx,imin;
1057 REAL norm_sh,norm_xhalf,normrmin,factor;
1062 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
1063 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
1064 REAL *t = sh+m, *xmin = t+m;
1067 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling BiCGstab solver (STR) ...\n");
1070 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1071 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1092 if ( normr <= tolb ) {
1115 for (iter=1;iter <= MaxIt;iter++) {
1120 if ((rho ==0.0 )|| (
ABS(rho) >= DBL_MAX )) {
1129 beta = (rho/rho1)*(alpha/omega);
1131 if ((beta == 0)||(
ABS(beta) > DBL_MAX )) {
1152 if (( rtv==0.0 )||(
ABS(rtv) > DBL_MAX )) {
1159 if (
ABS(alpha) > DBL_MAX ) {
1167 if (
ABS(alpha)*normph < DBL_EPSILON*normx )
1181 factor = absres/absres0;
1182 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
1185 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
1191 if (normr_act <= tolb){
1198 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1199 flag,stag,imin,half_step);
1203 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1205 moresteps = moresteps + 1;
1206 if (moresteps >= maxmsteps){
1215 if ( stag >= maxstagsteps ) {
1220 if ( normr_act < normrmin )
1222 normrmin = normr_act;
1227 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1228 flag,stag,imin,half_step);
1242 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
1249 if (
ABS(omega) > DBL_MAX ) {
1257 if (
ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
1268 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
1273 if ( normr_act <= tolb ) {
1278 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1280 moresteps = moresteps + 1;
1281 if ( moresteps >= maxmsteps ) {
1289 if ( normr_act < normrmin ) {
1290 normrmin = normr_act;
1295 if ( stag >= maxstagsteps )
1301 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
1309 relres = normr_act / n2b;
1315 if ( normr <= normr_act ) {
1321 relres = normr_act/n2b;
1325 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1328 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1329 flag,stag,imin,half_step);
1335 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1371 const SHORT StopType,
1378 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
1379 INT flag, maxstagsteps, half_step=0;
1381 REAL alpha,beta,omega,rho,rho1,rtv,tt;
1382 REAL normr,normr_act,normph,normx,imin;
1383 REAL norm_sh,norm_xhalf,normrmin,factor;
1388 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
1389 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
1390 REAL *t = sh+m, *xmin = t+m;
1393 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling BiCGstab solver (MatFree) ...\n");
1396 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1397 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1420 if (normr <= tolb) {
1444 for (iter=1;iter <= MaxIt;iter++) {
1449 if ((rho ==0.0 )|| (
ABS(rho) >= DBL_MAX )) {
1458 beta = (rho/rho1)*(alpha/omega);
1460 if ((beta == 0)||(
ABS(beta) > DBL_MAX )) {
1481 if (( rtv==0.0 )||(
ABS(rtv) > DBL_MAX )) {
1488 if (
ABS(alpha) > DBL_MAX ) {
1496 if (
ABS(alpha)*normph < DBL_EPSILON*normx )
1510 factor = absres/absres0;
1511 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
1514 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
1521 if (normr_act <= tolb){
1528 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1529 flag,stag,imin,half_step);
1533 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1535 moresteps = moresteps + 1;
1536 if (moresteps >= maxmsteps){
1545 if ( stag >= maxstagsteps ) {
1550 if ( normr_act < normrmin )
1552 normrmin = normr_act;
1557 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1558 flag,stag,imin,half_step);
1573 if ((tt == 0) ||( tt >= DBL_MAX )) {
1580 if (
ABS(omega) > DBL_MAX )
1589 if (
ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
1600 if ( (normr <= tolb)||(stag >= maxstagsteps)||moresteps )
1607 if (normr_act <= tolb)
1614 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1616 moresteps = moresteps + 1;
1617 if (moresteps >= maxmsteps)
1625 if (normr_act < normrmin)
1627 normrmin = normr_act;
1632 if (stag >= maxstagsteps)
1638 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
1646 relres = normr_act / n2b;
1653 if ( normr <= normr_act ) {
1659 relres = normr_act/n2b;
1663 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1666 printf(
"Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1667 flag,stag,imin,half_step);
1673 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
void fasp_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.
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
void fasp_blas_darray_axpyz(const INT n, const REAL a, const REAL *x, const REAL *y, REAL *z)
z = a*x + y
void fasp_blas_darray_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
REAL fasp_blas_darray_norm2(const INT n, const REAL *x)
L2 norm of array x.
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
void fasp_blas_dblc_mxv(const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
void fasp_blas_dblc_aAxpy(const REAL alpha, const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
void fasp_blas_dbsr_aAxpy(const REAL alpha, const dBSRmat *A, const REAL *x, REAL *y)
Compute y := alpha*A*x + y.
void fasp_blas_dbsr_mxv(const dBSRmat *A, const REAL *x, REAL *y)
Compute y := A*x.
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
void fasp_blas_dstr_aAxpy(const REAL alpha, const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
void fasp_blas_dstr_mxv(const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
INT fasp_solver_dstr_pbcgs(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b for STR matrix.
INT fasp_solver_dblc_pbcgs(dBLCmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b for BLC matrix.
INT fasp_solver_pbcgs(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b.
INT fasp_solver_dbsr_pbcgs(dBSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b for BSR matrix.
INT fasp_solver_dcsr_pbcgs(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b for CSR matrix.
Main header file for the FASP project.
#define SHORT
FASP integer and floating point numbers.
#define ERROR_SOLVER_MAXIT
#define BIGREAL
Some global constants.
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Block REAL CSR matrix format.
Block sparse row storage matrix of REAL type.
Sparse matrix of REAL type in CSR format.
Structure matrix of REAL type.
Vector with n entries of REAL type.
REAL * val
actual vector entries
Matrix-vector multiplication, replace the actual matrix.
void * data
data for MxV, can be a Matrix or something else
void(* fct)(const void *, const REAL *, REAL *)
action for MxV, void function pointer
Preconditioner data and action.
void * data
data for preconditioner, void pointer
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer