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