Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
KryPcg.c
Go to the documentation of this file.
1
62#include <math.h>
63
64#include "fasp.h"
65#include "fasp_functs.h"
66
67/*---------------------------------*/
68/*-- Declare Private Functions --*/
69/*---------------------------------*/
70
71#include "KryUtil.inl"
72
73/*---------------------------------*/
74/*-- Public Functions --*/
75/*---------------------------------*/
76
99 dvector *b,
100 dvector *u,
101 precond *pc,
102 const REAL tol,
103 const INT MaxIt,
104 const SHORT StopType,
105 const SHORT PrtLvl)
106{
107 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
108 const INT m = b->row;
109 const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
110 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
111
112 // local variables
113 INT iter = 0, stag = 1, more_step = 1;
114 REAL absres0 = BIGREAL, absres = BIGREAL;
115 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
116 REAL reldiff, factor, normuinf;
117 REAL alpha, beta, temp1, temp2;
118
119 // allocate temp memory (need 4*m REAL numbers)
120 REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
121 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
122
123 // Output some info for debugging
124 if ( PrtLvl > PRINT_NONE ) printf("\nCalling CG solver (CSR) ...\n");
125
126#if DEBUG_MODE > 0
127 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
128 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
129#endif
130
131 // r = b-A*u
132 fasp_darray_cp(m,b->val,r);
133 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
134
135 if ( pc != NULL )
136 pc->fct(r,z,pc->data); /* Apply preconditioner */
137 else
138 fasp_darray_cp(m,r,z); /* No preconditioner */
139
140 // compute initial residuals
141 switch ( StopType ) {
142 case STOP_REL_RES:
143 absres0 = fasp_blas_darray_norm2(m,r);
144 normr0 = MAX(SMALLREAL,absres0);
145 relres = absres0/normr0;
146 break;
147 case STOP_REL_PRECRES:
148 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
149 normr0 = MAX(SMALLREAL,absres0);
150 relres = absres0/normr0;
151 break;
152 case STOP_MOD_REL_RES:
153 absres0 = fasp_blas_darray_norm2(m,r);
155 relres = absres0/normu;
156 break;
157 default:
158 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
159 goto FINISHED;
160 }
161
162 // if initial residual is small, no need to iterate!
163 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
164
165 // output iteration information if needed
166 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
167
168 fasp_darray_cp(m,z,p);
169 temp1 = fasp_blas_darray_dotprod(m,z,r);
170
171 // main PCG loop
172 while ( iter++ < MaxIt ) {
173
174 // t = A*p
175 fasp_blas_dcsr_mxv(A,p,t);
176
177 // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
178 temp2 = fasp_blas_darray_dotprod(m,t,p);
179 if ( ABS(temp2) > SMALLREAL2 ) {
180 alpha = temp1/temp2;
181 }
182 else { // Possible breakdown
183 ITS_DIVZERO; goto FINISHED;
184 }
185
186 // u_k = u_{k-1} + alpha_k*p_{k-1}
187 fasp_blas_darray_axpy(m,alpha,p,u->val);
188
189 // r_k = r_{k-1} - alpha_k*A*p_{k-1}
190 fasp_blas_darray_axpy(m,-alpha,t,r);
191
192 // compute norm of residual
193 switch ( StopType ) {
194 case STOP_REL_RES:
195 absres = fasp_blas_darray_norm2(m,r);
196 relres = absres/normr0;
197 break;
198 case STOP_REL_PRECRES:
199 // z = B(r)
200 if ( pc != NULL )
201 pc->fct(r,z,pc->data); /* Apply preconditioner */
202 else
203 fasp_darray_cp(m,r,z); /* No preconditioner */
204 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
205 relres = absres/normr0;
206 break;
207 case STOP_MOD_REL_RES:
208 absres = fasp_blas_darray_norm2(m,r);
209 relres = absres/normu;
210 break;
211 }
212
213 // compute reduction factor of residual ||r||
214 factor = absres/absres0;
215
216 // output iteration information if needed
217 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
218
219 if ( factor > 0.9 ) { // Only check when converge slowly
220
221 // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
222 normuinf = fasp_blas_darray_norminf(m, u->val);
223 if ( normuinf <= sol_inf_tol ) {
224 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
226 break;
227 }
228
229 // Check II: if stagnated, try to restart
230 normu = fasp_blas_darray_norm2(m, u->val);
231
232 // compute relative difference
233 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
234 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
235
236 if ( PrtLvl >= PRINT_MORE ) {
237 ITS_DIFFRES(reldiff,relres);
238 ITS_RESTART;
239 }
240
241 fasp_darray_cp(m,b->val,r);
242 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
243
244 // compute residual norms
245 switch ( StopType ) {
246 case STOP_REL_RES:
247 absres = fasp_blas_darray_norm2(m,r);
248 relres = absres/normr0;
249 break;
250 case STOP_REL_PRECRES:
251 // z = B(r)
252 if ( pc != NULL )
253 pc->fct(r,z,pc->data); /* Apply preconditioner */
254 else
255 fasp_darray_cp(m,r,z); /* No preconditioner */
256 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
257 relres = absres/normr0;
258 break;
259 case STOP_MOD_REL_RES:
260 absres = fasp_blas_darray_norm2(m,r);
261 relres = absres/normu;
262 break;
263 }
264
265 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
266
267 if ( relres < tol )
268 break;
269 else {
270 if ( stag >= MaxStag ) {
271 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
272 iter = ERROR_SOLVER_STAG;
273 break;
274 }
275 fasp_darray_set(m,p,0.0);
276 ++stag;
277 }
278
279 } // end of stagnation check!
280
281 } // end of check I and II
282
283 // Check III: prevent false convergence
284 if ( relres < tol ) {
285
286 REAL updated_relres = relres;
287
288 // compute true residual r = b - Ax and update residual
289 fasp_darray_cp(m,b->val,r);
290 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
291
292 // compute residual norms
293 switch ( StopType ) {
294 case STOP_REL_RES:
295 absres = fasp_blas_darray_norm2(m,r);
296 relres = absres/normr0;
297 break;
298 case STOP_REL_PRECRES:
299 // z = B(r)
300 if ( pc != NULL )
301 pc->fct(r,z,pc->data); /* Apply preconditioner */
302 else
303 fasp_darray_cp(m,r,z); /* No preconditioner */
304 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
305 relres = absres/normr0;
306 break;
307 case STOP_MOD_REL_RES:
308 absres = fasp_blas_darray_norm2(m,r);
309 relres = absres/normu;
310 break;
311 }
312
313 // check convergence
314 if ( relres < tol ) break;
315
316 if ( PrtLvl >= PRINT_MORE ) {
317 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
318 }
319
320 if ( more_step >= MaxRestartStep ) {
321 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
323 break;
324 }
325
326 // prepare for restarting method
327 fasp_darray_set(m,p,0.0);
328 ++more_step;
329
330 } // end of safe-guard check!
331
332 // save residual for next iteration
333 absres0 = absres;
334
335 // compute z_k = B(r_k)
336 if ( StopType != STOP_REL_PRECRES ) {
337 if ( pc != NULL )
338 pc->fct(r,z,pc->data); /* Apply preconditioner */
339 else
340 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
341 }
342
343 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
344 temp2 = fasp_blas_darray_dotprod(m,z,r);
345 beta = temp2/temp1;
346 temp1 = temp2;
347
348 // compute p_k = z_k + beta_k*p_{k-1}
349 fasp_blas_darray_axpby(m,1.0,z,beta,p);
350
351 } // end of main PCG loop.
352
353FINISHED: // finish iterative method
354 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
355
356 // clean up temp memory
357 fasp_mem_free(work); work = NULL;
358
359#if DEBUG_MODE > 0
360 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
361#endif
362
363 if ( iter > MaxIt )
364 return ERROR_SOLVER_MAXIT;
365 else
366 return iter;
367}
368
391 dvector *b,
392 dvector *u,
393 precond *pc,
394 const REAL tol,
395 const INT MaxIt,
396 const SHORT StopType,
397 const SHORT PrtLvl)
398{
399 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
400 const INT m = b->row;
401 const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
402 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
403
404 // local variables
405 INT iter = 0, stag = 1, more_step = 1;
406 REAL absres0 = BIGREAL, absres = BIGREAL;
407 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
408 REAL reldiff, factor, normuinf;
409 REAL alpha, beta, temp1, temp2;
410
411 // allocate temp memory (need 4*m REAL numbers)
412 REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
413 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
414
415 // Output some info for debuging
416 if ( PrtLvl > PRINT_NONE ) printf("\nCalling CG solver (BSR) ...\n");
417
418#if DEBUG_MODE > 0
419 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
420 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
421#endif
422
423 // r = b-A*u
424 fasp_darray_cp(m,b->val,r);
425 fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
426
427 if ( pc != NULL )
428 pc->fct(r,z,pc->data); /* Apply preconditioner */
429 else
430 fasp_darray_cp(m,r,z); /* No preconditioner */
431
432 // compute initial residuals
433 switch ( StopType ) {
434 case STOP_REL_RES:
435 absres0 = fasp_blas_darray_norm2(m,r);
436 normr0 = MAX(SMALLREAL,absres0);
437 relres = absres0/normr0;
438 break;
439 case STOP_REL_PRECRES:
440 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
441 normr0 = MAX(SMALLREAL,absres0);
442 relres = absres0/normr0;
443 break;
444 case STOP_MOD_REL_RES:
445 absres0 = fasp_blas_darray_norm2(m,r);
447 relres = absres0/normu;
448 break;
449 default:
450 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
451 goto FINISHED;
452 }
453
454 // if initial residual is small, no need to iterate!
455 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
456
457 // output iteration information if needed
458 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
459
460 fasp_darray_cp(m,z,p);
461 temp1 = fasp_blas_darray_dotprod(m,z,r);
462
463 // main PCG loop
464 while ( iter++ < MaxIt ) {
465
466 // t = A*p
467 fasp_blas_dbsr_mxv(A,p,t);
468
469 // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
470 temp2 = fasp_blas_darray_dotprod(m,t,p);
471 if ( ABS(temp2) > SMALLREAL2 ) {
472 alpha = temp1/temp2;
473 }
474 else { // Possible breakdown
475 ITS_DIVZERO; goto FINISHED;
476 }
477
478 // u_k = u_{k-1} + alpha_k*p_{k-1}
479 fasp_blas_darray_axpy(m,alpha,p,u->val);
480
481 // r_k = r_{k-1} - alpha_k*A*p_{k-1}
482 fasp_blas_darray_axpy(m,-alpha,t,r);
483
484 // compute norm of residual
485 switch ( StopType ) {
486 case STOP_REL_RES:
487 absres = fasp_blas_darray_norm2(m,r);
488 relres = absres/normr0;
489 break;
490 case STOP_REL_PRECRES:
491 // z = B(r)
492 if ( pc != NULL )
493 pc->fct(r,z,pc->data); /* Apply preconditioner */
494 else
495 fasp_darray_cp(m,r,z); /* No preconditioner */
496 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
497 relres = absres/normr0;
498 break;
499 case STOP_MOD_REL_RES:
500 absres = fasp_blas_darray_norm2(m,r);
501 relres = absres/normu;
502 break;
503 }
504
505 // compute reduction factor of residual ||r||
506 factor = absres/absres0;
507
508 // output iteration information if needed
509 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
510
511 if ( factor > 0.9 ) { // Only check when converge slowly
512
513 // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
514 normuinf = fasp_blas_darray_norminf(m, u->val);
515 if ( normuinf <= sol_inf_tol ) {
516 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
518 break;
519 }
520
521 // Check II: if stagnated, try to restart
522 normu = fasp_blas_darray_norm2(m, u->val);
523
524 // compute relative difference
525 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
526 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
527
528 if ( PrtLvl >= PRINT_MORE ) {
529 ITS_DIFFRES(reldiff,relres);
530 ITS_RESTART;
531 }
532
533 fasp_darray_cp(m,b->val,r);
534 fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
535
536 // compute residual norms
537 switch ( StopType ) {
538 case STOP_REL_RES:
539 absres = fasp_blas_darray_norm2(m,r);
540 relres = absres/normr0;
541 break;
542 case STOP_REL_PRECRES:
543 // z = B(r)
544 if ( pc != NULL )
545 pc->fct(r,z,pc->data); /* Apply preconditioner */
546 else
547 fasp_darray_cp(m,r,z); /* No preconditioner */
548 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
549 relres = absres/normr0;
550 break;
551 case STOP_MOD_REL_RES:
552 absres = fasp_blas_darray_norm2(m,r);
553 relres = absres/normu;
554 break;
555 }
556
557 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
558
559 if ( relres < tol )
560 break;
561 else {
562 if ( stag >= MaxStag ) {
563 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
564 iter = ERROR_SOLVER_STAG;
565 break;
566 }
567 fasp_darray_set(m,p,0.0);
568 ++stag;
569 }
570
571 } // end of stagnation check!
572
573 } // end of check I and II
574
575 // Check III: prevent false convergence
576 if ( relres < tol ) {
577
578 REAL updated_relres = relres;
579
580 // compute true residual r = b - Ax and update residual
581 fasp_darray_cp(m,b->val,r);
582 fasp_blas_dbsr_aAxpy(-1.0,A,u->val,r);
583
584 // compute residual norms
585 switch ( StopType ) {
586 case STOP_REL_RES:
587 absres = fasp_blas_darray_norm2(m,r);
588 relres = absres/normr0;
589 break;
590 case STOP_REL_PRECRES:
591 // z = B(r)
592 if ( pc != NULL )
593 pc->fct(r,z,pc->data); /* Apply preconditioner */
594 else
595 fasp_darray_cp(m,r,z); /* No preconditioner */
596 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
597 relres = absres/normr0;
598 break;
599 case STOP_MOD_REL_RES:
600 absres = fasp_blas_darray_norm2(m,r);
601 relres = absres/normu;
602 break;
603 }
604
605 // check convergence
606 if ( relres < tol ) break;
607
608 if ( PrtLvl >= PRINT_MORE ) {
609 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
610 }
611
612 if ( more_step >= MaxRestartStep ) {
613 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
615 break;
616 }
617
618 // prepare for restarting method
619 fasp_darray_set(m,p,0.0);
620 ++more_step;
621
622 } // end of safe-guard check!
623
624 // save residual for next iteration
625 absres0 = absres;
626
627 // compute z_k = B(r_k)
628 if ( StopType != STOP_REL_PRECRES ) {
629 if ( pc != NULL )
630 pc->fct(r,z,pc->data); /* Apply preconditioner */
631 else
632 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
633 }
634
635 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
636 temp2 = fasp_blas_darray_dotprod(m,z,r);
637 beta = temp2/temp1;
638 temp1 = temp2;
639
640 // compute p_k = z_k + beta_k*p_{k-1}
641 fasp_blas_darray_axpby(m,1.0,z,beta,p);
642
643 } // end of main PCG loop.
644
645FINISHED: // finish iterative method
646 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
647
648 // clean up temp memory
649 fasp_mem_free(work); work = NULL;
650
651#if DEBUG_MODE > 0
652 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
653#endif
654
655 if ( iter > MaxIt )
656 return ERROR_SOLVER_MAXIT;
657 else
658 return iter;
659}
660
685 dvector *b,
686 dvector *u,
687 precond *pc,
688 const REAL tol,
689 const INT MaxIt,
690 const SHORT StopType,
691 const SHORT PrtLvl)
692{
693 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
694 const INT m = b->row;
695 const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
696 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
697
698 // local variables
699 INT iter = 0, stag = 1, more_step = 1;
700 REAL absres0 = BIGREAL, absres = BIGREAL;
701 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
702 REAL reldiff, factor, normuinf;
703 REAL alpha, beta, temp1, temp2;
704
705 // allocate temp memory (need 4*m REAL numbers)
706 REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
707 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
708
709 // Output some info for debuging
710 if ( PrtLvl > PRINT_NONE ) printf("\nCalling CG solver (BLC) ...\n");
711
712#if DEBUG_MODE > 0
713 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
714 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
715#endif
716
717 // r = b-A*u
718 fasp_darray_cp(m,b->val,r);
719 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
720
721 if ( pc != NULL )
722 pc->fct(r,z,pc->data); /* Apply preconditioner */
723 else
724 fasp_darray_cp(m,r,z); /* No preconditioner */
725
726 // compute initial residuals
727 switch ( StopType ) {
728 case STOP_REL_RES:
729 absres0 = fasp_blas_darray_norm2(m,r);
730 normr0 = MAX(SMALLREAL,absres0);
731 relres = absres0/normr0;
732 break;
733 case STOP_REL_PRECRES:
734 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
735 normr0 = MAX(SMALLREAL,absres0);
736 relres = absres0/normr0;
737 break;
738 case STOP_MOD_REL_RES:
739 absres0 = fasp_blas_darray_norm2(m,r);
741 relres = absres0/normu;
742 break;
743 default:
744 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
745 goto FINISHED;
746 }
747
748 // if initial residual is small, no need to iterate!
749 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
750
751 // output iteration information if needed
752 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
753
754 fasp_darray_cp(m,z,p);
755 temp1 = fasp_blas_darray_dotprod(m,z,r);
756
757 // main PCG loop
758 while ( iter++ < MaxIt ) {
759
760 // t = A*p
761 fasp_blas_dblc_mxv(A,p,t);
762
763 // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
764 temp2 = fasp_blas_darray_dotprod(m,t,p);
765 if ( ABS(temp2) > SMALLREAL2 ) {
766 alpha = temp1/temp2;
767 }
768 else { // Possible breakdown
769 ITS_DIVZERO; goto FINISHED;
770 }
771
772 // u_k = u_{k-1} + alpha_k*p_{k-1}
773 fasp_blas_darray_axpy(m,alpha,p,u->val);
774
775 // r_k = r_{k-1} - alpha_k*A*p_{k-1}
776 fasp_blas_darray_axpy(m,-alpha,t,r);
777
778 // compute norm of residual
779 switch ( StopType ) {
780 case STOP_REL_RES:
781 absres = fasp_blas_darray_norm2(m,r);
782 relres = absres/normr0;
783 break;
784 case STOP_REL_PRECRES:
785 // z = B(r)
786 if ( pc != NULL )
787 pc->fct(r,z,pc->data); /* Apply preconditioner */
788 else
789 fasp_darray_cp(m,r,z); /* No preconditioner */
790 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
791 relres = absres/normr0;
792 break;
793 case STOP_MOD_REL_RES:
794 absres = fasp_blas_darray_norm2(m,r);
795 relres = absres/normu;
796 break;
797 }
798
799 // compute reduction factor of residual ||r||
800 factor = absres/absres0;
801
802 // output iteration information if needed
803 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
804
805 if ( factor > 0.9 ) { // Only check when converge slowly
806
807 // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
808 normuinf = fasp_blas_darray_norminf(m, u->val);
809 if ( normuinf <= sol_inf_tol ) {
810 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
812 break;
813 }
814
815 // Check II: if stagnated, try to restart
816 normu = fasp_blas_darray_norm2(m, u->val);
817
818 // compute relative difference
819 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
820 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
821
822 if ( PrtLvl >= PRINT_MORE ) {
823 ITS_DIFFRES(reldiff,relres);
824 ITS_RESTART;
825 }
826
827 fasp_darray_cp(m,b->val,r);
828 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
829
830 // compute residual norms
831 switch ( StopType ) {
832 case STOP_REL_RES:
833 absres = fasp_blas_darray_norm2(m,r);
834 relres = absres/normr0;
835 break;
836 case STOP_REL_PRECRES:
837 // z = B(r)
838 if ( pc != NULL )
839 pc->fct(r,z,pc->data); /* Apply preconditioner */
840 else
841 fasp_darray_cp(m,r,z); /* No preconditioner */
842 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
843 relres = absres/normr0;
844 break;
845 case STOP_MOD_REL_RES:
846 absres = fasp_blas_darray_norm2(m,r);
847 relres = absres/normu;
848 break;
849 }
850
851 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
852
853 if ( relres < tol )
854 break;
855 else {
856 if ( stag >= MaxStag ) {
857 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
858 iter = ERROR_SOLVER_STAG;
859 break;
860 }
861 fasp_darray_set(m,p,0.0);
862 ++stag;
863 }
864
865 } // end of stagnation check!
866
867 } // end of check I and II
868
869 // Check III: prevent false convergence
870 if ( relres < tol ) {
871
872 REAL updated_relres = relres;
873
874 // compute true residual r = b - Ax and update residual
875 fasp_darray_cp(m,b->val,r);
876 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
877
878 // compute residual norms
879 switch ( StopType ) {
880 case STOP_REL_RES:
881 absres = fasp_blas_darray_norm2(m,r);
882 relres = absres/normr0;
883 break;
884 case STOP_REL_PRECRES:
885 // z = B(r)
886 if ( pc != NULL )
887 pc->fct(r,z,pc->data); /* Apply preconditioner */
888 else
889 fasp_darray_cp(m,r,z); /* No preconditioner */
890 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
891 relres = absres/normr0;
892 break;
893 case STOP_MOD_REL_RES:
894 absres = fasp_blas_darray_norm2(m,r);
895 relres = absres/normu;
896 break;
897 }
898
899 // check convergence
900 if ( relres < tol ) break;
901
902 if ( PrtLvl >= PRINT_MORE ) {
903 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
904 }
905
906 if ( more_step >= MaxRestartStep ) {
907 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
909 break;
910 }
911
912 // prepare for restarting method
913 fasp_darray_set(m,p,0.0);
914 ++more_step;
915
916 } // end of safe-guard check!
917
918 // save residual for next iteration
919 absres0 = absres;
920
921 // compute z_k = B(r_k)
922 if ( StopType != STOP_REL_PRECRES ) {
923 if ( pc != NULL )
924 pc->fct(r,z,pc->data); /* Apply preconditioner */
925 else
926 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
927 }
928
929 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
930 temp2 = fasp_blas_darray_dotprod(m,z,r);
931 beta = temp2/temp1;
932 temp1 = temp2;
933
934 // compute p_k = z_k + beta_k*p_{k-1}
935 fasp_blas_darray_axpby(m,1.0,z,beta,p);
936
937 } // end of main PCG loop.
938
939FINISHED: // finish iterative method
940 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
941
942 // clean up temp memory
943 fasp_mem_free(work); work = NULL;
944
945#if DEBUG_MODE > 0
946 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
947#endif
948
949 if ( iter > MaxIt )
950 return ERROR_SOLVER_MAXIT;
951 else
952 return iter;
953}
954
979 dvector *b,
980 dvector *u,
981 precond *pc,
982 const REAL tol,
983 const INT MaxIt,
984 const SHORT StopType,
985 const SHORT PrtLvl)
986{
987 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
988 const INT m = b->row;
989 const REAL maxdiff = tol*STAG_RATIO; // stagnation tolerance
990 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
991
992 // local variables
993 INT iter = 0, stag = 1, more_step = 1;
994 REAL absres0 = BIGREAL, absres = BIGREAL;
995 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
996 REAL reldiff, factor, normuinf;
997 REAL alpha, beta, temp1, temp2;
998
999 // allocate temp memory (need 4*m REAL numbers)
1000 REAL *work = (REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
1001 REAL *p = work, *z = work+m, *r = z+m, *t = r+m;
1002
1003 // Output some info for debuging
1004 if ( PrtLvl > PRINT_NONE ) printf("\nCalling CG solver (STR) ...\n");
1005
1006#if DEBUG_MODE > 0
1007 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1008 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1009#endif
1010
1011 // r = b-A*u
1012 fasp_darray_cp(m,b->val,r);
1013 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1014
1015 if ( pc != NULL )
1016 pc->fct(r,z,pc->data); /* Apply preconditioner */
1017 else
1018 fasp_darray_cp(m,r,z); /* No preconditioner */
1019
1020 // compute initial residuals
1021 switch ( StopType ) {
1022 case STOP_REL_RES:
1023 absres0 = fasp_blas_darray_norm2(m,r);
1024 normr0 = MAX(SMALLREAL,absres0);
1025 relres = absres0/normr0;
1026 break;
1027 case STOP_REL_PRECRES:
1028 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,z));
1029 normr0 = MAX(SMALLREAL,absres0);
1030 relres = absres0/normr0;
1031 break;
1032 case STOP_MOD_REL_RES:
1033 absres0 = fasp_blas_darray_norm2(m,r);
1034 normu = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
1035 relres = absres0/normu;
1036 break;
1037 default:
1038 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1039 goto FINISHED;
1040 }
1041
1042 // if initial residual is small, no need to iterate!
1043 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
1044
1045 // output iteration information if needed
1046 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
1047
1048 fasp_darray_cp(m,z,p);
1049 temp1 = fasp_blas_darray_dotprod(m,z,r);
1050
1051 // main PCG loop
1052 while ( iter++ < MaxIt ) {
1053
1054 // t = A*p
1055 fasp_blas_dstr_mxv(A,p,t);
1056
1057 // alpha_k = (z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
1058 temp2 = fasp_blas_darray_dotprod(m,t,p);
1059 if ( ABS(temp2) > SMALLREAL2 ) {
1060 alpha = temp1/temp2;
1061 }
1062 else { // Possible breakdown
1063 ITS_DIVZERO; goto FINISHED;
1064 }
1065
1066 // u_k = u_{k-1} + alpha_k*p_{k-1}
1067 fasp_blas_darray_axpy(m,alpha,p,u->val);
1068
1069 // r_k = r_{k-1} - alpha_k*A*p_{k-1}
1070 fasp_blas_darray_axpy(m,-alpha,t,r);
1071
1072 // compute norm of residual
1073 switch ( StopType ) {
1074 case STOP_REL_RES:
1075 absres = fasp_blas_darray_norm2(m,r);
1076 relres = absres/normr0;
1077 break;
1078 case STOP_REL_PRECRES:
1079 // z = B(r)
1080 if ( pc != NULL )
1081 pc->fct(r,z,pc->data); /* Apply preconditioner */
1082 else
1083 fasp_darray_cp(m,r,z); /* No preconditioner */
1084 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1085 relres = absres/normr0;
1086 break;
1087 case STOP_MOD_REL_RES:
1088 absres = fasp_blas_darray_norm2(m,r);
1089 relres = absres/normu;
1090 break;
1091 }
1092
1093 // compute reduction factor of residual ||r||
1094 factor = absres/absres0;
1095
1096 // output iteration information if needed
1097 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1098
1099 if ( factor > 0.9 ) { // Only check when converge slowly
1100
1101 // Check I: if solution is close to zero, return ERROR_SOLVER_SOLSTAG
1102 normuinf = fasp_blas_darray_norminf(m, u->val);
1103 if ( normuinf <= sol_inf_tol ) {
1104 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
1105 iter = ERROR_SOLVER_SOLSTAG;
1106 break;
1107 }
1108
1109 // Check II: if stagnated, try to restart
1110 normu = fasp_blas_darray_norm2(m, u->val);
1111
1112 // compute relative difference
1113 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m,p)/normu;
1114 if ( (stag <= MaxStag) & (reldiff < maxdiff) ) {
1115
1116 if ( PrtLvl >= PRINT_MORE ) {
1117 ITS_DIFFRES(reldiff,relres);
1118 ITS_RESTART;
1119 }
1120
1121 fasp_darray_cp(m,b->val,r);
1122 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1123
1124 // compute residual norms
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 // z = B(r)
1132 if ( pc != NULL )
1133 pc->fct(r,z,pc->data); /* Apply preconditioner */
1134 else
1135 fasp_darray_cp(m,r,z); /* No preconditioner */
1136 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1137 relres = absres/normr0;
1138 break;
1139 case STOP_MOD_REL_RES:
1140 absres = fasp_blas_darray_norm2(m,r);
1141 relres = absres/normu;
1142 break;
1143 }
1144
1145 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1146
1147 if ( relres < tol )
1148 break;
1149 else {
1150 if ( stag >= MaxStag ) {
1151 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
1152 iter = ERROR_SOLVER_STAG;
1153 break;
1154 }
1155 fasp_darray_set(m,p,0.0);
1156 ++stag;
1157 }
1158
1159 } // end of stagnation check!
1160
1161 } // end of check I and II
1162
1163 // Check III: prevent false convergence
1164 if ( relres < tol ) {
1165
1166 REAL updated_relres = relres;
1167
1168 // compute true residual r = b - Ax and update residual
1169 fasp_darray_cp(m,b->val,r);
1170 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1171
1172 // compute residual norms
1173 switch ( StopType ) {
1174 case STOP_REL_RES:
1175 absres = fasp_blas_darray_norm2(m,r);
1176 relres = absres/normr0;
1177 break;
1178 case STOP_REL_PRECRES:
1179 // z = B(r)
1180 if ( pc != NULL )
1181 pc->fct(r,z,pc->data); /* Apply preconditioner */
1182 else
1183 fasp_darray_cp(m,r,z); /* No preconditioner */
1184 absres = sqrt(ABS(fasp_blas_darray_dotprod(m,z,r)));
1185 relres = absres/normr0;
1186 break;
1187 case STOP_MOD_REL_RES:
1188 absres = fasp_blas_darray_norm2(m,r);
1189 relres = absres/normu;
1190 break;
1191 }
1192
1193 // check convergence
1194 if ( relres < tol ) break;
1195
1196 if ( PrtLvl >= PRINT_MORE ) {
1197 ITS_COMPRES(updated_relres); ITS_REALRES(relres);
1198 }
1199
1200 if ( more_step >= MaxRestartStep ) {
1201 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
1202 iter = ERROR_SOLVER_TOLSMALL;
1203 break;
1204 }
1205
1206 // prepare for restarting method
1207 fasp_darray_set(m,p,0.0);
1208 ++more_step;
1209
1210 } // end of safe-guard check!
1211
1212 // save residual for next iteration
1213 absres0 = absres;
1214
1215 // compute z_k = B(r_k)
1216 if ( StopType != STOP_REL_PRECRES ) {
1217 if ( pc != NULL )
1218 pc->fct(r,z,pc->data); /* Apply preconditioner */
1219 else
1220 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
1221 }
1222
1223 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
1224 temp2 = fasp_blas_darray_dotprod(m,z,r);
1225 beta = temp2/temp1;
1226 temp1 = temp2;
1227
1228 // compute p_k = z_k + beta_k*p_{k-1}
1229 fasp_blas_darray_axpby(m,1.0,z,beta,p);
1230
1231 } // end of main PCG loop.
1232
1233FINISHED: // finish iterative method
1234 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1235
1236 // clean up temp memory
1237 fasp_mem_free(work); work = NULL;
1238
1239#if DEBUG_MODE > 0
1240 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1241#endif
1242
1243 if ( iter > MaxIt )
1244 return ERROR_SOLVER_MAXIT;
1245 else
1246 return iter;
1247}
1248
1273 dvector *b,
1274 dvector *u,
1275 precond *pc,
1276 const REAL tol,
1277 const INT MaxIt,
1278 const SHORT StopType,
1279 const SHORT PrtLvl)
1280{
1281 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
1282 const INT m=b->row;
1283 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1284 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
1285
1286 // local variables
1287 INT iter = 0, stag, more_step, restart_step;
1288 REAL absres0 = BIGREAL, absres = BIGREAL;
1289 REAL relres = BIGREAL, normu = BIGREAL, normr0 = BIGREAL;
1290 REAL reldiff, factor, infnormu;
1291 REAL alpha, beta, temp1, temp2;
1292
1293 // allocate temp memory (need 4*m REAL numbers)
1294 REAL *work=(REAL *)fasp_mem_calloc(4*m,sizeof(REAL));
1295 REAL *p=work, *z=work+m, *r=z+m, *t=r+m;
1296
1297 // Output some info for debuging
1298 if ( PrtLvl > PRINT_NONE ) printf("\nCalling CG solver (MatFree) ...\n");
1299
1300#if DEBUG_MODE > 0
1301 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
1302 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1303#endif
1304
1305 // initialize counters
1306 stag=1; more_step=1; restart_step=1;
1307
1308 // r = b-A*u
1309 mf->fct(mf->data, u->val, r);
1310 fasp_blas_darray_axpby(m, 1.0, b->val, -1.0, r);
1311
1312 if (pc != NULL)
1313 pc->fct(r,z,pc->data); /* Apply preconditioner */
1314 else
1315 fasp_darray_cp(m,r,z); /* No preconditioner */
1316
1317 // compute initial relative residual
1318 switch (StopType) {
1319 case STOP_REL_PRECRES:
1320 absres0=sqrt(fasp_blas_darray_dotprod(m,r,z));
1321 normr0=MAX(SMALLREAL,absres0);
1322 relres=absres0/normr0;
1323 break;
1324 case STOP_MOD_REL_RES:
1325 absres0=fasp_blas_darray_norm2(m,r);
1327 relres=absres0/normu;
1328 break;
1329 default:
1330 absres0=fasp_blas_darray_norm2(m,r);
1331 normr0=MAX(SMALLREAL,absres0);
1332 relres=absres0/normr0;
1333 break;
1334 }
1335
1336 // if initial residual is small, no need to iterate!
1337 if ( relres < tol || absres0 < 1e-12*tol ) goto FINISHED;
1338
1339 fasp_darray_cp(m,z,p);
1340 temp1=fasp_blas_darray_dotprod(m,z,r);
1341
1342 while ( iter++ < MaxIt ) {
1343
1344 // t=A*p
1345 mf->fct(mf->data, p, t);
1346
1347 // alpha_k=(z_{k-1},r_{k-1})/(A*p_{k-1},p_{k-1})
1348 temp2=fasp_blas_darray_dotprod(m,t,p);
1349 alpha=temp1/temp2;
1350
1351 // u_k=u_{k-1} + alpha_k*p_{k-1}
1352 fasp_blas_darray_axpy(m,alpha,p,u->val);
1353
1354 // r_k=r_{k-1} - alpha_k*A*p_{k-1}
1355 fasp_blas_darray_axpy(m,-alpha,t,r);
1356 absres=fasp_blas_darray_norm2(m,r);
1357
1358 // compute reducation factor of residual ||r||
1359 factor=absres/absres0;
1360
1361 // compute relative residual
1362 switch (StopType) {
1363 case STOP_REL_PRECRES:
1364 // z = B(r)
1365 if (pc != NULL)
1366 pc->fct(r,z,pc->data); /* Apply preconditioner */
1367 else
1368 fasp_darray_cp(m,r,z); /* No preconditioner */
1369 temp2=fasp_blas_darray_dotprod(m,z,r);
1370 relres=sqrt(ABS(temp2))/normr0;
1371 break;
1372 case STOP_MOD_REL_RES:
1373 relres=absres/normu;
1374 break;
1375 default:
1376 relres=absres/normr0;
1377 break;
1378 }
1379
1380 // output iteration information if needed
1381 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1382
1383 // solution check, if soultion is too small, return ERROR_SOLVER_SOLSTAG.
1384 infnormu = fasp_blas_darray_norminf(m, u->val);
1385 if ( infnormu <= sol_inf_tol ) {
1386 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
1387 iter = ERROR_SOLVER_SOLSTAG;
1388 break;
1389 }
1390
1391 // compute relative difference
1392 normu = fasp_blas_darray_norm2(m, u->val);
1393 reldiff = ABS(alpha)*fasp_blas_darray_norm2(m, p)/normu;
1394
1395 // stagnation check
1396 if ( (stag<=MaxStag) & (reldiff<maxdiff) ) {
1397
1398 if ( PrtLvl >= PRINT_MORE ) {
1399 ITS_DIFFRES(reldiff,relres);
1400 ITS_RESTART;
1401 }
1402
1403 mf->fct(mf->data, u->val, r);
1404 fasp_blas_darray_axpby(m, 1.0, b->val, -1.0, r);
1405 absres = fasp_blas_darray_norm2(m,r);
1406
1407 // relative residual
1408 switch (StopType) {
1409 case STOP_REL_PRECRES:
1410 // z = B(r)
1411 if (pc != NULL)
1412 pc->fct(r,z,pc->data); /* Apply preconditioner */
1413 else
1414 fasp_darray_cp(m,r,z); /* No preconditioner */
1415 temp2=fasp_blas_darray_dotprod(m,z,r);
1416 relres=sqrt(ABS(temp2))/normr0;
1417 break;
1418 case STOP_MOD_REL_RES:
1419 relres=absres/normu;
1420 break;
1421 default:
1422 relres=absres/normr0;
1423 break;
1424 }
1425
1426 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1427
1428 if ( relres < tol )
1429 break;
1430 else {
1431 if ( stag >= MaxStag ) {
1432 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
1433 iter = ERROR_SOLVER_STAG;
1434 break;
1435 }
1436 fasp_darray_set(m,p,0.0);
1437 ++stag;
1438 ++restart_step;
1439 }
1440 } // end of stagnation check!
1441
1442 // safe-guard check
1443 if ( relres < tol ) {
1444 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
1445
1446 mf->fct(mf->data, u->val, r);
1447 fasp_blas_darray_axpby(m, 1.0, b->val, -1.0, r);
1448
1449 // relative residual
1450 switch (StopType) {
1451 case STOP_REL_PRECRES:
1452 // z = B(r)
1453 if (pc != NULL)
1454 pc->fct(r,z,pc->data); /* Apply preconditioner */
1455 else
1456 fasp_darray_cp(m,r,z); /* No preconditioner */
1457 temp2=fasp_blas_darray_dotprod(m,z,r);
1458 relres=sqrt(ABS(temp2))/normr0;
1459 break;
1460 case STOP_MOD_REL_RES:
1461 absres=fasp_blas_darray_norm2(m,r);
1462 relres=absres/normu;
1463 break;
1464 default:
1465 absres=fasp_blas_darray_norm2(m,r);
1466 relres=absres/normr0;
1467 break;
1468 }
1469
1470 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1471
1472 // check convergence
1473 if ( relres < tol ) break;
1474
1475 if ( more_step >= MaxRestartStep ) {
1476 if ( PrtLvl > PRINT_MIN) ITS_ZEROTOL;
1477 iter = ERROR_SOLVER_TOLSMALL;
1478 break;
1479 }
1480
1481 // prepare for restarting method
1482 fasp_darray_set(m,p,0.0);
1483 ++more_step;
1484 ++restart_step;
1485
1486 } // end of safe-guard check!
1487
1488 // update relative residual here
1489 absres0 = absres;
1490
1491 // compute z_k = B(r_k)
1492 if ( StopType != STOP_REL_PRECRES ) {
1493 if ( pc != NULL )
1494 pc->fct(r,z,pc->data); /* Apply preconditioner */
1495 else
1496 fasp_darray_cp(m,r,z); /* No preconditioner, B=I */
1497 }
1498
1499 // compute beta_k = (z_k, r_k)/(z_{k-1}, r_{k-1})
1500 temp2=fasp_blas_darray_dotprod(m,z,r);
1501 beta=temp2/temp1;
1502 temp1=temp2;
1503
1504 // compute p_k = z_k + beta_k*p_{k-1}
1505 fasp_blas_darray_axpby(m,1.0,z,beta,p);
1506
1507 } // end of main PCG loop.
1508
1509FINISHED: // finish iterative method
1510 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1511
1512 // clean up temp memory
1513 fasp_mem_free(work); work = NULL;
1514
1515#if DEBUG_MODE > 0
1516 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1517#endif
1518
1519 if (iter>MaxIt)
1520 return ERROR_SOLVER_MAXIT;
1521 else
1522 return iter;
1523}
1524
1525/*---------------------------------*/
1526/*-- End of File --*/
1527/*---------------------------------*/
void fasp_darray_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: AuxArray.c:41
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
REAL fasp_blas_darray_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: BlaArray.c:686
void fasp_blas_darray_axpby(const INT n, const REAL a, const REAL *x, const REAL b, REAL *y)
y = a*x + b*y
Definition: BlaArray.c:580
REAL fasp_blas_darray_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: BlaArray.c:657
void fasp_blas_darray_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_dbsr_pcg(dBSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:390
INT fasp_solver_dstr_pcg(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:978
INT fasp_solver_dcsr_pcg(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:98
INT fasp_solver_dblc_pcg(dBLCmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient method for solving Au=b.
Definition: KryPcg.c:684
INT fasp_solver_pcg(mxv_matfree *mf, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT StopType, const SHORT PrtLvl)
Preconditioned conjugate gradient (CG) method for solving Au=b.
Definition: KryPcg.c:1272
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 MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:74
#define INT
Definition: fasp.h:64
#define STAG_RATIO
Definition: fasp_const.h:258
#define ERROR_SOLVER_STAG
Definition: fasp_const.h:43
#define SMALLREAL2
Definition: fasp_const.h:250
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:48
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:131
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
#define ERROR_SOLVER_SOLSTAG
Definition: fasp_const.h:44
#define MAX_STAG
Definition: fasp_const.h:257
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
#define BIGREAL
Some global constants.
Definition: fasp_const.h:248
#define MAX_RESTART
Definition: fasp_const.h:256
#define PRINT_MORE
Definition: fasp_const.h:76
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:73
#define ERROR_SOLVER_TOLSMALL
Definition: fasp_const.h:45
#define SMALLREAL
Definition: fasp_const.h:249
#define PRINT_MIN
Definition: fasp_const.h:74
Block REAL CSR matrix format.
Definition: fasp_block.h:74
Block sparse row storage matrix of REAL type.
Definition: fasp_block.h:34
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:143
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