Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPgmres.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 // local variables
81 INT iter = 0;
82 int i, j, k; // must be signed! -zcs
83
84 REAL r_norm, r_normb, gamma, t;
85 REAL absres0 = BIGREAL, absres = BIGREAL;
86 REAL relres = BIGREAL, normu = BIGREAL;
87
88 // allocate temp memory (need about (restart+4)*n REAL numbers)
89 REAL *c = NULL, *s = NULL, *rs = NULL;
90 REAL *norms = NULL, *r = NULL, *w = NULL;
91 REAL *work = NULL;
92 REAL **p = NULL, **hh = NULL;
93
94 INT Restart = MIN(restart, MaxIt);
95 INT Restart1 = Restart + 1;
96 LONG worksize = (Restart+4)*(Restart+n)+1-n;
97
98 /* allocate memory and setup temp work space */
99 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
100
101 // Output some info for debugging
102 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GMRes solver (CSR) ...\n");
103
104#if DEBUG_MODE > 0
105 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
106 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
107#endif
108
109 /* check whether memory is enough for GMRES */
110 while ( (work == NULL) && (Restart > 5) ) {
111 Restart = Restart - 5;
112 Restart1 = Restart + 1;
113 worksize = (Restart+4)*(Restart+n)+1-n;
114 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
115 }
116
117 if ( work == NULL ) {
118 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
119 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
120 }
121
122 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
123 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
124 }
125
126 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
127 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
128 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
129
130 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
131
132 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
133
134 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
135
136 // compute initial residual: r = b-A*x
137 fasp_darray_cp(n, b->val, p[0]);
138 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, p[0]);
139 r_norm = fasp_blas_darray_norm2(n,p[0]);
140
141 // compute stopping criteria
142 switch (StopType) {
143 case STOP_REL_RES:
144 absres0 = MAX(SMALLREAL,r_norm);
145 relres = r_norm/absres0;
146 break;
147 case STOP_REL_PRECRES:
148 if ( pc == NULL )
149 fasp_darray_cp(n, p[0], r);
150 else
151 pc->fct(p[0], r, pc->data);
152 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
153 absres0 = MAX(SMALLREAL,r_normb);
154 relres = r_normb/absres0;
155 break;
156 case STOP_MOD_REL_RES:
158 absres0 = r_norm;
159 relres = absres0/normu;
160 break;
161 default:
162 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
163 goto FINISHED;
164 }
165
166 // if initial residual is small, no need to iterate!
167 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
168
169 // output iteration information if needed
170 fasp_itinfo(PrtLvl,StopType,0,relres,absres0,0.0);
171
172 // store initial residual
173 norms[0] = relres;
174
175 /* GMRES(M) outer iteration */
176 while ( iter < MaxIt && relres > tol ) {
177
178 rs[0] = r_norm;
179
180 t = 1.0 / r_norm;
181
182 fasp_blas_darray_ax(n, t, p[0]);
183
184 /* RESTART CYCLE (right-preconditioning) */
185 i = 0;
186 while ( i < Restart && iter < MaxIt ) {
187
188 i++; iter++;
189
190 /* apply preconditioner */
191 if ( pc == NULL )
192 fasp_darray_cp(n, p[i-1], r);
193 else
194 pc->fct(p[i-1], r, pc->data);
195
196 fasp_blas_dcsr_mxv(A, r, p[i]);
197
198 /* Modified Gram_Schmidt orthogonalization */
199 for ( j = 0; j < i; j++ ) {
200 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
201 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
202 }
203 t = fasp_blas_darray_norm2(n, p[i]);
204 hh[i][i-1] = t;
205
206 if ( ABS(t) > SMALLREAL ) { // If t=0, we get solution subspace
207 t = 1.0/t;
208 fasp_blas_darray_ax(n, t, p[i]);
209 }
210
211 for ( j = 1; j < i; ++j ) {
212 t = hh[j-1][i-1];
213 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
214 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
215 }
216 t = hh[i][i-1]*hh[i][i-1];
217 t += hh[i-1][i-1]*hh[i-1][i-1];
218
219 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
220 c[i-1] = hh[i-1][i-1] / gamma;
221 s[i-1] = hh[i][i-1] / gamma;
222 rs[i] = -s[i-1]*rs[i-1];
223 rs[i-1] = c[i-1]*rs[i-1];
224 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
225
226 absres = r_norm = fabs(rs[i]);
227
228 relres = absres/absres0;
229
230 norms[iter] = relres;
231
232 // output iteration information if needed
233 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
234 norms[iter]/norms[iter-1]);
235
236 // exit restart cycle if reaches tolerance
237 if ( relres < tol && iter >= MIN_ITER ) break;
238
239 } /* end of restart cycle */
240
241 /* compute solution, first solve upper triangular system */
242 rs[i-1] = rs[i-1] / hh[i-1][i-1];
243 for ( k = i-2; k >= 0; k-- ) {
244 t = 0.0;
245 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
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 // Check: prevent false convergence
265 if ( relres < tol && iter >= MIN_ITER ) {
266
267 REAL computed_relres = relres;
268
269 // compute residual
270 fasp_darray_cp(n, b->val, r);
271 fasp_blas_dcsr_aAxpy(-1.0, A, x->val, r);
272 r_norm = fasp_blas_darray_norm2(n, r);
273
274 switch ( StopType ) {
275 case STOP_REL_RES:
276 absres = r_norm;
277 relres = absres/absres0;
278 break;
279 case STOP_REL_PRECRES:
280 if ( pc == NULL )
281 fasp_darray_cp(n, r, w);
282 else
283 pc->fct(r, w, pc->data);
284 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
285 relres = absres/absres0;
286 break;
287 case STOP_MOD_REL_RES:
288 absres = r_norm;
290 relres = absres/normu;
291 break;
292 }
293
294 norms[iter] = relres;
295
296 if ( relres < tol ) {
297 break;
298 }
299 else { // Need to restart
300 fasp_darray_cp(n, r, p[0]); i = 0;
301 }
302
303 if ( PrtLvl >= PRINT_MORE ) {
304 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
305 }
306
307 } /* end of convergence check */
308
309 /* compute residual vector and continue loop */
310 for ( j = i; j > 0; j-- ) {
311 rs[j-1] = -s[j-1]*rs[j];
312 rs[j] = c[j-1]*rs[j];
313 }
314
315 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
316
317 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
318
319 if ( i ) {
320 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
321 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
322 }
323
324 } /* end of main while loop */
325FINISHED:
326 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
327
328 /*-------------------------------------------
329 * Clean up workspace
330 *------------------------------------------*/
331 fasp_mem_free(work); work = NULL;
332 fasp_mem_free(p); p = NULL;
333 fasp_mem_free(hh); hh = NULL;
334 fasp_mem_free(norms); norms = NULL;
335
336#if DEBUG_MODE > 0
337 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
338#endif
339
340 if ( iter >= MaxIt )
341 return ERROR_SOLVER_MAXIT;
342 else
343 return iter;
344}
345
371 dvector *b,
372 dvector *x,
373 precond *pc,
374 const REAL tol,
375 const INT MaxIt,
376 const SHORT restart,
377 const SHORT StopType,
378 const SHORT PrtLvl)
379{
380 const INT n = b->row;
381 const INT MIN_ITER = 0;
382
383 // local variables
384 INT iter = 0;
385 int i, j, k; // must be signed! -zcs
386
387 REAL r_norm, r_normb, gamma, t;
388 REAL absres0 = BIGREAL, absres = BIGREAL;
389 REAL relres = BIGREAL, normu = BIGREAL;
390
391 // allocate temp memory (need about (restart+4)*n REAL numbers)
392 REAL *c = NULL, *s = NULL, *rs = NULL;
393 REAL *norms = NULL, *r = NULL, *w = NULL;
394 REAL *work = NULL;
395 REAL **p = NULL, **hh = NULL;
396
397 INT Restart = MIN(restart, MaxIt);
398 INT Restart1 = Restart + 1;
399 LONG worksize = (Restart+4)*(Restart+n)+1-n;
400
401 /* allocate memory and setup temp work space */
402 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
403
404 // Output some info for debugging
405 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GMRes solver (BSR) ...\n");
406
407#if DEBUG_MODE > 0
408 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
409 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
410#endif
411
412 /* check whether memory is enough for GMRES */
413 while ( (work == NULL) && (Restart > 5) ) {
414 Restart = Restart - 5;
415 Restart1 = Restart + 1;
416 worksize = (Restart+4)*(Restart+n)+1-n;
417 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
418 }
419
420 if ( work == NULL ) {
421 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
422 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
423 }
424
425 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
426 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
427 }
428
429 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
430 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
431 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
432
433 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
434
435 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
436
437 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
438
439 // compute initial residual: r = b-A*x
440 fasp_darray_cp(n, b->val, p[0]);
441 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, p[0]);
442 r_norm = fasp_blas_darray_norm2(n,p[0]);
443
444 // compute stopping criteria
445 switch (StopType) {
446 case STOP_REL_RES:
447 absres0 = MAX(SMALLREAL,r_norm);
448 relres = r_norm/absres0;
449 break;
450 case STOP_REL_PRECRES:
451 if ( pc == NULL )
452 fasp_darray_cp(n, p[0], r);
453 else
454 pc->fct(p[0], r, pc->data);
455 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
456 absres0 = MAX(SMALLREAL,r_normb);
457 relres = r_normb/absres0;
458 break;
459 case STOP_MOD_REL_RES:
461 absres0 = r_norm;
462 relres = absres0/normu;
463 break;
464 default:
465 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
466 goto FINISHED;
467 }
468
469 // if initial residual is small, no need to iterate!
470 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
471
472 // output iteration information if needed
473 fasp_itinfo(PrtLvl,StopType,0,relres,absres0,0.0);
474
475 // store initial residual
476 norms[0] = relres;
477
478 /* GMRES(M) outer iteration */
479 while ( iter < MaxIt && relres > tol ) {
480
481 rs[0] = r_norm;
482
483 t = 1.0 / r_norm;
484
485 fasp_blas_darray_ax(n, t, p[0]);
486
487 /* RESTART CYCLE (right-preconditioning) */
488 i = 0;
489 while ( i < Restart && iter < MaxIt ) {
490
491 i++; iter++;
492
493 /* apply preconditioner */
494 if ( pc == NULL )
495 fasp_darray_cp(n, p[i-1], r);
496 else
497 pc->fct(p[i-1], r, pc->data);
498
499 fasp_blas_dbsr_mxv(A, r, p[i]);
500
501 /* Modified Gram_Schmidt orthogonalization */
502 for ( j = 0; j < i; j++ ) {
503 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
504 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
505 }
506 t = fasp_blas_darray_norm2(n, p[i]);
507 hh[i][i-1] = t;
508
509 if ( ABS(t) > SMALLREAL ) { // If t=0, we get solution subspace
510 t = 1.0/t;
511 fasp_blas_darray_ax(n, t, p[i]);
512 }
513
514 for ( j = 1; j < i; ++j ) {
515 t = hh[j-1][i-1];
516 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
517 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
518 }
519 t = hh[i][i-1]*hh[i][i-1];
520 t += hh[i-1][i-1]*hh[i-1][i-1];
521
522 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
523 c[i-1] = hh[i-1][i-1] / gamma;
524 s[i-1] = hh[i][i-1] / gamma;
525 rs[i] = -s[i-1]*rs[i-1];
526 rs[i-1] = c[i-1]*rs[i-1];
527 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
528
529 absres = r_norm = fabs(rs[i]);
530
531 relres = absres/absres0;
532
533 norms[iter] = relres;
534
535 // output iteration information if needed
536 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
537 norms[iter]/norms[iter-1]);
538
539 // exit restart cycle if reaches tolerance
540 if ( relres < tol && iter >= MIN_ITER ) break;
541
542 } /* end of restart cycle */
543
544 /* compute solution, first solve upper triangular system */
545 rs[i-1] = rs[i-1] / hh[i-1][i-1];
546 for ( k = i-2; k >= 0; k-- ) {
547 t = 0.0;
548 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
549 t += rs[k];
550 rs[k] = t / hh[k][k];
551 }
552
553 fasp_darray_cp(n, p[i-1], w);
554
555 fasp_blas_darray_ax(n, rs[i-1], w);
556
557 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
558
559 /* apply preconditioner */
560 if ( pc == NULL )
561 fasp_darray_cp(n, w, r);
562 else
563 pc->fct(w, r, pc->data);
564
565 fasp_blas_darray_axpy(n, 1.0, r, x->val);
566
567 // Check: prevent false convergence
568 if ( relres < tol && iter >= MIN_ITER ) {
569
570 REAL computed_relres = relres;
571
572 // compute residual
573 fasp_darray_cp(n, b->val, r);
574 fasp_blas_dbsr_aAxpy(-1.0, A, x->val, r);
575 r_norm = fasp_blas_darray_norm2(n, r);
576
577 switch ( StopType ) {
578 case STOP_REL_RES:
579 absres = r_norm;
580 relres = absres/absres0;
581 break;
582 case STOP_REL_PRECRES:
583 if ( pc == NULL )
584 fasp_darray_cp(n, r, w);
585 else
586 pc->fct(r, w, pc->data);
587 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
588 relres = absres/absres0;
589 break;
590 case STOP_MOD_REL_RES:
591 absres = r_norm;
593 relres = absres/normu;
594 break;
595 }
596
597 norms[iter] = relres;
598
599 if ( relres < tol ) {
600 break;
601 }
602 else { // Need to restart
603 fasp_darray_cp(n, r, p[0]); i = 0;
604 }
605
606 if ( PrtLvl >= PRINT_MORE ) {
607 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
608 }
609
610
611 } /* end of convergence check */
612
613 /* compute residual vector and continue loop */
614 for ( j = i; j > 0; j-- ) {
615 rs[j-1] = -s[j-1]*rs[j];
616 rs[j] = c[j-1]*rs[j];
617 }
618
619 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
620
621 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
622
623 if ( i ) {
624 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
625 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
626 }
627
628 } /* end of main while loop */
629
630FINISHED:
631 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
632
633 /*-------------------------------------------
634 * Clean up workspace
635 *------------------------------------------*/
636 fasp_mem_free(work); work = NULL;
637 fasp_mem_free(p); p = NULL;
638 fasp_mem_free(hh); hh = NULL;
639 fasp_mem_free(norms); norms = NULL;
640
641#if DEBUG_MODE > 0
642 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
643#endif
644
645 if ( iter >= MaxIt )
646 return ERROR_SOLVER_MAXIT;
647 else
648 return iter;
649}
650
676 dvector *b,
677 dvector *x,
678 precond *pc,
679 const REAL tol,
680 const INT MaxIt,
681 const SHORT restart,
682 const SHORT StopType,
683 const SHORT PrtLvl)
684{
685 const INT n = b->row;
686 const INT MIN_ITER = 0;
687
688 // local variables
689 INT iter = 0;
690 int i, j, k; // must be signed! -zcs
691
692 REAL r_norm, r_normb, gamma, t;
693 REAL absres0 = BIGREAL, absres = BIGREAL;
694 REAL relres = BIGREAL, normu = BIGREAL;
695
696 // allocate temp memory (need about (restart+4)*n REAL numbers)
697 REAL *c = NULL, *s = NULL, *rs = NULL;
698 REAL *norms = NULL, *r = NULL, *w = NULL;
699 REAL *work = NULL;
700 REAL **p = NULL, **hh = NULL;
701
702 INT Restart = MIN(restart, MaxIt);
703 INT Restart1 = Restart + 1;
704 LONG worksize = (Restart+4)*(Restart+n)+1-n;
705
706 /* allocate memory and setup temp work space */
707 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
708
709 // Output some info for debugging
710 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GMRes solver (BLC) ...\n");
711
712#if DEBUG_MODE > 0
713 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
714 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
715#endif
716
717 /* check whether memory is enough for GMRES */
718 while ( (work == NULL) && (Restart > 5) ) {
719 Restart = Restart - 5;
720 Restart1 = Restart + 1;
721 worksize = (Restart+4)*(Restart+n)+1-n;
722 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
723 }
724
725 if ( work == NULL ) {
726 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
727 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
728 }
729
730 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
731 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
732 }
733
734 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
735 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
736 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
737
738 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
739
740 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
741
742 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
743
744 // compute initial residual: r = b-A*x
745 fasp_darray_cp(n, b->val, p[0]);
746 fasp_blas_dblc_aAxpy(-1.0, A, x->val, p[0]);
747 r_norm = fasp_blas_darray_norm2(n,p[0]);
748
749 // compute stopping criteria
750 switch (StopType) {
751 case STOP_REL_RES:
752 absres0 = MAX(SMALLREAL,r_norm);
753 relres = r_norm/absres0;
754 break;
755 case STOP_REL_PRECRES:
756 if ( pc == NULL )
757 fasp_darray_cp(n, p[0], r);
758 else
759 pc->fct(p[0], r, pc->data);
760 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
761 absres0 = MAX(SMALLREAL,r_normb);
762 relres = r_normb/absres0;
763 break;
764 case STOP_MOD_REL_RES:
766 absres0 = r_norm;
767 relres = absres0/normu;
768 break;
769 default:
770 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
771 goto FINISHED;
772 }
773
774 // if initial residual is small, no need to iterate!
775 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
776
777 // output iteration information if needed
778 fasp_itinfo(PrtLvl,StopType,0,relres,absres0,0.0);
779
780 // store initial residual
781 norms[0] = relres;
782
783 /* GMRES(M) outer iteration */
784 while ( iter < MaxIt && relres > tol ) {
785
786 rs[0] = r_norm;
787
788 t = 1.0 / r_norm;
789
790 fasp_blas_darray_ax(n, t, p[0]);
791
792 /* RESTART CYCLE (right-preconditioning) */
793 i = 0;
794 while ( i < Restart && iter < MaxIt ) {
795
796 i++; iter++;
797
798 /* apply preconditioner */
799 if ( pc == NULL )
800 fasp_darray_cp(n, p[i-1], r);
801 else
802 pc->fct(p[i-1], r, pc->data);
803
804 fasp_blas_dblc_mxv(A, r, p[i]);
805
806 /* Modified Gram_Schmidt orthogonalization */
807 for ( j = 0; j < i; j++ ) {
808 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
809 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
810 }
811 t = fasp_blas_darray_norm2(n, p[i]);
812 hh[i][i-1] = t;
813
814 if ( ABS(t) > SMALLREAL ) { // If t=0, we get solution subspace
815 t = 1.0/t;
816 fasp_blas_darray_ax(n, t, p[i]);
817 }
818
819 for ( j = 1; j < i; ++j ) {
820 t = hh[j-1][i-1];
821 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
822 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
823 }
824 t = hh[i][i-1]*hh[i][i-1];
825 t += hh[i-1][i-1]*hh[i-1][i-1];
826
827 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
828 c[i-1] = hh[i-1][i-1] / gamma;
829 s[i-1] = hh[i][i-1] / gamma;
830 rs[i] = -s[i-1]*rs[i-1];
831 rs[i-1] = c[i-1]*rs[i-1];
832 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
833
834 absres = r_norm = fabs(rs[i]);
835
836 relres = absres/absres0;
837
838 norms[iter] = relres;
839
840 // output iteration information if needed
841 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
842 norms[iter]/norms[iter-1]);
843
844 // exit restart cycle if reaches tolerance
845 if ( relres < tol && iter >= MIN_ITER ) break;
846
847 } /* end of restart cycle */
848
849 /* compute solution, first solve upper triangular system */
850 rs[i-1] = rs[i-1] / hh[i-1][i-1];
851 for ( k = i-2; k >= 0; k-- ) {
852 t = 0.0;
853 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
854 t += rs[k];
855 rs[k] = t / hh[k][k];
856 }
857
858 fasp_darray_cp(n, p[i-1], w);
859
860 fasp_blas_darray_ax(n, rs[i-1], w);
861
862 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
863
864 /* apply preconditioner */
865 if ( pc == NULL )
866 fasp_darray_cp(n, w, r);
867 else
868 pc->fct(w, r, pc->data);
869
870 fasp_blas_darray_axpy(n, 1.0, r, x->val);
871
872 // Check: prevent false convergence
873 if ( relres < tol && iter >= MIN_ITER ) {
874
875 REAL computed_relres = relres;
876
877 // compute residual
878 fasp_darray_cp(n, b->val, r);
879 fasp_blas_dblc_aAxpy(-1.0, A, x->val, r);
880 r_norm = fasp_blas_darray_norm2(n, r);
881
882 switch ( StopType ) {
883 case STOP_REL_RES:
884 absres = r_norm;
885 relres = absres/absres0;
886 break;
887 case STOP_REL_PRECRES:
888 if ( pc == NULL )
889 fasp_darray_cp(n, r, w);
890 else
891 pc->fct(r, w, pc->data);
892 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
893 relres = absres/absres0;
894 break;
895 case STOP_MOD_REL_RES:
896 absres = r_norm;
898 relres = absres/normu;
899 break;
900 }
901
902 norms[iter] = relres;
903
904 if ( relres < tol ) {
905 break;
906 }
907 else { // Need to restart
908 fasp_darray_cp(n, r, p[0]); i = 0;
909 }
910
911 if ( PrtLvl >= PRINT_MORE ) {
912 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
913 }
914
915 } /* end of convergence check */
916
917 /* compute residual vector and continue loop */
918 for ( j = i; j > 0; j-- ) {
919 rs[j-1] = -s[j-1]*rs[j];
920 rs[j] = c[j-1]*rs[j];
921 }
922
923 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
924
925 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
926
927 if ( i ) {
928 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
929 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
930 }
931
932 } /* end of main while loop */
933
934FINISHED:
935 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
936
937 /*-------------------------------------------
938 * Clean up workspace
939 *------------------------------------------*/
940 fasp_mem_free(work); work = NULL;
941 fasp_mem_free(p); p = NULL;
942 fasp_mem_free(hh); hh = NULL;
943 fasp_mem_free(norms); norms = NULL;
944
945#if DEBUG_MODE > 0
946 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
947#endif
948
949 if ( iter >= MaxIt )
950 return ERROR_SOLVER_MAXIT;
951 else
952 return iter;
953}
954
980 dvector *b,
981 dvector *x,
982 precond *pc,
983 const REAL tol,
984 const INT MaxIt,
985 const SHORT restart,
986 const SHORT StopType,
987 const SHORT PrtLvl)
988{
989 const INT n = b->row;
990 const INT MIN_ITER = 0;
991
992 // local variables
993 INT iter = 0;
994 int i, j, k; // must be signed! -zcs
995
996 REAL r_norm, r_normb, gamma, t;
997 REAL absres0 = BIGREAL, absres = BIGREAL;
998 REAL relres = BIGREAL, normu = BIGREAL;
999
1000 // allocate temp memory (need about (restart+4)*n REAL numbers)
1001 REAL *c = NULL, *s = NULL, *rs = NULL;
1002 REAL *norms = NULL, *r = NULL, *w = NULL;
1003 REAL *work = NULL;
1004 REAL **p = NULL, **hh = NULL;
1005
1006 INT Restart = MIN(restart, MaxIt);
1007 INT Restart1 = Restart + 1;
1008 LONG worksize = (Restart+4)*(Restart+n)+1-n;
1009
1010 /* allocate memory and setup temp work space */
1011 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1012
1013 // Output some info for debugging
1014 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GMRes solver (STR) ...\n");
1015
1016#if DEBUG_MODE > 0
1017 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1018 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1019#endif
1020
1021 /* check whether memory is enough for GMRES */
1022 while ( (work == NULL) && (Restart > 5) ) {
1023 Restart = Restart - 5;
1024 Restart1 = Restart + 1;
1025 worksize = (Restart+4)*(Restart+n)+1-n;
1026 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1027 }
1028
1029 if ( work == NULL ) {
1030 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
1031 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1032 }
1033
1034 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
1035 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
1036 }
1037
1038 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1039 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1040 norms = (REAL *) fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1041
1042 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
1043
1044 for ( i = 0; i < Restart1; i++ ) p[i] = s + Restart + i*n;
1045
1046 for ( i = 0; i < Restart1; i++ ) hh[i] = p[Restart] + n + i*Restart;
1047
1048 // compute initial residual: r = b-A*x
1049 fasp_darray_cp(n, b->val, p[0]);
1050 fasp_blas_dstr_aAxpy(-1.0, A, x->val, p[0]);
1051 r_norm = fasp_blas_darray_norm2(n,p[0]);
1052
1053 // compute stopping criteria
1054 switch (StopType) {
1055 case STOP_REL_RES:
1056 absres0 = MAX(SMALLREAL,r_norm);
1057 relres = r_norm/absres0;
1058 break;
1059 case STOP_REL_PRECRES:
1060 if ( pc == NULL )
1061 fasp_darray_cp(n, p[0], r);
1062 else
1063 pc->fct(p[0], r, pc->data);
1064 r_normb = sqrt(fasp_blas_darray_dotprod(n,p[0],r));
1065 absres0 = MAX(SMALLREAL,r_normb);
1066 relres = r_normb/absres0;
1067 break;
1068 case STOP_MOD_REL_RES:
1069 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1070 absres0 = r_norm;
1071 relres = absres0/normu;
1072 break;
1073 default:
1074 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1075 goto FINISHED;
1076 }
1077
1078 // if initial residual is small, no need to iterate!
1079 if ( relres < tol || absres0 < 1e-312*tol ) goto FINISHED;
1080
1081 // output iteration information if needed
1082 fasp_itinfo(PrtLvl,StopType,0,relres,absres0,0.0);
1083
1084 // store initial residual
1085 norms[0] = relres;
1086
1087 /* GMRES(M) outer iteration */
1088 while ( iter < MaxIt && relres > tol ) {
1089
1090 rs[0] = r_norm;
1091
1092 t = 1.0 / r_norm;
1093
1094 fasp_blas_darray_ax(n, t, p[0]);
1095
1096 /* RESTART CYCLE (right-preconditioning) */
1097 i = 0;
1098 while ( i < Restart && iter < MaxIt ) {
1099
1100 i++; iter++;
1101
1102 /* apply preconditioner */
1103 if ( pc == NULL )
1104 fasp_darray_cp(n, p[i-1], r);
1105 else
1106 pc->fct(p[i-1], r, pc->data);
1107
1108 fasp_blas_dstr_mxv(A, r, p[i]);
1109
1110 /* Modified Gram_Schmidt orthogonalization */
1111 for ( j = 0; j < i; j++ ) {
1112 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1113 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
1114 }
1115 t = fasp_blas_darray_norm2(n, p[i]);
1116 hh[i][i-1] = t;
1117
1118 if ( ABS(t) > SMALLREAL ) { // If t=0, we get solution subspace
1119 t = 1.0/t;
1120 fasp_blas_darray_ax(n, t, p[i]);
1121 }
1122
1123 for ( j = 1; j < i; ++j ) {
1124 t = hh[j-1][i-1];
1125 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1126 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1127 }
1128 t = hh[i][i-1]*hh[i][i-1];
1129 t += hh[i-1][i-1]*hh[i-1][i-1];
1130
1131 gamma = MAX(sqrt(t), SMALLREAL); // Possible breakdown?
1132 c[i-1] = hh[i-1][i-1] / gamma;
1133 s[i-1] = hh[i][i-1] / gamma;
1134 rs[i] = -s[i-1]*rs[i-1];
1135 rs[i-1] = c[i-1]*rs[i-1];
1136 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1137
1138 absres = r_norm = fabs(rs[i]);
1139
1140 relres = absres/absres0;
1141
1142 norms[iter] = relres;
1143
1144 // output iteration information if needed
1145 fasp_itinfo(PrtLvl, StopType, iter, relres, absres,
1146 norms[iter]/norms[iter-1]);
1147
1148 // exit restart cycle if reaches tolerance
1149 if ( relres < tol && iter >= MIN_ITER ) break;
1150
1151 } /* end of restart cycle */
1152
1153 /* compute solution, first solve upper triangular system */
1154 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1155 for ( k = i-2; k >= 0; k-- ) {
1156 t = 0.0;
1157 for ( j = k+1; j < i; j++ ) t -= hh[k][j]*rs[j];
1158 t += rs[k];
1159 rs[k] = t / hh[k][k];
1160 }
1161
1162 fasp_darray_cp(n, p[i-1], w);
1163
1164 fasp_blas_darray_ax(n, rs[i-1], w);
1165
1166 for ( j = i-2; j >= 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1167
1168 /* apply preconditioner */
1169 if ( pc == NULL )
1170 fasp_darray_cp(n, w, r);
1171 else
1172 pc->fct(w, r, pc->data);
1173
1174 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1175
1176 // Check: prevent false convergence
1177 if ( relres < tol && iter >= MIN_ITER ) {
1178
1179 REAL computed_relres = relres;
1180
1181 // compute residual
1182 fasp_darray_cp(n, b->val, r);
1183 fasp_blas_dstr_aAxpy(-1.0, A, x->val, r);
1184 r_norm = fasp_blas_darray_norm2(n, r);
1185
1186 switch ( StopType ) {
1187 case STOP_REL_RES:
1188 absres = r_norm;
1189 relres = absres/absres0;
1190 break;
1191 case STOP_REL_PRECRES:
1192 if ( pc == NULL )
1193 fasp_darray_cp(n, r, w);
1194 else
1195 pc->fct(r, w, pc->data);
1196 absres = sqrt(fasp_blas_darray_dotprod(n,w,r));
1197 relres = absres/absres0;
1198 break;
1199 case STOP_MOD_REL_RES:
1200 absres = r_norm;
1201 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(n,x->val));
1202 relres = absres/normu;
1203 break;
1204 }
1205
1206 norms[iter] = relres;
1207
1208 if ( relres < tol ) {
1209 break;
1210 }
1211 else { // Need to restart
1212 fasp_darray_cp(n, r, p[0]); i = 0;
1213 }
1214
1215 if ( PrtLvl >= PRINT_MORE ) {
1216 ITS_COMPRES(computed_relres); ITS_REALRES(relres);
1217 }
1218
1219 } /* end of convergence check */
1220
1221 /* compute residual vector and continue loop */
1222 for ( j = i; j > 0; j-- ) {
1223 rs[j-1] = -s[j-1]*rs[j];
1224 rs[j] = c[j-1]*rs[j];
1225 }
1226
1227 if ( i ) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1228
1229 for ( j = i-1 ; j > 0; j-- ) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1230
1231 if ( i ) {
1232 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1233 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1234 }
1235
1236 } /* end of main while loop */
1237
1238FINISHED:
1239 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1240
1241 /*-------------------------------------------
1242 * Clean up workspace
1243 *------------------------------------------*/
1244 fasp_mem_free(work); work = NULL;
1245 fasp_mem_free(p); p = NULL;
1246 fasp_mem_free(hh); hh = NULL;
1247 fasp_mem_free(norms); norms = NULL;
1248
1249#if DEBUG_MODE > 0
1250 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1251#endif
1252
1253 if ( iter >= MaxIt )
1254 return ERROR_SOLVER_MAXIT;
1255 else
1256 return iter;
1257}
1258
1284 dvector *b,
1285 dvector *x,
1286 precond *pc,
1287 const REAL tol,
1288 const INT MaxIt,
1289 const SHORT restart,
1290 const SHORT StopType,
1291 const SHORT PrtLvl)
1292{
1293 const INT n = b->row;
1294 const INT min_iter = 0;
1295
1296 // local variables
1297 INT iter = 0;
1298 int i, j, k; // must be signed! -zcs
1299
1300 REAL epsmac = SMALLREAL;
1301 REAL r_norm, b_norm, den_norm;
1302 REAL epsilon, gamma, t;
1303
1304 // allocate temp memory (need about (restart+4)*n REAL numbers)
1305 REAL *c = NULL, *s = NULL, *rs = NULL;
1306 REAL *norms = NULL, *r = NULL, *w = NULL;
1307 REAL *work = NULL;
1308 REAL **p = NULL, **hh = NULL;
1309
1310 INT Restart = restart;
1311 INT Restart1 = Restart + 1;
1312 LONG worksize = (Restart+4)*(Restart+n)+1-n;
1313
1314 // Output some info for debugging
1315 if ( PrtLvl > PRINT_NONE ) printf("\nCalling GMRes solver (MatFree) ...\n");
1316
1317#if DEBUG_MODE > 0
1318 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1319 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1320#endif
1321
1322 /* allocate memory and setup temp work space */
1323 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1324
1325 /* check whether memory is enough for GMRES */
1326 while ( (work == NULL) && (Restart > 5) ) {
1327 Restart = Restart - 5;
1328 worksize = (Restart+4)*(Restart+n)+1-n;
1329 work = (REAL *) fasp_mem_calloc(worksize, sizeof(REAL));
1330 Restart1 = Restart + 1;
1331 }
1332
1333 if ( work == NULL ) {
1334 printf("### ERROR: No enough memory! [%s:%d]\n", __FILE__, __LINE__ );
1335 fasp_chkerr(ERROR_ALLOC_MEM, __FUNCTION__);
1336 }
1337
1338 if ( PrtLvl > PRINT_MIN && Restart < restart ) {
1339 printf("### WARNING: GMRES restart number set to %d!\n", Restart);
1340 }
1341
1342 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1343 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1344 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1345
1346 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + Restart;
1347
1348 for (i = 0; i < Restart1; i ++) p[i] = s + Restart + i*n;
1349
1350 for (i = 0; i < Restart1; i ++) hh[i] = p[Restart] + n + i*Restart;
1351
1352 /* initialization */
1353 mf->fct(mf->data, x->val, p[0]);
1354 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, p[0]);
1355
1356 b_norm = fasp_blas_darray_norm2(n, b->val);
1357 r_norm = fasp_blas_darray_norm2(n, p[0]);
1358
1359 if ( PrtLvl > PRINT_NONE) {
1360 norms[0] = r_norm;
1361 if ( PrtLvl >= PRINT_SOME ) {
1362 ITS_PUTNORM("right-hand side", b_norm);
1363 ITS_PUTNORM("residual", r_norm);
1364 }
1365 }
1366
1367 if (b_norm > 0.0) den_norm = b_norm;
1368 else den_norm = r_norm;
1369
1370 epsilon = tol*den_norm;
1371
1372 /* outer iteration cycle */
1373 while (iter < MaxIt) {
1374
1375 rs[0] = r_norm;
1376 if (r_norm == 0.0) {
1377 fasp_mem_free(work); work = NULL;
1378 fasp_mem_free(p); p = NULL;
1379 fasp_mem_free(hh); hh = NULL;
1380 fasp_mem_free(norms); norms = NULL;
1381 return iter;
1382 }
1383
1384 if (r_norm <= epsilon && iter >= min_iter) {
1385 mf->fct(mf->data, x->val, r);
1386 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1387 r_norm = fasp_blas_darray_norm2(n, r);
1388
1389 if (r_norm <= epsilon) {
1390 break;
1391 }
1392 else {
1393 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1394 }
1395 }
1396
1397 t = 1.0 / r_norm;
1398 //for (j = 0; j < n; j ++) p[0][j] *= t;
1399 fasp_blas_darray_ax(n, t, p[0]);
1400
1401 /* RESTART CYCLE (right-preconditioning) */
1402 i = 0;
1403 while (i < Restart && iter < MaxIt) {
1404
1405 i ++; iter ++;
1406
1407 /* apply preconditioner */
1408 if (pc == NULL)
1409 fasp_darray_cp(n, p[i-1], r);
1410 else
1411 pc->fct(p[i-1], r, pc->data);
1412
1413 mf->fct(mf->data, r, p[i]);
1414
1415 /* modified Gram_Schmidt */
1416 for (j = 0; j < i; j ++) {
1417 hh[j][i-1] = fasp_blas_darray_dotprod(n, p[j], p[i]);
1418 fasp_blas_darray_axpy(n, -hh[j][i-1], p[j], p[i]);
1419 }
1420 t = fasp_blas_darray_norm2(n, p[i]);
1421 hh[i][i-1] = t;
1422 if (t != 0.0) {
1423 t = 1.0/t;
1424 //for (j = 0; j < n; j ++) p[i][j] *= t;
1425 fasp_blas_darray_ax(n, t, p[i]);
1426 }
1427
1428 for (j = 1; j < i; ++j) {
1429 t = hh[j-1][i-1];
1430 hh[j-1][i-1] = s[j-1]*hh[j][i-1] + c[j-1]*t;
1431 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
1432 }
1433 t= hh[i][i-1]*hh[i][i-1];
1434 t+= hh[i-1][i-1]*hh[i-1][i-1];
1435 gamma = sqrt(t);
1436 if (gamma == 0.0) gamma = epsmac;
1437 c[i-1] = hh[i-1][i-1] / gamma;
1438 s[i-1] = hh[i][i-1] / gamma;
1439 rs[i] = -s[i-1]*rs[i-1];
1440 rs[i-1] = c[i-1]*rs[i-1];
1441 hh[i-1][i-1] = s[i-1]*hh[i][i-1] + c[i-1]*hh[i-1][i-1];
1442 r_norm = fabs(rs[i]);
1443
1444 norms[iter] = r_norm;
1445
1446 if (b_norm > 0 ) {
1447 fasp_itinfo(PrtLvl,StopType,iter,norms[iter]/b_norm,
1448 norms[iter],norms[iter]/norms[iter-1]);
1449 }
1450 else {
1451 fasp_itinfo(PrtLvl,StopType,iter,norms[iter],norms[iter],
1452 norms[iter]/norms[iter-1]);
1453 }
1454
1455 /* should we exit restart cycle? */
1456 if (r_norm <= epsilon && iter >= min_iter) {
1457 break;
1458 }
1459 } /* end of restart cycle */
1460
1461 /* now compute solution, first solve upper triangular system */
1462 rs[i-1] = rs[i-1] / hh[i-1][i-1];
1463 for (k = i-2; k >= 0; k --) {
1464 t = 0.0;
1465 for (j = k+1; j < i; j ++) t -= hh[k][j]*rs[j];
1466
1467 t += rs[k];
1468 rs[k] = t / hh[k][k];
1469 }
1470 fasp_darray_cp(n, p[i-1], w);
1471 //for (j = 0; j < n; j ++) w[j] *= rs[i-1];
1472 fasp_blas_darray_ax(n, rs[i-1], w);
1473 for (j = i-2; j >= 0; j --) fasp_blas_darray_axpy(n, rs[j], p[j], w);
1474
1475 /* apply preconditioner */
1476 if (pc == NULL)
1477 fasp_darray_cp(n, w, r);
1478 else
1479 pc->fct(w, r, pc->data);
1480
1481 fasp_blas_darray_axpy(n, 1.0, r, x->val);
1482
1483 if (r_norm <= epsilon && iter >= min_iter) {
1484 mf->fct(mf->data, x->val, r);
1485 fasp_blas_darray_axpby(n, 1.0, b->val, -1.0, r);
1486 r_norm = fasp_blas_darray_norm2(n, r);
1487
1488 if (r_norm <= epsilon) {
1489 break;
1490 }
1491 else {
1492 if (PrtLvl >= PRINT_SOME) ITS_FACONV;
1493 fasp_darray_cp(n, r, p[0]); i = 0;
1494 }
1495 } /* end of convergence check */
1496
1497 /* compute residual vector and continue loop */
1498 for (j = i; j > 0; j--) {
1499 rs[j-1] = -s[j-1]*rs[j];
1500 rs[j] = c[j-1]*rs[j];
1501 }
1502
1503 if (i) fasp_blas_darray_axpy(n, rs[i]-1.0, p[i], p[i]);
1504
1505 for (j = i-1 ; j > 0; j --) fasp_blas_darray_axpy(n, rs[j], p[j], p[i]);
1506
1507 if (i) {
1508 fasp_blas_darray_axpy(n, rs[0]-1.0, p[0], p[0]);
1509 fasp_blas_darray_axpy(n, 1.0, p[i], p[0]);
1510 }
1511 } /* end of iteration while loop */
1512
1513 if (PrtLvl > PRINT_NONE) ITS_FINAL(iter,MaxIt,r_norm);
1514
1515 /*-------------------------------------------
1516 * Clean up workspace
1517 *------------------------------------------*/
1518 fasp_mem_free(work); work = NULL;
1519 fasp_mem_free(p); p = NULL;
1520 fasp_mem_free(hh); hh = NULL;
1521 fasp_mem_free(norms); norms = NULL;
1522
1523#if DEBUG_MODE > 0
1524 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1525#endif
1526
1527 if (iter>=MaxIt)
1528 return ERROR_SOLVER_MAXIT;
1529 else
1530 return iter;
1531}
1532
1533#if 0
1534static double estimate_spectral_radius (const double **A, int n, size_t k = 20)
1535{
1536 double *x = (double *)malloc(n* sizeof(double));
1537 double *y = (double *)malloc(n* sizeof(double));
1538 double *z = (double *)malloc(n* sizeof(double));
1539 double t;
1540 int i1,j1;
1541
1542 // initialize x to random values in [0,1)
1543 // cusp::copy(cusp::detail::random_reals<ValueType>(N), x);
1544 dvector px;
1545 px.row = n;
1546 px.val = x;
1547
1548 fasp_dvec_rand(n, &px);
1549
1550 for(size_t i = 0; i < k; i++)
1551 {
1552 //cusp::blas::scal(x, ValueType(1.0) / cusp::blas::nrmmax(x));
1553 t= 1.0/ fasp_blas_darray_norminf(n, px);
1554 for(i1= 0; i1 <n; i1++) x[i1] *= t;
1555
1556 //cusp::multiply(A, x, y);
1557
1558 for(i1= 0; i1 <n; i1++) {
1559 t= 0.0
1560 for(j1= 0; j1 <n; j1++) t += A[i1][j1] * x[j1];
1561 y[i1] = t; }
1562 // x.swap(y);
1563 for(i1= 0; i1 <n; i1++) z[i1] = x[i1];
1564 for(i1= 0; i1 <n; i1++) x[i1] = y[i1];
1565 for(i1= 0; i1 <n; i1++) y[i1] = z[i1];
1566 }
1567
1568 free(x);
1569 free(y);
1570 free(z);
1571
1572 if (k == 0)
1573 return 0;
1574 else
1575 //return cusp::blas::nrm2(x) / cusp::blas::nrm2(y);
1577}
1578
1579static double spectral_radius (dCSRmat *A,
1580 const SHORT restart)
1581{
1582 const INT n = A->row;
1583 const INT MIN_ITER = 0;
1584
1585 // local variables
1586 INT iter = 0;
1587 INT Restart1 = restart + 1;
1588 INT i, j, k;
1589
1590 REAL r_norm, den_norm;
1591 REAL epsilon, gamma, t;
1592
1593 REAL *c = NULL, *s = NULL, *rs = NULL;
1594 REAL *norms = NULL, *r = NULL, *w = NULL;
1595 REAL **p = NULL, **hh = NULL;
1596 REAL *work = NULL;
1597
1598 /* allocate memory */
1599 work = (REAL *)fasp_mem_calloc((restart+4)*(restart+n)+1-n, sizeof(REAL));
1600 p = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1601 hh = (REAL **)fasp_mem_calloc(Restart1, sizeof(REAL *));
1602
1603 norms = (REAL *)fasp_mem_calloc(MaxIt+1, sizeof(REAL));
1604
1605 r = work; w = r + n; rs = w + n; c = rs + Restart1; s = c + restart;
1606
1607 for (i = 0; i < Restart1; i ++) p[i] = s + restart + i*n;
1608 for (i = 0; i < Restart1; i ++) hh[i] = p[restart] + n + i*restart;
1609
1610 /* initialization */
1611 dvector p0;
1612 p0.row = n;
1613 p0.val = p[0];
1614 fasp_dvec_rand(n, &p0);
1615
1616 r_norm = fasp_blas_darray_norm2(n, p[0]);
1617 t = 1.0 / r_norm;
1618 for (j = 0; j < n; j ++) p[0][j] *= t;
1619
1620 int maxiter = MIN(n, restart) ;
1621 for ( j = 0; j < maxiter; j++ ) {
1622 fasp_blas_bdbsr_mxv(A, p[j], p[j+1]);
1623
1624 for( i = 0; i <= j; i++ ) {
1625 hh[i][j] = fasp_blas_darray_dotprod(n, p[i], p[j+1]);
1626 fasp_blas_darray_axpy(n, -hh[i][j], p[i], p[ j+1 ]);
1627 }
1628
1629 hh[j+1][j] = fasp_blas_darray_norm2 (n, p[j+1]);
1630 if ( hh[j+1][j] < 1e-10) break;
1631 t = 1.0/hh[j+1][j];
1632 for (k = 0; k < n; k ++) p[j+1][k] *= t;
1633 }
1634
1635 H = (REAL **)fasp_mem_calloc(j, sizeof(REAL *));
1636 H[0] = (REAL *)fasp_mem_calloc(j*j, sizeof(REAL));
1637 for (i = 1; i < j; i ++) H[i] = H[i-1] + j;
1638
1639
1640 for( size_t row = 0; row < j; row++ )
1641 for( size_t col = 0; col < j; col++ )
1642 H[row][col] = hh[row][col];
1643
1644 double spectral_radius = estimate_spectral_radius( H, j, 20);
1645
1646 /*-------------------------------------------
1647 * Clean up workspace
1648 *------------------------------------------*/
1649 fasp_mem_free(work); work = NULL;
1650 fasp_mem_free(p); p = NULL;
1651 fasp_mem_free(hh); hh = NULL;
1652 fasp_mem_free(norms); norms = NULL;
1653 fasp_mem_free(H[0]); H[0] = NULL;
1654 fasp_mem_free(H); H = NULL;
1655
1656 return spectral_radius;
1657}
1658#endif
1659
1660/*---------------------------------*/
1661/*-- End of File --*/
1662/*---------------------------------*/
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
void fasp_dvec_rand(const INT n, dvector *x)
Generate fake random REAL vector in the range from 0 to 1.
Definition: AuxVector.c:192
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
REAL fasp_blas_darray_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: BlaArray.c:686
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
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
void fasp_blas_dstr_mxv(const dSTRmat *A, const REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: BlaSpmvSTR.c:131
INT fasp_solver_dcsr_pgmres(dCSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Right preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:67
INT fasp_solver_dblc_pgmres(dBLCmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:675
INT fasp_solver_pgmres(mxv_matfree *mf, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Solve "Ax=b" using PGMRES (right preconditioned) iterative method.
Definition: KryPgmres.c:1283
INT fasp_solver_dbsr_pgmres(dBSRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:370
INT fasp_solver_dstr_pgmres(dSTRmat *A, dvector *b, dvector *x, precond *pc, const REAL tol, const INT MaxIt, const SHORT restart, const SHORT StopType, const SHORT PrtLvl)
Preconditioned GMRES method for solving Au=b.
Definition: KryPgmres.c:979
Main header file for the FASP project.
#define MIN(a, b)
Definition: fasp.h:75
#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 ABS(a)
Definition: fasp.h:76
#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 BIGREAL
Some global constants.
Definition: fasp_const.h:248
#define PRINT_SOME
Definition: fasp_const.h:75
#define ERROR_ALLOC_MEM
Definition: fasp_const.h:30
#define PRINT_MORE
Definition: fasp_const.h:76
#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
INT row
row number of matrix A, m
Definition: fasp.h:146
Structure matrix of REAL type.
Definition: fasp.h:308
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