20#include "fasp_functs.h"
69 int i, j, k, rst = -1;
71 REAL gamma, alpha, beta, checktol;
76 REAL *c = NULL, *z = NULL, *alp = NULL, *tmpx = NULL;
77 REAL *norms = NULL, *r = NULL, *work = NULL;
80 INT Restart =
MIN(restart, MaxIt);
81 LONG worksize = n+2*Restart*n+Restart+Restart;
84 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GCR solver (CSR) ...\n");
87 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
88 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
94 while ( (work == NULL) && (Restart > 5) ) {
95 Restart = Restart - 5;
96 worksize = n+2*Restart*n+Restart+Restart;
100 if ( work == NULL ) {
101 printf(
"### ERROR: No enough memory for GCR! [%s:%d]\n",
102 __FILE__, __LINE__ );
106 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
107 printf(
"### WARNING: GCR restart number set to %d!\n", Restart);
110 r = work; z = r+n; c = z + Restart*n; alp = c + Restart*n; tmpx = alp + Restart;
125 relres = absres/absres0;
128 fasp_itinfo(PrtLvl,StopType,0,relres,sqrt(absres0),0.0);
133 checktol =
MAX(tol*tol*absres0, absres*1.0e-4);
135 while ( iter < MaxIt && sqrt(relres) > tol ) {
139 while ( i < Restart-1 && iter < MaxIt ) {
153 for ( j = 0; j < i; j++ ) {
155 h[i][j] = gamma/h[j][j];
174 absres = absres - alpha*alpha/gamma;
176 if (absres < checktol) {
178 checktol =
MAX(tol*tol*absres0, absres*1.0e-4);
181 relres = absres / absres0;
183 norms[iter] = relres;
185 fasp_itinfo(PrtLvl, StopType, iter, sqrt(relres), sqrt(absres),
186 sqrt(norms[iter]/norms[iter-1]));
188 if (sqrt(relres) < tol)
break;
191 for ( k = i; k >=0; k-- ) {
193 for (j=0; j<k; ++j) {
194 alp[j] -= h[k][j]*tmpx[k];
198 if (rst==0) dense_aAtxpby(n, i+1, z, 1.0, tmpx, 0.0, x->
val);
199 else dense_aAtxpby(n, i+1, z, 1.0, tmpx, 1.0, x->
val);
203 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,sqrt(relres));
206 for (i = 0; i < Restart; i++) {
215 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
256 const SHORT StopType,
263 int i, j, k, rst = -1;
265 REAL gamma, alpha, beta, checktol;
270 REAL *c = NULL, *z = NULL, *alp = NULL, *tmpx = NULL;
271 REAL *norms = NULL, *r = NULL, *work = NULL;
274 INT Restart =
MIN(restart, MaxIt);
275 LONG worksize = n+2*Restart*n+Restart+Restart;
278 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GCR solver (BLC) ...\n");
281 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
282 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
288 while ( (work == NULL) && (Restart > 5) ) {
289 Restart = Restart - 5;
290 worksize = n+2*Restart*n+Restart+Restart;
294 if ( work == NULL ) {
295 printf(
"### ERROR: No enough memory for GCR! [%s:%d]\n",
296 __FILE__, __LINE__ );
300 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
301 printf(
"### WARNING: GCR restart number set to %d!\n", Restart);
304 r = work; z = r+n; c = z + Restart*n; alp = c + Restart*n; tmpx = alp + Restart;
319 relres = absres/absres0;
322 fasp_itinfo(PrtLvl,StopType,0,relres,sqrt(absres0),0.0);
327 checktol =
MAX(tol*tol*absres0, absres*1.0e-4);
329 while ( iter < MaxIt && sqrt(relres) > tol ) {
332 while ( i < Restart && iter < MaxIt ) {
346 for ( j = 0; j < i; j++ ) {
348 h[i][j] = gamma/h[j][j];
367 absres = absres - alpha*alpha/gamma;
369 if (absres < checktol) {
371 checktol =
MAX(tol*tol*absres0, absres*1.0e-4);
374 relres = absres / absres0;
376 norms[iter] = relres;
378 fasp_itinfo(PrtLvl, StopType, iter, sqrt(relres), sqrt(absres),
379 sqrt(norms[iter]/norms[iter-1]));
381 if (sqrt(relres) < tol)
break;
386 for ( k = i; k >=0; k-- ) {
388 for (j=0; j<k; ++j) {
389 alp[j] -= h[k][j]*tmpx[k];
393 if (rst==0) dense_aAtxpby(n, i+1, z, 1.0, tmpx, 0.0, x->
val);
394 else dense_aAtxpby(n, i+1, z, 1.0, tmpx, 1.0, x->
val);
397 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,sqrt(relres));
400 for (i = 0; i < Restart; i++) {
409 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
441static void dense_aAtxpby (
INT n,
453 for (j=1; j<m; j++) {
454 for (i=0; i<n; i++) {
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.
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
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_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
void fasp_blas_darray_ax(const INT n, const REAL a, REAL *x)
x = a*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_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.
INT fasp_solver_dcsr_pgcr(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
A preconditioned GCR method for solving Au=b.
INT fasp_solver_dblc_pgcr(dBLCmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
A preconditioned GCR 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_MAXIT
#define BIGREAL
Some global constants.
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Block REAL CSR matrix format.
Sparse matrix of REAL type in CSR format.
Vector with n entries of REAL type.
REAL * val
actual vector entries
Preconditioner data and action.
void * data
data for preconditioner, void pointer
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer