28#include "fasp_functs.h"
78 const INT MIN_ITER = 0;
84 REAL r_norm, r_normb, gamma, t;
89 REAL *c = NULL, *s = NULL, *rs = NULL;
90 REAL *norms = NULL, *r = NULL, *w = NULL;
92 REAL **p = NULL, **hh = NULL;
94 INT Restart =
MIN(restart, MaxIt);
95 INT Restart1 = Restart + 1;
96 LONG worksize = (Restart+4)*(Restart+n)+1-n;
102 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GMRes solver (CSR) ...\n");
105 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
106 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
110 while ( (work == NULL) && (Restart > 5) ) {
111 Restart = Restart - 5;
112 Restart1 = Restart + 1;
113 worksize = (Restart+4)*(Restart+n)+1-n;
117 if ( work == NULL ) {
118 printf(
"### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
122 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
123 printf(
"### WARNING: GMRES restart number set to %d!\n", Restart);
130 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
132 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
134 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
145 relres = r_norm/absres0;
154 relres = r_normb/absres0;
159 relres = absres0/normu;
162 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
167 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
176 while ( iter < MaxIt && relres > tol ) {
186 while ( i < Restart && iter < MaxIt ) {
199 for ( j = 0; j < i; j++ ) {
211 for ( j = 1; j < i; ++j ) {
213 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
214 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
216 t = hh[i][i-1]*hh[i][i-1];
217 t += hh[i-1][i-1]*hh[i-1][i-1];
220 c[i-1] = hh[i-1][i-1] / gamma;
221 s[i-1] = hh[i][i-1] / gamma;
222 rs[i] = -s[i-1]*rs[i-1];
223 rs[i-1] = c[i-1]*rs[i-1];
224 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
226 absres = r_norm = fabs(rs[i]);
228 relres = absres/absres0;
230 norms[iter] = relres;
233 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
234 norms[iter]/norms[iter-1]);
237 if ( relres < tol && iter >= MIN_ITER )
break;
242 rs[i-1] = rs[i-1] / hh[i-1][i-1];
243 for ( k = i-2; k >= 0; k-- ) {
245 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
247 rs[k] = t / hh[k][k];
265 if ( relres < tol && iter >= MIN_ITER ) {
267 REAL computed_relres = relres;
274 switch ( StopType ) {
277 relres = absres/absres0;
285 relres = absres/absres0;
290 relres = absres/normu;
294 norms[iter] = relres;
296 if ( relres < tol ) {
304 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
310 for ( j = i; j > 0; j-- ) {
311 rs[j-1] = -s[j-1]*rs[j];
312 rs[j] = c[j-1]*rs[j];
326 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
337 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
377 const SHORT StopType,
381 const INT MIN_ITER = 0;
387 REAL r_norm, r_normb, gamma, t;
392 REAL *c = NULL, *s = NULL, *rs = NULL;
393 REAL *norms = NULL, *r = NULL, *w = NULL;
395 REAL **p = NULL, **hh = NULL;
397 INT Restart =
MIN(restart, MaxIt);
398 INT Restart1 = Restart + 1;
399 LONG worksize = (Restart+4)*(Restart+n)+1-n;
405 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GMRes solver (BSR) ...\n");
408 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
409 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
413 while ( (work == NULL) && (Restart > 5) ) {
414 Restart = Restart - 5;
415 Restart1 = Restart + 1;
416 worksize = (Restart+4)*(Restart+n)+1-n;
420 if ( work == NULL ) {
421 printf(
"### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
425 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
426 printf(
"### WARNING: GMRES restart number set to %d!\n", Restart);
433 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
435 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
437 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
448 relres = r_norm/absres0;
457 relres = r_normb/absres0;
462 relres = absres0/normu;
465 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
470 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
479 while ( iter < MaxIt && relres > tol ) {
489 while ( i < Restart && iter < MaxIt ) {
502 for ( j = 0; j < i; j++ ) {
514 for ( j = 1; j < i; ++j ) {
516 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
517 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
519 t = hh[i][i-1]*hh[i][i-1];
520 t += hh[i-1][i-1]*hh[i-1][i-1];
523 c[i-1] = hh[i-1][i-1] / gamma;
524 s[i-1] = hh[i][i-1] / gamma;
525 rs[i] = -s[i-1]*rs[i-1];
526 rs[i-1] = c[i-1]*rs[i-1];
527 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
529 absres = r_norm = fabs(rs[i]);
531 relres = absres/absres0;
533 norms[iter] = relres;
536 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
537 norms[iter]/norms[iter-1]);
540 if ( relres < tol && iter >= MIN_ITER )
break;
545 rs[i-1] = rs[i-1] / hh[i-1][i-1];
546 for ( k = i-2; k >= 0; k-- ) {
548 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
550 rs[k] = t / hh[k][k];
568 if ( relres < tol && iter >= MIN_ITER ) {
570 REAL computed_relres = relres;
577 switch ( StopType ) {
580 relres = absres/absres0;
588 relres = absres/absres0;
593 relres = absres/normu;
597 norms[iter] = relres;
599 if ( relres < tol ) {
607 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
614 for ( j = i; j > 0; j-- ) {
615 rs[j-1] = -s[j-1]*rs[j];
616 rs[j] = c[j-1]*rs[j];
631 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
642 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
682 const SHORT StopType,
686 const INT MIN_ITER = 0;
692 REAL r_norm, r_normb, gamma, t;
697 REAL *c = NULL, *s = NULL, *rs = NULL;
698 REAL *norms = NULL, *r = NULL, *w = NULL;
700 REAL **p = NULL, **hh = NULL;
702 INT Restart =
MIN(restart, MaxIt);
703 INT Restart1 = Restart + 1;
704 LONG worksize = (Restart+4)*(Restart+n)+1-n;
710 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GMRes solver (BLC) ...\n");
713 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
714 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
718 while ( (work == NULL) && (Restart > 5) ) {
719 Restart = Restart - 5;
720 Restart1 = Restart + 1;
721 worksize = (Restart+4)*(Restart+n)+1-n;
725 if ( work == NULL ) {
726 printf(
"### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
730 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
731 printf(
"### WARNING: GMRES restart number set to %d!\n", Restart);
738 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
740 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
742 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
753 relres = r_norm/absres0;
762 relres = r_normb/absres0;
767 relres = absres0/normu;
770 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
775 if ( relres < tol || absres0 < 1e-12*tol )
goto FINISHED;
784 while ( iter < MaxIt && relres > tol ) {
794 while ( i < Restart && iter < MaxIt ) {
807 for ( j = 0; j < i; j++ ) {
819 for ( j = 1; j < i; ++j ) {
821 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
822 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
824 t = hh[i][i-1]*hh[i][i-1];
825 t += hh[i-1][i-1]*hh[i-1][i-1];
828 c[i-1] = hh[i-1][i-1] / gamma;
829 s[i-1] = hh[i][i-1] / gamma;
830 rs[i] = -s[i-1]*rs[i-1];
831 rs[i-1] = c[i-1]*rs[i-1];
832 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
834 absres = r_norm = fabs(rs[i]);
836 relres = absres/absres0;
838 norms[iter] = relres;
841 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
842 norms[iter]/norms[iter-1]);
845 if ( relres < tol && iter >= MIN_ITER )
break;
850 rs[i-1] = rs[i-1] / hh[i-1][i-1];
851 for ( k = i-2; k >= 0; k-- ) {
853 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
855 rs[k] = t / hh[k][k];
873 if ( relres < tol && iter >= MIN_ITER ) {
875 REAL computed_relres = relres;
882 switch ( StopType ) {
885 relres = absres/absres0;
893 relres = absres/absres0;
898 relres = absres/normu;
902 norms[iter] = relres;
904 if ( relres < tol ) {
912 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
918 for ( j = i; j > 0; j-- ) {
919 rs[j-1] = -s[j-1]*rs[j];
920 rs[j] = c[j-1]*rs[j];
935 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
946 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
986 const SHORT StopType,
990 const INT MIN_ITER = 0;
996 REAL r_norm, r_normb, gamma, t;
1001 REAL *c = NULL, *s = NULL, *rs = NULL;
1002 REAL *norms = NULL, *r = NULL, *w = NULL;
1004 REAL **p = NULL, **hh = NULL;
1006 INT Restart =
MIN(restart, MaxIt);
1007 INT Restart1 = Restart + 1;
1008 LONG worksize = (Restart+4)*(Restart+n)+1-n;
1014 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GMRes solver (STR) ...\n");
1017 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1018 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1022 while ( (work == NULL) && (Restart > 5) ) {
1023 Restart = Restart - 5;
1024 Restart1 = Restart + 1;
1025 worksize = (Restart+4)*(Restart+n)+1-n;
1029 if ( work == NULL ) {
1030 printf(
"### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
1034 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
1035 printf(
"### WARNING: GMRES restart number set to %d!\n", Restart);
1042 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
1044 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
1046 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
1057 relres = r_norm/absres0;
1066 relres = r_normb/absres0;
1071 relres = absres0/normu;
1074 printf(
"### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1079 if ( relres < tol || absres0 < 1e-312*tol )
goto FINISHED;
1082 fasp_itinfo(PrtLvl,StopType,0,relres,absres0,0.0);
1088 while ( iter < MaxIt && relres > tol ) {
1098 while ( i < Restart && iter < MaxIt ) {
1111 for ( j = 0; j < i; j++ ) {
1123 for ( j = 1; j < i; ++j ) {
1125 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1126 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1128 t = hh[i][i-1]*hh[i][i-1];
1129 t += hh[i-1][i-1]*hh[i-1][i-1];
1132 c[i-1] = hh[i-1][i-1] / gamma;
1133 s[i-1] = hh[i][i-1] / gamma;
1134 rs[i] = -s[i-1]*rs[i-1];
1135 rs[i-1] = c[i-1]*rs[i-1];
1136 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1138 absres = r_norm = fabs(rs[i]);
1140 relres = absres/absres0;
1142 norms[iter] = relres;
1145 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1146 norms[iter]/norms[iter-1]);
1149 if ( relres < tol && iter >= MIN_ITER )
break;
1154 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1155 for ( k = i-2; k >= 0; k-- ) {
1157 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
1159 rs[k] = t / hh[k][k];
1177 if ( relres < tol && iter >= MIN_ITER ) {
1179 REAL computed_relres = relres;
1186 switch ( StopType ) {
1189 relres = absres/absres0;
1197 relres = absres/absres0;
1202 relres = absres/normu;
1206 norms[iter] = relres;
1208 if ( relres < tol ) {
1216 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
1222 for ( j = i; j > 0; j-- ) {
1223 rs[j-1] = -s[j-1]*rs[j];
1224 rs[j] = c[j-1]*rs[j];
1239 if ( PrtLvl >
PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1250 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1253 if ( iter >= MaxIt )
1289 const SHORT restart,
1290 const SHORT StopType,
1294 const INT min_iter = 0;
1301 REAL r_norm, b_norm, den_norm;
1302 REAL epsilon, gamma, t;
1305 REAL *c = NULL, *s = NULL, *rs = NULL;
1306 REAL *norms = NULL, *r = NULL, *w = NULL;
1308 REAL **p = NULL, **hh = NULL;
1310 INT Restart = restart;
1311 INT Restart1 = Restart + 1;
1312 LONG worksize = (Restart+4)*(Restart+n)+1-n;
1315 if ( PrtLvl >
PRINT_NONE ) printf(
"\nCalling GMRes solver (MatFree) ...\n");
1318 printf(
"### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1319 printf(
"### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1326 while ( (work == NULL) && (Restart > 5) ) {
1327 Restart = Restart - 5;
1328 worksize = (Restart+4)*(Restart+n)+1-n;
1330 Restart1 = Restart + 1;
1333 if ( work == NULL ) {
1334 printf(
"### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
1338 if ( PrtLvl >
PRINT_MIN && Restart < restart ) {
1339 printf(
"### WARNING: GMRES restart number set to %d!\n", Restart);
1346 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
1348 for (i = 0; i < Restart1; i ++) p[i] = s + Restart + i*n;
1350 for (i = 0; i < Restart1; i ++) hh[i] = p[Restart] + n + i*Restart;
1362 ITS_PUTNORM(
"right-hand side", b_norm);
1363 ITS_PUTNORM(
"residual", r_norm);
1367 if (b_norm > 0.0) den_norm = b_norm;
1368 else den_norm = r_norm;
1370 epsilon = tol*den_norm;
1373 while (iter < MaxIt) {
1376 if (r_norm == 0.0) {
1384 if (r_norm <= epsilon && iter >= min_iter) {
1389 if (r_norm <= epsilon) {
1403 while (i < Restart && iter < MaxIt) {
1416 for (j = 0; j < i; j ++) {
1428 for (j = 1; j < i; ++j) {
1430 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1431 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1433 t= hh[i][i-1]*hh[i][i-1];
1434 t+= hh[i-1][i-1]*hh[i-1][i-1];
1436 if (gamma == 0.0) gamma = epsmac;
1437 c[i-1] = hh[i-1][i-1] / gamma;
1438 s[i-1] = hh[i][i-1] / gamma;
1439 rs[i] = -s[i-1]*rs[i-1];
1440 rs[i-1] = c[i-1]*rs[i-1];
1441 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1442 r_norm = fabs(rs[i]);
1444 norms[iter] = r_norm;
1447 fasp_itinfo(PrtLvl,StopType,iter,norms[iter]/b_norm,
1448 norms[iter],norms[iter]/norms[iter-1]);
1451 fasp_itinfo(PrtLvl,StopType,iter,norms[iter],norms[iter],
1452 norms[iter]/norms[iter-1]);
1456 if (r_norm <= epsilon && iter >= min_iter) {
1462 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1463 for (k = i-2; k >= 0; k --) {
1465 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
1468 rs[k] = t / hh[k][k];
1483 if (r_norm <= epsilon && iter >= min_iter) {
1488 if (r_norm <= epsilon) {
1498 for (j = i; j > 0; j--) {
1499 rs[j-1] = -s[j-1]*rs[j];
1500 rs[j] = c[j-1]*rs[j];
1513 if (PrtLvl >
PRINT_NONE) ITS_FINAL(iter,MaxIt,r_norm);
1524 printf(
"### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1534static double estimate_spectral_radius (
const double **A,
int n,
size_t k = 20)
1536 double *x = (
double *)malloc(n*
sizeof(
double));
1537 double *y = (
double *)malloc(n*
sizeof(
double));
1538 double *z = (
double *)malloc(n*
sizeof(
double));
1550 for(
size_t i = 0; i < k; i++)
1554 for(i1= 0; i1 <n; i1++) x[i1] *= t;
1558 for(i1= 0; i1 <n; i1++) {
1560 for(j1= 0; j1 <n; j1++) t += A[i1][j1] * x[j1];
1563 for(i1= 0; i1 <n; i1++) z[i1] = x[i1];
1564 for(i1= 0; i1 <n; i1++) x[i1] = y[i1];
1565 for(i1= 0; i1 <n; i1++) y[i1] = z[i1];
1579static double spectral_radius (
dCSRmat *A,
1580 const SHORT restart)
1583 const INT MIN_ITER = 0;
1587 INT Restart1 = restart + 1;
1590 REAL r_norm, den_norm;
1591 REAL epsilon, gamma, t;
1593 REAL *c = NULL, *s = NULL, *rs = NULL;
1594 REAL *norms = NULL, *r = NULL, *w = NULL;
1595 REAL **p = NULL, **hh = NULL;
1605 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + restart;
1607 for (i = 0; i < Restart1; i ++) p[i] = s + restart + i*n;
1608 for (i = 0; i < Restart1; i ++) hh[i] = p[restart] + n + i*restart;
1618 for (j = 0; j < n; j ++) p[0][j] *= t;
1620 int maxiter =
MIN(n, restart) ;
1621 for ( j = 0; j < maxiter; j++ ) {
1622 fasp_blas_bdbsr_mxv(A, p[j], p[j+1]);
1624 for( i = 0; i <= j; i++ ) {
1630 if ( hh[j+1][j] < 1e-10)
break;
1632 for (k = 0; k < n; k ++) p[j+1][k] *= t;
1637 for (i = 1; i < j; i ++) H[i] = H[i-1] + j;
1640 for(
size_t row = 0; row < j; row++ )
1641 for(
size_t col = 0; col < j; col++ )
1642 H[row][col] = hh[row][col];
1644 double spectral_radius = estimate_spectral_radius( H, j, 20);
1656 return spectral_radius;
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.
void fasp_dvec_rand(const INT n, dvector *x)
Generate fake random REAL vector in the range from 0 to 1.
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_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_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_dcsr_pgmres(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Right preconditioned GMRES method for solving Au=b.
INT fasp_solver_dblc_pgmres(dBLCmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
INT fasp_solver_pgmres(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Solve "Ax=b" using PGMRES (right preconditioned) iterative method.
INT fasp_solver_dbsr_pgmres(dBSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
INT fasp_solver_dstr_pgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES 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 STOP_REL_RES
Definition of iterative solver stopping criteria types.
#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.
INT row
row number of matrix A, m
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