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