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