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