65#include "fasp_functs.h"
104 const SHORT StopType,
113 INT iter = 0, stag = 1, more_step = 1;
116 REAL reldiff, factor, normuinf;
117 REAL alpha, beta, temp1, temp2;
121 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
124 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling CG solver (CSR) ...\n");
127 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
128 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
141 switch ( StopType ) {
145 relres = absres0/normr0;
150 relres = absres0/normr0;
155 relres = absres0/normu;
158 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
163 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
166 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
172 while ( iter++ < MaxIt ) {
183 ITS_DIVZERO;
goto FINISHED;
193 switch ( StopType ) {
196 relres = absres/normr0;
205 relres = absres/normr0;
209 relres = absres/normu;
214 factor = absres/absres0;
217 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
219 if ( factor > 0.9 ) {
223 if ( normuinf <= sol_inf_tol ) {
234 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
237 ITS_DIFFRES(reldiff,relres);
245 switch ( StopType ) {
248 relres = absres/normr0;
257 relres = absres/normr0;
261 relres = absres/normu;
265 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
270 if ( stag >= MaxStag ) {
284 if ( relres < tol ) {
286 REAL updated_relres = relres;
293 switch ( StopType ) {
296 relres = absres/normr0;
305 relres = absres/normr0;
309 relres = absres/normu;
314 if ( relres < tol )
break;
317 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
320 if ( more_step >= MaxRestartStep ) {
354 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
360 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
396 const SHORT StopType,
405 INT iter = 0, stag = 1, more_step = 1;
408 REAL reldiff, factor, normuinf;
409 REAL alpha, beta, temp1, temp2;
413 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
416 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling CG solver (BSR) ...\n");
419 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
420 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
433 switch ( StopType ) {
437 relres = absres0/normr0;
442 relres = absres0/normr0;
447 relres = absres0/normu;
450 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
455 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
458 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
464 while ( iter++ < MaxIt ) {
475 ITS_DIVZERO;
goto FINISHED;
485 switch ( StopType ) {
488 relres = absres/normr0;
497 relres = absres/normr0;
501 relres = absres/normu;
506 factor = absres/absres0;
509 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
511 if ( factor > 0.9 ) {
515 if ( normuinf <= sol_inf_tol ) {
526 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
529 ITS_DIFFRES(reldiff,relres);
537 switch ( StopType ) {
540 relres = absres/normr0;
549 relres = absres/normr0;
553 relres = absres/normu;
557 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
562 if ( stag >= MaxStag ) {
576 if ( relres < tol ) {
578 REAL updated_relres = relres;
585 switch ( StopType ) {
588 relres = absres/normr0;
597 relres = absres/normr0;
601 relres = absres/normu;
606 if ( relres < tol )
break;
609 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
612 if ( more_step >= MaxRestartStep ) {
646 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
652 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
690 const SHORT StopType,
699 INT iter = 0, stag = 1, more_step = 1;
702 REAL reldiff, factor, normuinf;
703 REAL alpha, beta, temp1, temp2;
707 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
710 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling CG solver (BLC) ...\n");
713 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
714 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
727 switch ( StopType ) {
731 relres = absres0/normr0;
736 relres = absres0/normr0;
741 relres = absres0/normu;
744 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
749 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
752 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
758 while ( iter++ < MaxIt ) {
769 ITS_DIVZERO;
goto FINISHED;
779 switch ( StopType ) {
782 relres = absres/normr0;
791 relres = absres/normr0;
795 relres = absres/normu;
800 factor = absres/absres0;
803 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
805 if ( factor > 0.9 ) {
809 if ( normuinf <= sol_inf_tol ) {
820 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
823 ITS_DIFFRES(reldiff,relres);
831 switch ( StopType ) {
834 relres = absres/normr0;
843 relres = absres/normr0;
847 relres = absres/normu;
851 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
856 if ( stag >= MaxStag ) {
870 if ( relres < tol ) {
872 REAL updated_relres = relres;
879 switch ( StopType ) {
882 relres = absres/normr0;
891 relres = absres/normr0;
895 relres = absres/normu;
900 if ( relres < tol )
break;
903 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
906 if ( more_step >= MaxRestartStep ) {
940 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
946 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
984 const SHORT StopType,
993 INT iter = 0, stag = 1, more_step = 1;
996 REAL reldiff, factor, normuinf;
997 REAL alpha, beta, temp1, temp2;
1001 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
1004 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling CG solver (STR) ...\n");
1007 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1008 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1021 switch ( StopType ) {
1025 relres = absres0/normr0;
1030 relres = absres0/normr0;
1035 relres = absres0/normu;
1038 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1043 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
1046 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
1052 while ( iter++ < MaxIt ) {
1060 alpha = temp1/temp2;
1063 ITS_DIVZERO;
goto FINISHED;
1073 switch ( StopType ) {
1076 relres = absres/normr0;
1085 relres = absres/normr0;
1089 relres = absres/normu;
1094 factor = absres/absres0;
1097 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1099 if ( factor > 0.9 ) {
1103 if ( normuinf <= sol_inf_tol ) {
1114 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
1117 ITS_DIFFRES(reldiff,relres);
1125 switch ( StopType ) {
1128 relres = absres/normr0;
1137 relres = absres/normr0;
1141 relres = absres/normu;
1145 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
1150 if ( stag >= MaxStag ) {
1164 if ( relres < tol ) {
1166 REAL updated_relres = relres;
1173 switch ( StopType ) {
1176 relres = absres/normr0;
1185 relres = absres/normr0;
1189 relres = absres/normu;
1194 if ( relres < tol )
break;
1197 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
1200 if ( more_step >= MaxRestartStep ) {
1234 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1240 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1278 const SHORT StopType,
1287 INT iter = 0, stag, more_step, restart_step;
1290 REAL reldiff, factor, infnormu;
1291 REAL alpha, beta, temp1, temp2;
1295 REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
1298 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling CG solver (MatFree) ...\n");
1301 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1302 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1306 stag=1; more_step=1; restart_step=1;
1322 relres=absres0/normr0;
1327 relres=absres0/normu;
1332 relres=absres0/normr0;
1337 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
1342 while ( iter++ < MaxIt ) {
1359 factor=absres/absres0;
1370 relres=sqrt(
ABS(temp2))/normr0;
1373 relres=absres/normu;
1376 relres=absres/normr0;
1381 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1385 if ( infnormu <= sol_inf_tol ) {
1396 if ( (stag<=MaxStag) & (reldiff<maxdiff) ) {
1399 ITS_DIFFRES(reldiff,relres);
1416 relres=sqrt(
ABS(temp2))/normr0;
1419 relres=absres/normu;
1422 relres=absres/normr0;
1426 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
1431 if ( stag >= MaxStag ) {
1443 if ( relres < tol ) {
1444 if ( PrtLvl >=
PRINT_MORE ) ITS_COMPRES(relres);
1458 relres=sqrt(
ABS(temp2))/normr0;
1462 relres=absres/normu;
1466 relres=absres/normr0;
1470 if ( PrtLvl >=
PRINT_MORE ) ITS_REALRES(relres);
1473 if ( relres < tol )
break;
1475 if ( more_step >= MaxRestartStep ) {
1510 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1516 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
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.
REAL fasp_blas_darray_norminf(const INT n, const REAL *x)
Linf norm of array x.
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_dbsr_pcg(dBSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
INT fasp_solver_dstr_pcg(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
INT fasp_solver_dcsr_pcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
INT fasp_solver_dblc_pcg(dBLCmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
INT fasp_solver_pcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient (CG) method for solving Au=b.
Main header file for the FASP project.
#define SHORT
FASP integer and floating point numbers.
#define MAX(a, b)
Definition of max, min, abs.
#define ERROR_SOLVER_STAG
#define ERROR_SOLVER_MAXIT
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
#define ERROR_SOLVER_SOLSTAG
#define BIGREAL
Some global constants.
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
#define ERROR_SOLVER_TOLSMALL
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