Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPvfgmres.c
Go to the documentation of this file.
1
25#include <math.h>
26
27#include "fasp.h"
28#include "fasp_functs.h"
29
30/*---------------------------------*/
31/*-- Declare Private Functions --*/
32/*---------------------------------*/
33
34#include "KryUtil.inl"
35
36/*---------------------------------*/
37/*-- Public Functions --*/
38/*---------------------------------*/
39
68 dvector *b,
69 dvector *x,
70 precond *pc,
71 const REAL tol,
72 const INT MaxIt,
73 const SHORT restart,
74 const SHORT StopType,
75 const SHORT PrtLvl)
76{
77 const INT n = b->row;
78 const INT min_iter = 0;
79
80 //--------------------------------------------//
81 // Newly added parameters to monitor when //
82 // to change the restart parameter //
83 //--------------------------------------------//
84 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
85 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
86
87 // local variables
88 INT iter = 0;
89 int i, j, k; // must be signed! -zcs
90
91 REAL epsmac = SMALLREAL;
92 REAL r_norm, b_norm, den_norm;
93 REAL epsilon, gamma, t;
94 REAL relres, normu, r_normb;
95
96 REAL *c = NULL, *s = NULL, *rs = NULL, *norms = NULL, *r = NULL;
97 REAL **p = NULL, **hh = NULL, **z=NULL;
98
99 REAL cr = 1.0; // convergence rate
100 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
101 INT d = 3; // reduction for restart parameter
102 INT restart_max = restart; // upper bound for restart in each restart cycle
103 INT restart_min = 3; // lower bound for restart in each restart cycle
104
105 INT Restart = restart; // real restart in some fixed restarted cycle
106 INT Restart1 = Restart + 1;
107 LONG worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
108
109 // Output some info for debuging
110 if ( PrtLvl > PRINT_NONE ) printf("\nCalling VFGMRes solver (CSR) ...\n");
111
112#if DEBUG_MODE > 0
113 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
114 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
115#endif
116
117 /* allocate memory and setup temp work space */
118 REAL *work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
119
120 /* check whether memory is enough for GMRES */
121 while ( (work == NULL) && (Restart > 5) ) {
122 Restart = Restart - 5;
123 worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
124 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
125 Restart1 = Restart + 1;
126 }
127
128 if ( work == NULL ) {
129 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
130 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
131 }
132
133 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
134 printf("### WARNING: vFGMRES restart number set to %d!\n", Restart);
135 }
136
137 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
138 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
139 z = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
140 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
141
142 r = work; rs = r + n; c = rs + Restart1; s = c + Restart;
143 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
144 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
145 for ( i = 0; i < Restart1; i++ ) z[i] = hh[Restart] + Restart + i*n;
146
147 /* initialization */
148 fasp_darray_cp(n, b->val, p[0]);
149 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
150
151 b_norm = fasp_blas_darray_norm2(n, b->val);
152 r_norm = fasp_blas_darray_norm2(n, p[0]);
153 norms[0] = r_norm;
154
155 if ( PrtLvl >= PRINT_SOME) {
156 ITS_PUTNORM("right-hand side", b_norm);
157 ITS_PUTNORM("residual", r_norm);
158 }
159
160 if ( b_norm > 0.0 ) den_norm = b_norm;
161 else den_norm = r_norm;
162
163 epsilon = tol*den_norm;
164
165 // if initial residual is small, no need to iterate!
166 if ( r_norm < epsilon || r_norm < 1e-12*tol ) goto FINISHED;
167
168 if ( b_norm > 0.0 ) {
169 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm, norms[iter], 0);
170 }
171 else {
172 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter], 0);
173 }
174
175 /* outer iteration cycle */
176 while ( iter < MaxIt ) {
177
178 rs[0] = r_norm;
179 r_norm_old = r_norm;
180 if ( r_norm == 0.0 ) {
181 fasp_mem_free(work); work = NULL;
182 fasp_mem_free(p); p = NULL;
183 fasp_mem_free(hh); hh = NULL;
184 fasp_mem_free(norms); norms = NULL;
185 fasp_mem_free(z); z = NULL;
186 return iter;
187 }
188
189 //-----------------------------------//
190 // adjust the restart parameter //
191 //-----------------------------------//
192
193 if ( cr > cr_max || iter == 0 ) {
194 Restart = restart_max;
195 }
196 else if ( cr < cr_min ) {
197 // Restart = Restart;
198 }
199 else {
200 if ( Restart - d > restart_min ) Restart -= d;
201 else Restart = restart_max;
202 }
203
204 // Enter the cycle at the first iteration for at least one iteration
205 t = 1.0 / r_norm;
206 fasp_blas_darray_ax(n, t, p[0]);
207 i = 0;
208
209 // RESTART CYCLE (right-preconditioning)
210 while ( i < Restart && iter < MaxIt ) {
211
212 i ++; iter ++;
213
214 /* apply preconditioner */
215 if ( pc == NULL )
216 fasp_darray_cp(n, p[i-1], z[i-1]);
217 else
218 pc->fct(p[i-1], z[i-1], pc->data);
219
220 fasp_blas_dcsr_mxv(A, z[i-1], p[i]);
221
222 /* modified Gram_Schmidt */
223 for ( j = 0; j < i; j++ ) {
224 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
225 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
226 }
227 t = fasp_blas_darray_norm2(n, p[i]);
228 hh[i][i-1] = t;
229 if ( t != 0.0 ) {
230 t = 1.0 / t;
231 fasp_blas_darray_ax(n, t, p[i]);
232 }
233
234 for ( j = 1; j < i; ++j ) {
235 t = hh[j-1][i-1];
236 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
237 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
238 }
239 t = hh[i][i-1] * hh[i][i-1];
240 t += hh[i-1][i-1] * hh[i-1][i-1];
241 gamma = sqrt(t);
242 if (gamma == 0.0) gamma = epsmac;
243 c[i-1] = hh[i-1][i-1] / gamma;
244 s[i-1] = hh[i][i-1] / gamma;
245 rs[i] = -s[i-1] * rs[i-1];
246 rs[i-1] = c[i-1] * rs[i-1];
247 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
248
249 r_norm = fabs(rs[i]);
250 norms[iter] = r_norm;
251
252 if ( b_norm > 0 ) {
253 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm,
254 norms[iter], norms[iter]/norms[iter-1]);
255 }
256 else {
257 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
258 norms[iter]/norms[iter-1]);
259 }
260
261 /* Check: Exit the restart cycle? */
262 if (r_norm <= epsilon && iter >= min_iter) break;
263
264 } /* end of restart cycle */
265
266 /* now compute solution, first solve upper triangular system */
267
268 rs[i-1] = rs[i-1] / hh[i-1][i-1];
269 for ( k = i-2; k >= 0; k-- ) {
270 t = 0.0;
271 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
272
273 t += rs[k];
274 rs[k] = t / hh[k][k];
275 }
276
277 fasp_darray_cp(n, z[i-1], r);
278 fasp_blas_darray_ax(n, rs[i-1], r);
279
280 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], z[j], r);
281
282 fasp_blas_darray_axpy(n, 1.0, r, x->val);
283
284 if ( r_norm <= epsilon && iter >= min_iter ) {
285 fasp_darray_cp(n, b->val, r);
286 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
287 r_norm = fasp_blas_darray_norm2(n, r);
288
289 switch (StopType) {
290 case STOP_REL_RES:
291 relres = r_norm/den_norm;
292 break;
293 case STOP_REL_PRECRES:
294 if ( pc == NULL ) fasp_darray_cp(n, r, p[0]);
295 else pc->fct(r, p[0], pc->data);
296 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
297 relres = r_normb/den_norm;
298 break;
299 case STOP_MOD_REL_RES:
301 relres = r_norm/normu;
302 break;
303 default:
304 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
305 goto FINISHED;
306 }
307
308 if ( relres <= tol ) {
309 break;
310 }
311 else {
312 if ( PrtLvl >= PRINT_SOME ) ITS_FACONV;
313 fasp_darray_cp(n, r, p[0]); i = 0;
314 }
315
316 } /* end of convergence check */
317
318 /* compute residual vector and continue loop */
319 for ( j = i; j > 0; j-- ) {
320 rs[j-1] = -s[j-1]*rs[j];
321 rs[j] = c[j-1]*rs[j];
322 }
323
324 if (i) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
325
326 for ( j = i-1; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
327
328 if (i) {
329 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
330 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
331 }
332
333 //-----------------------------------//
334 // compute the convergence rate //
335 //-----------------------------------//
336 cr = r_norm / r_norm_old;
337
338 } /* end of iteration while loop */
339
340 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,r_norm/den_norm);
341
342FINISHED:
343 /*-------------------------------------------
344 * Free some stuff
345 *------------------------------------------*/
346 fasp_mem_free(work); work = NULL;
347 fasp_mem_free(p); p = NULL;
348 fasp_mem_free(hh); hh = NULL;
349 fasp_mem_free(norms); norms = NULL;
350 fasp_mem_free(z); z = NULL;
351
352#if DEBUG_MODE > 0
353 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
354#endif
355
356 if ( iter >= MaxIt )
357 return ERROR_SOLVER_MAXIT;
358 else
359 return iter;
360}
361
390 dvector *b,
391 dvector *x,
392 precond *pc,
393 const REAL tol,
394 const INT MaxIt,
395 const SHORT restart,
396 const SHORT StopType,
397 const SHORT PrtLvl)
398{
399 const INT n = b->row;
400 const INT min_iter = 0;
401
402 //--------------------------------------------//
403 // Newly added parameters to monitor when //
404 // to change the restart parameter //
405 //--------------------------------------------//
406 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
407 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
408
409 // local variables
410 INT iter = 0;
411 int i, j, k; // must be signed! -zcs
412
413 REAL epsmac = SMALLREAL;
414 REAL r_norm, b_norm, den_norm;
415 REAL epsilon, gamma, t;
416 REAL relres, normu, r_normb;
417
418 REAL *c = NULL, *s = NULL, *rs = NULL, *norms = NULL, *r = NULL;
419 REAL **p = NULL, **hh = NULL, **z=NULL;
420
421 REAL cr = 1.0; // convergence rate
422 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
423 INT d = 3; // reduction for restart parameter
424 INT restart_max = restart; // upper bound for restart in each restart cycle
425 INT restart_min = 3; // lower bound for restart in each restart cycle
426
427 INT Restart = restart; // real restart in some fixed restarted cycle
428 INT Restart1 = Restart + 1;
429 LONG worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
430
431 // Output some info for debuging
432 if ( PrtLvl > PRINT_NONE ) printf("\nCalling VFGMRes solver (BSR) ...\n");
433
434#if DEBUG_MODE > 0
435 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
436 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
437#endif
438
439 /* allocate memory and setup temp work space */
440 REAL *work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
441
442 /* check whether memory is enough for GMRES */
443 while ( (work == NULL) && (Restart > 5) ) {
444 Restart = Restart - 5;
445 worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
446 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
447 Restart1 = Restart + 1;
448 }
449
450 if ( work == NULL ) {
451 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
452 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
453 }
454
455 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
456 printf("### WARNING: vFGMRES restart number set to %d!\n", Restart);
457 }
458
459 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
460 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
461 z = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
462 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
463
464 r = work; rs = r + n; c = rs + Restart1; s = c + Restart;
465 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
466 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
467 for ( i = 0; i < Restart1; i++ ) z[i] = hh[Restart] + Restart + i*n;
468
469 /* initialization */
470 fasp_darray_cp(n, b->val, p[0]);
471 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
472
473 b_norm = fasp_blas_darray_norm2(n, b->val);
474 r_norm = fasp_blas_darray_norm2(n, p[0]);
475 norms[0] = r_norm;
476
477 if ( PrtLvl >= PRINT_SOME) {
478 ITS_PUTNORM("right-hand side", b_norm);
479 ITS_PUTNORM("residual", r_norm);
480 }
481
482 if ( b_norm > 0.0 ) den_norm = b_norm;
483 else den_norm = r_norm;
484
485 epsilon = tol*den_norm;
486
487 // if initial residual is small, no need to iterate!
488 if ( r_norm < epsilon || r_norm < 1e-12*tol ) goto FINISHED;
489
490 if ( b_norm > 0.0 ) {
491 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm, norms[iter], 0);
492 }
493 else {
494 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter], 0);
495 }
496
497 /* outer iteration cycle */
498 while ( iter < MaxIt ) {
499
500 rs[0] = r_norm;
501 r_norm_old = r_norm;
502 if ( r_norm == 0.0 ) {
503 fasp_mem_free(work); work = NULL;
504 fasp_mem_free(p); p = NULL;
505 fasp_mem_free(hh); hh = NULL;
506 fasp_mem_free(norms); norms = NULL;
507 fasp_mem_free(z); z = NULL;
508 return iter;
509 }
510
511 //-----------------------------------//
512 // adjust the restart parameter //
513 //-----------------------------------//
514
515 if ( cr > cr_max || iter == 0 ) {
516 Restart = restart_max;
517 }
518 else if ( cr < cr_min ) {
519 // Restart = Restart;
520 }
521 else {
522 if ( Restart - d > restart_min ) Restart -= d;
523 else Restart = restart_max;
524 }
525
526 // Enter the cycle at the first iteration for at least one iteration
527 t = 1.0 / r_norm;
528 fasp_blas_darray_ax(n, t, p[0]);
529 i = 0;
530
531 // RESTART CYCLE (right-preconditioning)
532 while ( i < Restart && iter < MaxIt ) {
533
534 i ++; iter ++;
535
536 /* apply preconditioner */
537 if ( pc == NULL )
538 fasp_darray_cp(n, p[i-1], z[i-1]);
539 else
540 pc->fct(p[i-1], z[i-1], pc->data);
541
542 fasp_blas_dbsr_mxv(A, z[i-1], p[i]);
543
544 /* modified Gram_Schmidt */
545 for ( j = 0; j < i; j++ ) {
546 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
547 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
548 }
549 t = fasp_blas_darray_norm2(n, p[i]);
550 hh[i][i-1] = t;
551 if ( t != 0.0 ) {
552 t = 1.0 / t;
553 fasp_blas_darray_ax(n, t, p[i]);
554 }
555
556 for ( j = 1; j < i; ++j ) {
557 t = hh[j-1][i-1];
558 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
559 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
560 }
561 t = hh[i][i-1] * hh[i][i-1];
562 t += hh[i-1][i-1] * hh[i-1][i-1];
563 gamma = sqrt(t);
564 if (gamma == 0.0) gamma = epsmac;
565 c[i-1] = hh[i-1][i-1] / gamma;
566 s[i-1] = hh[i][i-1] / gamma;
567 rs[i] = -s[i-1] * rs[i-1];
568 rs[i-1] = c[i-1] * rs[i-1];
569 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
570
571 r_norm = fabs(rs[i]);
572 norms[iter] = r_norm;
573
574 if ( b_norm > 0 ) {
575 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm,
576 norms[iter], norms[iter]/norms[iter-1]);
577 }
578 else {
579 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
580 norms[iter]/norms[iter-1]);
581 }
582
583 /* Check: Exit the restart cycle? */
584 if (r_norm <= epsilon && iter >= min_iter) break;
585
586 } /* end of restart cycle */
587
588 /* now compute solution, first solve upper triangular system */
589
590 rs[i-1] = rs[i-1] / hh[i-1][i-1];
591 for ( k = i-2; k >= 0; k-- ) {
592 t = 0.0;
593 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
594
595 t += rs[k];
596 rs[k] = t / hh[k][k];
597 }
598
599 fasp_darray_cp(n, z[i-1], r);
600 fasp_blas_darray_ax(n, rs[i-1], r);
601
602 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], z[j], r);
603
604 fasp_blas_darray_axpy(n, 1.0, r, x->val);
605
606 if ( r_norm <= epsilon && iter >= min_iter ) {
607 fasp_darray_cp(n, b->val, r);
608 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
609 r_norm = fasp_blas_darray_norm2(n, r);
610
611 switch (StopType) {
612 case STOP_REL_RES:
613 relres = r_norm/den_norm;
614 break;
615 case STOP_REL_PRECRES:
616 if ( pc == NULL ) fasp_darray_cp(n, r, p[0]);
617 else pc->fct(r, p[0], pc->data);
618 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
619 relres = r_normb/den_norm;
620 break;
621 case STOP_MOD_REL_RES:
623 relres = r_norm/normu;
624 break;
625 default:
626 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
627 goto FINISHED;
628 }
629
630 if ( relres <= tol ) {
631 break;
632 }
633 else {
634 if ( PrtLvl >= PRINT_SOME ) ITS_FACONV;
635 fasp_darray_cp(n, r, p[0]); i = 0;
636 }
637
638 } /* end of convergence check */
639
640 /* compute residual vector and continue loop */
641 for ( j = i; j > 0; j-- ) {
642 rs[j-1] = -s[j-1]*rs[j];
643 rs[j] = c[j-1]*rs[j];
644 }
645
646 if (i) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
647
648 for ( j = i-1; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
649
650 if (i) {
651 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
652 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
653 }
654
655 //-----------------------------------//
656 // compute the convergence rate //
657 //-----------------------------------//
658 cr = r_norm / r_norm_old;
659
660 } /* end of iteration while loop */
661
662 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,r_norm/den_norm);
663
664FINISHED:
665 /*-------------------------------------------
666 * Free some stuff
667 *------------------------------------------*/
668 fasp_mem_free(work); work = NULL;
669 fasp_mem_free(p); p = NULL;
670 fasp_mem_free(hh); hh = NULL;
671 fasp_mem_free(norms); norms = NULL;
672 fasp_mem_free(z); z = NULL;
673
674#if DEBUG_MODE > 0
675 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
676#endif
677
678 if ( iter >= MaxIt )
679 return ERROR_SOLVER_MAXIT;
680 else
681 return iter;
682}
683
715 dvector *b,
716 dvector *x,
717 precond *pc,
718 const REAL tol,
719 const INT MaxIt,
720 const SHORT restart,
721 const SHORT StopType,
722 const SHORT PrtLvl)
723{
724 const INT n = b->row;
725 const INT min_iter = 0;
726
727 //--------------------------------------------//
728 // Newly added parameters to monitor when //
729 // to change the restart parameter //
730 //--------------------------------------------//
731 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
732 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
733
734 // local variables
735 INT iter = 0;
736 int i, j, k; // must be signed! -zcs
737
738 REAL epsmac = SMALLREAL;
739 REAL r_norm, b_norm, den_norm;
740 REAL epsilon, gamma, t;
741 REAL relres, normu, r_normb;
742
743 REAL *c = NULL, *s = NULL, *rs = NULL, *norms = NULL, *r = NULL;
744 REAL **p = NULL, **hh = NULL, **z=NULL;
745
746 REAL cr = 1.0; // convergence rate
747 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
748 INT d = 3; // reduction for restart parameter
749 INT restart_max = restart; // upper bound for restart in each restart cycle
750 INT restart_min = 3; // lower bound for restart in each restart cycle
751
752 INT Restart = restart; // real restart in some fixed restarted cycle
753 INT Restart1 = Restart + 1;
754 LONG worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
755
756 // Output some info for debuging
757 if ( PrtLvl > PRINT_NONE ) printf("\nCalling VFGMRes solver (BLC) ...\n");
758
759#if DEBUG_MODE > 0
760 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
761 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
762#endif
763
764 /* allocate memory and setup temp work space */
765 REAL *work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
766
767 /* check whether memory is enough for GMRES */
768 while ( (work == NULL) && (Restart > 5) ) {
769 Restart = Restart - 5;
770 worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
771 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
772 Restart1 = Restart + 1;
773 }
774
775 if ( work == NULL ) {
776 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
777 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
778 }
779
780 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
781 printf("### WARNING: vFGMRES restart number set to %d!\n", Restart);
782 }
783
784 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
785 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
786 z = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
787 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
788
789 r = work; rs = r + n; c = rs + Restart1; s = c + Restart;
790 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
791 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
792 for ( i = 0; i < Restart1; i++ ) z[i] = hh[Restart] + Restart + i*n;
793
794 /* initialization */
795 fasp_darray_cp(n, b->val, p[0]);
796 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
797
798 b_norm = fasp_blas_darray_norm2(n, b->val);
799 r_norm = fasp_blas_darray_norm2(n, p[0]);
800 norms[0] = r_norm;
801
802 if ( PrtLvl >= PRINT_SOME) {
803 ITS_PUTNORM("right-hand side", b_norm);
804 ITS_PUTNORM("residual", r_norm);
805 }
806
807 if ( b_norm > 0.0 ) den_norm = b_norm;
808 else den_norm = r_norm;
809
810 epsilon = tol*den_norm;
811
812 // if initial residual is small, no need to iterate!
813 if ( r_norm < epsilon || r_norm < 1e-12 * tol ) goto FINISHED;
814
815 if ( b_norm > 0.0 ) {
816 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm, norms[iter], 0);
817 }
818 else {
819 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter], 0);
820 }
821
822 /* outer iteration cycle */
823 while ( iter < MaxIt ) {
824
825 rs[0] = r_norm;
826 r_norm_old = r_norm;
827 if ( r_norm == 0.0 ) {
828 fasp_mem_free(work); work = NULL;
829 fasp_mem_free(p); p = NULL;
830 fasp_mem_free(hh); hh = NULL;
831 fasp_mem_free(norms); norms = NULL;
832 fasp_mem_free(z); z = NULL;
833 return iter;
834 }
835
836 //-----------------------------------//
837 // adjust the restart parameter //
838 //-----------------------------------//
839
840 if ( cr > cr_max || iter == 0 ) {
841 Restart = restart_max;
842 }
843 else if ( cr < cr_min ) {
844 // Restart = Restart;
845 }
846 else {
847 if ( Restart - d > restart_min ) Restart -= d;
848 else Restart = restart_max;
849 }
850
851 // Enter the cycle at the first iteration for at least one iteration
852 t = 1.0 / r_norm;
853 fasp_blas_darray_ax(n, t, p[0]);
854 i = 0;
855
856 // RESTART CYCLE (right-preconditioning)
857 while ( i < Restart && iter < MaxIt ) {
858
859 i ++; iter ++;
860
861 /* apply preconditioner */
862 if ( pc == NULL )
863 fasp_darray_cp(n, p[i-1], z[i-1]);
864 else
865 pc->fct(p[i-1], z[i-1], pc->data);
866
867 fasp_blas_dblc_mxv(A, z[i-1], p[i]);
868
869 /* modified Gram_Schmidt */
870 for ( j = 0; j < i; j++ ) {
871 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
872 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
873 }
874 t = fasp_blas_darray_norm2(n, p[i]);
875 hh[i][i-1] = t;
876 if ( t != 0.0 ) {
877 t = 1.0 / t;
878 fasp_blas_darray_ax(n, t, p[i]);
879 }
880
881 for ( j = 1; j < i; ++j ) {
882 t = hh[j-1][i-1];
883 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
884 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
885 }
886 t = hh[i][i-1] * hh[i][i-1];
887 t += hh[i-1][i-1] * hh[i-1][i-1];
888 gamma = sqrt(t);
889 if (gamma == 0.0) gamma = epsmac;
890 c[i-1] = hh[i-1][i-1] / gamma;
891 s[i-1] = hh[i][i-1] / gamma;
892 rs[i] = -s[i-1] * rs[i-1];
893 rs[i-1] = c[i-1] * rs[i-1];
894 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
895
896 r_norm = fabs(rs[i]);
897 norms[iter] = r_norm;
898
899 if ( b_norm > 0 ) {
900 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm,
901 norms[iter], norms[iter]/norms[iter-1]);
902 }
903 else {
904 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
905 norms[iter]/norms[iter-1]);
906 }
907
908 /* Check: Exit the restart cycle? */
909 if (r_norm <= epsilon && iter >= min_iter) break;
910
911 } /* end of restart cycle */
912
913 /* now compute solution, first solve upper triangular system */
914
915 rs[i-1] = rs[i-1] / hh[i-1][i-1];
916 for ( k = i-2; k >= 0; k-- ) {
917 t = 0.0;
918 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
919
920 t += rs[k];
921 rs[k] = t / hh[k][k];
922 }
923
924 fasp_darray_cp(n, z[i-1], r);
925 fasp_blas_darray_ax(n, rs[i-1], r);
926
927 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], z[j], r);
928
929 fasp_blas_darray_axpy(n, 1.0, r, x->val);
930
931 if ( r_norm <= epsilon && iter >= min_iter ) {
932 fasp_darray_cp(n, b->val, r);
933 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
934 r_norm = fasp_blas_darray_norm2(n, r);
935
936 switch (StopType) {
937 case STOP_REL_RES:
938 relres = r_norm/den_norm;
939 break;
940 case STOP_REL_PRECRES:
941 if ( pc == NULL ) fasp_darray_cp(n, r, p[0]);
942 else pc->fct(r, p[0], pc->data);
943 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
944 relres = r_normb/den_norm;
945 break;
946 case STOP_MOD_REL_RES:
948 relres = r_norm/normu;
949 break;
950 default:
951 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
952 goto FINISHED;
953 }
954
955 if ( relres <= tol ) {
956 break;
957 }
958 else {
959 if ( PrtLvl >= PRINT_SOME ) ITS_FACONV;
960 fasp_darray_cp(n, r, p[0]); i = 0;
961 }
962
963 } /* end of convergence check */
964
965 /* compute residual vector and continue loop */
966 for ( j = i; j > 0; j-- ) {
967 rs[j-1] = -s[j-1]*rs[j];
968 rs[j] = c[j-1]*rs[j];
969 }
970
971 if (i) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
972
973 for ( j = i-1; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
974
975 if (i) {
976 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
977 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
978 }
979
980 //-----------------------------------//
981 // compute the convergence rate //
982 //-----------------------------------//
983 cr = r_norm / r_norm_old;
984
985 } /* end of iteration while loop */
986
987 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,r_norm/den_norm);
988
989FINISHED:
990 /*-------------------------------------------
991 * Free some stuff
992 *------------------------------------------*/
993 fasp_mem_free(work); work = NULL;
994 fasp_mem_free(p); p = NULL;
995 fasp_mem_free(hh); hh = NULL;
996 fasp_mem_free(norms); norms = NULL;
997 fasp_mem_free(z); z = NULL;
998
999#if DEBUG_MODE > 0
1000 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1001#endif
1002
1003 if ( iter >= MaxIt )
1004 return ERROR_SOLVER_MAXIT;
1005 else
1006 return iter;
1007}
1008
1037 dvector *b,
1038 dvector *x,
1039 precond *pc,
1040 const REAL tol,
1041 const INT MaxIt,
1042 const SHORT restart,
1043 const SHORT StopType,
1044 const SHORT PrtLvl)
1045{
1046 const INT n = b->row;
1047 const INT min_iter = 0;
1048
1049 //--------------------------------------------//
1050 // Newly added parameters to monitor when //
1051 // to change the restart parameter //
1052 //--------------------------------------------//
1053 const REAL cr_max = 0.99; // = cos(8^o) (experimental)
1054 const REAL cr_min = 0.174; // = cos(80^o) (experimental)
1055
1056 // local variables
1057 INT iter = 0;
1058 int i, j, k; // must be signed! -zcs
1059
1060 REAL epsmac = SMALLREAL;
1061 REAL r_norm, b_norm, den_norm;
1062 REAL epsilon, gamma, t;
1063
1064 REAL *c = NULL, *s = NULL, *rs = NULL;
1065 REAL *norms = NULL, *r = NULL;
1066 REAL **p = NULL, **hh = NULL, **z=NULL;
1067 REAL *work = NULL;
1068
1069 REAL cr = 1.0; // convergence rate
1070 REAL r_norm_old = 0.0; // save residual norm of previous restart cycle
1071 INT d = 3; // reduction for restart parameter
1072 INT restart_max = restart; // upper bound for restart in each restart cycle
1073 INT restart_min = 3; // lower bound for restart in each restart cycle
1074
1075 INT Restart = restart; // real restart in some fixed restarted cycle
1076 INT Restart1 = Restart + 1;
1077 LONG worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
1078
1079 // Output some info for debuging
1080 if ( PrtLvl > PRINT_NONE ) printf("\nCalling VFGMRes solver (MatFree) ...\n");
1081
1082#if DEBUG_MODE > 0
1083 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1084 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1085#endif
1086
1087 /* allocate memory and setup temp work space */
1088 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1089
1090 /* check whether memory is enough for GMRES */
1091 while ( (work == NULL) && (Restart > 5) ) {
1092 Restart = Restart - 5;
1093 worksize = (Restart+4)*(Restart+n)+1-n+Restart*n;
1094 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1095 Restart1 = Restart + 1;
1096 }
1097
1098 if ( work == NULL ) {
1099 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
1100 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1101 }
1102
1103 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
1104 printf("### WARNING: vFGMRES restart number set to %d!\n", Restart);
1105 }
1106
1107 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1108 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1109 z = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1110 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1111
1112 r = work; rs = r + n; c = rs + Restart1; s = c + Restart;
1113 for (i = 0; i < Restart1; i ++) p[i] = s + Restart + i*n;
1114 for (i = 0; i < Restart1; i ++) hh[i] = p[Restart] + n + i*Restart;
1115 for (i = 0; i < Restart1; i ++) z[i] = hh[Restart] + Restart + i*n;
1116
1117 /* initialization */
1118 mf->fct(mf->data, x->val, p[0]);
1119 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, p[0]);
1120
1121 b_norm = fasp_blas_darray_norm2(n, b->val);
1122 r_norm = fasp_blas_darray_norm2(n, p[0]);
1123 norms[0] = r_norm;
1124
1125 if ( PrtLvl >= PRINT_SOME) {
1126 ITS_PUTNORM("right-hand side", b_norm);
1127 ITS_PUTNORM("residual", r_norm);
1128 }
1129
1130 if (b_norm > 0.0) den_norm = b_norm;
1131 else den_norm = r_norm;
1132
1133 epsilon = tol*den_norm;
1134
1135 /* outer iteration cycle */
1136 while (iter < MaxIt) {
1137 rs[0] = r_norm;
1138 r_norm_old = r_norm;
1139 if (r_norm == 0.0) {
1140 fasp_mem_free(work); work = NULL;
1141 fasp_mem_free(p); p = NULL;
1142 fasp_mem_free(hh); hh = NULL;
1143 fasp_mem_free(norms); norms = NULL;
1144 fasp_mem_free(z); z = NULL;
1145 return iter;
1146 }
1147
1148 //-----------------------------------//
1149 // adjust the restart parameter //
1150 //-----------------------------------//
1151
1152 if (cr > cr_max || iter == 0) {
1153 Restart = restart_max;
1154 }
1155 else if (cr < cr_min) {
1156 // Restart = Restart;
1157 }
1158 else {
1159 if ( Restart - d > restart_min ) Restart -= d;
1160 else Restart = restart_max;
1161 }
1162
1163 if (r_norm <= epsilon && iter >= min_iter) {
1164 mf->fct(mf->data, x->val, r);
1165 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1166 r_norm = fasp_blas_darray_norm2(n, r);
1167
1168 if (r_norm <= epsilon) {
1169 break;
1170 }
1171 else {
1172 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1173 }
1174 }
1175
1176 t = 1.0 / r_norm;
1177 fasp_blas_darray_ax(n, t, p[0]);
1178
1179 /* RESTART CYCLE (right-preconditioning) */
1180 i = 0;
1181 while (i < Restart && iter < MaxIt) {
1182
1183 i ++; iter ++;
1184
1185 /* apply preconditioner */
1186 if (pc == NULL) fasp_darray_cp(n, p[i-1], z[i-1]);
1187 else pc->fct(p[i-1], z[i-1], pc->data);
1188
1189 mf->fct(mf->data, z[i-1], p[i]);
1190
1191 /* modified Gram_Schmidt */
1192 for (j = 0; j < i; j ++) {
1193 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1194 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
1195 }
1196 t = fasp_blas_darray_norm2(n, p[i]);
1197 hh[i][i-1] = t;
1198 if (t != 0.0) {
1199 t = 1.0/t;
1200 fasp_blas_darray_ax(n, t, p[i]);
1201 }
1202
1203 for (j = 1; j < i; ++j) {
1204 t = hh[j-1][i-1];
1205 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1206 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1207 }
1208 t= hh[i][i-1]*hh[i][i-1];
1209 t+= hh[i-1][i-1]*hh[i-1][i-1];
1210 gamma = sqrt(t);
1211 if (gamma == 0.0) gamma = epsmac;
1212 c[i-1] = hh[i-1][i-1] / gamma;
1213 s[i-1] = hh[i][i-1] / gamma;
1214 rs[i] = -s[i-1]*rs[i-1];
1215 rs[i-1] = c[i-1]*rs[i-1];
1216 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1217
1218 r_norm = fabs(rs[i]);
1219 norms[iter] = r_norm;
1220
1221 if (b_norm > 0 ) {
1222 fasp_itinfo(PrtLvl, StopType, iter, norms[iter]/b_norm,
1223 norms[iter], norms[iter]/norms[iter-1]);
1224 }
1225 else {
1226 fasp_itinfo(PrtLvl, StopType, iter, norms[iter], norms[iter],
1227 norms[iter]/norms[iter-1]);
1228 }
1229
1230 /* Check: Exit restart cycle? */
1231 if (r_norm <= epsilon && iter >= min_iter) break;
1232
1233 } /* end of restart cycle */
1234
1235 /* now compute solution, first solve upper triangular system */
1236
1237 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1238 for (k = i-2; k >= 0; k --) {
1239 t = 0.0;
1240 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
1241
1242 t += rs[k];
1243 rs[k] = t / hh[k][k];
1244 }
1245
1246 fasp_darray_cp(n, z[i-1], r);
1247 fasp_blas_darray_ax(n, rs[i-1], r);
1248 for (j = i-2; j >= 0; j --) fasp_blas_darray_axpy(n, rs[j], z[j], r);
1249
1250 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1251
1252 if (r_norm <= epsilon && iter >= min_iter) {
1253 mf->fct(mf->data, x->val, r);
1254 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1255 r_norm = fasp_blas_darray_norm2(n, r);
1256
1257 if (r_norm <= epsilon) {
1258 break;
1259 }
1260 else {
1261 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1262 fasp_darray_cp(n, r, p[0]); i = 0;
1263 }
1264 } /* end of convergence check */
1265
1266
1267 /* compute residual vector and continue loop */
1268 for (j = i; j > 0; j--) {
1269 rs[j-1] = -s[j-1]*rs[j];
1270 rs[j] = c[j-1]*rs[j];
1271 }
1272
1273 if (i) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1274
1275 for (j = i-1 ; j > 0; j --) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1276
1277 if (i) {
1278 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1279 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1280 }
1281
1282 //-----------------------------------//
1283 // compute the convergence rate //
1284 //-----------------------------------//
1285 cr = r_norm / r_norm_old;
1286
1287 } /* end of iteration while loop */
1288
1289 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter,MaxIt,r_norm);
1290
1291 /*-------------------------------------------
1292 * Free some stuff
1293 *------------------------------------------*/
1294 fasp_mem_free(work); work = NULL;
1295 fasp_mem_free(p); p = NULL;
1296 fasp_mem_free(hh); hh = NULL;
1297 fasp_mem_free(norms); norms = NULL;
1298 fasp_mem_free(z); z = NULL;
1299
1300#if DEBUG_MODE > 0
1301 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1302#endif
1303
1304 if (iter>=MaxIt)
1305 return ERROR_SOLVER_MAXIT;
1306 else
1307 return iter;
1308}
1309
1310/*---------------------------------*/
1311/*-- End of File --*/
1312/*---------------------------------*/
void fasp_darray_cp(const INT n, const REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: AuxArray.c:164
void fasp_mem_free(void *mem)
Free up previous allocated memory body and set pointer to NULL.
Definition: AuxMemory.c:155
void * fasp_mem_calloc(const unsigned int size, const unsigned int type)
Allocate, initiate, and check memory.
Definition: AuxMemory.c:65
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.
Definition: AuxMessage.c:41
void fasp_chkerr(const SHORT status, const char *fctname)
Check error status and print out error messages before quit.
Definition: AuxMessage.c:213
REAL fasp_blas_darray_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies x and y.
Definition: BlaArray.c:741
void fasp_blas_darray_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: BlaArray.c:580
REAL fasp_blas_darray_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: BlaArray.c:657
void fasp_blas_darray_ax(const INT n, const REAL a, REAL *x)
x = a*x
Definition: BlaArray.c:43
void fasp_blas_darray_axpy(const INT n, const REAL a, const REAL *x, REAL *y)
y = a*x + y
Definition: BlaArray.c:93
void fasp_blas_dblc_mxv(const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvBLC.c:164
void fasp_blas_dblc_aAxpy(const REAL alpha, const dBLCmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvBLC.c:38
void fasp_blas_dbsr_aAxpy(const REAL alpha, const dBSRmat *A, const REAL *x, REAL *y)
Compute y := alpha*A*x + y.
Definition: BlaSpmvBSR.c:348
void fasp_blas_dbsr_mxv(const dBSRmat *A, const REAL *x, REAL *y)
Compute y := A*x.
Definition: BlaSpmvBSR.c:910
void fasp_blas_dcsr_mxv(const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvCSR.c:241
void fasp_blas_dcsr_aAxpy(const REAL alpha, const dCSRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvCSR.c:493
INT fasp_solver_dblc_pvfgmres(dBLCmat *A, 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 PFGMRES (right preconditioned) iterative method in which the restart parameter can...
Definition: KryPvfgmres.c:714
INT fasp_solver_dcsr_pvfgmres(dCSRmat *A, 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 PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: KryPvfgmres.c:67
INT fasp_solver_dbsr_pvfgmres(dBSRmat *A, 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 PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: KryPvfgmres.c:389
INT fasp_solver_pvfgmres(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 PFGMRES(right preconditioned) iterative method in which the restart parameter can ...
Definition: KryPvfgmres.c:1036
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define LONG
Definition: fasp.h:65
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:74
#define INT
Definition: fasp.h:64
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:131
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
#define PRINT_SOME
Definition: fasp_const.h:75
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
#define SMALLREAL
Definition: fasp_const.h:249
#define PRINT_MIN
Definition: fasp_const.h:74
Block REAL CSR matrix format.
Definition: fasp_block.h:74
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
Vector with n entries of REAL type.
Definition: fasp.h:346
REAL * val
actual vector entries
Definition: fasp.h:352
INT row
number of rows
Definition: fasp.h:349
Matrix-vector multiplication, replace the actual matrix.
Definition: fasp.h:1095
void * data
data for MxV, can be a Matrix or something else
Definition: fasp.h:1098
void(* fct)(const void *, const REAL *, REAL *)
action for MxV, void function pointer
Definition: fasp.h:1101
Preconditioner data and action.
Definition: fasp.h:1081
void * data
data for preconditioner, void pointer
Definition: fasp.h:1084
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1087