Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPbcgs.c
Go to the documentation of this file.
1
25#include <math.h>
26#include <float.h>
27
28#include "fasp.h"
29#include "fasp_functs.h"
30
31/*---------------------------------*/
32/*-- Declare Private Functions --*/
33/*---------------------------------*/
34
35#include "KryUtil.inl"
36
37/*---------------------------------*/
38/*-- Public Functions --*/
39/*---------------------------------*/
40
63 dvector *b,
64 dvector *u,
65 precond *pc,
66 const REAL tol,
67 const INT MaxIt,
68 const SHORT StopType,
69 const SHORT PrtLvl)
70{
71 const INT m = b->row;
72
73 // local variables
74 REAL n2b,tolb;
75 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
76 INT flag, maxstagsteps, half_step=0;
77 REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
78 REAL alpha,beta,omega,rho,rho1,rtv,tt;
79 REAL normr,normr_act,normph,normx,imin;
80 REAL norm_sh,norm_xhalf,normrmin,factor;
81 REAL *x = u->val, *bval=b->val;
82
83 // allocate temp memory (need 10*m REAL)
84 REAL *work=(REAL *)fasp_mem_calloc(10*m,sizeof(REAL));
85 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
86 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
87 REAL *t = sh+m, *xmin = t+m;
88
89 // Output some info for debuging
90 if ( PrtLvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (CSR) ...\n");
91
92#if DEBUG_MODE > 0
93 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
94 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
95#endif
96
97 // r = b-A*u
98 fasp_darray_cp(m,bval,r);
99 n2b = fasp_blas_darray_norm2(m,r);
100
101 flag = 1;
102 fasp_darray_cp(m,x,xmin);
103 imin = 0;
104
105 iter = 0;
106
107 tolb = n2b*tol;
108
109 fasp_blas_dcsr_aAxpy(-1.0, A, x, r);
110 normr = fasp_blas_darray_norm2(m,r);
111 normr_act = normr;
112 relres = normr/n2b;
113
114 // if initial residual is small, no need to iterate!
115 if ( normr <= tolb ) {
116 flag = 0;
117 iter = 0;
118 goto FINISHED;
119 }
120
121 // output iteration information if needed
122 fasp_itinfo(PrtLvl,StopType,iter,relres,n2b,0.0);
123
124 // shadow residual rt = r* := r
125 fasp_darray_cp(m,r,rt);
126 normrmin = normr;
127
128 rho = 1.0;
129 omega = 1.0;
130 stag = 0;
131 alpha = 0.0;
132
133 moresteps = 0;
134 maxmsteps = 10;
135 maxstagsteps = 3;
136
137 // loop over maxit iterations (unless convergence or failure)
138 for (iter=1;iter <= MaxIt;iter++) {
139
140 rho1 = rho;
141 rho = fasp_blas_darray_dotprod(m,rt,r);
142
143 if ((rho ==0.0 )|| (ABS(rho) >= DBL_MAX )) {
144 flag = 4;
145 goto FINISHED;
146 }
147
148 if (iter==1) {
149 fasp_darray_cp(m,r,p);
150 }
151 else {
152 beta = (rho/rho1)*(alpha/omega);
153
154 if ((beta == 0)||( ABS(beta) > DBL_MAX )) {
155 flag = 4;
156 goto FINISHED;
157 }
158
159 // p = r + beta * (p - omega * v);
160 fasp_blas_darray_axpy(m,-omega,v,p); //p=p - omega*v
161 fasp_blas_darray_axpby(m,1.0, r, beta, p); //p = 1.0*r +beta*p
162 }
163
164 // pp = precond(p) ,ph
165 if ( pc != NULL )
166 pc->fct(p,ph,pc->data); /* Apply preconditioner */
167 // if ph all is infinite then exit need add
168 else
169 fasp_darray_cp(m,p,ph); /* No preconditioner */
170
171 // v = A*ph
172 fasp_blas_dcsr_mxv(A,ph,v);
173 rtv = fasp_blas_darray_dotprod(m,rt,v);
174
175 if (( rtv==0.0 )||( ABS(rtv) > DBL_MAX )){
176 flag = 4;
177 goto FINISHED;
178 }
179
180 alpha = rho/rtv;
181
182 if ( ABS(alpha) > DBL_MAX ){
183 flag = 4;
184 ITS_DIVZERO;
185 goto FINISHED;
186 }
187
188 normx = fasp_blas_darray_norm2(m,x);
189 normph = fasp_blas_darray_norm2(m,ph);
190 if (ABS(alpha)*normph < DBL_EPSILON*normx )
191 stag = stag + 1;
192 else
193 stag = 0;
194
195 // xhalf = x + alpha * ph; // form the "half" iterate
196 // s = r - alpha * v; // residual associated with xhalf
197 fasp_blas_darray_axpyz(m, alpha, ph, x , xhalf); // z= ax + y
198 fasp_blas_darray_axpyz(m, -alpha, v, r, s);
199 normr = fasp_blas_darray_norm2(m,s); // normr = norm(s);
200 normr_act = normr;
201
202 // compute reduction factor of residual ||r||
203 absres = normr_act;
204 factor = absres/absres0;
205 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
206
207 // check for convergence
208 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
209 {
210 fasp_darray_cp(m,bval,s);
211 fasp_blas_dcsr_aAxpy(-1.0,A,xhalf,s);
212 normr_act = fasp_blas_darray_norm2(m,s);
213
214 if (normr_act <= tolb) {
215 // x = xhalf;
216 fasp_darray_cp(m,xhalf,x); // x = xhalf;
217 flag = 0;
218 imin = iter - 0.5;
219 half_step++;
220 if ( PrtLvl >= PRINT_MORE )
221 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
222 flag,stag,imin,half_step);
223 goto FINISHED;
224 }
225 else {
226 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
227
228 moresteps = moresteps + 1;
229 if (moresteps >= maxmsteps) {
230 // if ~warned
231 flag = 3;
232 fasp_darray_cp(m,xhalf,x);
233 goto FINISHED;
234 }
235 }
236 }
237
238 if ( stag >= maxstagsteps ) {
239 flag = 3;
240 goto FINISHED;
241 }
242
243 if ( normr_act < normrmin ) // update minimal norm quantities
244 {
245 normrmin = normr_act;
246 fasp_darray_cp(m,xhalf,xmin);
247 imin = iter - 0.5;
248 half_step++;
249 if ( PrtLvl >= PRINT_MORE )
250 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
251 flag,stag,imin,half_step);
252 }
253
254 // sh = precond(s)
255 if ( pc != NULL ) {
256 pc->fct(s,sh,pc->data); /* Apply preconditioner */
257 }
258 else
259 fasp_darray_cp(m,s,sh); /* No preconditioner */
260
261 // t = A*sh;
262 fasp_blas_dcsr_mxv(A,sh,t);
263 // tt = t' * t;
264 tt = fasp_blas_darray_dotprod(m,t,t);
265 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
266 flag = 4;
267 goto FINISHED;
268 }
269
270 // omega = (t' * s) / tt;
271 omega = fasp_blas_darray_dotprod(m,s,t)/tt;
272 if ( ABS(omega) > DBL_MAX ) {
273 flag = 4;
274 goto FINISHED;
275 }
276
277 norm_sh = fasp_blas_darray_norm2(m,sh);
278 norm_xhalf = fasp_blas_darray_norm2(m,xhalf);
279
280 if ( ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
281 stag = stag + 1;
282 else
283 stag = 0;
284
285 fasp_blas_darray_axpyz(m, omega,sh,xhalf, x); // x = xhalf + omega * sh;
286 fasp_blas_darray_axpyz(m, -omega, t, s, r); // r = s - omega * t;
287 normr = fasp_blas_darray_norm2(m,r); // normr = norm(r);
288 normr_act = normr;
289
290 // check for convergence
291 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
292 {
293 fasp_darray_cp(m,bval,r);
294 fasp_blas_dcsr_aAxpy(-1.0,A,x,r);
295 normr_act = fasp_blas_darray_norm2(m,r);
296 if ( normr_act <= tolb ) {
297 flag = 0;
298 goto FINISHED;
299 }
300 else {
301 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
302
303 moresteps = moresteps + 1;
304 if ( moresteps >= maxmsteps ) {
305 flag = 3;
306 goto FINISHED;
307 }
308 }
309 }
310
311 // update minimal norm quantities
312 if ( normr_act < normrmin ) {
313 normrmin = normr_act;
314 fasp_darray_cp(m,x,xmin);
315 imin = iter;
316 }
317
318 if ( stag >= maxstagsteps ) {
319 flag = 3;
320 goto FINISHED;
321 }
322
323 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
324
325 absres0 = absres;
326 } // for iter = 1 : maxit
327
328FINISHED: // finish iterative method
329 // returned solution is first with minimal residual
330 if (flag == 0)
331 relres = normr_act / n2b;
332 else {
333 fasp_darray_cp(m, bval,r);
334 fasp_blas_dcsr_aAxpy(-1.0,A,xmin,r);
335 normr = fasp_blas_darray_norm2(m,r);
336
337 if ( normr <= normr_act) {
338 fasp_darray_cp(m, xmin, x);
339 iter = imin;
340 relres = normr/n2b;
341 }
342 else {
343 relres = normr_act/n2b;
344 }
345 }
346
347 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
348
349 if ( PrtLvl >= PRINT_MORE )
350 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
351 flag,stag,imin,half_step);
352
353 // clean up temp memory
354 fasp_mem_free(work); work = NULL;
355
356#if DEBUG_MODE > 0
357 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
358#endif
359
360 if ( iter > MaxIt )
361 return ERROR_SOLVER_MAXIT;
362 else
363 return iter;
364}
365
388 dvector *b,
389 dvector *u,
390 precond *pc,
391 const REAL tol,
392 const INT MaxIt,
393 const SHORT StopType,
394 const SHORT PrtLvl)
395{
396 const INT m = b->row;
397
398 // local variables
399 REAL n2b,tolb;
400 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
401 INT flag, maxstagsteps, half_step=0;
402 REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
403 REAL alpha,beta,omega,rho,rho1,rtv,tt;
404 REAL normr,normr_act,normph,normx,imin;
405 REAL norm_sh,norm_xhalf,normrmin,factor;
406 REAL *x = u->val, *bval=b->val;
407
408 // allocate temp memory (need 10*m REAL)
409 REAL *work=(REAL *)fasp_mem_calloc(10*m,sizeof(REAL));
410 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
411 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
412 REAL *t = sh+m, *xmin = t+m;
413
414 // Output some info for debuging
415 if ( PrtLvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (BSR) ...\n");
416
417#if DEBUG_MODE > 0
418 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
419 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
420#endif
421
422 // r = b-A*u
423 fasp_darray_cp(m,bval,r);
424 n2b = fasp_blas_darray_norm2(m,r);
425
426 flag = 1;
427 fasp_darray_cp(m,x,xmin);
428 imin = 0;
429
430 iter = 0;
431
432 tolb = n2b*tol;
433
434 fasp_blas_dbsr_aAxpy(-1.0, A, x, r);
435 normr = fasp_blas_darray_norm2(m,r);
436 normr_act = normr;
437 relres = normr/n2b;
438
439 // if initial residual is small, no need to iterate!
440 if ( normr <= tolb ) {
441 flag = 0;
442 iter = 0;
443 goto FINISHED;
444 }
445
446 // output iteration information if needed
447 fasp_itinfo(PrtLvl,StopType,iter,relres,n2b,0.0);
448
449 // shadow residual rt = r* := r
450 fasp_darray_cp(m,r,rt);
451 normrmin = normr;
452
453 rho = 1.0;
454 omega = 1.0;
455 stag = 0;
456 alpha = 0.0;
457
458 moresteps = 0;
459 maxmsteps = 10;
460 maxstagsteps = 3;
461
462 // loop over maxit iterations (unless convergence or failure)
463 for (iter=1;iter <= MaxIt;iter++) {
464
465 rho1 = rho;
466 rho = fasp_blas_darray_dotprod(m,rt,r);
467
468 if ((rho ==0.0 )|| (ABS(rho) >= DBL_MAX )) {
469 flag = 4;
470 goto FINISHED;
471 }
472
473 if (iter==1) {
474 fasp_darray_cp(m,r,p);
475 }
476 else {
477 beta = (rho/rho1)*(alpha/omega);
478
479 if ((beta == 0)||( ABS(beta) > DBL_MAX )) {
480 flag = 4;
481 goto FINISHED;
482 }
483
484 // p = r + beta * (p - omega * v);
485 fasp_blas_darray_axpy(m,-omega,v,p); //p=p - omega*v
486 fasp_blas_darray_axpby(m,1.0, r, beta, p); //p = 1.0*r +beta*p
487 }
488
489 // pp = precond(p) ,ph
490 if ( pc != NULL )
491 pc->fct(p,ph,pc->data); /* Apply preconditioner */
492 // if ph all is infinite then exit need add
493 else
494 fasp_darray_cp(m,p,ph); /* No preconditioner */
495
496 // v = A*ph
497 fasp_blas_dbsr_mxv(A,ph,v);
498 rtv = fasp_blas_darray_dotprod(m,rt,v);
499
500 if (( rtv==0.0 )||( ABS(rtv) > DBL_MAX )) {
501 flag = 4;
502 goto FINISHED;
503 }
504
505 alpha = rho/rtv;
506
507 if ( ABS(alpha) > DBL_MAX ) {
508 flag = 4;
509 ITS_DIVZERO;
510 goto FINISHED;
511 }
512
513 normx = fasp_blas_darray_norm2(m,x);
514 normph = fasp_blas_darray_norm2(m,ph);
515 if (ABS(alpha)*normph < DBL_EPSILON*normx )
516 stag = stag + 1;
517 else
518 stag = 0;
519
520 // xhalf = x + alpha * ph; // form the "half" iterate
521 // s = r - alpha * v; // residual associated with xhalf
522 fasp_blas_darray_axpyz(m, alpha, ph, x , xhalf); // z= ax + y
523 fasp_blas_darray_axpyz(m, -alpha, v, r, s);
524 normr = fasp_blas_darray_norm2(m,s); // normr = norm(s);
525 normr_act = normr;
526
527 // compute reduction factor of residual ||r||
528 absres = normr_act;
529 factor = absres/absres0;
530 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
531
532 // check for convergence
533 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
534 {
535 fasp_darray_cp(m,bval,s);
536 fasp_blas_dbsr_aAxpy(-1.0,A,xhalf,s);
537 normr_act = fasp_blas_darray_norm2(m,s);
538
539 if (normr_act <= tolb) {
540 // x = xhalf;
541 fasp_darray_cp(m,xhalf,x); // x = xhalf;
542 flag = 0;
543 imin = iter - 0.5;
544 half_step++;
545 if ( PrtLvl >= PRINT_MORE )
546 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
547 flag,stag,imin,half_step);
548 goto FINISHED;
549 }
550 else {
551 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
552
553 moresteps = moresteps + 1;
554 if (moresteps >= maxmsteps){
555 // if ~warned
556 flag = 3;
557 fasp_darray_cp(m,xhalf,x);
558 goto FINISHED;
559 }
560 }
561 }
562
563 if ( stag >= maxstagsteps ) {
564 flag = 3;
565 goto FINISHED;
566 }
567
568 if ( normr_act < normrmin ) // update minimal norm quantities
569 {
570 normrmin = normr_act;
571 fasp_darray_cp(m,xhalf,xmin);
572 imin = iter - 0.5;
573 half_step++;
574 if ( PrtLvl >= PRINT_MORE )
575 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
576 flag,stag,imin,half_step);
577 }
578
579 // sh = precond(s)
580 if ( pc != NULL ) {
581 pc->fct(s,sh,pc->data); /* Apply preconditioner */
582 }
583 else
584 fasp_darray_cp(m,s,sh); /* No preconditioner */
585
586 // t = A*sh;
587 fasp_blas_dbsr_mxv(A,sh,t);
588 // tt = t' * t;
589 tt = fasp_blas_darray_dotprod(m,t,t);
590 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
591 flag = 4;
592 goto FINISHED;
593 }
594
595 // omega = (t' * s) / tt;
596 omega = fasp_blas_darray_dotprod(m,s,t)/tt;
597 if ( ABS(omega) > DBL_MAX ) {
598 flag = 4;
599 goto FINISHED;
600 }
601
602 norm_sh = fasp_blas_darray_norm2(m,sh);
603 norm_xhalf = fasp_blas_darray_norm2(m,xhalf);
604
605 if ( ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
606 stag = stag + 1;
607 else
608 stag = 0;
609
610 fasp_blas_darray_axpyz(m, omega,sh,xhalf, x); // x = xhalf + omega * sh;
611 fasp_blas_darray_axpyz(m, -omega, t, s, r); // r = s - omega * t;
612 normr = fasp_blas_darray_norm2(m,r); // normr = norm(r);
613 normr_act = normr;
614
615 // check for convergence
616 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
617 {
618 fasp_darray_cp(m,bval,r);
619 fasp_blas_dbsr_aAxpy(-1.0,A,x,r);
620 normr_act = fasp_blas_darray_norm2(m,r);
621 if ( normr_act <= tolb ) {
622 flag = 0;
623 goto FINISHED;
624 }
625 else {
626 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
627
628 moresteps = moresteps + 1;
629 if ( moresteps >= maxmsteps ) {
630 flag = 3;
631 goto FINISHED;
632 }
633 }
634 }
635
636 // update minimal norm quantities
637 if ( normr_act < normrmin ) {
638 normrmin = normr_act;
639 fasp_darray_cp(m,x,xmin);
640 imin = iter;
641 }
642
643 if ( stag >= maxstagsteps )
644 {
645 flag = 3;
646 goto FINISHED;
647 }
648
649 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
650
651 absres0 = absres;
652 } // for iter = 1 : maxit
653
654FINISHED: // finish iterative method
655 // returned solution is first with minimal residual
656 if (flag == 0)
657 relres = normr_act / n2b;
658 else {
659 fasp_darray_cp(m, bval,r);
660 fasp_blas_dbsr_aAxpy(-1.0,A,xmin,r);
661 normr = fasp_blas_darray_norm2(m,r);
662
663 if ( normr <= normr_act) {
664 fasp_darray_cp(m, xmin,x);
665 iter = imin;
666 relres = normr/n2b;
667 }
668 else {
669 relres = normr_act/n2b;
670 }
671 }
672
673 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
674
675 if ( PrtLvl >= PRINT_MORE )
676 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
677 flag,stag,imin,half_step);
678
679 // clean up temp memory
680 fasp_mem_free(work); work = NULL;
681
682#if DEBUG_MODE > 0
683 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
684#endif
685
686 if ( iter > MaxIt )
687 return ERROR_SOLVER_MAXIT;
688 else
689 return iter;
690}
691
714 dvector *b,
715 dvector *u,
716 precond *pc,
717 const REAL tol,
718 const INT MaxIt,
719 const SHORT StopType,
720 const SHORT PrtLvl)
721{
722 const INT m = b->row;
723
724 // local variables
725 REAL n2b,tolb;
726 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
727 INT flag, maxstagsteps, half_step=0;
728 REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
729 REAL alpha,beta,omega,rho,rho1,rtv,tt;
730 REAL normr,normr_act,normph,normx,imin;
731 REAL norm_sh,norm_xhalf,normrmin,factor;
732 REAL *x = u->val, *bval=b->val;
733
734 // allocate temp memory (need 10*m REAL)
735 REAL *work=(REAL *)fasp_mem_calloc(10*m,sizeof(REAL));
736 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
737 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
738 REAL *t = sh+m, *xmin = t+m;
739
740 // Output some info for debuging
741 if ( PrtLvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (BLC) ...\n");
742
743#if DEBUG_MODE > 0
744 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
745 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
746#endif
747
748 // r = b-A*u
749 fasp_darray_cp(m,bval,r);
750 n2b = fasp_blas_darray_norm2(m,r);
751
752 flag = 1;
753 fasp_darray_cp(m,x,xmin);
754 imin = 0;
755
756 iter = 0;
757
758 tolb = n2b*tol;
759
760 fasp_blas_dblc_aAxpy(-1.0, A, x, r);
761 normr = fasp_blas_darray_norm2(m,r);
762 normr_act = normr;
763 relres = normr/n2b;
764
765 // if initial residual is small, no need to iterate!
766 if ( normr <= tolb ) {
767 flag = 0;
768 iter = 0;
769 goto FINISHED;
770 }
771
772 // output iteration information if needed
773 fasp_itinfo(PrtLvl,StopType,iter,relres,n2b,0.0);
774
775 // shadow residual rt = r* := r
776 fasp_darray_cp(m,r,rt);
777 normrmin = normr;
778
779 rho = 1.0;
780 omega = 1.0;
781 stag = 0;
782 alpha = 0.0;
783
784 moresteps = 0;
785 maxmsteps = 10;
786 maxstagsteps = 3;
787
788 // loop over maxit iterations (unless convergence or failure)
789 for (iter=1;iter <= MaxIt;iter++) {
790
791 rho1 = rho;
792 rho = fasp_blas_darray_dotprod(m,rt,r);
793
794 if ((rho ==0.0 )|| (ABS(rho) >= DBL_MAX )) {
795 flag = 4;
796 goto FINISHED;
797 }
798
799 if (iter==1) {
800 fasp_darray_cp(m,r,p);
801 }
802 else {
803 beta = (rho/rho1)*(alpha/omega);
804
805 if ((beta == 0)||( ABS(beta) > DBL_MAX )) {
806 flag = 4;
807 goto FINISHED;
808 }
809
810 // p = r + beta * (p - omega * v);
811 fasp_blas_darray_axpy(m,-omega,v,p); //p=p - omega*v
812 fasp_blas_darray_axpby(m,1.0, r, beta, p); //p = 1.0*r +beta*p
813 }
814
815 // pp = precond(p) ,ph
816 if ( pc != NULL )
817 pc->fct(p,ph,pc->data); /* Apply preconditioner */
818 // if ph all is infinite then exit need add
819 else
820 fasp_darray_cp(m,p,ph); /* No preconditioner */
821
822 // v = A*ph
823 fasp_blas_dblc_mxv(A,ph,v);
824 rtv = fasp_blas_darray_dotprod(m,rt,v);
825
826 if (( rtv==0.0 )||( ABS(rtv) > DBL_MAX )) {
827 flag = 4;
828 goto FINISHED;
829 }
830
831 alpha = rho/rtv;
832
833 if ( ABS(alpha) > DBL_MAX ) {
834 flag = 4;
835 ITS_DIVZERO;
836 goto FINISHED;
837 }
838
839 normx = fasp_blas_darray_norm2(m,x);
840 normph = fasp_blas_darray_norm2(m,ph);
841 if (ABS(alpha)*normph < DBL_EPSILON*normx )
842 stag = stag + 1;
843 else
844 stag = 0;
845
846 // xhalf = x + alpha * ph; // form the "half" iterate
847 // s = r - alpha * v; // residual associated with xhalf
848 fasp_blas_darray_axpyz(m, alpha, ph, x , xhalf); // z= ax + y
849 fasp_blas_darray_axpyz(m, -alpha, v, r, s);
850 normr = fasp_blas_darray_norm2(m,s); // normr = norm(s);
851 normr_act = normr;
852
853 // compute reduction factor of residual ||r||
854 absres = normr_act;
855 factor = absres/absres0;
856 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
857
858 // check for convergence
859 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
860 {
861 fasp_darray_cp(m,bval,s);
862 fasp_blas_dblc_aAxpy(-1.0,A,xhalf,s);
863 normr_act = fasp_blas_darray_norm2(m,s);
864
865 if (normr_act <= tolb) {
866 // x = xhalf;
867 fasp_darray_cp(m,xhalf,x); // x = xhalf;
868 flag = 0;
869 imin = iter - 0.5;
870 half_step++;
871 if ( PrtLvl >= PRINT_MORE )
872 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
873 flag,stag,imin,half_step);
874 goto FINISHED;
875 }
876 else {
877 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
878
879 moresteps = moresteps + 1;
880 if (moresteps >= maxmsteps) {
881 // if ~warned
882 flag = 3;
883 fasp_darray_cp(m,xhalf,x);
884 goto FINISHED;
885 }
886 }
887 }
888
889 if ( stag >= maxstagsteps ) {
890 flag = 3;
891 goto FINISHED;
892 }
893
894 if ( normr_act < normrmin ) // update minimal norm quantities
895 {
896 normrmin = normr_act;
897 fasp_darray_cp(m,xhalf,xmin);
898 imin = iter - 0.5;
899 half_step++;
900 if ( PrtLvl >= PRINT_MORE )
901 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
902 flag,stag,imin,half_step);
903 }
904
905 // sh = precond(s)
906 if ( pc != NULL ) {
907 pc->fct(s,sh,pc->data); /* Apply preconditioner */
908 }
909 else
910 fasp_darray_cp(m,s,sh); /* No preconditioner */
911
912 // t = A*sh;
913 fasp_blas_dblc_mxv(A,sh,t);
914 // tt = t' * t;
915 tt = fasp_blas_darray_dotprod(m,t,t);
916 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
917 flag = 4;
918 goto FINISHED;
919 }
920
921 // omega = (t' * s) / tt;
922 omega = fasp_blas_darray_dotprod(m,s,t)/tt;
923 if ( ABS(omega) > DBL_MAX ) {
924 flag = 4;
925 goto FINISHED;
926 }
927
928 norm_sh = fasp_blas_darray_norm2(m,sh);
929 norm_xhalf = fasp_blas_darray_norm2(m,xhalf);
930
931 if ( ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
932 stag = stag + 1;
933 else
934 stag = 0;
935
936 fasp_blas_darray_axpyz(m, omega,sh,xhalf, x); // x = xhalf + omega * sh;
937 fasp_blas_darray_axpyz(m, -omega, t, s, r); // r = s - omega * t;
938 normr = fasp_blas_darray_norm2(m,r); // normr = norm(r);
939 normr_act = normr;
940
941 // check for convergence
942 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
943 {
944 fasp_darray_cp(m,bval,r);
945 fasp_blas_dblc_aAxpy(-1.0,A,x,r);
946 normr_act = fasp_blas_darray_norm2(m,r);
947 if ( normr_act <= tolb ) {
948 flag = 0;
949 goto FINISHED;
950 }
951 else {
952 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
953
954 moresteps = moresteps + 1;
955 if ( moresteps >= maxmsteps ) {
956 flag = 3;
957 goto FINISHED;
958 }
959 }
960 }
961
962 // update minimal norm quantities
963 if ( normr_act < normrmin ) {
964 normrmin = normr_act;
965 fasp_darray_cp(m,x,xmin);
966 imin = iter;
967 }
968
969 if ( stag >= maxstagsteps )
970 {
971 flag = 3;
972 goto FINISHED;
973 }
974
975 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
976
977 absres0 = absres;
978 } // for iter = 1 : maxit
979
980FINISHED: // finish iterative method
981 // returned solution is first with minimal residual
982 if (flag == 0)
983 relres = normr_act / n2b;
984 else {
985 fasp_darray_cp(m, bval,r);
986 fasp_blas_dblc_aAxpy(-1.0,A,xmin,r);
987 normr = fasp_blas_darray_norm2(m,r);
988
989 if ( normr <= normr_act) {
990 fasp_darray_cp(m, xmin,x);
991 iter = imin;
992 relres = normr/n2b;
993 }
994 else {
995 relres = normr_act/n2b;
996 }
997 }
998
999 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1000
1001 if ( PrtLvl >= PRINT_MORE )
1002 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1003 flag,stag,imin,half_step);
1004
1005 // clean up temp memory
1006 fasp_mem_free(work); work = NULL;
1007
1008#if DEBUG_MODE > 0
1009 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1010#endif
1011
1012 if ( iter > MaxIt )
1013 return ERROR_SOLVER_MAXIT;
1014 else
1015 return iter;
1016}
1017
1040 dvector *b,
1041 dvector *u,
1042 precond *pc,
1043 const REAL tol,
1044 const INT MaxIt,
1045 const SHORT StopType,
1046 const SHORT PrtLvl)
1047{
1048 const INT m = b->row;
1049
1050 // local variables
1051 REAL n2b,tolb;
1052 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
1053 INT flag, maxstagsteps, half_step=0;
1054 REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
1055 REAL alpha,beta,omega,rho,rho1,rtv,tt;
1056 REAL normr,normr_act,normph,normx,imin;
1057 REAL norm_sh,norm_xhalf,normrmin,factor;
1058 REAL *x = u->val, *bval=b->val;
1059
1060 // allocate temp memory (need 10*m REAL)
1061 REAL *work=(REAL *)fasp_mem_calloc(10*m,sizeof(REAL));
1062 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
1063 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
1064 REAL *t = sh+m, *xmin = t+m;
1065
1066 // Output some info for debuging
1067 if ( PrtLvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (STR) ...\n");
1068
1069#if DEBUG_MODE > 0
1070 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1071 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1072#endif
1073
1074 // r = b-A*u
1075 fasp_darray_cp(m,bval,r);
1076 n2b = fasp_blas_darray_norm2(m,r);
1077
1078 flag = 1;
1079 fasp_darray_cp(m,x,xmin);
1080 imin = 0;
1081
1082 iter = 0;
1083
1084 tolb = n2b*tol;
1085
1086 fasp_blas_dstr_aAxpy(-1.0, A, x, r);
1087 normr = fasp_blas_darray_norm2(m,r);
1088 normr_act = normr;
1089 relres = normr/n2b;
1090
1091 // if initial residual is small, no need to iterate!
1092 if ( normr <= tolb ) {
1093 flag = 0;
1094 iter = 0;
1095 goto FINISHED;
1096 }
1097
1098 // output iteration information if needed
1099 fasp_itinfo(PrtLvl,StopType,iter,relres,n2b,0.0);
1100
1101 // shadow residual rt = r* := r
1102 fasp_darray_cp(m,r,rt);
1103 normrmin = normr;
1104
1105 rho = 1.0;
1106 omega = 1.0;
1107 stag = 0;
1108 alpha = 0.0;
1109
1110 moresteps = 0;
1111 maxmsteps = 10;
1112 maxstagsteps = 3;
1113
1114 // loop over maxit iterations (unless convergence or failure)
1115 for (iter=1;iter <= MaxIt;iter++) {
1116
1117 rho1 = rho;
1118 rho = fasp_blas_darray_dotprod(m,rt,r);
1119
1120 if ((rho ==0.0 )|| (ABS(rho) >= DBL_MAX )) {
1121 flag = 4;
1122 goto FINISHED;
1123 }
1124
1125 if (iter==1) {
1126 fasp_darray_cp(m,r,p);
1127 }
1128 else {
1129 beta = (rho/rho1)*(alpha/omega);
1130
1131 if ((beta == 0)||( ABS(beta) > DBL_MAX )) {
1132 flag = 4;
1133 goto FINISHED;
1134 }
1135
1136 // p = r + beta * (p - omega * v);
1137 fasp_blas_darray_axpy(m,-omega,v,p); //p=p - omega*v
1138 fasp_blas_darray_axpby(m,1.0, r, beta, p); //p = 1.0*r +beta*p
1139 }
1140
1141 // pp = precond(p) ,ph
1142 if ( pc != NULL )
1143 pc->fct(p,ph,pc->data); /* Apply preconditioner */
1144 // if ph all is infinite then exit need add
1145 else
1146 fasp_darray_cp(m,p,ph); /* No preconditioner */
1147
1148 // v = A*ph
1149 fasp_blas_dstr_mxv(A,ph,v);
1150 rtv = fasp_blas_darray_dotprod(m,rt,v);
1151
1152 if (( rtv==0.0 )||( ABS(rtv) > DBL_MAX )) {
1153 flag = 4;
1154 goto FINISHED;
1155 }
1156
1157 alpha = rho/rtv;
1158
1159 if ( ABS(alpha) > DBL_MAX ) {
1160 flag = 4;
1161 ITS_DIVZERO;
1162 goto FINISHED;
1163 }
1164
1165 normx = fasp_blas_darray_norm2(m,x);
1166 normph = fasp_blas_darray_norm2(m,ph);
1167 if (ABS(alpha)*normph < DBL_EPSILON*normx )
1168 stag = stag + 1;
1169 else
1170 stag = 0;
1171
1172 // xhalf = x + alpha * ph; // form the "half" iterate
1173 // s = r - alpha * v; // residual associated with xhalf
1174 fasp_blas_darray_axpyz(m, alpha, ph, x , xhalf); // z= ax + y
1175 fasp_blas_darray_axpyz(m, -alpha, v, r, s);
1176 normr = fasp_blas_darray_norm2(m,s); // normr = norm(s);
1177 normr_act = normr;
1178
1179 // compute reduction factor of residual ||r||
1180 absres = normr_act;
1181 factor = absres/absres0;
1182 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
1183
1184 // check for convergence
1185 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
1186 {
1187 fasp_darray_cp(m,bval,s);
1188 fasp_blas_dstr_aAxpy(-1.0,A,xhalf,s);
1189 normr_act = fasp_blas_darray_norm2(m,s);
1190
1191 if (normr_act <= tolb){
1192 // x = xhalf;
1193 fasp_darray_cp(m,xhalf,x); // x = xhalf;
1194 flag = 0;
1195 imin = iter - 0.5;
1196 half_step++;
1197 if ( PrtLvl >= PRINT_MORE )
1198 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1199 flag,stag,imin,half_step);
1200 goto FINISHED;
1201 }
1202 else {
1203 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1204
1205 moresteps = moresteps + 1;
1206 if (moresteps >= maxmsteps){
1207 // if ~warned
1208 flag = 3;
1209 fasp_darray_cp(m,xhalf,x);
1210 goto FINISHED;
1211 }
1212 }
1213 }
1214
1215 if ( stag >= maxstagsteps ) {
1216 flag = 3;
1217 goto FINISHED;
1218 }
1219
1220 if ( normr_act < normrmin ) // update minimal norm quantities
1221 {
1222 normrmin = normr_act;
1223 fasp_darray_cp(m,xhalf,xmin);
1224 imin = iter - 0.5;
1225 half_step++;
1226 if ( PrtLvl >= PRINT_MORE )
1227 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1228 flag,stag,imin,half_step);
1229 }
1230
1231 // sh = precond(s)
1232 if ( pc != NULL ) {
1233 pc->fct(s,sh,pc->data); /* Apply preconditioner */
1234 }
1235 else
1236 fasp_darray_cp(m,s,sh); /* No preconditioner */
1237
1238 // t = A*sh;
1239 fasp_blas_dstr_mxv(A,sh,t);
1240 // tt = t' * t;
1241 tt = fasp_blas_darray_dotprod(m,t,t);
1242 if ( (tt == 0) ||( tt >= DBL_MAX ) ) {
1243 flag = 4;
1244 goto FINISHED;
1245 }
1246
1247 // omega = (t' * s) / tt;
1248 omega = fasp_blas_darray_dotprod(m,s,t)/tt;
1249 if ( ABS(omega) > DBL_MAX ) {
1250 flag = 4;
1251 goto FINISHED;
1252 }
1253
1254 norm_sh = fasp_blas_darray_norm2(m,sh);
1255 norm_xhalf = fasp_blas_darray_norm2(m,xhalf);
1256
1257 if ( ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
1258 stag = stag + 1;
1259 else
1260 stag = 0;
1261
1262 fasp_blas_darray_axpyz(m, omega,sh,xhalf, x); // x = xhalf + omega * sh;
1263 fasp_blas_darray_axpyz(m, -omega, t, s, r); // r = s - omega * t;
1264 normr = fasp_blas_darray_norm2(m,r); // normr = norm(r);
1265 normr_act = normr;
1266
1267 // check for convergence
1268 if ( (normr <= tolb) || (stag >= maxstagsteps) || moresteps )
1269 {
1270 fasp_darray_cp(m,bval,r);
1271 fasp_blas_dstr_aAxpy(-1.0,A,x,r);
1272 normr_act = fasp_blas_darray_norm2(m,r);
1273 if ( normr_act <= tolb ) {
1274 flag = 0;
1275 goto FINISHED;
1276 }
1277 else {
1278 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1279
1280 moresteps = moresteps + 1;
1281 if ( moresteps >= maxmsteps ) {
1282 flag = 3;
1283 goto FINISHED;
1284 }
1285 }
1286 }
1287
1288 // update minimal norm quantities
1289 if ( normr_act < normrmin ) {
1290 normrmin = normr_act;
1291 fasp_darray_cp(m,x,xmin);
1292 imin = iter;
1293 }
1294
1295 if ( stag >= maxstagsteps )
1296 {
1297 flag = 3;
1298 goto FINISHED;
1299 }
1300
1301 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1302
1303 absres0 = absres;
1304 } // for iter = 1 : maxit
1305
1306FINISHED: // finish iterative method
1307 // returned solution is first with minimal residual
1308 if (flag == 0)
1309 relres = normr_act / n2b;
1310 else {
1311 fasp_darray_cp(m, bval,r);
1312 fasp_blas_dstr_aAxpy(-1.0,A,xmin,r);
1313 normr = fasp_blas_darray_norm2(m,r);
1314
1315 if ( normr <= normr_act ) {
1316 fasp_darray_cp(m, xmin,x);
1317 iter = imin;
1318 relres = normr/n2b;
1319 }
1320 else {
1321 relres = normr_act/n2b;
1322 }
1323 }
1324
1325 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1326
1327 if ( PrtLvl >= PRINT_MORE )
1328 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1329 flag,stag,imin,half_step);
1330
1331 // clean up temp memory
1332 fasp_mem_free(work); work = NULL;
1333
1334#if DEBUG_MODE > 0
1335 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1336#endif
1337
1338 if ( iter > MaxIt )
1339 return ERROR_SOLVER_MAXIT;
1340 else
1341 return iter;
1342}
1343
1366 dvector *b,
1367 dvector *u,
1368 precond *pc,
1369 const REAL tol,
1370 const INT MaxIt,
1371 const SHORT StopType,
1372 const SHORT PrtLvl)
1373{
1374 const INT m = b->row;
1375
1376 // local variables
1377 REAL n2b,tolb;
1378 INT iter=0, stag = 1, moresteps = 1, maxmsteps=1;
1379 INT flag, maxstagsteps, half_step=0;
1380 REAL absres0 = BIGREAL, absres = BIGREAL, relres = BIGREAL;
1381 REAL alpha,beta,omega,rho,rho1,rtv,tt;
1382 REAL normr,normr_act,normph,normx,imin;
1383 REAL norm_sh,norm_xhalf,normrmin,factor;
1384 REAL *x = u->val, *bval=b->val;
1385
1386 // allocate temp memory (need 10*m REAL)
1387 REAL *work=(REAL *)fasp_mem_calloc(10*m,sizeof(REAL));
1388 REAL *r=work, *rt=r+m, *p=rt+m, *v=p+m;
1389 REAL *ph=v+m, *xhalf=ph+m, *s=xhalf+m, *sh=s+m;
1390 REAL *t = sh+m, *xmin = t+m;
1391
1392 // Output some info for debuging
1393 if ( PrtLvl > PRINT_NONE ) printf("\nCalling BiCGstab solver (MatFree) ...\n");
1394
1395#if DEBUG_MODE > 0
1396 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1397 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1398#endif
1399
1400 // r = b-A*u
1401 fasp_darray_cp(m,bval,r);
1402 n2b = fasp_blas_darray_norm2(m,r);
1403
1404 flag = 1;
1405 fasp_darray_cp(m,x,xmin);
1406 imin = 0;
1407
1408 iter = 0;
1409
1410 tolb = n2b*tol;
1411
1412 // r = b-A*x
1413 mf->fct(mf->data, x, r);
1414 fasp_blas_darray_axpby(m, 1.0, bval, -1.0, r);
1415 normr = fasp_blas_darray_norm2(m,r);
1416 normr_act = normr;
1417
1418 relres = normr/n2b;
1419 // if initial residual is small, no need to iterate!
1420 if (normr <= tolb) {
1421 flag =0;
1422 iter =0;
1423 goto FINISHED;
1424 }
1425
1426 // output iteration information if needed
1427
1428 fasp_itinfo(PrtLvl,StopType,iter,relres,n2b,0.0);
1429
1430 // shadow residual rt = r* := r
1431 fasp_darray_cp(m,r,rt);
1432 normrmin = normr;
1433
1434 rho = 1.0;
1435 omega = 1.0;
1436 stag = 0;
1437 alpha =0.0;
1438
1439 moresteps = 0;
1440 maxmsteps = 10;
1441 maxstagsteps = 3;
1442
1443 // loop over maxit iterations (unless convergence or failure)
1444 for (iter=1;iter <= MaxIt;iter++) {
1445
1446 rho1 = rho;
1447 rho = fasp_blas_darray_dotprod(m,rt,r);
1448
1449 if ((rho ==0.0 )|| (ABS(rho) >= DBL_MAX )) {
1450 flag = 4;
1451 goto FINISHED;
1452 }
1453
1454 if (iter==1) {
1455 fasp_darray_cp(m,r,p);
1456 }
1457 else {
1458 beta = (rho/rho1)*(alpha/omega);
1459
1460 if ((beta == 0)||( ABS(beta) > DBL_MAX )) {
1461 flag = 4;
1462 goto FINISHED;
1463 }
1464
1465 // p = r + beta * (p - omega * v);
1466 fasp_blas_darray_axpy(m,-omega,v,p); //p=p - omega*v
1467 fasp_blas_darray_axpby(m,1.0, r, beta, p); //p = 1.0*r +beta*p
1468 }
1469
1470 // pp = precond(p) ,ph
1471 if ( pc != NULL )
1472 pc->fct(p,ph,pc->data); /* Apply preconditioner */
1473 // if ph all is infinite then exit need add
1474 else
1475 fasp_darray_cp(m,p,ph); /* No preconditioner */
1476
1477 // v = A*ph
1478 mf->fct(mf->data, ph, v);
1479 rtv = fasp_blas_darray_dotprod(m,rt,v);
1480
1481 if (( rtv==0.0 )||( ABS(rtv) > DBL_MAX )) {
1482 flag = 4;
1483 goto FINISHED;
1484 }
1485
1486 alpha = rho/rtv;
1487
1488 if ( ABS(alpha) > DBL_MAX ) {
1489 flag = 4;
1490 ITS_DIVZERO;
1491 goto FINISHED;
1492 }
1493
1494 normx = fasp_blas_darray_norm2(m,x);
1495 normph = fasp_blas_darray_norm2(m,ph);
1496 if (ABS(alpha)*normph < DBL_EPSILON*normx )
1497 stag = stag + 1;
1498 else
1499 stag = 0;
1500
1501 // xhalf = x + alpha * ph; // form the "half" iterate
1502 // s = r - alpha * v; // residual associated with xhalf
1503 fasp_blas_darray_axpyz(m, alpha, ph, x , xhalf); // z= ax + y
1504 fasp_blas_darray_axpyz(m, -alpha, v, r, s);
1505 normr = fasp_blas_darray_norm2(m,s); // normr = norm(s);
1506 normr_act = normr;
1507
1508 // compute reduction factor of residual ||r||
1509 absres = normr_act;
1510 factor = absres/absres0;
1511 fasp_itinfo(PrtLvl,StopType,iter,normr_act/n2b,absres,factor);
1512
1513 // check for convergence
1514 if ((normr <= tolb)||(stag >= maxstagsteps)||moresteps)
1515 {
1516 // s = b-A*xhalf
1517 mf->fct(mf->data, xhalf, s);
1518 fasp_blas_darray_axpby(m, 1.0, bval, -1.0, s);
1519 normr_act = fasp_blas_darray_norm2(m,s);
1520
1521 if (normr_act <= tolb){
1522 // x = xhalf;
1523 fasp_darray_cp(m,xhalf,x); // x = xhalf;
1524 flag = 0;
1525 imin = iter - 0.5;
1526 half_step++;
1527 if ( PrtLvl >= PRINT_MORE )
1528 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1529 flag,stag,imin,half_step);
1530 goto FINISHED;
1531 }
1532 else {
1533 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1534
1535 moresteps = moresteps + 1;
1536 if (moresteps >= maxmsteps){
1537 // if ~warned
1538 flag = 3;
1539 fasp_darray_cp(m,xhalf,x);
1540 goto FINISHED;
1541 }
1542 }
1543 }
1544
1545 if ( stag >= maxstagsteps ) {
1546 flag = 3;
1547 goto FINISHED;
1548 }
1549
1550 if ( normr_act < normrmin ) // update minimal norm quantities
1551 {
1552 normrmin = normr_act;
1553 fasp_darray_cp(m,xhalf,xmin);
1554 imin = iter - 0.5;
1555 half_step++;
1556 if ( PrtLvl >= PRINT_MORE )
1557 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1558 flag,stag,imin,half_step);
1559 }
1560
1561 // sh = precond(s)
1562 if ( pc != NULL ){
1563 pc->fct(s,sh,pc->data); /* Apply preconditioner */
1564 //if all is finite
1565 }
1566 else
1567 fasp_darray_cp(m,s,sh); /* No preconditioner */
1568
1569 // t = A*sh;
1570 mf->fct(mf->data, sh, t);
1571 // tt = t' * t;
1572 tt = fasp_blas_darray_dotprod(m,t,t);
1573 if ((tt == 0) ||( tt >= DBL_MAX )) {
1574 flag = 4;
1575 goto FINISHED;
1576 }
1577
1578 // omega = (t' * s) / tt;
1579 omega = fasp_blas_darray_dotprod(m,s,t)/tt;
1580 if (ABS(omega) > DBL_MAX )
1581 {
1582 flag = 4;
1583 goto FINISHED;
1584 }
1585
1586 norm_sh = fasp_blas_darray_norm2(m,sh);
1587 norm_xhalf = fasp_blas_darray_norm2(m,xhalf);
1588
1589 if (ABS(omega)*norm_sh < DBL_EPSILON*norm_xhalf )
1590 stag = stag + 1;
1591 else
1592 stag = 0;
1593
1594 fasp_blas_darray_axpyz(m, omega,sh,xhalf, x); // x = xhalf + omega * sh;
1595 fasp_blas_darray_axpyz(m, -omega, t, s, r); // r = s - omega * t;
1596 normr = fasp_blas_darray_norm2(m,r); //normr = norm(r);
1597 normr_act = normr;
1598
1599 // check for convergence
1600 if ( (normr <= tolb)||(stag >= maxstagsteps)||moresteps )
1601 {
1602 // normr_act = norm(r);
1603 // r = b-A*x
1604 mf->fct(mf->data, x, r);
1605 fasp_blas_darray_axpby(m, 1.0, bval, -1.0, r);
1606 normr_act = fasp_blas_darray_norm2(m,r);
1607 if (normr_act <= tolb)
1608 {
1609 flag = 0;
1610 goto FINISHED;
1611 }
1612 else
1613 {
1614 if ((stag >= maxstagsteps) && (moresteps == 0)) stag = 0;
1615
1616 moresteps = moresteps + 1;
1617 if (moresteps >= maxmsteps)
1618 {
1619 flag = 3;
1620 goto FINISHED;
1621 }
1622 }
1623 }
1624
1625 if (normr_act < normrmin) // update minimal norm quantities
1626 {
1627 normrmin = normr_act;
1628 fasp_darray_cp(m,x,xmin);
1629 imin = iter;
1630 }
1631
1632 if (stag >= maxstagsteps)
1633 {
1634 flag = 3;
1635 goto FINISHED;
1636 }
1637
1638 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1639
1640 absres0 = absres;
1641 } // for iter = 1 : maxit
1642
1643FINISHED: // finish iterative method
1644 // returned solution is first with minimal residual
1645 if (flag == 0)
1646 relres = normr_act / n2b;
1647 else {
1648 // r = b-A*xmin
1649 mf->fct(mf->data, xmin, r);
1650 fasp_blas_darray_axpby(m, 1.0, bval, -1.0, r);
1651 normr = fasp_blas_darray_norm2(m,r);
1652
1653 if ( normr <= normr_act ) {
1654 fasp_darray_cp(m, xmin,x);
1655 iter = imin;
1656 relres = normr/n2b;
1657 }
1658 else {
1659 relres = normr_act/n2b;
1660 }
1661 }
1662
1663 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1664
1665 if ( PrtLvl >= PRINT_MORE )
1666 printf("Flag = %d Stag = %d Itermin = %.1f Half_step = %d\n",
1667 flag,stag,imin,half_step);
1668
1669 // clean up temp memory
1670 fasp_mem_free(work); work = NULL;
1671
1672#if DEBUG_MODE > 0
1673 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1674#endif
1675
1676 if ( iter > MaxIt )
1677 return ERROR_SOLVER_MAXIT;
1678 else
1679 return iter;
1680}
1681
1682/*---------------------------------*/
1683/*-- End of File --*/
1684/*---------------------------------*/
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
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
void fasp_blas_darray_axpyz(const INT n, const REAL a, const REAL *x, const REAL *y, REAL *z)
z = a*x + y
Definition: BlaArray.c:347
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_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_dstr_pbcgs(dSTRmat *A, 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 for STR matrix.
Definition: KryPbcgs.c:1039
INT fasp_solver_dblc_pbcgs(dBLCmat *A, 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 for BLC matrix.
Definition: KryPbcgs.c:713
INT fasp_solver_pbcgs(mxv_matfree *mf, 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.
Definition: KryPbcgs.c:1365
INT fasp_solver_dbsr_pbcgs(dBSRmat *A, 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 for BSR matrix.
Definition: KryPbcgs.c:387
INT fasp_solver_dcsr_pbcgs(dCSRmat *A, 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 for CSR matrix.
Definition: KryPbcgs.c:62
Main header file for the FASP project.
#define REAL
Definition: fasp.h:67
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
#define ABS(a)
Definition: fasp.h:76
#define INT
Definition: fasp.h:64
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
#define BIGREAL
Some global constants.
Definition: fasp_const.h:248
#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
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
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