Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KrySPgmres.c
Go to the documentation of this file.
1
26#include <math.h>
27
28#include "fasp.h"
29#include "fasp_functs.h"
30
31/*---------------------------------*/
32/*-- Declare Private Functions --*/
33/*---------------------------------*/
34
35#include "KryUtil.inl"
36
37/*---------------------------------*/
38/*-- Public Functions --*/
39/*---------------------------------*/
40
66INT fasp_solver_dcsr_spgmres (const dCSRmat *A,
67 const dvector *b,
68 dvector *x,
69 precond *pc,
70 const REAL tol,
71 const INT MaxIt,
72 SHORT restart,
73 const SHORT StopType,
74 const SHORT PrtLvl)
75{
76 const INT n = b->row;
77 const INT MIN_ITER = 0;
78 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
79 const REAL epsmac = SMALLREAL;
80
81 // local variables
82 INT iter = 0;
83 INT restart1 = restart + 1;
84 int i, j, k; // must be signed! -zcs
85
86 REAL r_norm, r_normb, gamma, t;
87 REAL normr0 = BIGREAL, absres = BIGREAL;
88 REAL relres = BIGREAL, normu = BIGREAL;
89
90 INT iter_best = 0; // initial best known iteration
91 REAL absres_best = BIGREAL; // initial best known residual
92
93 // allocate temp memory (need about (restart+4)*n REAL numbers)
94 REAL *c = NULL, *s = NULL, *rs = NULL;
95 REAL *norms = NULL, *r = NULL, *w = NULL;
96 REAL *work = NULL, *x_best = NULL;
97 REAL **p = NULL, **hh = NULL;
98
99 // Output some info for debuging
100 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe GMRes solver (CSR) ...\n");
101
102#if DEBUG_MODE > 0
103 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
104 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
105#endif
106
107 /* allocate memory and setup temp work space */
108 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
109
110 /* check whether memory is enough for GMRES */
111 while ( (work == NULL) && (restart > 5 ) ) {
112 restart = restart - 5 ;
113 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
114 printf("### WARNING: GMRES restart number set to %d!\n", restart);
115 restart1 = restart + 1;
116 }
117
118 if ( work == NULL ) {
119 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
120 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
121 }
122
123 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
124 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
125 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
126
127 r = work; w = r + n; rs = w + n; c = rs + restart1;
128 x_best = c + restart; s = x_best + n;
129
130 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
131
132 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
133
134 // r = b-A*x
135 fasp_darray_cp(n, b->val, p[0]);
136 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
137
138 r_norm = fasp_blas_darray_norm2(n,p[0]);
139
140 // compute initial residuals
141 switch (StopType) {
142 case STOP_REL_RES:
143 normr0 = MAX(SMALLREAL,r_norm);
144 relres = r_norm/normr0;
145 break;
146 case STOP_REL_PRECRES:
147 if ( pc == NULL )
148 fasp_darray_cp(n, p[0], r);
149 else
150 pc->fct(p[0], r, pc->data);
151 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
152 normr0 = MAX(SMALLREAL,r_normb);
153 relres = r_normb/normr0;
154 break;
155 case STOP_MOD_REL_RES:
156 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
157 normr0 = r_norm;
158 relres = normr0/normu;
159 break;
160 default:
161 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
162 goto FINISHED;
163 }
164
165 // if initial residual is small, no need to iterate!
166 if ( relres < tol ) goto FINISHED;
167
168 // output iteration information if needed
169 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
170
171 // store initial residual
172 norms[0] = relres;
173
174 /* outer iteration cycle */
175 while ( iter < MaxIt ) {
176
177 rs[0] = r_norm;
178
179 t = 1.0 / r_norm;
180
181 fasp_blas_darray_ax(n, t, p[0]);
182
183 /* RESTART CYCLE (right-preconditioning) */
184 i = 0;
185 while ( i < restart && iter < MaxIt ) {
186
187 i++; iter++;
188
189 /* apply preconditioner */
190 if ( pc == NULL )
191 fasp_darray_cp(n, p[i-1], r);
192 else
193 pc->fct(p[i-1], r, pc->data);
194
195 fasp_blas_dcsr_mxv(A, r, p[i]);
196
197 /* modified Gram_Schmidt */
198 for ( j = 0; j < i; j++ ) {
199 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
200 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
201 }
202 t = fasp_blas_darray_norm2(n, p[i]);
203 hh[i][i-1] = t;
204 if (t != 0.0) {
205 t = 1.0/t;
206 fasp_blas_darray_ax(n, t, p[i]);
207 }
208
209 for (j = 1; j < i; ++j) {
210 t = hh[j-1][i-1];
211 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
212 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
213 }
214 t= hh[i][i-1]*hh[i][i-1];
215 t+= hh[i-1][i-1]*hh[i-1][i-1];
216
217 gamma = sqrt(t);
218 if (gamma == 0.0) gamma = epsmac;
219 c[i-1] = hh[i-1][i-1] / gamma;
220 s[i-1] = hh[i][i-1] / gamma;
221 rs[i] = -s[i-1]*rs[i-1];
222 rs[i-1] = c[i-1]*rs[i-1];
223 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
224
225 absres = r_norm = fabs(rs[i]);
226
227 relres = absres/normr0;
228
229 norms[iter] = relres;
230
231 // output iteration information if needed
232 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
233 norms[iter]/norms[iter-1]);
234
235 // should we exit restart cycle
236 if ( relres <= tol && iter >= MIN_ITER ) break;
237
238 } /* end of restart cycle */
239
240 /* compute solution, first solve upper triangular system */
241 rs[i-1] = rs[i-1] / hh[i-1][i-1];
242 for ( k = i-2; k >= 0; k-- ) {
243 t = 0.0;
244 for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
245
246 t += rs[k];
247 rs[k] = t / hh[k][k];
248 }
249
250 fasp_darray_cp(n, p[i-1], w);
251
252 fasp_blas_darray_ax(n, rs[i-1], w);
253
254 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
255
256 /* apply preconditioner */
257 if ( pc == NULL )
258 fasp_darray_cp(n, w, r);
259 else
260 pc->fct(w, r, pc->data);
261
262 fasp_blas_darray_axpy(n, 1.0, r, x->val);
263
264 // safety net check: save the best-so-far solution
265 if ( fasp_dvec_isnan(x) ) {
266 // If the solution is NAN, restore the best solution
267 absres = BIGREAL;
268 goto RESTORE_BESTSOL;
269 }
270
271 if ( absres < absres_best - maxdiff) {
272 absres_best = absres;
273 iter_best = iter;
274 fasp_darray_cp(n,x->val,x_best);
275 }
276
277 // Check: prevent false convergence
278 if ( relres <= tol && iter >= MIN_ITER ) {
279
280 fasp_darray_cp(n, b->val, r);
281 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
282
283 r_norm = fasp_blas_darray_norm2(n, r);
284
285 switch ( StopType ) {
286 case STOP_REL_RES:
287 absres = r_norm;
288 relres = absres/normr0;
289 break;
290 case STOP_REL_PRECRES:
291 if ( pc == NULL )
292 fasp_darray_cp(n, r, w);
293 else
294 pc->fct(r, w, pc->data);
295 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
296 relres = absres/normr0;
297 break;
298 case STOP_MOD_REL_RES:
299 absres = r_norm;
300 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
301 relres = absres/normu;
302 break;
303 }
304
305 norms[iter] = relres;
306
307 if ( relres <= tol ) {
308 break;
309 }
310 else {
311 // Need to restart
312 fasp_darray_cp(n, r, p[0]); i = 0;
313 }
314
315 } /* end of convergence check */
316
317 /* compute residual vector and continue loop */
318 for (j = i; j > 0; j--) {
319 rs[j-1] = -s[j-1]*rs[j];
320 rs[j] = c[j-1]*rs[j];
321 }
322
323 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
324
325 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
326
327 if ( i ) {
328 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
329 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
330 }
331
332 } /* end of main while loop */
333
334RESTORE_BESTSOL: // restore the best-so-far solution if necessary
335 if ( iter != iter_best ) {
336
337 // compute best residual
338 fasp_darray_cp(n,b->val,r);
339 fasp_blas_dcsr_aAxpy(-1.0,A,x_best,r);
340
341 switch ( StopType ) {
342 case STOP_REL_RES:
343 absres_best = fasp_blas_darray_norm2(n,r);
344 break;
345 case STOP_REL_PRECRES:
346 // z = B(r)
347 if ( pc != NULL )
348 pc->fct(r,w,pc->data); /* Apply preconditioner */
349 else
350 fasp_darray_cp(n,r,w); /* No preconditioner */
351 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
352 break;
353 case STOP_MOD_REL_RES:
354 absres_best = fasp_blas_darray_norm2(n,r);
355 break;
356 }
357
358 if ( absres > absres_best + maxdiff || isnan(absres) ) {
359 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
360 fasp_darray_cp(n,x_best,x->val);
361 relres = absres_best / normr0;
362 }
363 }
364
365FINISHED:
366 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
367
368 /*-------------------------------------------
369 * Clean up workspace
370 *------------------------------------------*/
371 fasp_mem_free(work); work = NULL;
372 fasp_mem_free(p); p = NULL;
373 fasp_mem_free(hh); hh = NULL;
374 fasp_mem_free(norms); norms = NULL;
375
376#if DEBUG_MODE > 0
377 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
378#endif
379
380 if ( iter >= MaxIt )
381 return ERROR_SOLVER_MAXIT;
382 else
383 return iter;
384}
385
409INT fasp_solver_dbsr_spgmres (const dBSRmat *A,
410 const dvector *b,
411 dvector *x,
412 precond *pc,
413 const REAL tol,
414 const INT MaxIt,
415 SHORT restart,
416 const SHORT StopType,
417 const SHORT PrtLvl)
418{
419 const INT n = b->row;
420 const INT MIN_ITER = 0;
421 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
422 const REAL epsmac = SMALLREAL;
423
424 // local variables
425 INT iter = 0;
426 INT restart1 = restart + 1;
427 int i, j, k; // must be signed! -zcs
428
429 REAL r_norm, r_normb, gamma, t;
430 REAL normr0 = BIGREAL, absres = BIGREAL;
431 REAL relres = BIGREAL, normu = BIGREAL;
432
433 INT iter_best = 0; // initial best known iteration
434 REAL absres_best = BIGREAL; // initial best known residual
435
436 // allocate temp memory (need about (restart+4)*n REAL numbers)
437 REAL *c = NULL, *s = NULL, *rs = NULL;
438 REAL *norms = NULL, *r = NULL, *w = NULL;
439 REAL *work = NULL, *x_best = NULL;
440 REAL **p = NULL, **hh = NULL;
441
442 // Output some info for debuging
443 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe GMRes solver (BSR) ...\n");
444
445#if DEBUG_MODE > 0
446 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
447 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
448#endif
449
450 /* allocate memory and setup temp work space */
451 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
452
453 /* check whether memory is enough for GMRES */
454 while ( (work == NULL) && (restart > 5 ) ) {
455 restart = restart - 5 ;
456 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
457 printf("### WARNING: GMRES restart number set to %d!\n", restart);
458 restart1 = restart + 1;
459 }
460
461 if ( work == NULL ) {
462 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
463 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
464 }
465
466 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
467 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
468 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
469
470 r = work; w = r + n; rs = w + n; c = rs + restart1;
471 x_best = c + restart; s = x_best + n;
472
473 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
474
475 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
476
477 // r = b-A*x
478 fasp_darray_cp(n, b->val, p[0]);
479 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
480
481 r_norm = fasp_blas_darray_norm2(n,p[0]);
482
483 // compute initial residuals
484 switch (StopType) {
485 case STOP_REL_RES:
486 normr0 = MAX(SMALLREAL,r_norm);
487 relres = r_norm/normr0;
488 break;
489 case STOP_REL_PRECRES:
490 if ( pc == NULL )
491 fasp_darray_cp(n, p[0], r);
492 else
493 pc->fct(p[0], r, pc->data);
494 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
495 normr0 = MAX(SMALLREAL,r_normb);
496 relres = r_normb/normr0;
497 break;
498 case STOP_MOD_REL_RES:
499 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
500 normr0 = r_norm;
501 relres = normr0/normu;
502 break;
503 default:
504 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
505 goto FINISHED;
506 }
507
508 // if initial residual is small, no need to iterate!
509 if ( relres < tol ) goto FINISHED;
510
511 // output iteration information if needed
512 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
513
514 // store initial residual
515 norms[0] = relres;
516
517 /* outer iteration cycle */
518 while ( iter < MaxIt ) {
519
520 rs[0] = r_norm;
521
522 t = 1.0 / r_norm;
523
524 fasp_blas_darray_ax(n, t, p[0]);
525
526 /* RESTART CYCLE (right-preconditioning) */
527 i = 0;
528 while ( i < restart && iter < MaxIt ) {
529
530 i++; iter++;
531
532 /* apply preconditioner */
533 if ( pc == NULL )
534 fasp_darray_cp(n, p[i-1], r);
535 else
536 pc->fct(p[i-1], r, pc->data);
537
538 fasp_blas_dbsr_mxv(A, r, p[i]);
539
540 /* modified Gram_Schmidt */
541 for ( j = 0; j < i; j++ ) {
542 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
543 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
544 }
545 t = fasp_blas_darray_norm2(n, p[i]);
546 hh[i][i-1] = t;
547 if (t != 0.0) {
548 t = 1.0/t;
549 fasp_blas_darray_ax(n, t, p[i]);
550 }
551
552 for (j = 1; j < i; ++j) {
553 t = hh[j-1][i-1];
554 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
555 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
556 }
557 t= hh[i][i-1]*hh[i][i-1];
558 t+= hh[i-1][i-1]*hh[i-1][i-1];
559
560 gamma = sqrt(t);
561 if (gamma == 0.0) gamma = epsmac;
562 c[i-1] = hh[i-1][i-1] / gamma;
563 s[i-1] = hh[i][i-1] / gamma;
564 rs[i] = -s[i-1]*rs[i-1];
565 rs[i-1] = c[i-1]*rs[i-1];
566 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
567
568 absres = r_norm = fabs(rs[i]);
569
570 relres = absres/normr0;
571
572 norms[iter] = relres;
573
574 // output iteration information if needed
575 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
576 norms[iter]/norms[iter-1]);
577
578 // should we exit restart cycle
579 if ( relres <= tol && iter >= MIN_ITER ) break;
580
581 } /* end of restart cycle */
582
583 /* compute solution, first solve upper triangular system */
584 rs[i-1] = rs[i-1] / hh[i-1][i-1];
585 for ( k = i-2; k >= 0; k-- ) {
586 t = 0.0;
587 for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
588
589 t += rs[k];
590 rs[k] = t / hh[k][k];
591 }
592
593 fasp_darray_cp(n, p[i-1], w);
594
595 fasp_blas_darray_ax(n, rs[i-1], w);
596
597 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
598
599 /* apply preconditioner */
600 if ( pc == NULL )
601 fasp_darray_cp(n, w, r);
602 else
603 pc->fct(w, r, pc->data);
604
605 fasp_blas_darray_axpy(n, 1.0, r, x->val);
606
607 // safety net check: save the best-so-far solution
608 if ( fasp_dvec_isnan(x) ) {
609 // If the solution is NAN, restore the best solution
610 absres = BIGREAL;
611 goto RESTORE_BESTSOL;
612 }
613
614 if ( absres < absres_best - maxdiff) {
615 absres_best = absres;
616 iter_best = iter;
617 fasp_darray_cp(n,x->val,x_best);
618 }
619
620 // Check: prevent false convergence
621 if ( relres <= tol && iter >= MIN_ITER ) {
622
623 fasp_darray_cp(n, b->val, r);
624 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
625
626 r_norm = fasp_blas_darray_norm2(n, r);
627
628 switch ( StopType ) {
629 case STOP_REL_RES:
630 absres = r_norm;
631 relres = absres/normr0;
632 break;
633 case STOP_REL_PRECRES:
634 if ( pc == NULL )
635 fasp_darray_cp(n, r, w);
636 else
637 pc->fct(r, w, pc->data);
638 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
639 relres = absres/normr0;
640 break;
641 case STOP_MOD_REL_RES:
642 absres = r_norm;
643 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
644 relres = absres/normu;
645 break;
646 }
647
648 norms[iter] = relres;
649
650 if ( relres <= tol ) {
651 break;
652 }
653 else {
654 // Need to restart
655 fasp_darray_cp(n, r, p[0]); i = 0;
656 }
657
658 } /* end of convergence check */
659
660 /* compute residual vector and continue loop */
661 for (j = i; j > 0; j--) {
662 rs[j-1] = -s[j-1]*rs[j];
663 rs[j] = c[j-1]*rs[j];
664 }
665
666 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
667
668 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
669
670 if ( i ) {
671 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
672 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
673 }
674
675 } /* end of main while loop */
676
677RESTORE_BESTSOL: // restore the best-so-far solution if necessary
678 if ( iter != iter_best ) {
679
680 // compute best residual
681 fasp_darray_cp(n,b->val,r);
682 fasp_blas_dbsr_aAxpy(-1.0,A,x_best,r);
683
684 switch ( StopType ) {
685 case STOP_REL_RES:
686 absres_best = fasp_blas_darray_norm2(n,r);
687 break;
688 case STOP_REL_PRECRES:
689 // z = B(r)
690 if ( pc != NULL )
691 pc->fct(r,w,pc->data); /* Apply preconditioner */
692 else
693 fasp_darray_cp(n,r,w); /* No preconditioner */
694 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
695 break;
696 case STOP_MOD_REL_RES:
697 absres_best = fasp_blas_darray_norm2(n,r);
698 break;
699 }
700
701 if ( absres > absres_best + maxdiff || isnan(absres) ) {
702 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
703 fasp_darray_cp(n,x_best,x->val);
704 relres = absres_best / normr0;
705 }
706 }
707
708FINISHED:
709 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
710
711 /*-------------------------------------------
712 * Clean up workspace
713 *------------------------------------------*/
714 fasp_mem_free(work); work = NULL;
715 fasp_mem_free(p); p = NULL;
716 fasp_mem_free(hh); hh = NULL;
717 fasp_mem_free(norms); norms = NULL;
718
719#if DEBUG_MODE > 0
720 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
721#endif
722
723 if ( iter >= MaxIt )
724 return ERROR_SOLVER_MAXIT;
725 else
726 return iter;
727}
728
752INT fasp_solver_dblc_spgmres (const dBLCmat *A,
753 const dvector *b,
754 dvector *x,
755 precond *pc,
756 const REAL tol,
757 const INT MaxIt,
758 SHORT restart,
759 const SHORT StopType,
760 const SHORT PrtLvl)
761{
762 const INT n = b->row;
763 const INT MIN_ITER = 0;
764 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
765 const REAL epsmac = SMALLREAL;
766
767 // local variables
768 INT iter = 0;
769 INT restart1 = restart + 1;
770 int i, j, k; // must be signed! -zcs
771
772 REAL r_norm, r_normb, gamma, t;
773 REAL normr0 = BIGREAL, absres = BIGREAL;
774 REAL relres = BIGREAL, normu = BIGREAL;
775
776 INT iter_best = 0; // initial best known iteration
777 REAL absres_best = BIGREAL; // initial best known residual
778
779 // allocate temp memory (need about (restart+4)*n REAL numbers)
780 REAL *c = NULL, *s = NULL, *rs = NULL;
781 REAL *norms = NULL, *r = NULL, *w = NULL;
782 REAL *work = NULL, *x_best = NULL;
783 REAL **p = NULL, **hh = NULL;
784
785 // Output some info for debuging
786 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe GMRes solver (BLC) ...\n");
787
788#if DEBUG_MODE > 0
789 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
790 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
791#endif
792
793 /* allocate memory and setup temp work space */
794 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
795
796 /* check whether memory is enough for GMRES */
797 while ( (work == NULL) && (restart > 5 ) ) {
798 restart = restart - 5 ;
799 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
800 printf("### WARNING: GMRES restart number set to %d!\n", restart);
801 restart1 = restart + 1;
802 }
803
804 if ( work == NULL ) {
805 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
806 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
807 }
808
809 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
810 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
811 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
812
813 r = work; w = r + n; rs = w + n; c = rs + restart1;
814 x_best = c + restart; s = x_best + n;
815
816 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
817
818 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
819
820 // r = b-A*x
821 fasp_darray_cp(n, b->val, p[0]);
822 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
823
824 r_norm = fasp_blas_darray_norm2(n,p[0]);
825
826 // compute initial residuals
827 switch (StopType) {
828 case STOP_REL_RES:
829 normr0 = MAX(SMALLREAL,r_norm);
830 relres = r_norm/normr0;
831 break;
832 case STOP_REL_PRECRES:
833 if ( pc == NULL )
834 fasp_darray_cp(n, p[0], r);
835 else
836 pc->fct(p[0], r, pc->data);
837 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
838 normr0 = MAX(SMALLREAL,r_normb);
839 relres = r_normb/normr0;
840 break;
841 case STOP_MOD_REL_RES:
842 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
843 normr0 = r_norm;
844 relres = normr0/normu;
845 break;
846 default:
847 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
848 goto FINISHED;
849 }
850
851 // if initial residual is small, no need to iterate!
852 if ( relres < tol ) goto FINISHED;
853
854 // output iteration information if needed
855 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
856
857 // store initial residual
858 norms[0] = relres;
859
860 /* outer iteration cycle */
861 while ( iter < MaxIt ) {
862
863 rs[0] = r_norm;
864
865 t = 1.0 / r_norm;
866
867 fasp_blas_darray_ax(n, t, p[0]);
868
869 /* RESTART CYCLE (right-preconditioning) */
870 i = 0;
871 while ( i < restart && iter < MaxIt ) {
872
873 i++; iter++;
874
875 /* apply preconditioner */
876 if ( pc == NULL )
877 fasp_darray_cp(n, p[i-1], r);
878 else
879 pc->fct(p[i-1], r, pc->data);
880
881 fasp_blas_dblc_mxv(A, r, p[i]);
882
883 /* modified Gram_Schmidt */
884 for ( j = 0; j < i; j++ ) {
885 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
886 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
887 }
888 t = fasp_blas_darray_norm2(n, p[i]);
889 hh[i][i-1] = t;
890 if (t != 0.0) {
891 t = 1.0/t;
892 fasp_blas_darray_ax(n, t, p[i]);
893 }
894
895 for (j = 1; j < i; ++j) {
896 t = hh[j-1][i-1];
897 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
898 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
899 }
900 t= hh[i][i-1]*hh[i][i-1];
901 t+= hh[i-1][i-1]*hh[i-1][i-1];
902
903 gamma = sqrt(t);
904 if (gamma == 0.0) gamma = epsmac;
905 c[i-1] = hh[i-1][i-1] / gamma;
906 s[i-1] = hh[i][i-1] / gamma;
907 rs[i] = -s[i-1]*rs[i-1];
908 rs[i-1] = c[i-1]*rs[i-1];
909 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
910
911 absres = r_norm = fabs(rs[i]);
912
913 relres = absres/normr0;
914
915 norms[iter] = relres;
916
917 // output iteration information if needed
918 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
919 norms[iter]/norms[iter-1]);
920
921 // should we exit restart cycle
922 if ( relres <= tol && iter >= MIN_ITER ) break;
923
924 } /* end of restart cycle */
925
926 /* compute solution, first solve upper triangular system */
927 rs[i-1] = rs[i-1] / hh[i-1][i-1];
928 for ( k = i-2; k >= 0; k-- ) {
929 t = 0.0;
930 for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
931
932 t += rs[k];
933 rs[k] = t / hh[k][k];
934 }
935
936 fasp_darray_cp(n, p[i-1], w);
937
938 fasp_blas_darray_ax(n, rs[i-1], w);
939
940 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
941
942 /* apply preconditioner */
943 if ( pc == NULL )
944 fasp_darray_cp(n, w, r);
945 else
946 pc->fct(w, r, pc->data);
947
948 fasp_blas_darray_axpy(n, 1.0, r, x->val);
949
950 // safety net check: save the best-so-far solution
951 if ( fasp_dvec_isnan(x) ) {
952 // If the solution is NAN, restore the best solution
953 absres = BIGREAL;
954 goto RESTORE_BESTSOL;
955 }
956
957 if ( absres < absres_best - maxdiff) {
958 absres_best = absres;
959 iter_best = iter;
960 fasp_darray_cp(n,x->val,x_best);
961 }
962
963 // Check: prevent false convergence
964 if ( relres <= tol && iter >= MIN_ITER ) {
965
966 fasp_darray_cp(n, b->val, r);
967 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
968
969 r_norm = fasp_blas_darray_norm2(n, r);
970
971 switch ( StopType ) {
972 case STOP_REL_RES:
973 absres = r_norm;
974 relres = absres/normr0;
975 break;
976 case STOP_REL_PRECRES:
977 if ( pc == NULL )
978 fasp_darray_cp(n, r, w);
979 else
980 pc->fct(r, w, pc->data);
981 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
982 relres = absres/normr0;
983 break;
984 case STOP_MOD_REL_RES:
985 absres = r_norm;
986 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
987 relres = absres/normu;
988 break;
989 }
990
991 norms[iter] = relres;
992
993 if ( relres <= tol ) {
994 break;
995 }
996 else {
997 // Need to restart
998 fasp_darray_cp(n, r, p[0]); i = 0;
999 }
1000
1001 } /* end of convergence check */
1002
1003 /* compute residual vector and continue loop */
1004 for (j = i; j > 0; j--) {
1005 rs[j-1] = -s[j-1]*rs[j];
1006 rs[j] = c[j-1]*rs[j];
1007 }
1008
1009 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1010
1011 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1012
1013 if ( i ) {
1014 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1015 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1016 }
1017
1018 } /* end of main while loop */
1019
1020RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1021 if ( iter != iter_best ) {
1022
1023 // compute best residual
1024 fasp_darray_cp(n,b->val,r);
1025 fasp_blas_dblc_aAxpy(-1.0,A,x_best,r);
1026
1027 switch ( StopType ) {
1028 case STOP_REL_RES:
1029 absres_best = fasp_blas_darray_norm2(n,r);
1030 break;
1031 case STOP_REL_PRECRES:
1032 // z = B(r)
1033 if ( pc != NULL )
1034 pc->fct(r,w,pc->data); /* Apply preconditioner */
1035 else
1036 fasp_darray_cp(n,r,w); /* No preconditioner */
1037 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
1038 break;
1039 case STOP_MOD_REL_RES:
1040 absres_best = fasp_blas_darray_norm2(n,r);
1041 break;
1042 }
1043
1044 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1045 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1046 fasp_darray_cp(n,x_best,x->val);
1047 relres = absres_best / normr0;
1048 }
1049 }
1050
1051FINISHED:
1052 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1053
1054 /*-------------------------------------------
1055 * Clean up workspace
1056 *------------------------------------------*/
1057 fasp_mem_free(work); work = NULL;
1058 fasp_mem_free(p); p = NULL;
1059 fasp_mem_free(hh); hh = NULL;
1060 fasp_mem_free(norms); norms = NULL;
1061
1062#if DEBUG_MODE > 0
1063 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1064#endif
1065
1066 if ( iter >= MaxIt )
1067 return ERROR_SOLVER_MAXIT;
1068 else
1069 return iter;
1070}
1071
1095INT fasp_solver_dstr_spgmres (const dSTRmat *A,
1096 const dvector *b,
1097 dvector *x,
1098 precond *pc,
1099 const REAL tol,
1100 const INT MaxIt,
1101 SHORT restart,
1102 const SHORT StopType,
1103 const SHORT PrtLvl)
1104{
1105 const INT n = b->row;
1106 const INT MIN_ITER = 0;
1107 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1108 const REAL epsmac = SMALLREAL;
1109
1110 // local variables
1111 INT iter = 0;
1112 INT restart1 = restart + 1;
1113 int i, j, k; // must be signed! -zcs
1114
1115 REAL r_norm, r_normb, gamma, t;
1116 REAL normr0 = BIGREAL, absres = BIGREAL;
1117 REAL relres = BIGREAL, normu = BIGREAL;
1118
1119 INT iter_best = 0; // initial best known iteration
1120 REAL absres_best = BIGREAL; // initial best known residual
1121
1122 // allocate temp memory (need about (restart+4)*n REAL numbers)
1123 REAL *c = NULL, *s = NULL, *rs = NULL;
1124 REAL *norms = NULL, *r = NULL, *w = NULL;
1125 REAL *work = NULL, *x_best = NULL;
1126 REAL **p = NULL, **hh = NULL;
1127
1128 // Output some info for debuging
1129 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe GMRes solver (STR) ...\n");
1130
1131#if DEBUG_MODE > 0
1132 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1133 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1134#endif
1135
1136 /* allocate memory and setup temp work space */
1137 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1138
1139 /* check whether memory is enough for GMRES */
1140 while ( (work == NULL) && (restart > 5 ) ) {
1141 restart = restart - 5 ;
1142 work = (REAL *) fasp_mem_calloc((restart+4)*(restart+n)+1, sizeof(REAL));
1143 printf("### WARNING: GMRES restart number set to %d!\n", restart);
1144 restart1 = restart + 1;
1145 }
1146
1147 if ( work == NULL ) {
1148 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__);
1149 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1150 }
1151
1152 p = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1153 hh = (REAL **)fasp_mem_calloc(restart1, sizeof(REAL *));
1154 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1155
1156 r = work; w = r + n; rs = w + n; c = rs + restart1;
1157 x_best = c + restart; s = x_best + n;
1158
1159 for ( i = 0; i < restart1; i++ ) p[i] = s + restart + i*n;
1160
1161 for ( i = 0; i < restart1; i++ ) hh[i] = p[restart] + n + i*restart;
1162
1163 // r = b-A*x
1164 fasp_darray_cp(n, b->val, p[0]);
1165 fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1166
1167 r_norm = fasp_blas_darray_norm2(n,p[0]);
1168
1169 // compute initial residuals
1170 switch (StopType) {
1171 case STOP_REL_RES:
1172 normr0 = MAX(SMALLREAL,r_norm);
1173 relres = r_norm/normr0;
1174 break;
1175 case STOP_REL_PRECRES:
1176 if ( pc == NULL )
1177 fasp_darray_cp(n, p[0], r);
1178 else
1179 pc->fct(p[0], r, pc->data);
1180 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
1181 normr0 = MAX(SMALLREAL,r_normb);
1182 relres = r_normb/normr0;
1183 break;
1184 case STOP_MOD_REL_RES:
1185 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1186 normr0 = r_norm;
1187 relres = normr0/normu;
1188 break;
1189 default:
1190 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1191 goto FINISHED;
1192 }
1193
1194 // if initial residual is small, no need to iterate!
1195 if ( relres < tol ) goto FINISHED;
1196
1197 // output iteration information if needed
1198 fasp_itinfo(PrtLvl,StopType,0,relres,normr0,0.0);
1199
1200 // store initial residual
1201 norms[0] = relres;
1202
1203 /* outer iteration cycle */
1204 while ( iter < MaxIt ) {
1205
1206 rs[0] = r_norm;
1207
1208 t = 1.0 / r_norm;
1209
1210 fasp_blas_darray_ax(n, t, p[0]);
1211
1212 /* RESTART CYCLE (right-preconditioning) */
1213 i = 0;
1214 while ( i < restart && iter < MaxIt ) {
1215
1216 i++; iter++;
1217
1218 /* apply preconditioner */
1219 if ( pc == NULL )
1220 fasp_darray_cp(n, p[i-1], r);
1221 else
1222 pc->fct(p[i-1], r, pc->data);
1223
1224 fasp_blas_dstr_mxv(A, r, p[i]);
1225
1226 /* modified Gram_Schmidt */
1227 for ( j = 0; j < i; j++ ) {
1228 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1229 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
1230 }
1231 t = fasp_blas_darray_norm2(n, p[i]);
1232 hh[i][i-1] = t;
1233 if (t != 0.0) {
1234 t = 1.0/t;
1235 fasp_blas_darray_ax(n, t, p[i]);
1236 }
1237
1238 for (j = 1; j < i; ++j) {
1239 t = hh[j-1][i-1];
1240 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1241 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1242 }
1243 t= hh[i][i-1]*hh[i][i-1];
1244 t+= hh[i-1][i-1]*hh[i-1][i-1];
1245
1246 gamma = sqrt(t);
1247 if (gamma == 0.0) gamma = epsmac;
1248 c[i-1] = hh[i-1][i-1] / gamma;
1249 s[i-1] = hh[i][i-1] / gamma;
1250 rs[i] = -s[i-1]*rs[i-1];
1251 rs[i-1] = c[i-1]*rs[i-1];
1252 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1253
1254 absres = r_norm = fabs(rs[i]);
1255
1256 relres = absres/normr0;
1257
1258 norms[iter] = relres;
1259
1260 // output iteration information if needed
1261 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1262 norms[iter]/norms[iter-1]);
1263
1264 // should we exit restart cycle
1265 if ( relres <= tol && iter >= MIN_ITER ) break;
1266
1267 } /* end of restart cycle */
1268
1269 /* compute solution, first solve upper triangular system */
1270 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1271 for ( k = i-2; k >= 0; k-- ) {
1272 t = 0.0;
1273 for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
1274
1275 t += rs[k];
1276 rs[k] = t / hh[k][k];
1277 }
1278
1279 fasp_darray_cp(n, p[i-1], w);
1280
1281 fasp_blas_darray_ax(n, rs[i-1], w);
1282
1283 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1284
1285 /* apply preconditioner */
1286 if ( pc == NULL )
1287 fasp_darray_cp(n, w, r);
1288 else
1289 pc->fct(w, r, pc->data);
1290
1291 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1292
1293 // safety net check: save the best-so-far solution
1294 if ( fasp_dvec_isnan(x) ) {
1295 // If the solution is NAN, restore the best solution
1296 absres = BIGREAL;
1297 goto RESTORE_BESTSOL;
1298 }
1299
1300 if ( absres < absres_best - maxdiff) {
1301 absres_best = absres;
1302 iter_best = iter;
1303 fasp_darray_cp(n,x->val,x_best);
1304 }
1305
1306 // Check: prevent false convergence
1307 if ( relres <= tol && iter >= MIN_ITER ) {
1308
1309 fasp_darray_cp(n, b->val, r);
1310 fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1311
1312 r_norm = fasp_blas_darray_norm2(n, r);
1313
1314 switch ( StopType ) {
1315 case STOP_REL_RES:
1316 absres = r_norm;
1317 relres = absres/normr0;
1318 break;
1319 case STOP_REL_PRECRES:
1320 if ( pc == NULL )
1321 fasp_darray_cp(n, r, w);
1322 else
1323 pc->fct(r, w, pc->data);
1324 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
1325 relres = absres/normr0;
1326 break;
1327 case STOP_MOD_REL_RES:
1328 absres = r_norm;
1329 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1330 relres = absres/normu;
1331 break;
1332 }
1333
1334 norms[iter] = relres;
1335
1336 if ( relres <= tol ) {
1337 break;
1338 }
1339 else {
1340 // Need to restart
1341 fasp_darray_cp(n, r, p[0]); i = 0;
1342 }
1343
1344 } /* end of convergence check */
1345
1346 /* compute residual vector and continue loop */
1347 for (j = i; j > 0; j--) {
1348 rs[j-1] = -s[j-1]*rs[j];
1349 rs[j] = c[j-1]*rs[j];
1350 }
1351
1352 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1353
1354 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1355
1356 if ( i ) {
1357 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1358 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1359 }
1360
1361 } /* end of main while loop */
1362
1363RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1364 if ( iter != iter_best ) {
1365
1366 // compute best residual
1367 fasp_darray_cp(n,b->val,r);
1368 fasp_blas_dstr_aAxpy(-1.0,A,x_best,r);
1369
1370 switch ( StopType ) {
1371 case STOP_REL_RES:
1372 absres_best = fasp_blas_darray_norm2(n,r);
1373 break;
1374 case STOP_REL_PRECRES:
1375 // z = B(r)
1376 if ( pc != NULL )
1377 pc->fct(r,w,pc->data); /* Apply preconditioner */
1378 else
1379 fasp_darray_cp(n,r,w); /* No preconditioner */
1380 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(n,w,r)));
1381 break;
1382 case STOP_MOD_REL_RES:
1383 absres_best = fasp_blas_darray_norm2(n,r);
1384 break;
1385 }
1386
1387 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1388 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1389 fasp_darray_cp(n,x_best,x->val);
1390 relres = absres_best / normr0;
1391 }
1392 }
1393
1394FINISHED:
1395 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1396
1397 /*-------------------------------------------
1398 * Clean up workspace
1399 *------------------------------------------*/
1400 fasp_mem_free(work); work = NULL;
1401 fasp_mem_free(p); p = NULL;
1402 fasp_mem_free(hh); hh = NULL;
1403 fasp_mem_free(norms); norms = NULL;
1404
1405#if DEBUG_MODE > 0
1406 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1407#endif
1408
1409 if ( iter >= MaxIt )
1410 return ERROR_SOLVER_MAXIT;
1411 else
1412 return iter;
1413}
1414
1415/*---------------------------------*/
1416/*-- End of File --*/
1417/*---------------------------------*/
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_dvec_isnan
SHORT fasp_dvec_isnan(const dvector *u)
Check a dvector whether there is NAN.
Definition: AuxVector.c:39
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_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_blas_dstr_aAxpy
void fasp_blas_dstr_aAxpy(const REAL alpha, const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: BlaSpmvSTR.c:61
fasp_blas_dstr_mxv
void fasp_blas_dstr_mxv(const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvSTR.c:131
fasp_solver_dbsr_spgmres
INT fasp_solver_dbsr_spgmres(const dBSRmat *A, const dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b with safe-guard.
Definition: KrySPgmres.c:409
fasp_solver_dcsr_spgmres
INT fasp_solver_dcsr_spgmres(const dCSRmat *A, const dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b with safe-guard.
Definition: KrySPgmres.c:66
fasp_solver_dstr_spgmres
INT fasp_solver_dstr_spgmres(const dSTRmat *A, const dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b with safe-guard.
Definition: KrySPgmres.c:1095
fasp_solver_dblc_spgmres
INT fasp_solver_dblc_spgmres(const dBLCmat *A, const dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b with safe-guard.
Definition: KrySPgmres.c:752
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
ABS
#define ABS(a)
Definition: fasp.h:84
MAX
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:82
INT
#define INT
Definition: fasp.h:72
STAG_RATIO
#define STAG_RATIO
Definition: fasp_const.h:265
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
BIGREAL
#define BIGREAL
Some global constants.
Definition: fasp_const.h:255
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
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
dSTRmat
Structure matrix of REAL type.
Definition: fasp.h:316
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
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