Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KrySPcg.c
Go to the documentation of this file.
1
24#include <math.h>
25
26#include "fasp.h"
27#include "fasp_functs.h"
28
29/*---------------------------------*/
30/*-- Declare Private Functions --*/
31/*---------------------------------*/
32
33#include "KryUtil.inl"
34
35/*---------------------------------*/
36/*-- Public Functions --*/
37/*---------------------------------*/
38
60INT fasp_solver_dcsr_spcg (const dCSRmat *A,
61 const dvector *b,
62 dvector *u,
63 precond *pc,
64 const REAL tol,
65 const INT MaxIt,
66 const SHORT StopType,
67 const SHORT PrtLvl)
68{
69 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
70 const INT m = b->row;
71 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
72 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
73
74 // local variables
75 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
76 REAL absres0 = BIGREAL, absres = BIGREAL;
77 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
78 REAL reldiff, factor, normuinf;
79 REAL alpha, beta, temp1, temp2;
80 INT iter_best = 0; // initial best known iteration
81 REAL absres_best = BIGREAL; // initial best known residual
82
83 // allocate temp memory (need 5*m REAL numbers)
84 REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
85 REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
86
87 // Output some info for debuging
88 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe CG solver (CSR) ...\n");
89
90#if DEBUG_MODE > 0
91 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
92 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
93#endif
94
95 // r = b-A*u
96 fasp_darray_cp(m,b->val,r);
97 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
98
99 if (pc != NULL)
100 pc->fct(r,z,pc->data); /* Apply preconditioner */
101 else
102 fasp_darray_cp(m,r,z); /* No preconditioner */
103
104 // compute initial residuals
105 switch (StopType) {
106 case STOP_REL_RES:
107 absres0 = fasp_blas_darray_norm2(m,r);
108 normr0 = MAX(SMALLREAL,absres0);
109 relres = absres0/normr0;
110 break;
111 case STOP_REL_PRECRES:
112 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
113 normr0 = MAX(SMALLREAL,absres0);
114 relres = absres0/normr0;
115 break;
116 case STOP_MOD_REL_RES:
117 absres0 = fasp_blas_darray_norm2(m,r);
118 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
119 relres = absres0/normu;
120 break;
121 default:
122 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
123 goto FINISHED;
124 }
125
126 // if initial residual is small, no need to iterate!
127 if ( relres < tol ) goto FINISHED;
128
129 // output iteration information if needed
130 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
131
132 fasp_darray_cp(m,z,p);
133 temp1 = fasp_blas_darray_dotprod(m,z,r);
134
135 // main PCG loop
136 while ( iter++ < MaxIt ) {
137
138 // t=A*p
139 fasp_blas_dcsr_mxv(A,p,t);
140
141 // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
142 temp2 = fasp_blas_darray_dotprod(m,t,p);
143 if ( ABS(temp2) > SMALLREAL2 ) {
144 alpha = temp1/temp2;
145 }
146 else { // Possible breakdown
147 goto RESTORE_BESTSOL;
148 }
149
150 // u_k=u_{k-1} + alpha_k*p_{k-1}
151 fasp_blas_darray_axpy(m,alpha,p,u->val);
152
153 // r_k=r_{k-1} - alpha_k*A*p_{k-1}
154 fasp_blas_darray_axpy(m,-alpha,t,r);
155
156 // compute residuals
157 switch ( StopType ) {
158 case STOP_REL_RES:
159 absres = fasp_blas_darray_norm2(m,r);
160 relres = absres/normr0;
161 break;
162 case STOP_REL_PRECRES:
163 // z = B(r)
164 if ( pc != NULL )
165 pc->fct(r,z,pc->data); /* Apply preconditioner */
166 else
167 fasp_darray_cp(m,r,z); /* No preconditioner */
168 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
169 relres = absres/normr0;
170 break;
171 case STOP_MOD_REL_RES:
172 absres = fasp_blas_darray_norm2(m,r);
173 relres = absres/normu;
174 break;
175 }
176
177 // compute reducation factor of residual ||r||
178 factor = absres/absres0;
179
180 // output iteration information if needed
181 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
182
183 // if the solution is NAN, restore the best solution
184 if ( fasp_dvec_isnan(u) ) {
185 absres = BIGREAL;
186 goto RESTORE_BESTSOL;
187 }
188
189 // safety net check: save the best-so-far solution
190 if ( absres < absres_best - maxdiff) {
191 absres_best = absres;
192 iter_best = iter;
193 fasp_darray_cp(m,u->val,u_best);
194 }
195
196 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
197 normuinf = fasp_blas_darray_norminf(m, u->val);
198 if ( normuinf <= sol_inf_tol ) {
199 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
200 iter = ERROR_SOLVER_SOLSTAG;
201 break;
202 }
203
204 // Check II: if staggenated, try to restart
205 normu = fasp_blas_dvec_norm2(u);
206
207 // compute relative difference
208 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
209 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
210
211 if ( PrtLvl >= PRINT_MORE ) {
212 ITS_DIFFRES(reldiff,relres);
213 ITS_RESTART;
214 }
215
216 fasp_darray_cp(m,b->val,r);
217 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
218
219 // compute residuals
220 switch ( StopType ) {
221 case STOP_REL_RES:
222 absres = fasp_blas_darray_norm2(m,r);
223 relres = absres/normr0;
224 break;
225 case STOP_REL_PRECRES:
226 // z = B(r)
227 if ( pc != NULL )
228 pc->fct(r,z,pc->data); /* Apply preconditioner */
229 else
230 fasp_darray_cp(m,r,z); /* No preconditioner */
231 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
232 relres = absres/normr0;
233 break;
234 case STOP_MOD_REL_RES:
235 absres = fasp_blas_darray_norm2(m,r);
236 relres = absres/normu;
237 break;
238 }
239
240 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
241
242 if ( relres < tol )
243 break;
244 else {
245 if ( stag >= MaxStag ) {
246 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
247 iter = ERROR_SOLVER_STAG;
248 break;
249 }
250 fasp_darray_set(m,p,0.0);
251 ++stag;
252 ++restart_step;
253 }
254 } // end of staggnation check!
255
256 // Check III: prevent false convergence
257 if ( relres < tol ) {
258
259 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
260
261 // compute residual r = b - Ax again
262 fasp_darray_cp(m,b->val,r);
263 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
264
265 // compute residuals
266 switch ( StopType ) {
267 case STOP_REL_RES:
268 absres = fasp_blas_darray_norm2(m,r);
269 relres = absres/normr0;
270 break;
271 case STOP_REL_PRECRES:
272 // z = B(r)
273 if ( pc != NULL )
274 pc->fct(r,z,pc->data); /* Apply preconditioner */
275 else
276 fasp_darray_cp(m,r,z); /* No preconditioner */
277 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
278 relres = absres/normr0;
279 break;
280 case STOP_MOD_REL_RES:
281 absres = fasp_blas_darray_norm2(m,r);
282 relres = absres/normu;
283 break;
284 }
285
286 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
287
288 // check convergence
289 if ( relres < tol ) break;
290
291 if ( more_step >= MaxRestartStep ) {
292 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
293 iter = ERROR_SOLVER_TOLSMALL;
294 break;
295 }
296
297 // prepare for restarting the method
298 fasp_darray_set(m,p,0.0);
299 ++more_step;
300 ++restart_step;
301
302 } // end of safe-guard check!
303
304 // save residual for next iteration
305 absres0 = absres;
306
307 // compute z_k = B(r_k)
308 if ( StopType != STOP_REL_PRECRES ) {
309 if ( pc != NULL )
310 pc->fct(r,z,pc->data); /* Apply preconditioner */
311 else
312 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
313 }
314
315 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
316 temp2 = fasp_blas_darray_dotprod(m,z,r);
317 beta = temp2/temp1;
318 temp1 = temp2;
319
320 // compute p_k = z_k + beta_k*p_{k-1}
321 fasp_blas_darray_axpby(m,1.0,z,beta,p);
322
323 } // end of main PCG loop.
324
325RESTORE_BESTSOL: // restore the best-so-far solution if necessary
326 if ( iter != iter_best ) {
327
328 // compute best residual
329 fasp_darray_cp(m,b->val,r);
330 fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
331
332 switch ( StopType ) {
333 case STOP_REL_RES:
334 absres_best = fasp_blas_darray_norm2(m,r);
335 break;
336 case STOP_REL_PRECRES:
337 // z = B(r)
338 if ( pc != NULL )
339 pc->fct(r,z,pc->data); /* Apply preconditioner */
340 else
341 fasp_darray_cp(m,r,z); /* No preconditioner */
342 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
343 break;
344 case STOP_MOD_REL_RES:
345 absres_best = fasp_blas_darray_norm2(m,r);
346 break;
347 }
348
349 if ( absres > absres_best + maxdiff || isnan(absres) ) {
350 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
351 fasp_darray_cp(m,u_best,u->val);
352 relres = absres_best / normr0;
353 }
354 }
355
356FINISHED: // finish the iterative method
357 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
358
359 // clean up temp memory
360 fasp_mem_free(work); work = NULL;
361
362#if DEBUG_MODE > 0
363 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
364#endif
365
366 if ( iter > MaxIt )
367 return ERROR_SOLVER_MAXIT;
368 else
369 return iter;
370}
371
393INT fasp_solver_dblc_spcg (const dBLCmat *A,
394 const dvector *b,
395 dvector *u,
396 precond *pc,
397 const REAL tol,
398 const INT MaxIt,
399 const SHORT StopType,
400 const SHORT PrtLvl)
401{
402 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
403 const INT m = b->row;
404 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
405 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
406
407 // local variables
408 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
409 REAL absres0 = BIGREAL, absres = BIGREAL;
410 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
411 REAL reldiff, factor, normuinf;
412 REAL alpha, beta, temp1, temp2;
413 INT iter_best = 0; // initial best known iteration
414 REAL absres_best = BIGREAL; // initial best known residual
415
416 // allocate temp memory (need 5*m REAL numbers)
417 REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
418 REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
419
420 // Output some info for debuging
421 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe CG solver (BLC) ...\n");
422
423#if DEBUG_MODE > 0
424 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
425 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
426#endif
427
428 // r = b-A*u
429 fasp_darray_cp(m,b->val,r);
430 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
431
432 if (pc != NULL)
433 pc->fct(r,z,pc->data); /* Apply preconditioner */
434 else
435 fasp_darray_cp(m,r,z); /* No preconditioner */
436
437 // compute initial residuals
438 switch (StopType) {
439 case STOP_REL_RES:
440 absres0 = fasp_blas_darray_norm2(m,r);
441 normr0 = MAX(SMALLREAL,absres0);
442 relres = absres0/normr0;
443 break;
444 case STOP_REL_PRECRES:
445 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
446 normr0 = MAX(SMALLREAL,absres0);
447 relres = absres0/normr0;
448 break;
449 case STOP_MOD_REL_RES:
450 absres0 = fasp_blas_darray_norm2(m,r);
451 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
452 relres = absres0/normu;
453 break;
454 default:
455 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
456 goto FINISHED;
457 }
458
459 // if initial residual is small, no need to iterate!
460 if ( relres < tol ) goto FINISHED;
461
462 // output iteration information if needed
463 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
464
465 fasp_darray_cp(m,z,p);
466 temp1 = fasp_blas_darray_dotprod(m,z,r);
467
468 // main PCG loop
469 while ( iter++ < MaxIt ) {
470
471 // t=A*p
472 fasp_blas_dblc_mxv(A,p,t);
473
474 // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
475 temp2 = fasp_blas_darray_dotprod(m,t,p);
476 if ( ABS(temp2) > SMALLREAL2 ) {
477 alpha = temp1/temp2;
478 }
479 else { // Possible breakdown
480 goto RESTORE_BESTSOL;
481 }
482
483 // u_k=u_{k-1} + alpha_k*p_{k-1}
484 fasp_blas_darray_axpy(m,alpha,p,u->val);
485
486 // r_k=r_{k-1} - alpha_k*A*p_{k-1}
487 fasp_blas_darray_axpy(m,-alpha,t,r);
488
489 // compute residuals
490 switch ( StopType ) {
491 case STOP_REL_RES:
492 absres = fasp_blas_darray_norm2(m,r);
493 relres = absres/normr0;
494 break;
495 case STOP_REL_PRECRES:
496 // z = B(r)
497 if ( pc != NULL )
498 pc->fct(r,z,pc->data); /* Apply preconditioner */
499 else
500 fasp_darray_cp(m,r,z); /* No preconditioner */
501 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
502 relres = absres/normr0;
503 break;
504 case STOP_MOD_REL_RES:
505 absres = fasp_blas_darray_norm2(m,r);
506 relres = absres/normu;
507 break;
508 }
509
510 // compute reducation factor of residual ||r||
511 factor = absres/absres0;
512
513 // output iteration information if needed
514 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
515
516 // if the solution is NAN, restore the best solution
517 if ( fasp_dvec_isnan(u) ) {
518 absres = BIGREAL;
519 goto RESTORE_BESTSOL;
520 }
521
522 // safety net check: save the best-so-far solution
523 if ( absres < absres_best - maxdiff) {
524 absres_best = absres;
525 iter_best = iter;
526 fasp_darray_cp(m,u->val,u_best);
527 }
528
529 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
530 normuinf = fasp_blas_darray_norminf(m, u->val);
531 if ( normuinf <= sol_inf_tol ) {
532 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
533 iter = ERROR_SOLVER_SOLSTAG;
534 break;
535 }
536
537 // Check II: if staggenated, try to restart
538 normu = fasp_blas_dvec_norm2(u);
539
540 // compute relative difference
541 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
542 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
543
544 if ( PrtLvl >= PRINT_MORE ) {
545 ITS_DIFFRES(reldiff,relres);
546 ITS_RESTART;
547 }
548
549 fasp_darray_cp(m,b->val,r);
550 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
551
552 // compute residuals
553 switch ( StopType ) {
554 case STOP_REL_RES:
555 absres = fasp_blas_darray_norm2(m,r);
556 relres = absres/normr0;
557 break;
558 case STOP_REL_PRECRES:
559 // z = B(r)
560 if ( pc != NULL )
561 pc->fct(r,z,pc->data); /* Apply preconditioner */
562 else
563 fasp_darray_cp(m,r,z); /* No preconditioner */
564 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
565 relres = absres/normr0;
566 break;
567 case STOP_MOD_REL_RES:
568 absres = fasp_blas_darray_norm2(m,r);
569 relres = absres/normu;
570 break;
571 }
572
573 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
574
575 if ( relres < tol )
576 break;
577 else {
578 if ( stag >= MaxStag ) {
579 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
580 iter = ERROR_SOLVER_STAG;
581 break;
582 }
583 fasp_darray_set(m,p,0.0);
584 ++stag;
585 ++restart_step;
586 }
587 } // end of staggnation check!
588
589 // Check III: prevent false convergence
590 if ( relres < tol ) {
591
592 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
593
594 // compute residual r = b - Ax again
595 fasp_darray_cp(m,b->val,r);
596 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
597
598 // compute residuals
599 switch ( StopType ) {
600 case STOP_REL_RES:
601 absres = fasp_blas_darray_norm2(m,r);
602 relres = absres/normr0;
603 break;
604 case STOP_REL_PRECRES:
605 // z = B(r)
606 if ( pc != NULL )
607 pc->fct(r,z,pc->data); /* Apply preconditioner */
608 else
609 fasp_darray_cp(m,r,z); /* No preconditioner */
610 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
611 relres = absres/normr0;
612 break;
613 case STOP_MOD_REL_RES:
614 absres = fasp_blas_darray_norm2(m,r);
615 relres = absres/normu;
616 break;
617 }
618
619 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
620
621 // check convergence
622 if ( relres < tol ) break;
623
624 if ( more_step >= MaxRestartStep ) {
625 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
626 iter = ERROR_SOLVER_TOLSMALL;
627 break;
628 }
629
630 // prepare for restarting the method
631 fasp_darray_set(m,p,0.0);
632 ++more_step;
633 ++restart_step;
634
635 } // end of safe-guard check!
636
637 // save residual for next iteration
638 absres0 = absres;
639
640 // compute z_k = B(r_k)
641 if ( StopType != STOP_REL_PRECRES ) {
642 if ( pc != NULL )
643 pc->fct(r,z,pc->data); /* Apply preconditioner */
644 else
645 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
646 }
647
648 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
649 temp2 = fasp_blas_darray_dotprod(m,z,r);
650 beta = temp2/temp1;
651 temp1 = temp2;
652
653 // compute p_k = z_k + beta_k*p_{k-1}
654 fasp_blas_darray_axpby(m,1.0,z,beta,p);
655
656 } // end of main PCG loop.
657
658RESTORE_BESTSOL: // restore the best-so-far solution if necessary
659 if ( iter != iter_best ) {
660
661 // compute best residual
662 fasp_darray_cp(m,b->val,r);
663 fasp_blas_dblc_aAxpy(-1.0,A,u_best,r);
664
665 switch ( StopType ) {
666 case STOP_REL_RES:
667 absres_best = fasp_blas_darray_norm2(m,r);
668 break;
669 case STOP_REL_PRECRES:
670 // z = B(r)
671 if ( pc != NULL )
672 pc->fct(r,z,pc->data); /* Apply preconditioner */
673 else
674 fasp_darray_cp(m,r,z); /* No preconditioner */
675 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
676 break;
677 case STOP_MOD_REL_RES:
678 absres_best = fasp_blas_darray_norm2(m,r);
679 break;
680 }
681
682 if ( absres > absres_best + maxdiff || isnan(absres) ) {
683 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
684 fasp_darray_cp(m,u_best,u->val);
685 relres = absres_best / normr0;
686 }
687 }
688
689FINISHED: // finish the iterative method
690 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
691
692 // clean up temp memory
693 fasp_mem_free(work); work = NULL;
694
695#if DEBUG_MODE > 0
696 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
697#endif
698
699 if ( iter > MaxIt )
700 return ERROR_SOLVER_MAXIT;
701 else
702 return iter;
703}
704
726INT fasp_solver_dstr_spcg (const dSTRmat *A,
727 const dvector *b,
728 dvector *u,
729 precond *pc,
730 const REAL tol,
731 const INT MaxIt,
732 const SHORT StopType,
733 const SHORT PrtLvl)
734{
735 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
736 const INT m = b->row;
737 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
738 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
739
740 // local variables
741 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
742 REAL absres0 = BIGREAL, absres = BIGREAL;
743 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
744 REAL reldiff, factor, normuinf;
745 REAL alpha, beta, temp1, temp2;
746 INT iter_best = 0; // initial best known iteration
747 REAL absres_best = BIGREAL; // initial best known residual
748
749 // allocate temp memory (need 5*m REAL numbers)
750 REAL *work = (REAL *)fasp_mem_calloc(5*m,sizeof(REAL));
751 REAL *p = work, *z = work+m, *r = z+m, *t = r+m, *u_best = t+m;
752
753 // Output some info for debuging
754 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe CG solver (STR) ...\n");
755
756#if DEBUG_MODE > 0
757 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
758 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
759#endif
760
761 // r = b-A*u
762 fasp_darray_cp(m,b->val,r);
763 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
764
765 if (pc != NULL)
766 pc->fct(r,z,pc->data); /* Apply preconditioner */
767 else
768 fasp_darray_cp(m,r,z); /* No preconditioner */
769
770 // compute initial residuals
771 switch (StopType) {
772 case STOP_REL_RES:
773 absres0 = fasp_blas_darray_norm2(m,r);
774 normr0 = MAX(SMALLREAL,absres0);
775 relres = absres0/normr0;
776 break;
777 case STOP_REL_PRECRES:
778 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
779 normr0 = MAX(SMALLREAL,absres0);
780 relres = absres0/normr0;
781 break;
782 case STOP_MOD_REL_RES:
783 absres0 = fasp_blas_darray_norm2(m,r);
784 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
785 relres = absres0/normu;
786 break;
787 default:
788 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
789 goto FINISHED;
790 }
791
792 // if initial residual is small, no need to iterate!
793 if ( relres < tol ) goto FINISHED;
794
795 // output iteration information if needed
796 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
797
798 fasp_darray_cp(m,z,p);
799 temp1 = fasp_blas_darray_dotprod(m,z,r);
800
801 // main PCG loop
802 while ( iter++ < MaxIt ) {
803
804 // t=A*p
805 fasp_blas_dstr_mxv(A,p,t);
806
807 // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
808 temp2 = fasp_blas_darray_dotprod(m,t,p);
809 if ( ABS(temp2) > SMALLREAL2 ) {
810 alpha = temp1/temp2;
811 }
812 else { // Possible breakdown
813 goto RESTORE_BESTSOL;
814 }
815
816 // u_k=u_{k-1} + alpha_k*p_{k-1}
817 fasp_blas_darray_axpy(m,alpha,p,u->val);
818
819 // r_k=r_{k-1} - alpha_k*A*p_{k-1}
820 fasp_blas_darray_axpy(m,-alpha,t,r);
821
822 // compute residuals
823 switch ( StopType ) {
824 case STOP_REL_RES:
825 absres = fasp_blas_darray_norm2(m,r);
826 relres = absres/normr0;
827 break;
828 case STOP_REL_PRECRES:
829 // z = B(r)
830 if ( pc != NULL )
831 pc->fct(r,z,pc->data); /* Apply preconditioner */
832 else
833 fasp_darray_cp(m,r,z); /* No preconditioner */
834 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
835 relres = absres/normr0;
836 break;
837 case STOP_MOD_REL_RES:
838 absres = fasp_blas_darray_norm2(m,r);
839 relres = absres/normu;
840 break;
841 }
842
843 // compute reducation factor of residual ||r||
844 factor = absres/absres0;
845
846 // output iteration information if needed
847 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
848
849 // if the solution is NAN, restore the best solution
850 if ( fasp_dvec_isnan(u) ) {
851 absres = BIGREAL;
852 goto RESTORE_BESTSOL;
853 }
854
855 // safety net check: save the best-so-far solution
856 if ( absres < absres_best - maxdiff) {
857 absres_best = absres;
858 iter_best = iter;
859 fasp_darray_cp(m,u->val,u_best);
860 }
861
862 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
863 normuinf = fasp_blas_darray_norminf(m, u->val);
864 if ( normuinf <= sol_inf_tol ) {
865 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
866 iter = ERROR_SOLVER_SOLSTAG;
867 break;
868 }
869
870 // Check II: if staggenated, try to restart
871 normu = fasp_blas_dvec_norm2(u);
872
873 // compute relative difference
874 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
875 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
876
877 if ( PrtLvl >= PRINT_MORE ) {
878 ITS_DIFFRES(reldiff,relres);
879 ITS_RESTART;
880 }
881
882 fasp_darray_cp(m,b->val,r);
883 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
884
885 // compute residuals
886 switch ( StopType ) {
887 case STOP_REL_RES:
888 absres = fasp_blas_darray_norm2(m,r);
889 relres = absres/normr0;
890 break;
891 case STOP_REL_PRECRES:
892 // z = B(r)
893 if ( pc != NULL )
894 pc->fct(r,z,pc->data); /* Apply preconditioner */
895 else
896 fasp_darray_cp(m,r,z); /* No preconditioner */
897 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
898 relres = absres/normr0;
899 break;
900 case STOP_MOD_REL_RES:
901 absres = fasp_blas_darray_norm2(m,r);
902 relres = absres/normu;
903 break;
904 }
905
906 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
907
908 if ( relres < tol )
909 break;
910 else {
911 if ( stag >= MaxStag ) {
912 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
913 iter = ERROR_SOLVER_STAG;
914 break;
915 }
916 fasp_darray_set(m,p,0.0);
917 ++stag;
918 ++restart_step;
919 }
920 } // end of staggnation check!
921
922 // Check III: prevent false convergence
923 if ( relres < tol ) {
924
925 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
926
927 // compute residual r = b - Ax again
928 fasp_darray_cp(m,b->val,r);
929 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
930
931 // compute residuals
932 switch ( StopType ) {
933 case STOP_REL_RES:
934 absres = fasp_blas_darray_norm2(m,r);
935 relres = absres/normr0;
936 break;
937 case STOP_REL_PRECRES:
938 // z = B(r)
939 if ( pc != NULL )
940 pc->fct(r,z,pc->data); /* Apply preconditioner */
941 else
942 fasp_darray_cp(m,r,z); /* No preconditioner */
943 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
944 relres = absres/normr0;
945 break;
946 case STOP_MOD_REL_RES:
947 absres = fasp_blas_darray_norm2(m,r);
948 relres = absres/normu;
949 break;
950 }
951
952 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
953
954 // check convergence
955 if ( relres < tol ) break;
956
957 if ( more_step >= MaxRestartStep ) {
958 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
959 iter = ERROR_SOLVER_TOLSMALL;
960 break;
961 }
962
963 // prepare for restarting the method
964 fasp_darray_set(m,p,0.0);
965 ++more_step;
966 ++restart_step;
967
968 } // end of safe-guard check!
969
970 // save residual for next iteration
971 absres0 = absres;
972
973 // compute z_k = B(r_k)
974 if ( StopType != STOP_REL_PRECRES ) {
975 if ( pc != NULL )
976 pc->fct(r,z,pc->data); /* Apply preconditioner */
977 else
978 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
979 }
980
981 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
982 temp2 = fasp_blas_darray_dotprod(m,z,r);
983 beta = temp2/temp1;
984 temp1 = temp2;
985
986 // compute p_k = z_k + beta_k*p_{k-1}
987 fasp_blas_darray_axpby(m,1.0,z,beta,p);
988
989 } // end of main PCG loop.
990
991RESTORE_BESTSOL: // restore the best-so-far solution if necessary
992 if ( iter != iter_best ) {
993
994 // compute best residual
995 fasp_darray_cp(m,b->val,r);
996 fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
997
998 switch ( StopType ) {
999 case STOP_REL_RES:
1000 absres_best = fasp_blas_darray_norm2(m,r);
1001 break;
1002 case STOP_REL_PRECRES:
1003 // z = B(r)
1004 if ( pc != NULL )
1005 pc->fct(r,z,pc->data); /* Apply preconditioner */
1006 else
1007 fasp_darray_cp(m,r,z); /* No preconditioner */
1008 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1009 break;
1010 case STOP_MOD_REL_RES:
1011 absres_best = fasp_blas_darray_norm2(m,r);
1012 break;
1013 }
1014
1015 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1016 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1017 fasp_darray_cp(m,u_best,u->val);
1018 relres = absres_best / normr0;
1019 }
1020 }
1021
1022FINISHED: // finish the iterative method
1023 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1024
1025 // clean up temp memory
1026 fasp_mem_free(work); work = NULL;
1027
1028#if DEBUG_MODE > 0
1029 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1030#endif
1031
1032 if ( iter > MaxIt )
1033 return ERROR_SOLVER_MAXIT;
1034 else
1035 return iter;
1036}
1037
1038/*---------------------------------*/
1039/*-- End of File --*/
1040/*---------------------------------*/
fasp_darray_set
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:41
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_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_norminf
REAL fasp_blas_darray_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: BlaArray.c:719
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_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_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_blas_dvec_norm2
REAL fasp_blas_dvec_norm2(const dvector *x)
L2 norm of dvector x.
Definition: BlaVector.c:170
fasp_solver_dstr_spcg
INT fasp_solver_dstr_spcg(const dSTRmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b with safety net.
Definition: KrySPcg.c:726
fasp_solver_dcsr_spcg
INT fasp_solver_dcsr_spcg(const dCSRmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b with safety net.
Definition: KrySPcg.c:60
fasp_solver_dblc_spcg
INT fasp_solver_dblc_spcg(const dBLCmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b with safety net.
Definition: KrySPcg.c:393
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_STAG
#define ERROR_SOLVER_STAG
Definition: fasp_const.h:43
SMALLREAL2
#define SMALLREAL2
Definition: fasp_const.h:257
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
ERROR_SOLVER_SOLSTAG
#define ERROR_SOLVER_SOLSTAG
Definition: fasp_const.h:44
MAX_STAG
#define MAX_STAG
Definition: fasp_const.h:264
STOP_REL_PRECRES
#define STOP_REL_PRECRES
Definition: fasp_const.h:133
BIGREAL
#define BIGREAL
Some global constants.
Definition: fasp_const.h:255
MAX_RESTART
#define MAX_RESTART
Definition: fasp_const.h:263
PRINT_MORE
#define PRINT_MORE
Definition: fasp_const.h:76
PRINT_NONE
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
ERROR_SOLVER_TOLSMALL
#define ERROR_SOLVER_TOLSMALL
Definition: fasp_const.h:45
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
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