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