Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KrySPbcgs.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
61INT fasp_solver_dcsr_spbcgs (const dCSRmat *A,
62 const dvector *b,
63 dvector *u,
64 precond *pc,
65 const REAL tol,
66 const INT MaxIt,
67 const SHORT StopType,
68 const SHORT PrtLvl)
69{
70 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
71 const INT m = b->row;
72 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
73 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
74 const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
75
76 // local variables
77 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
78 REAL alpha, beta, omega, temp1, temp2;
79 REAL absres0 = BIGREAL, absres = BIGREAL;
80 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
81 REAL reldiff, factor, normd, tempr, normuinf;
82 REAL *uval = u->val, *bval = b->val;
83 INT iter_best = 0; // initial best known iteration
84 REAL absres_best = BIGREAL; // initial best known residual
85
86 // allocate temp memory (need 8*m REAL)
87 REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
88 REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
89 REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
90
91 // Output some info for debuging
92 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe BiCGstab solver (CSR) ...\n");
93
94#if DEBUG_MODE > 0
95 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
96 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
97#endif
98
99 // r = b-A*u
100 fasp_darray_cp(m,bval,r);
101 fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
102 absres0 = fasp_blas_darray_norm2(m,r);
103
104 // compute initial relative residual
105 switch (StopType) {
106 case STOP_REL_RES:
107 normr0 = MAX(SMALLREAL,absres0);
108 relres = absres0/normr0;
109 break;
110 case STOP_REL_PRECRES:
111 normr0 = MAX(SMALLREAL,absres0);
112 relres = absres0/normr0;
113 break;
114 case STOP_MOD_REL_RES:
115 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,uval));
116 relres = absres0/normu;
117 break;
118 default:
119 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
120 goto FINISHED;
121 }
122
123 // if initial residual is small, no need to iterate!
124 if (relres<tol) goto FINISHED;
125
126 // output iteration information if needed
127 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
128
129 // rho = r* := r
130 fasp_darray_cp(m,r,rho);
131 temp1 = fasp_blas_darray_dotprod(m,r,rho);
132
133 // p = r
134 fasp_darray_cp(m,r,p);
135
136 // main BiCGstab loop
137 while ( iter++ < MaxIt ) {
138
139 // pp = precond(p)
140 if ( pc != NULL )
141 pc->fct(p,pp,pc->data); /* Apply preconditioner */
142 else
143 fasp_darray_cp(m,p,pp); /* No preconditioner */
144
145 // z = A*pp
146 fasp_blas_dcsr_mxv(A,pp,z);
147
148 // alpha = (r,rho)/(A*p,rho)
149 temp2 = fasp_blas_darray_dotprod(m,z,rho);
150 if ( ABS(temp2) > SMALLREAL ) {
151 alpha = temp1/temp2;
152 }
153 else {
154 ITS_DIVZERO; goto FINISHED;
155 }
156
157 // s = r - alpha z
158 fasp_darray_cp(m,r,s);
159 fasp_blas_darray_axpy(m,-alpha,z,s);
160
161 // sp = precond(s)
162 if ( pc != NULL )
163 pc->fct(s,sp,pc->data); /* Apply preconditioner */
164 else
165 fasp_darray_cp(m,s,sp); /* No preconditioner */
166
167 // t = A*sp;
168 fasp_blas_dcsr_mxv(A,sp,t);
169
170 // omega = (t,s)/(t,t)
171 tempr = fasp_blas_darray_dotprod(m,t,t);
172
173 if ( ABS(tempr) > SMALLREAL ) {
174 omega = fasp_blas_darray_dotprod(m,s,t)/tempr;
175 }
176 else {
177 omega = 0.0;
178 if ( PrtLvl >= PRINT_SOME ) ITS_DIVZERO;
179 }
180
181 // delu = alpha pp + omega sp
182 fasp_blas_darray_axpby(m,alpha,pp,omega,sp);
183
184 // u = u + delu
185 fasp_blas_darray_axpy(m,1.0,sp,uval);
186
187 // r = s - omega t
188 fasp_blas_darray_axpy(m,-omega,t,s);
189 fasp_darray_cp(m,s,r);
190
191 // beta = (r,rho)/(rp,rho)
192 temp2 = temp1;
193 temp1 = fasp_blas_darray_dotprod(m,r,rho);
194
195 if ( ABS(temp2) > SMALLREAL ) {
196 beta = (temp1*alpha)/(temp2*omega);
197 }
198 else {
199 ITS_DIVZERO; goto RESTORE_BESTSOL;
200 }
201
202 // p = p - omega z
203 fasp_blas_darray_axpy(m,-omega,z,p);
204
205 // p = r + beta p
206 fasp_blas_darray_axpby(m,1.0,r,beta,p);
207
208 // compute difference
209 normd = fasp_blas_darray_norm2(m,sp);
210 normu = fasp_blas_darray_norm2(m,uval);
211 reldiff = normd/normu;
212
213 if ( normd < TOL_s ) {
214 ITS_SMALLSP; goto FINISHED;
215 }
216
217 // compute residuals
218 switch (StopType) {
219 case STOP_REL_RES:
220 absres = fasp_blas_darray_norm2(m,r);
221 relres = absres/normr0;
222 break;
223 case STOP_REL_PRECRES:
224 if ( pc == NULL )
225 fasp_darray_cp(m,r,z);
226 else
227 pc->fct(r,z,pc->data);
228 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
229 relres = absres/normr0;
230 break;
231 case STOP_MOD_REL_RES:
232 absres = fasp_blas_darray_norm2(m,r);
233 relres = absres/normu;
234 break;
235 }
236
237 // safety net check: save the best-so-far solution
238 if ( fasp_dvec_isnan(u) ) {
239 // If the solution is NAN, restrore the best solution
240 absres = BIGREAL;
241 goto RESTORE_BESTSOL;
242 }
243
244 if ( absres < absres_best - maxdiff) {
245 absres_best = absres;
246 iter_best = iter;
247 fasp_darray_cp(m,uval,u_best);
248 }
249
250 // compute reducation factor of residual ||r||
251 factor = absres/absres0;
252
253 // output iteration information if needed
254 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
255
256 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
257 normuinf = fasp_blas_darray_norminf(m, uval);
258 if ( normuinf <= sol_inf_tol ) {
259 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
260 iter = ERROR_SOLVER_SOLSTAG;
261 goto FINISHED;
262 }
263
264 // Check II: if staggenated, try to restart
265 if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
266
267 if ( PrtLvl >= PRINT_MORE ) {
268 ITS_DIFFRES(reldiff,relres);
269 ITS_RESTART;
270 }
271
272 // re-init iteration param
273 fasp_darray_cp(m,bval,r);
274 fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
275
276 // pp = precond(p)
277 fasp_darray_cp(m,r,p);
278 if ( pc != NULL )
279 pc->fct(p,pp,pc->data); /* Apply preconditioner */
280 else
281 fasp_darray_cp(m,p,pp); /* No preconditioner */
282
283 // rho = r* := r
284 fasp_darray_cp(m,r,rho);
285 temp1 = fasp_blas_darray_dotprod(m,r,rho);
286
287 // compute residuals
288 switch (StopType) {
289 case STOP_REL_RES:
290 absres = fasp_blas_darray_norm2(m,r);
291 relres = absres/normr0;
292 break;
293 case STOP_REL_PRECRES:
294 if ( pc != NULL )
295 pc->fct(r,z,pc->data);
296 else
297 fasp_darray_cp(m,r,z);
298 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
299 relres = absres/normr0;
300 break;
301 case STOP_MOD_REL_RES:
302 absres = fasp_blas_darray_norm2(m,r);
303 relres = absres/normu;
304 break;
305 }
306
307 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
308
309 if ( relres < tol )
310 break;
311 else {
312 if ( stag >= MaxStag ) {
313 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
314 iter = ERROR_SOLVER_STAG;
315 goto FINISHED;
316 }
317 ++stag;
318 ++restart_step;
319 }
320
321 } // end of stagnation check!
322
323 // Check III: prevent false convergence
324 if ( relres < tol ) {
325 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
326
327 // re-init iteration param
328 fasp_darray_cp(m,bval,r);
329 fasp_blas_dcsr_aAxpy(-1.0,A,uval,r);
330
331 // pp = precond(p)
332 fasp_darray_cp(m,r,p);
333 if ( pc != NULL )
334 pc->fct(p,pp,pc->data); /* Apply preconditioner */
335 else
336 fasp_darray_cp(m,p,pp); /* No preconditioner */
337
338 // rho = r* := r
339 fasp_darray_cp(m,r,rho);
340 temp1 = fasp_blas_darray_dotprod(m,r,rho);
341
342 // compute residuals
343 switch (StopType) {
344 case STOP_REL_RES:
345 absres = fasp_blas_darray_norm2(m,r);
346 relres = absres/normr0;
347 break;
348 case STOP_REL_PRECRES:
349 if ( pc != NULL )
350 pc->fct(r,z,pc->data);
351 else
352 fasp_darray_cp(m,r,z);
353 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
354 relres = tempr/normr0;
355 break;
356 case STOP_MOD_REL_RES:
357 absres = fasp_blas_darray_norm2(m,r);
358 relres = absres/normu;
359 break;
360 }
361
362 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
363
364 // check convergence
365 if ( relres < tol ) break;
366
367 if ( more_step >= MaxRestartStep ) {
368 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
369 iter = ERROR_SOLVER_TOLSMALL;
370 goto FINISHED;
371 }
372 else {
373 if ( PrtLvl > PRINT_NONE ) ITS_RESTART;
374 }
375
376 ++more_step;
377 ++restart_step;
378 } // end if safe guard
379
380 absres0 = absres;
381
382 } // end of main BiCGstab loop
383
384RESTORE_BESTSOL: // restore the best-so-far solution if necessary
385 if ( iter != iter_best ) {
386
387 // compute best residual
388 fasp_darray_cp(m,b->val,r);
389 fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
390
391 switch ( StopType ) {
392 case STOP_REL_RES:
393 absres_best = fasp_blas_darray_norm2(m,r);
394 break;
395 case STOP_REL_PRECRES:
396 // z = B(r)
397 if ( pc != NULL )
398 pc->fct(r,z,pc->data); /* Apply preconditioner */
399 else
400 fasp_darray_cp(m,r,z); /* No preconditioner */
401 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
402 break;
403 case STOP_MOD_REL_RES:
404 absres_best = fasp_blas_darray_norm2(m,r);
405 break;
406 }
407
408 if ( absres > absres_best + maxdiff || isnan(absres) ) {
409 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
410 fasp_darray_cp(m,u_best,u->val);
411 relres = absres_best / normr0;
412 }
413 }
414
415FINISHED: // finish the iterative method
416 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
417
418 // clean up temp memory
419 fasp_mem_free(work); work = NULL;
420
421#if DEBUG_MODE > 0
422 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
423#endif
424
425 if ( iter > MaxIt )
426 return ERROR_SOLVER_MAXIT;
427 else
428 return iter;
429}
430
452INT fasp_solver_dbsr_spbcgs (const dBSRmat *A,
453 const dvector *b,
454 dvector *u,
455 precond *pc,
456 const REAL tol,
457 const INT MaxIt,
458 const SHORT StopType,
459 const SHORT PrtLvl)
460{
461 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
462 const INT m = b->row;
463 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
464 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
465 const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
466
467 // local variables
468 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
469 REAL alpha, beta, omega, temp1, temp2;
470 REAL absres0 = BIGREAL, absres = BIGREAL;
471 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
472 REAL reldiff, factor, normd, tempr, normuinf;
473 REAL *uval = u->val, *bval = b->val;
474 INT iter_best = 0; // initial best known iteration
475 REAL absres_best = BIGREAL; // initial best known residual
476
477 // allocate temp memory (need 8*m REAL)
478 REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
479 REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
480 REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
481
482 // Output some info for debuging
483 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe BiCGstab solver (BSR) ...\n");
484
485#if DEBUG_MODE > 0
486 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
487 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
488#endif
489
490 // r = b-A*u
491 fasp_darray_cp(m,bval,r);
492 fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
493 absres0 = fasp_blas_darray_norm2(m,r);
494
495 // compute initial relative residual
496 switch (StopType) {
497 case STOP_REL_RES:
498 normr0 = MAX(SMALLREAL,absres0);
499 relres = absres0/normr0;
500 break;
501 case STOP_REL_PRECRES:
502 normr0 = MAX(SMALLREAL,absres0);
503 relres = absres0/normr0;
504 break;
505 case STOP_MOD_REL_RES:
506 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,uval));
507 relres = absres0/normu;
508 break;
509 default:
510 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
511 goto FINISHED;
512 }
513
514 // if initial residual is small, no need to iterate!
515 if (relres<tol) goto FINISHED;
516
517 // output iteration information if needed
518 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
519
520 // rho = r* := r
521 fasp_darray_cp(m,r,rho);
522 temp1 = fasp_blas_darray_dotprod(m,r,rho);
523
524 // p = r
525 fasp_darray_cp(m,r,p);
526
527 // main BiCGstab loop
528 while ( iter++ < MaxIt ) {
529
530 // pp = precond(p)
531 if ( pc != NULL )
532 pc->fct(p,pp,pc->data); /* Apply preconditioner */
533 else
534 fasp_darray_cp(m,p,pp); /* No preconditioner */
535
536 // z = A*pp
537 fasp_blas_dbsr_mxv(A,pp,z);
538
539 // alpha = (r,rho)/(A*p,rho)
540 temp2 = fasp_blas_darray_dotprod(m,z,rho);
541 if ( ABS(temp2) > SMALLREAL ) {
542 alpha = temp1/temp2;
543 }
544 else {
545 ITS_DIVZERO; goto FINISHED;
546 }
547
548 // s = r - alpha z
549 fasp_darray_cp(m,r,s);
550 fasp_blas_darray_axpy(m,-alpha,z,s);
551
552 // sp = precond(s)
553 if ( pc != NULL )
554 pc->fct(s,sp,pc->data); /* Apply preconditioner */
555 else
556 fasp_darray_cp(m,s,sp); /* No preconditioner */
557
558 // t = A*sp;
559 fasp_blas_dbsr_mxv(A,sp,t);
560
561 // omega = (t,s)/(t,t)
562 tempr = fasp_blas_darray_dotprod(m,t,t);
563
564 if ( ABS(tempr) > SMALLREAL ) {
565 omega = fasp_blas_darray_dotprod(m,s,t)/tempr;
566 }
567 else {
568 omega = 0.0;
569 if ( PrtLvl >= PRINT_SOME ) ITS_DIVZERO;
570 }
571
572 // delu = alpha pp + omega sp
573 fasp_blas_darray_axpby(m,alpha,pp,omega,sp);
574
575 // u = u + delu
576 fasp_blas_darray_axpy(m,1.0,sp,uval);
577
578 // r = s - omega t
579 fasp_blas_darray_axpy(m,-omega,t,s);
580 fasp_darray_cp(m,s,r);
581
582 // beta = (r,rho)/(rp,rho)
583 temp2 = temp1;
584 temp1 = fasp_blas_darray_dotprod(m,r,rho);
585
586 if ( ABS(temp2) > SMALLREAL ) {
587 beta = (temp1*alpha)/(temp2*omega);
588 }
589 else {
590 ITS_DIVZERO; goto RESTORE_BESTSOL;
591 }
592
593 // p = p - omega z
594 fasp_blas_darray_axpy(m,-omega,z,p);
595
596 // p = r + beta p
597 fasp_blas_darray_axpby(m,1.0,r,beta,p);
598
599 // compute difference
600 normd = fasp_blas_darray_norm2(m,sp);
601 normu = fasp_blas_darray_norm2(m,uval);
602 reldiff = normd/normu;
603
604 if ( normd < TOL_s ) {
605 ITS_SMALLSP; goto FINISHED;
606 }
607
608 // compute residuals
609 switch (StopType) {
610 case STOP_REL_RES:
611 absres = fasp_blas_darray_norm2(m,r);
612 relres = absres/normr0;
613 break;
614 case STOP_REL_PRECRES:
615 if ( pc == NULL )
616 fasp_darray_cp(m,r,z);
617 else
618 pc->fct(r,z,pc->data);
619 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
620 relres = absres/normr0;
621 break;
622 case STOP_MOD_REL_RES:
623 absres = fasp_blas_darray_norm2(m,r);
624 relres = absres/normu;
625 break;
626 }
627
628 // safety net check: save the best-so-far solution
629 if ( fasp_dvec_isnan(u) ) {
630 // If the solution is NAN, restrore the best solution
631 absres = BIGREAL;
632 goto RESTORE_BESTSOL;
633 }
634
635 if ( absres < absres_best - maxdiff) {
636 absres_best = absres;
637 iter_best = iter;
638 fasp_darray_cp(m,uval,u_best);
639 }
640
641 // compute reducation factor of residual ||r||
642 factor = absres/absres0;
643
644 // output iteration information if needed
645 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
646
647 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
648 normuinf = fasp_blas_darray_norminf(m, uval);
649 if ( normuinf <= sol_inf_tol ) {
650 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
651 iter = ERROR_SOLVER_SOLSTAG;
652 goto FINISHED;
653 }
654
655 // Check II: if staggenated, try to restart
656 if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
657
658 if ( PrtLvl >= PRINT_MORE ) {
659 ITS_DIFFRES(reldiff,relres);
660 ITS_RESTART;
661 }
662
663 // re-init iteration param
664 fasp_darray_cp(m,bval,r);
665 fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
666
667 // pp = precond(p)
668 fasp_darray_cp(m,r,p);
669 if ( pc != NULL )
670 pc->fct(p,pp,pc->data); /* Apply preconditioner */
671 else
672 fasp_darray_cp(m,p,pp); /* No preconditioner */
673
674 // rho = r* := r
675 fasp_darray_cp(m,r,rho);
676 temp1 = fasp_blas_darray_dotprod(m,r,rho);
677
678 // compute residuals
679 switch (StopType) {
680 case STOP_REL_RES:
681 absres = fasp_blas_darray_norm2(m,r);
682 relres = absres/normr0;
683 break;
684 case STOP_REL_PRECRES:
685 if ( pc != NULL )
686 pc->fct(r,z,pc->data);
687 else
688 fasp_darray_cp(m,r,z);
689 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
690 relres = absres/normr0;
691 break;
692 case STOP_MOD_REL_RES:
693 absres = fasp_blas_darray_norm2(m,r);
694 relres = absres/normu;
695 break;
696 }
697
698 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
699
700 if ( relres < tol )
701 break;
702 else {
703 if ( stag >= MaxStag ) {
704 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
705 iter = ERROR_SOLVER_STAG;
706 goto FINISHED;
707 }
708 ++stag;
709 ++restart_step;
710 }
711
712 } // end of stagnation check!
713
714 // Check III: prevent false convergence
715 if ( relres < tol ) {
716 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
717
718 // re-init iteration param
719 fasp_darray_cp(m,bval,r);
720 fasp_blas_dbsr_aAxpy(-1.0,A,uval,r);
721
722 // pp = precond(p)
723 fasp_darray_cp(m,r,p);
724 if ( pc != NULL )
725 pc->fct(p,pp,pc->data); /* Apply preconditioner */
726 else
727 fasp_darray_cp(m,p,pp); /* No preconditioner */
728
729 // rho = r* := r
730 fasp_darray_cp(m,r,rho);
731 temp1 = fasp_blas_darray_dotprod(m,r,rho);
732
733 // compute residuals
734 switch (StopType) {
735 case STOP_REL_RES:
736 absres = fasp_blas_darray_norm2(m,r);
737 relres = absres/normr0;
738 break;
739 case STOP_REL_PRECRES:
740 if ( pc != NULL )
741 pc->fct(r,z,pc->data);
742 else
743 fasp_darray_cp(m,r,z);
744 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
745 relres = tempr/normr0;
746 break;
747 case STOP_MOD_REL_RES:
748 absres = fasp_blas_darray_norm2(m,r);
749 relres = absres/normu;
750 break;
751 }
752
753 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
754
755 // check convergence
756 if ( relres < tol ) break;
757
758 if ( more_step >= MaxRestartStep ) {
759 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
760 iter = ERROR_SOLVER_TOLSMALL;
761 goto FINISHED;
762 }
763 else {
764 if ( PrtLvl > PRINT_NONE ) ITS_RESTART;
765 }
766
767 ++more_step;
768 ++restart_step;
769 } // end if safe guard
770
771 absres0 = absres;
772
773 } // end of main BiCGstab loop
774
775RESTORE_BESTSOL: // restore the best-so-far solution if necessary
776 if ( iter != iter_best ) {
777
778 // compute best residual
779 fasp_darray_cp(m,b->val,r);
780 fasp_blas_dbsr_aAxpy(-1.0,A,u_best,r);
781
782 switch ( StopType ) {
783 case STOP_REL_RES:
784 absres_best = fasp_blas_darray_norm2(m,r);
785 break;
786 case STOP_REL_PRECRES:
787 // z = B(r)
788 if ( pc != NULL )
789 pc->fct(r,z,pc->data); /* Apply preconditioner */
790 else
791 fasp_darray_cp(m,r,z); /* No preconditioner */
792 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
793 break;
794 case STOP_MOD_REL_RES:
795 absres_best = fasp_blas_darray_norm2(m,r);
796 break;
797 }
798
799 if ( absres > absres_best + maxdiff || isnan(absres) ) {
800 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
801 fasp_darray_cp(m,u_best,u->val);
802 relres = absres_best / normr0;
803 }
804 }
805
806FINISHED: // finish the iterative method
807 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
808
809 // clean up temp memory
810 fasp_mem_free(work); work = NULL;
811
812#if DEBUG_MODE > 0
813 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
814#endif
815
816 if ( iter > MaxIt )
817 return ERROR_SOLVER_MAXIT;
818 else
819 return iter;
820}
821
843INT fasp_solver_dblc_spbcgs (const dBLCmat *A,
844 const dvector *b,
845 dvector *u,
846 precond *pc,
847 const REAL tol,
848 const INT MaxIt,
849 const SHORT StopType,
850 const SHORT PrtLvl)
851{
852 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
853 const INT m = b->row;
854 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
855 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
856 const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
857
858 // local variables
859 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
860 REAL alpha, beta, omega, temp1, temp2;
861 REAL absres0 = BIGREAL, absres = BIGREAL;
862 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
863 REAL reldiff, factor, normd, tempr, normuinf;
864 REAL *uval = u->val, *bval = b->val;
865 INT iter_best = 0; // initial best known iteration
866 REAL absres_best = BIGREAL; // initial best known residual
867
868 // allocate temp memory (need 8*m REAL)
869 REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
870 REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
871 REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
872
873 // Output some info for debuging
874 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe BiCGstab solver (BLC) ...\n");
875
876#if DEBUG_MODE > 0
877 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
878 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
879#endif
880
881 // r = b-A*u
882 fasp_darray_cp(m,bval,r);
883 fasp_blas_dblc_aAxpy(-1.0,A,uval,r);
884 absres0 = fasp_blas_darray_norm2(m,r);
885
886 // compute initial relative residual
887 switch (StopType) {
888 case STOP_REL_RES:
889 normr0 = MAX(SMALLREAL,absres0);
890 relres = absres0/normr0;
891 break;
892 case STOP_REL_PRECRES:
893 normr0 = MAX(SMALLREAL,absres0);
894 relres = absres0/normr0;
895 break;
896 case STOP_MOD_REL_RES:
897 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,uval));
898 relres = absres0/normu;
899 break;
900 default:
901 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
902 goto FINISHED;
903 }
904
905 // if initial residual is small, no need to iterate!
906 if (relres<tol) goto FINISHED;
907
908 // output iteration information if needed
909 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
910
911 // rho = r* := r
912 fasp_darray_cp(m,r,rho);
913 temp1 = fasp_blas_darray_dotprod(m,r,rho);
914
915 // p = r
916 fasp_darray_cp(m,r,p);
917
918 // main BiCGstab loop
919 while ( iter++ < MaxIt ) {
920
921 // pp = precond(p)
922 if ( pc != NULL )
923 pc->fct(p,pp,pc->data); /* Apply preconditioner */
924 else
925 fasp_darray_cp(m,p,pp); /* No preconditioner */
926
927 // z = A*pp
928 fasp_blas_dblc_mxv(A,pp,z);
929
930 // alpha = (r,rho)/(A*p,rho)
931 temp2 = fasp_blas_darray_dotprod(m,z,rho);
932 if ( ABS(temp2) > SMALLREAL ) {
933 alpha = temp1/temp2;
934 }
935 else {
936 ITS_DIVZERO; goto FINISHED;
937 }
938
939 // s = r - alpha z
940 fasp_darray_cp(m,r,s);
941 fasp_blas_darray_axpy(m,-alpha,z,s);
942
943 // sp = precond(s)
944 if ( pc != NULL )
945 pc->fct(s,sp,pc->data); /* Apply preconditioner */
946 else
947 fasp_darray_cp(m,s,sp); /* No preconditioner */
948
949 // t = A*sp;
950 fasp_blas_dblc_mxv(A,sp,t);
951
952 // omega = (t,s)/(t,t)
953 tempr = fasp_blas_darray_dotprod(m,t,t);
954
955 if ( ABS(tempr) > SMALLREAL ) {
956 omega = fasp_blas_darray_dotprod(m,s,t)/tempr;
957 }
958 else {
959 omega = 0.0;
960 if ( PrtLvl >= PRINT_SOME ) ITS_DIVZERO;
961 }
962
963 // delu = alpha pp + omega sp
964 fasp_blas_darray_axpby(m,alpha,pp,omega,sp);
965
966 // u = u + delu
967 fasp_blas_darray_axpy(m,1.0,sp,uval);
968
969 // r = s - omega t
970 fasp_blas_darray_axpy(m,-omega,t,s);
971 fasp_darray_cp(m,s,r);
972
973 // beta = (r,rho)/(rp,rho)
974 temp2 = temp1;
975 temp1 = fasp_blas_darray_dotprod(m,r,rho);
976
977 if ( ABS(temp2) > SMALLREAL ) {
978 beta = (temp1*alpha)/(temp2*omega);
979 }
980 else {
981 ITS_DIVZERO; goto RESTORE_BESTSOL;
982 }
983
984 // p = p - omega z
985 fasp_blas_darray_axpy(m,-omega,z,p);
986
987 // p = r + beta p
988 fasp_blas_darray_axpby(m,1.0,r,beta,p);
989
990 // compute difference
991 normd = fasp_blas_darray_norm2(m,sp);
992 normu = fasp_blas_darray_norm2(m,uval);
993 reldiff = normd/normu;
994
995 if ( normd < TOL_s ) {
996 ITS_SMALLSP; goto FINISHED;
997 }
998
999 // compute residuals
1000 switch (StopType) {
1001 case STOP_REL_RES:
1002 absres = fasp_blas_darray_norm2(m,r);
1003 relres = absres/normr0;
1004 break;
1005 case STOP_REL_PRECRES:
1006 if ( pc == NULL )
1007 fasp_darray_cp(m,r,z);
1008 else
1009 pc->fct(r,z,pc->data);
1010 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1011 relres = absres/normr0;
1012 break;
1013 case STOP_MOD_REL_RES:
1014 absres = fasp_blas_darray_norm2(m,r);
1015 relres = absres/normu;
1016 break;
1017 }
1018
1019 // safety net check: save the best-so-far solution
1020 if ( fasp_dvec_isnan(u) ) {
1021 // If the solution is NAN, restrore the best solution
1022 absres = BIGREAL;
1023 goto RESTORE_BESTSOL;
1024 }
1025
1026 if ( absres < absres_best - maxdiff) {
1027 absres_best = absres;
1028 iter_best = iter;
1029 fasp_darray_cp(m,uval,u_best);
1030 }
1031
1032 // compute reducation factor of residual ||r||
1033 factor = absres/absres0;
1034
1035 // output iteration information if needed
1036 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1037
1038 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1039 normuinf = fasp_blas_darray_norminf(m, uval);
1040 if ( normuinf <= sol_inf_tol ) {
1041 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
1042 iter = ERROR_SOLVER_SOLSTAG;
1043 goto FINISHED;
1044 }
1045
1046 // Check II: if staggenated, try to restart
1047 if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
1048
1049 if ( PrtLvl >= PRINT_MORE ) {
1050 ITS_DIFFRES(reldiff,relres);
1051 ITS_RESTART;
1052 }
1053
1054 // re-init iteration param
1055 fasp_darray_cp(m,bval,r);
1056 fasp_blas_dblc_aAxpy(-1.0,A,uval,r);
1057
1058 // pp = precond(p)
1059 fasp_darray_cp(m,r,p);
1060 if ( pc != NULL )
1061 pc->fct(p,pp,pc->data); /* Apply preconditioner */
1062 else
1063 fasp_darray_cp(m,p,pp); /* No preconditioner */
1064
1065 // rho = r* := r
1066 fasp_darray_cp(m,r,rho);
1067 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1068
1069 // compute residuals
1070 switch (StopType) {
1071 case STOP_REL_RES:
1072 absres = fasp_blas_darray_norm2(m,r);
1073 relres = absres/normr0;
1074 break;
1075 case STOP_REL_PRECRES:
1076 if ( pc != NULL )
1077 pc->fct(r,z,pc->data);
1078 else
1079 fasp_darray_cp(m,r,z);
1080 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1081 relres = absres/normr0;
1082 break;
1083 case STOP_MOD_REL_RES:
1084 absres = fasp_blas_darray_norm2(m,r);
1085 relres = absres/normu;
1086 break;
1087 }
1088
1089 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1090
1091 if ( relres < tol )
1092 break;
1093 else {
1094 if ( stag >= MaxStag ) {
1095 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
1096 iter = ERROR_SOLVER_STAG;
1097 goto FINISHED;
1098 }
1099 ++stag;
1100 ++restart_step;
1101 }
1102
1103 } // end of stagnation check!
1104
1105 // Check III: prevent false convergence
1106 if ( relres < tol ) {
1107 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
1108
1109 // re-init iteration param
1110 fasp_darray_cp(m,bval,r);
1111 fasp_blas_dblc_aAxpy(-1.0,A,uval,r);
1112
1113 // pp = precond(p)
1114 fasp_darray_cp(m,r,p);
1115 if ( pc != NULL )
1116 pc->fct(p,pp,pc->data); /* Apply preconditioner */
1117 else
1118 fasp_darray_cp(m,p,pp); /* No preconditioner */
1119
1120 // rho = r* := r
1121 fasp_darray_cp(m,r,rho);
1122 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1123
1124 // compute residuals
1125 switch (StopType) {
1126 case STOP_REL_RES:
1127 absres = fasp_blas_darray_norm2(m,r);
1128 relres = absres/normr0;
1129 break;
1130 case STOP_REL_PRECRES:
1131 if ( pc != NULL )
1132 pc->fct(r,z,pc->data);
1133 else
1134 fasp_darray_cp(m,r,z);
1135 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1136 relres = tempr/normr0;
1137 break;
1138 case STOP_MOD_REL_RES:
1139 absres = fasp_blas_darray_norm2(m,r);
1140 relres = absres/normu;
1141 break;
1142 }
1143
1144 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1145
1146 // check convergence
1147 if ( relres < tol ) break;
1148
1149 if ( more_step >= MaxRestartStep ) {
1150 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
1151 iter = ERROR_SOLVER_TOLSMALL;
1152 goto FINISHED;
1153 }
1154 else {
1155 if ( PrtLvl > PRINT_NONE ) ITS_RESTART;
1156 }
1157
1158 ++more_step;
1159 ++restart_step;
1160 } // end if safe guard
1161
1162 absres0 = absres;
1163
1164 } // end of main BiCGstab loop
1165
1166RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1167 if ( iter != iter_best ) {
1168
1169 // compute best residual
1170 fasp_darray_cp(m,b->val,r);
1171 fasp_blas_dblc_aAxpy(-1.0,A,u_best,r);
1172
1173 switch ( StopType ) {
1174 case STOP_REL_RES:
1175 absres_best = fasp_blas_darray_norm2(m,r);
1176 break;
1177 case STOP_REL_PRECRES:
1178 // z = B(r)
1179 if ( pc != NULL )
1180 pc->fct(r,z,pc->data); /* Apply preconditioner */
1181 else
1182 fasp_darray_cp(m,r,z); /* No preconditioner */
1183 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1184 break;
1185 case STOP_MOD_REL_RES:
1186 absres_best = fasp_blas_darray_norm2(m,r);
1187 break;
1188 }
1189
1190 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1191 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1192 fasp_darray_cp(m,u_best,u->val);
1193 relres = absres_best / normr0;
1194 }
1195 }
1196
1197FINISHED: // finish the iterative method
1198 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1199
1200 // clean up temp memory
1201 fasp_mem_free(work); work = NULL;
1202
1203#if DEBUG_MODE > 0
1204 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1205#endif
1206
1207 if ( iter > MaxIt )
1208 return ERROR_SOLVER_MAXIT;
1209 else
1210 return iter;
1211}
1212
1234INT fasp_solver_dstr_spbcgs (const dSTRmat *A,
1235 const dvector *b,
1236 dvector *u,
1237 precond *pc,
1238 const REAL tol,
1239 const INT MaxIt,
1240 const SHORT StopType,
1241 const SHORT PrtLvl)
1242{
1243 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
1244 const INT m = b->row;
1245 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1246 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
1247 const REAL TOL_s = tol*1e-2; // tolerance for norm(p)
1248
1249 // local variables
1250 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
1251 REAL alpha, beta, omega, temp1, temp2;
1252 REAL absres0 = BIGREAL, absres = BIGREAL;
1253 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
1254 REAL reldiff, factor, normd, tempr, normuinf;
1255 REAL *uval = u->val, *bval = b->val;
1256 INT iter_best = 0; // initial best known iteration
1257 REAL absres_best = BIGREAL; // initial best known residual
1258
1259 // allocate temp memory (need 8*m REAL)
1260 REAL *work = (REAL *)fasp_mem_calloc(9*m,sizeof(REAL));
1261 REAL *p = work, *z = work + m, *r = z + m, *t = r + m;
1262 REAL *rho = t + m, *pp = rho + m, *s = pp + m, *sp = s + m, *u_best = sp + m;
1263
1264 // Output some info for debuging
1265 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe BiCGstab solver (STR) ...\n");
1266
1267#if DEBUG_MODE > 0
1268 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1269 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1270#endif
1271
1272 // r = b-A*u
1273 fasp_darray_cp(m,bval,r);
1274 fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1275 absres0 = fasp_blas_darray_norm2(m,r);
1276
1277 // compute initial relative residual
1278 switch (StopType) {
1279 case STOP_REL_RES:
1280 normr0 = MAX(SMALLREAL,absres0);
1281 relres = absres0/normr0;
1282 break;
1283 case STOP_REL_PRECRES:
1284 normr0 = MAX(SMALLREAL,absres0);
1285 relres = absres0/normr0;
1286 break;
1287 case STOP_MOD_REL_RES:
1288 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,uval));
1289 relres = absres0/normu;
1290 break;
1291 default:
1292 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1293 goto FINISHED;
1294 }
1295
1296 // if initial residual is small, no need to iterate!
1297 if (relres<tol) goto FINISHED;
1298
1299 // output iteration information if needed
1300 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
1301
1302 // rho = r* := r
1303 fasp_darray_cp(m,r,rho);
1304 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1305
1306 // p = r
1307 fasp_darray_cp(m,r,p);
1308
1309 // main BiCGstab loop
1310 while ( iter++ < MaxIt ) {
1311
1312 // pp = precond(p)
1313 if ( pc != NULL )
1314 pc->fct(p,pp,pc->data); /* Apply preconditioner */
1315 else
1316 fasp_darray_cp(m,p,pp); /* No preconditioner */
1317
1318 // z = A*pp
1319 fasp_blas_dstr_mxv(A,pp,z);
1320
1321 // alpha = (r,rho)/(A*p,rho)
1322 temp2 = fasp_blas_darray_dotprod(m,z,rho);
1323 if ( ABS(temp2) > SMALLREAL ) {
1324 alpha = temp1/temp2;
1325 }
1326 else {
1327 ITS_DIVZERO; goto FINISHED;
1328 }
1329
1330 // s = r - alpha z
1331 fasp_darray_cp(m,r,s);
1332 fasp_blas_darray_axpy(m,-alpha,z,s);
1333
1334 // sp = precond(s)
1335 if ( pc != NULL )
1336 pc->fct(s,sp,pc->data); /* Apply preconditioner */
1337 else
1338 fasp_darray_cp(m,s,sp); /* No preconditioner */
1339
1340 // t = A*sp;
1341 fasp_blas_dstr_mxv(A,sp,t);
1342
1343 // omega = (t,s)/(t,t)
1344 tempr = fasp_blas_darray_dotprod(m,t,t);
1345
1346 if ( ABS(tempr) > SMALLREAL ) {
1347 omega = fasp_blas_darray_dotprod(m,s,t)/tempr;
1348 }
1349 else {
1350 omega = 0.0;
1351 if ( PrtLvl >= PRINT_SOME ) ITS_DIVZERO;
1352 }
1353
1354 // delu = alpha pp + omega sp
1355 fasp_blas_darray_axpby(m,alpha,pp,omega,sp);
1356
1357 // u = u + delu
1358 fasp_blas_darray_axpy(m,1.0,sp,uval);
1359
1360 // r = s - omega t
1361 fasp_blas_darray_axpy(m,-omega,t,s);
1362 fasp_darray_cp(m,s,r);
1363
1364 // beta = (r,rho)/(rp,rho)
1365 temp2 = temp1;
1366 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1367
1368 if ( ABS(temp2) > SMALLREAL ) {
1369 beta = (temp1*alpha)/(temp2*omega);
1370 }
1371 else {
1372 ITS_DIVZERO; goto RESTORE_BESTSOL;
1373 }
1374
1375 // p = p - omega z
1376 fasp_blas_darray_axpy(m,-omega,z,p);
1377
1378 // p = r + beta p
1379 fasp_blas_darray_axpby(m,1.0,r,beta,p);
1380
1381 // compute difference
1382 normd = fasp_blas_darray_norm2(m,sp);
1383 normu = fasp_blas_darray_norm2(m,uval);
1384 reldiff = normd/normu;
1385
1386 if ( normd < TOL_s ) {
1387 ITS_SMALLSP; goto FINISHED;
1388 }
1389
1390 // compute residuals
1391 switch (StopType) {
1392 case STOP_REL_RES:
1393 absres = fasp_blas_darray_norm2(m,r);
1394 relres = absres/normr0;
1395 break;
1396 case STOP_REL_PRECRES:
1397 if ( pc == NULL )
1398 fasp_darray_cp(m,r,z);
1399 else
1400 pc->fct(r,z,pc->data);
1401 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1402 relres = absres/normr0;
1403 break;
1404 case STOP_MOD_REL_RES:
1405 absres = fasp_blas_darray_norm2(m,r);
1406 relres = absres/normu;
1407 break;
1408 }
1409
1410 // safety net check: save the best-so-far solution
1411 if ( fasp_dvec_isnan(u) ) {
1412 // If the solution is NAN, restrore the best solution
1413 absres = BIGREAL;
1414 goto RESTORE_BESTSOL;
1415 }
1416
1417 if ( absres < absres_best - maxdiff) {
1418 absres_best = absres;
1419 iter_best = iter;
1420 fasp_darray_cp(m,uval,u_best);
1421 }
1422
1423 // compute reducation factor of residual ||r||
1424 factor = absres/absres0;
1425
1426 // output iteration information if needed
1427 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1428
1429 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1430 normuinf = fasp_blas_darray_norminf(m, uval);
1431 if ( normuinf <= sol_inf_tol ) {
1432 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
1433 iter = ERROR_SOLVER_SOLSTAG;
1434 goto FINISHED;
1435 }
1436
1437 // Check II: if staggenated, try to restart
1438 if ( (stag <= MaxStag) && (reldiff < maxdiff) ) {
1439
1440 if ( PrtLvl >= PRINT_MORE ) {
1441 ITS_DIFFRES(reldiff,relres);
1442 ITS_RESTART;
1443 }
1444
1445 // re-init iteration param
1446 fasp_darray_cp(m,bval,r);
1447 fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1448
1449 // pp = precond(p)
1450 fasp_darray_cp(m,r,p);
1451 if ( pc != NULL )
1452 pc->fct(p,pp,pc->data); /* Apply preconditioner */
1453 else
1454 fasp_darray_cp(m,p,pp); /* No preconditioner */
1455
1456 // rho = r* := r
1457 fasp_darray_cp(m,r,rho);
1458 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1459
1460 // compute residuals
1461 switch (StopType) {
1462 case STOP_REL_RES:
1463 absres = fasp_blas_darray_norm2(m,r);
1464 relres = absres/normr0;
1465 break;
1466 case STOP_REL_PRECRES:
1467 if ( pc != NULL )
1468 pc->fct(r,z,pc->data);
1469 else
1470 fasp_darray_cp(m,r,z);
1471 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1472 relres = absres/normr0;
1473 break;
1474 case STOP_MOD_REL_RES:
1475 absres = fasp_blas_darray_norm2(m,r);
1476 relres = absres/normu;
1477 break;
1478 }
1479
1480 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1481
1482 if ( relres < tol )
1483 break;
1484 else {
1485 if ( stag >= MaxStag ) {
1486 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
1487 iter = ERROR_SOLVER_STAG;
1488 goto FINISHED;
1489 }
1490 ++stag;
1491 ++restart_step;
1492 }
1493
1494 } // end of stagnation check!
1495
1496 // Check III: prevent false convergence
1497 if ( relres < tol ) {
1498 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
1499
1500 // re-init iteration param
1501 fasp_darray_cp(m,bval,r);
1502 fasp_blas_dstr_aAxpy(-1.0,A,uval,r);
1503
1504 // pp = precond(p)
1505 fasp_darray_cp(m,r,p);
1506 if ( pc != NULL )
1507 pc->fct(p,pp,pc->data); /* Apply preconditioner */
1508 else
1509 fasp_darray_cp(m,p,pp); /* No preconditioner */
1510
1511 // rho = r* := r
1512 fasp_darray_cp(m,r,rho);
1513 temp1 = fasp_blas_darray_dotprod(m,r,rho);
1514
1515 // compute residuals
1516 switch (StopType) {
1517 case STOP_REL_RES:
1518 absres = fasp_blas_darray_norm2(m,r);
1519 relres = absres/normr0;
1520 break;
1521 case STOP_REL_PRECRES:
1522 if ( pc != NULL )
1523 pc->fct(r,z,pc->data);
1524 else
1525 fasp_darray_cp(m,r,z);
1526 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,r,z)));
1527 relres = tempr/normr0;
1528 break;
1529 case STOP_MOD_REL_RES:
1530 absres = fasp_blas_darray_norm2(m,r);
1531 relres = absres/normu;
1532 break;
1533 }
1534
1535 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1536
1537 // check convergence
1538 if ( relres < tol ) break;
1539
1540 if ( more_step >= MaxRestartStep ) {
1541 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
1542 iter = ERROR_SOLVER_TOLSMALL;
1543 goto FINISHED;
1544 }
1545 else {
1546 if ( PrtLvl > PRINT_NONE ) ITS_RESTART;
1547 }
1548
1549 ++more_step;
1550 ++restart_step;
1551 } // end if safe guard
1552
1553 absres0 = absres;
1554
1555 } // end of main BiCGstab loop
1556
1557RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1558 if ( iter != iter_best ) {
1559
1560 // compute best residual
1561 fasp_darray_cp(m,b->val,r);
1562 fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
1563
1564 switch ( StopType ) {
1565 case STOP_REL_RES:
1566 absres_best = fasp_blas_darray_norm2(m,r);
1567 break;
1568 case STOP_REL_PRECRES:
1569 // z = B(r)
1570 if ( pc != NULL )
1571 pc->fct(r,z,pc->data); /* Apply preconditioner */
1572 else
1573 fasp_darray_cp(m,r,z); /* No preconditioner */
1574 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1575 break;
1576 case STOP_MOD_REL_RES:
1577 absres_best = fasp_blas_darray_norm2(m,r);
1578 break;
1579 }
1580
1581 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1582 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1583 fasp_darray_cp(m,u_best,u->val);
1584 relres = absres_best / normr0;
1585 }
1586 }
1587
1588FINISHED: // finish the iterative method
1589 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1590
1591 // clean up temp memory
1592 fasp_mem_free(work); work = NULL;
1593
1594#if DEBUG_MODE > 0
1595 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1596#endif
1597
1598 if ( iter > MaxIt )
1599 return ERROR_SOLVER_MAXIT;
1600 else
1601 return iter;
1602}
1603
1604/*---------------------------------*/
1605/*-- End of File --*/
1606/*---------------------------------*/
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_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_spbcgs
INT fasp_solver_dbsr_spbcgs(const dBSRmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: KrySPbcgs.c:452
fasp_solver_dstr_spbcgs
INT fasp_solver_dstr_spbcgs(const dSTRmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: KrySPbcgs.c:1234
fasp_solver_dcsr_spbcgs
INT fasp_solver_dcsr_spbcgs(const dCSRmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: KrySPbcgs.c:61
fasp_solver_dblc_spbcgs
INT fasp_solver_dblc_spbcgs(const dBLCmat *A, const dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned BiCGstab method for solving Au=b with safety net.
Definition: KrySPbcgs.c:843
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
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
PRINT_SOME
#define PRINT_SOME
Definition: fasp_const.h:75
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
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