Fast Auxiliary Space Preconditioning 2.7.7 Aug/28/2022
Loading...
Searching...
No Matches
KrySPminres.c
Go to the documentation of this file.
1
24#include <math.h>
25
26#include "fasp.h"
27#include "fasp_functs.h"
28
29/*---------------------------------*/
30/*-- Declare Private Functions --*/
31/*---------------------------------*/
32
33#include "KryUtil.inl"
34
35/*---------------------------------*/
36/*-- Public Functions --*/
37/*---------------------------------*/
38
60INT fasp_solver_dcsr_spminres (const dCSRmat *A,
61 const dvector *b,
62 dvector *u,
63 precond *pc,
64 const REAL tol,
65 const INT MaxIt,
66 const SHORT StopType,
67 const SHORT PrtLvl)
68{
69 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
70 const INT m = b->row;
71 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
72 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
73
74 // local variables
75 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
76 REAL absres0 = BIGREAL, absres = BIGREAL;
77 REAL normr0 = BIGREAL, relres = BIGREAL;
78 REAL normu2, normuu, normp, normuinf, factor;
79 REAL alpha, alpha0, alpha1, temp2;
80 INT iter_best = 0; // initial best known iteration
81 REAL absres_best = BIGREAL; // initial best known residual
82
83 // allocate temp memory (need 12*m REAL)
84 REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
85 REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
86 REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
87
88 // Output some info for debuging
89 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe 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 ) 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 // safety net check: save the best-so-far solution
248 if ( fasp_dvec_isnan(u) ) {
249 // If the solution is NAN, restore the best solution
250 absres = BIGREAL;
251 goto RESTORE_BESTSOL;
252 }
253
254 if ( absres < absres_best - maxdiff) {
255 absres_best = absres;
256 iter_best = iter;
257 fasp_darray_cp(m,u->val,u_best);
258 }
259
260 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
261 normuinf = fasp_blas_darray_norminf(m, u->val);
262 if (normuinf <= sol_inf_tol) {
263 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
264 iter = ERROR_SOLVER_SOLSTAG;
265 break;
266 }
267
268 // Check II: if staggenated, try to restart
269 normuu = fasp_blas_darray_norm2(m,p1);
270 normuu = ABS(alpha)*(normuu/normu2);
271
272 if ( normuu < maxdiff ) {
273
274 if ( stag < MaxStag ) {
275 if ( PrtLvl >= PRINT_MORE ) {
276 ITS_DIFFRES(normuu,relres);
277 ITS_RESTART;
278 }
279 }
280
281 fasp_darray_cp(m,b->val,r);
282 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
283
284 // compute residuals
285 switch (StopType) {
286 case STOP_REL_RES:
287 temp2 = fasp_blas_darray_dotprod(m,r,r);
288 absres = sqrt(temp2);
289 relres = absres/normr0;
290 break;
291 case STOP_REL_PRECRES:
292 if (pc == NULL)
293 fasp_darray_cp(m,r,t);
294 else
295 pc->fct(r,t,pc->data);
296 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
297 absres = sqrt(temp2);
298 relres = absres/normr0;
299 break;
300 case STOP_MOD_REL_RES:
301 temp2 = fasp_blas_darray_dotprod(m,r,r);
302 absres = sqrt(temp2);
303 relres = absres/normu2;
304 break;
305 }
306
307 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
308
309 if ( relres < tol )
310 break;
311 else {
312 if ( stag >= MaxStag ) {
313 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
314 iter = ERROR_SOLVER_STAG;
315 break;
316 }
317 fasp_darray_set(m,p0,0.0);
318 ++stag;
319 ++restart_step;
320
321 // p1 = B(r)
322 if ( pc != NULL )
323 pc->fct(r,p1,pc->data); /* Apply preconditioner */
324 else
325 fasp_darray_cp(m,r,p1); /* No preconditioner */
326
327 // tp = A*p1
328 fasp_blas_dcsr_mxv(A,p1,tp);
329
330 // tz = B(tp)
331 if ( pc != NULL )
332 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
333 else
334 fasp_darray_cp(m,tp,tz); /* No preconditioner */
335
336 // p1 = p1/normp
337 normp = fasp_blas_darray_dotprod(m,tz,tp);
338 normp = sqrt(normp);
339 fasp_darray_cp(m,p1,t);
340
341 // t0 = A*p0=0
342 fasp_darray_set(m,t0,0.0);
343 fasp_darray_cp(m,t0,z0);
344 fasp_darray_cp(m,t0,t1);
345 fasp_darray_cp(m,t0,z1);
346 fasp_darray_cp(m,t0,p1);
347
348 fasp_blas_darray_axpy(m,1/normp,t,p1);
349
350 // t1 = tp/normp, z1 = tz/normp
351 fasp_blas_darray_axpy(m,1/normp,tp,t1);
352 fasp_blas_darray_axpy(m,1/normp,tz,z1);
353 }
354 }
355
356 // Check III: prevent false convergence
357 if ( relres < tol ) {
358
359 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
360
361 // compute residual r = b - Ax again
362 fasp_darray_cp(m,b->val,r);
363 fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
364
365 // compute residuals
366 switch (StopType) {
367 case STOP_REL_RES:
368 temp2 = fasp_blas_darray_dotprod(m,r,r);
369 absres = sqrt(temp2);
370 relres = absres/normr0;
371 break;
372 case STOP_REL_PRECRES:
373 if (pc == NULL)
374 fasp_darray_cp(m,r,t);
375 else
376 pc->fct(r,t,pc->data);
377 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
378 absres = sqrt(temp2);
379 relres = absres/normr0;
380 break;
381 case STOP_MOD_REL_RES:
382 temp2 = fasp_blas_darray_dotprod(m,r,r);
383 absres = sqrt(temp2);
384 relres = absres/normu2;
385 break;
386 }
387
388 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
389
390 // check convergence
391 if ( relres < tol ) break;
392
393 if ( more_step >= MaxRestartStep ) {
394 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
395 iter = ERROR_SOLVER_TOLSMALL;
396 break;
397 }
398
399 // prepare for restarting method
400 fasp_darray_set(m,p0,0.0);
401 ++more_step;
402 ++restart_step;
403
404 // p1 = B(r)
405 if ( pc != NULL )
406 pc->fct(r,p1,pc->data); /* Apply preconditioner */
407 else
408 fasp_darray_cp(m,r,p1); /* No preconditioner */
409
410 // tp = A*p1
411 fasp_blas_dcsr_mxv(A,p1,tp);
412
413 // tz = B(tp)
414 if ( pc != NULL )
415 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
416 else
417 fasp_darray_cp(m,tp,tz); /* No preconditioner */
418
419 // p1 = p1/normp
420 normp = fasp_blas_darray_dotprod(m,tz,tp);
421 normp = sqrt(normp);
422 fasp_darray_cp(m,p1,t);
423
424 // t0 = A*p0 = 0
425 fasp_darray_set(m,t0,0.0);
426 fasp_darray_cp(m,t0,z0);
427 fasp_darray_cp(m,t0,t1);
428 fasp_darray_cp(m,t0,z1);
429 fasp_darray_cp(m,t0,p1);
430
431 fasp_blas_darray_axpy(m,1/normp,t,p1);
432
433 // t1=tp/normp,z1=tz/normp
434 fasp_blas_darray_axpy(m,1/normp,tp,t1);
435 fasp_blas_darray_axpy(m,1/normp,tz,z1);
436
437 } // end of convergence check
438
439 // update relative residual here
440 absres0 = absres;
441
442 } // end of the main loop
443
444RESTORE_BESTSOL: // restore the best-so-far solution if necessary
445 if ( iter != iter_best ) {
446
447 // compute best residual
448 fasp_darray_cp(m,b->val,r);
449 fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
450
451 switch ( StopType ) {
452 case STOP_REL_RES:
453 absres_best = fasp_blas_darray_norm2(m,r);
454 break;
455 case STOP_REL_PRECRES:
456 if ( pc != NULL )
457 pc->fct(r,t,pc->data); /* Apply preconditioner */
458 else
459 fasp_darray_cp(m,r,t); /* No preconditioner */
460 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,t,r)));
461 break;
462 case STOP_MOD_REL_RES:
463 absres_best = fasp_blas_darray_norm2(m,r);
464 break;
465 }
466
467 if ( absres > absres_best + maxdiff || isnan(absres) ) {
468 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
469 fasp_darray_cp(m,u_best,u->val);
470 relres = absres_best / normr0;
471 }
472 }
473
474FINISHED: // finish iterative method
475 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
476
477 // clean up temp memory
478 fasp_mem_free(work); work = NULL;
479
480#if DEBUG_MODE > 0
481 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
482#endif
483
484 if ( iter > MaxIt )
485 return ERROR_SOLVER_MAXIT;
486 else
487 return iter;
488}
489
511INT fasp_solver_dblc_spminres (const dBLCmat *A,
512 const dvector *b,
513 dvector *u,
514 precond *pc,
515 const REAL tol,
516 const INT MaxIt,
517 const SHORT StopType,
518 const SHORT PrtLvl)
519{
520 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
521 const INT m = b->row;
522 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
523 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
524
525 // local variables
526 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
527 REAL absres0 = BIGREAL, absres = BIGREAL;
528 REAL normr0 = BIGREAL, relres = BIGREAL;
529 REAL normu2, normuu, normp, normuinf, factor;
530 REAL alpha, alpha0, alpha1, temp2;
531 INT iter_best = 0; // initial best known iteration
532 REAL absres_best = BIGREAL; // initial best known residual
533
534 // allocate temp memory (need 12*m REAL)
535 REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
536 REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
537 REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
538
539 // Output some info for debuging
540 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe MinRes solver (BLC) ...\n");
541
542#if DEBUG_MODE > 0
543 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
544 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
545#endif
546
547 // p0 = 0
548 fasp_darray_set(m,p0,0.0);
549
550 // r = b-A*u
551 fasp_darray_cp(m,b->val,r);
552 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
553
554 // p1 = B(r)
555 if ( pc != NULL )
556 pc->fct(r,p1,pc->data); /* Apply preconditioner */
557 else
558 fasp_darray_cp(m,r,p1); /* No preconditioner */
559
560 // compute initial residuals
561 switch ( StopType ) {
562 case STOP_REL_RES:
563 absres0 = fasp_blas_darray_norm2(m,r);
564 normr0 = MAX(SMALLREAL,absres0);
565 relres = absres0/normr0;
566 break;
567 case STOP_REL_PRECRES:
568 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,p1));
569 normr0 = MAX(SMALLREAL,absres0);
570 relres = absres0/normr0;
571 break;
572 case STOP_MOD_REL_RES:
573 absres0 = fasp_blas_darray_norm2(m,r);
574 normu2 = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
575 relres = absres0/normu2;
576 break;
577 default:
578 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
579 goto FINISHED;
580 }
581
582 // if initial residual is small, no need to iterate!
583 if ( relres < tol ) goto FINISHED;
584
585 // output iteration information if needed
586 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
587
588 // tp = A*p1
589 fasp_blas_dblc_mxv(A,p1,tp);
590
591 // tz = B(tp)
592 if ( pc != NULL )
593 pc->fct(tp,tz,pc->data); /* Apply preconditioner */
594 else
595 fasp_darray_cp(m,tp,tz); /* No preconditioner */
596
597 // p1 = p1/normp
598 normp = ABS(fasp_blas_darray_dotprod(m,tz,tp));
599 normp = sqrt(normp);
600 fasp_darray_cp(m,p1,t);
601 fasp_darray_set(m,p1,0.0);
602 fasp_blas_darray_axpy(m,1/normp,t,p1);
603
604 // t0 = A*p0 = 0
605 fasp_darray_set(m,t0,0.0);
606 fasp_darray_cp(m,t0,z0);
607 fasp_darray_cp(m,t0,t1);
608 fasp_darray_cp(m,t0,z1);
609
610 // t1 = tp/normp, z1 = tz/normp
611 fasp_blas_darray_axpy(m,1.0/normp,tp,t1);
612 fasp_blas_darray_axpy(m,1.0/normp,tz,z1);
613
614 // main MinRes loop
615 while ( iter++ < MaxIt ) {
616
617 // alpha = <r,z1>
618 alpha=fasp_blas_darray_dotprod(m,r,z1);
619
620 // u = u+alpha*p1
621 fasp_blas_darray_axpy(m,alpha,p1,u->val);
622
623 // r = r-alpha*Ap1
624 fasp_blas_darray_axpy(m,-alpha,t1,r);
625
626 // compute t = A*z1 alpha1 = <z1,t>
627 fasp_blas_dblc_mxv(A,z1,t);
628 alpha1=fasp_blas_darray_dotprod(m,z1,t);
629
630 // compute t = A*z0 alpha0 = <z1,t>
631 fasp_blas_dblc_mxv(A,z0,t);
632 alpha0=fasp_blas_darray_dotprod(m,z1,t);
633
634 // p2 = z1-alpha1*p1-alpha0*p0
635 fasp_darray_cp(m,z1,p2);
636 fasp_blas_darray_axpy(m,-alpha1,p1,p2);
637 fasp_blas_darray_axpy(m,-alpha0,p0,p2);
638
639 // tp = A*p2
640 fasp_blas_dblc_mxv(A,p2,tp);
641
642 // tz = B(tp)
643 if ( pc != NULL )
644 pc->fct(tp,tz,pc->data); /* Apply preconditioner */
645 else
646 fasp_darray_cp(m,tp,tz); /* No preconditioner */
647
648 // p2 = p2/normp
649 normp = ABS(fasp_blas_darray_dotprod(m,tz,tp));
650 normp = sqrt(normp);
651 fasp_darray_cp(m,p2,t);
652 fasp_darray_set(m,p2,0.0);
653 fasp_blas_darray_axpy(m,1/normp,t,p2);
654
655 // prepare for next iteration
656 fasp_darray_cp(m,p1,p0);
657 fasp_darray_cp(m,p2,p1);
658 fasp_darray_cp(m,t1,t0);
659 fasp_darray_cp(m,z1,z0);
660
661 // t1=tp/normp,z1=tz/normp
662 fasp_darray_set(m,t1,0.0);
663 fasp_darray_cp(m,t1,z1);
664 fasp_blas_darray_axpy(m,1/normp,tp,t1);
665 fasp_blas_darray_axpy(m,1/normp,tz,z1);
666
667 normu2 = fasp_blas_darray_norm2(m,u->val);
668
669 // compute residuals
670 switch ( StopType ) {
671 case STOP_REL_RES:
672 temp2 = fasp_blas_darray_dotprod(m,r,r);
673 absres = sqrt(temp2);
674 relres = absres/normr0;
675 break;
676 case STOP_REL_PRECRES:
677 if (pc == NULL)
678 fasp_darray_cp(m,r,t);
679 else
680 pc->fct(r,t,pc->data);
681 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
682 absres = sqrt(temp2);
683 relres = absres/normr0;
684 break;
685 case STOP_MOD_REL_RES:
686 temp2 = fasp_blas_darray_dotprod(m,r,r);
687 absres = sqrt(temp2);
688 relres = absres/normu2;
689 break;
690 }
691
692 // compute reducation factor of residual ||r||
693 factor = absres/absres0;
694
695 // output iteration information if needed
696 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
697
698 // safety net check: save the best-so-far solution
699 if ( fasp_dvec_isnan(u) ) {
700 // If the solution is NAN, restore the best solution
701 absres = BIGREAL;
702 goto RESTORE_BESTSOL;
703 }
704
705 if ( absres < absres_best - maxdiff) {
706 absres_best = absres;
707 iter_best = iter;
708 fasp_darray_cp(m,u->val,u_best);
709 }
710
711 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
712 normuinf = fasp_blas_darray_norminf(m, u->val);
713 if (normuinf <= sol_inf_tol) {
714 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
715 iter = ERROR_SOLVER_SOLSTAG;
716 break;
717 }
718
719 // Check II: if staggenated, try to restart
720 normuu = fasp_blas_darray_norm2(m,p1);
721 normuu = ABS(alpha)*(normuu/normu2);
722
723 if ( normuu < maxdiff ) {
724
725 if ( stag < MaxStag ) {
726 if ( PrtLvl >= PRINT_MORE ) {
727 ITS_DIFFRES(normuu,relres);
728 ITS_RESTART;
729 }
730 }
731
732 fasp_darray_cp(m,b->val,r);
733 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
734
735 // compute residuals
736 switch (StopType) {
737 case STOP_REL_RES:
738 temp2 = fasp_blas_darray_dotprod(m,r,r);
739 absres = sqrt(temp2);
740 relres = absres/normr0;
741 break;
742 case STOP_REL_PRECRES:
743 if (pc == NULL)
744 fasp_darray_cp(m,r,t);
745 else
746 pc->fct(r,t,pc->data);
747 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
748 absres = sqrt(temp2);
749 relres = absres/normr0;
750 break;
751 case STOP_MOD_REL_RES:
752 temp2 = fasp_blas_darray_dotprod(m,r,r);
753 absres = sqrt(temp2);
754 relres = absres/normu2;
755 break;
756 }
757
758 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
759
760 if ( relres < tol )
761 break;
762 else {
763 if ( stag >= MaxStag ) {
764 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
765 iter = ERROR_SOLVER_STAG;
766 break;
767 }
768 fasp_darray_set(m,p0,0.0);
769 ++stag;
770 ++restart_step;
771
772 // p1 = B(r)
773 if ( pc != NULL )
774 pc->fct(r,p1,pc->data); /* Apply preconditioner */
775 else
776 fasp_darray_cp(m,r,p1); /* No preconditioner */
777
778 // tp = A*p1
779 fasp_blas_dblc_mxv(A,p1,tp);
780
781 // tz = B(tp)
782 if ( pc != NULL )
783 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
784 else
785 fasp_darray_cp(m,tp,tz); /* No preconditioner */
786
787 // p1 = p1/normp
788 normp = fasp_blas_darray_dotprod(m,tz,tp);
789 normp = sqrt(normp);
790 fasp_darray_cp(m,p1,t);
791
792 // t0 = A*p0=0
793 fasp_darray_set(m,t0,0.0);
794 fasp_darray_cp(m,t0,z0);
795 fasp_darray_cp(m,t0,t1);
796 fasp_darray_cp(m,t0,z1);
797 fasp_darray_cp(m,t0,p1);
798
799 fasp_blas_darray_axpy(m,1/normp,t,p1);
800
801 // t1 = tp/normp, z1 = tz/normp
802 fasp_blas_darray_axpy(m,1/normp,tp,t1);
803 fasp_blas_darray_axpy(m,1/normp,tz,z1);
804 }
805 }
806
807 // Check III: prevent false convergence
808 if ( relres < tol ) {
809
810 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
811
812 // compute residual r = b - Ax again
813 fasp_darray_cp(m,b->val,r);
814 fasp_blas_dblc_aAxpy(-1.0,A,u->val,r);
815
816 // compute residuals
817 switch (StopType) {
818 case STOP_REL_RES:
819 temp2 = fasp_blas_darray_dotprod(m,r,r);
820 absres = sqrt(temp2);
821 relres = absres/normr0;
822 break;
823 case STOP_REL_PRECRES:
824 if (pc == NULL)
825 fasp_darray_cp(m,r,t);
826 else
827 pc->fct(r,t,pc->data);
828 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
829 absres = sqrt(temp2);
830 relres = absres/normr0;
831 break;
832 case STOP_MOD_REL_RES:
833 temp2 = fasp_blas_darray_dotprod(m,r,r);
834 absres = sqrt(temp2);
835 relres = absres/normu2;
836 break;
837 }
838
839 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
840
841 // check convergence
842 if ( relres < tol ) break;
843
844 if ( more_step >= MaxRestartStep ) {
845 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
846 iter = ERROR_SOLVER_TOLSMALL;
847 break;
848 }
849
850 // prepare for restarting method
851 fasp_darray_set(m,p0,0.0);
852 ++more_step;
853 ++restart_step;
854
855 // p1 = B(r)
856 if ( pc != NULL )
857 pc->fct(r,p1,pc->data); /* Apply preconditioner */
858 else
859 fasp_darray_cp(m,r,p1); /* No preconditioner */
860
861 // tp = A*p1
862 fasp_blas_dblc_mxv(A,p1,tp);
863
864 // tz = B(tp)
865 if ( pc != NULL )
866 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
867 else
868 fasp_darray_cp(m,tp,tz); /* No preconditioner */
869
870 // p1 = p1/normp
871 normp = fasp_blas_darray_dotprod(m,tz,tp);
872 normp = sqrt(normp);
873 fasp_darray_cp(m,p1,t);
874
875 // t0 = A*p0 = 0
876 fasp_darray_set(m,t0,0.0);
877 fasp_darray_cp(m,t0,z0);
878 fasp_darray_cp(m,t0,t1);
879 fasp_darray_cp(m,t0,z1);
880 fasp_darray_cp(m,t0,p1);
881
882 fasp_blas_darray_axpy(m,1/normp,t,p1);
883
884 // t1=tp/normp,z1=tz/normp
885 fasp_blas_darray_axpy(m,1/normp,tp,t1);
886 fasp_blas_darray_axpy(m,1/normp,tz,z1);
887
888 } // end of convergence check
889
890 // update relative residual here
891 absres0 = absres;
892
893 } // end of the main loop
894
895RESTORE_BESTSOL: // restore the best-so-far solution if necessary
896 if ( iter != iter_best ) {
897
898 // compute best residual
899 fasp_darray_cp(m,b->val,r);
900 fasp_blas_dblc_aAxpy(-1.0,A,u_best,r);
901
902 switch ( StopType ) {
903 case STOP_REL_RES:
904 absres_best = fasp_blas_darray_norm2(m,r);
905 break;
906 case STOP_REL_PRECRES:
907 if ( pc != NULL )
908 pc->fct(r,t,pc->data); /* Apply preconditioner */
909 else
910 fasp_darray_cp(m,r,t); /* No preconditioner */
911 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,t,r)));
912 break;
913 case STOP_MOD_REL_RES:
914 absres_best = fasp_blas_darray_norm2(m,r);
915 break;
916 }
917
918 if ( absres > absres_best + maxdiff || isnan(absres) ) {
919 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
920 fasp_darray_cp(m,u_best,u->val);
921 relres = absres_best / normr0;
922 }
923 }
924
925FINISHED: // finish iterative method
926 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
927
928 // clean up temp memory
929 fasp_mem_free(work); work = NULL;
930
931#if DEBUG_MODE > 0
932 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
933#endif
934
935 if ( iter > MaxIt )
936 return ERROR_SOLVER_MAXIT;
937 else
938 return iter;
939}
940
962INT fasp_solver_dstr_spminres (const dSTRmat *A,
963 const dvector *b,
964 dvector *u,
965 precond *pc,
966 const REAL tol,
967 const INT MaxIt,
968 const SHORT StopType,
969 const SHORT PrtLvl)
970{
971 const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
972 const INT m = b->row;
973 const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
974 const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
975
976 // local variables
977 INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
978 REAL absres0 = BIGREAL, absres = BIGREAL;
979 REAL normr0 = BIGREAL, relres = BIGREAL;
980 REAL normu2, normuu, normp, normuinf, factor;
981 REAL alpha, alpha0, alpha1, temp2;
982 INT iter_best = 0; // initial best known iteration
983 REAL absres_best = BIGREAL; // initial best known residual
984
985 // allocate temp memory (need 12*m REAL)
986 REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
987 REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
988 REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
989
990 // Output some info for debuging
991 if ( PrtLvl > PRINT_NONE ) printf("\nCalling Safe MinRes solver (STR) ...\n");
992
993#if DEBUG_MODE > 0
994 printf("### DEBUG: [-Begin-] %s ...\n", __FUNCTION__);
995 printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
996#endif
997
998 // p0 = 0
999 fasp_darray_set(m,p0,0.0);
1000
1001 // r = b-A*u
1002 fasp_darray_cp(m,b->val,r);
1003 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1004
1005 // p1 = B(r)
1006 if ( pc != NULL )
1007 pc->fct(r,p1,pc->data); /* Apply preconditioner */
1008 else
1009 fasp_darray_cp(m,r,p1); /* No preconditioner */
1010
1011 // compute initial residuals
1012 switch ( StopType ) {
1013 case STOP_REL_RES:
1014 absres0 = fasp_blas_darray_norm2(m,r);
1015 normr0 = MAX(SMALLREAL,absres0);
1016 relres = absres0/normr0;
1017 break;
1018 case STOP_REL_PRECRES:
1019 absres0 = sqrt(fasp_blas_darray_dotprod(m,r,p1));
1020 normr0 = MAX(SMALLREAL,absres0);
1021 relres = absres0/normr0;
1022 break;
1023 case STOP_MOD_REL_RES:
1024 absres0 = fasp_blas_darray_norm2(m,r);
1025 normu2 = MAX(SMALLREAL,fasp_blas_darray_norm2(m,u->val));
1026 relres = absres0/normu2;
1027 break;
1028 default:
1029 printf("### ERROR: Unknown stopping type! [%s]\n", __FUNCTION__);
1030 goto FINISHED;
1031 }
1032
1033 // if initial residual is small, no need to iterate!
1034 if ( relres < tol ) goto FINISHED;
1035
1036 // output iteration information if needed
1037 fasp_itinfo(PrtLvl,StopType,iter,relres,absres0,0.0);
1038
1039 // tp = A*p1
1040 fasp_blas_dstr_mxv(A,p1,tp);
1041
1042 // tz = B(tp)
1043 if ( pc != NULL )
1044 pc->fct(tp,tz,pc->data); /* Apply preconditioner */
1045 else
1046 fasp_darray_cp(m,tp,tz); /* No preconditioner */
1047
1048 // p1 = p1/normp
1049 normp = ABS(fasp_blas_darray_dotprod(m,tz,tp));
1050 normp = sqrt(normp);
1051 fasp_darray_cp(m,p1,t);
1052 fasp_darray_set(m,p1,0.0);
1053 fasp_blas_darray_axpy(m,1/normp,t,p1);
1054
1055 // t0 = A*p0 = 0
1056 fasp_darray_set(m,t0,0.0);
1057 fasp_darray_cp(m,t0,z0);
1058 fasp_darray_cp(m,t0,t1);
1059 fasp_darray_cp(m,t0,z1);
1060
1061 // t1 = tp/normp, z1 = tz/normp
1062 fasp_blas_darray_axpy(m,1.0/normp,tp,t1);
1063 fasp_blas_darray_axpy(m,1.0/normp,tz,z1);
1064
1065 // main MinRes loop
1066 while ( iter++ < MaxIt ) {
1067
1068 // alpha = <r,z1>
1069 alpha=fasp_blas_darray_dotprod(m,r,z1);
1070
1071 // u = u+alpha*p1
1072 fasp_blas_darray_axpy(m,alpha,p1,u->val);
1073
1074 // r = r-alpha*Ap1
1075 fasp_blas_darray_axpy(m,-alpha,t1,r);
1076
1077 // compute t = A*z1 alpha1 = <z1,t>
1078 fasp_blas_dstr_mxv(A,z1,t);
1079 alpha1=fasp_blas_darray_dotprod(m,z1,t);
1080
1081 // compute t = A*z0 alpha0 = <z1,t>
1082 fasp_blas_dstr_mxv(A,z0,t);
1083 alpha0=fasp_blas_darray_dotprod(m,z1,t);
1084
1085 // p2 = z1-alpha1*p1-alpha0*p0
1086 fasp_darray_cp(m,z1,p2);
1087 fasp_blas_darray_axpy(m,-alpha1,p1,p2);
1088 fasp_blas_darray_axpy(m,-alpha0,p0,p2);
1089
1090 // tp = A*p2
1091 fasp_blas_dstr_mxv(A,p2,tp);
1092
1093 // tz = B(tp)
1094 if ( pc != NULL )
1095 pc->fct(tp,tz,pc->data); /* Apply preconditioner */
1096 else
1097 fasp_darray_cp(m,tp,tz); /* No preconditioner */
1098
1099 // p2 = p2/normp
1100 normp = ABS(fasp_blas_darray_dotprod(m,tz,tp));
1101 normp = sqrt(normp);
1102 fasp_darray_cp(m,p2,t);
1103 fasp_darray_set(m,p2,0.0);
1104 fasp_blas_darray_axpy(m,1/normp,t,p2);
1105
1106 // prepare for next iteration
1107 fasp_darray_cp(m,p1,p0);
1108 fasp_darray_cp(m,p2,p1);
1109 fasp_darray_cp(m,t1,t0);
1110 fasp_darray_cp(m,z1,z0);
1111
1112 // t1=tp/normp,z1=tz/normp
1113 fasp_darray_set(m,t1,0.0);
1114 fasp_darray_cp(m,t1,z1);
1115 fasp_blas_darray_axpy(m,1/normp,tp,t1);
1116 fasp_blas_darray_axpy(m,1/normp,tz,z1);
1117
1118 normu2 = fasp_blas_darray_norm2(m,u->val);
1119
1120 // compute residuals
1121 switch ( StopType ) {
1122 case STOP_REL_RES:
1123 temp2 = fasp_blas_darray_dotprod(m,r,r);
1124 absres = sqrt(temp2);
1125 relres = absres/normr0;
1126 break;
1127 case STOP_REL_PRECRES:
1128 if (pc == NULL)
1129 fasp_darray_cp(m,r,t);
1130 else
1131 pc->fct(r,t,pc->data);
1132 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
1133 absres = sqrt(temp2);
1134 relres = absres/normr0;
1135 break;
1136 case STOP_MOD_REL_RES:
1137 temp2 = fasp_blas_darray_dotprod(m,r,r);
1138 absres = sqrt(temp2);
1139 relres = absres/normu2;
1140 break;
1141 }
1142
1143 // compute reducation factor of residual ||r||
1144 factor = absres/absres0;
1145
1146 // output iteration information if needed
1147 fasp_itinfo(PrtLvl,StopType,iter,relres,absres,factor);
1148
1149 // safety net check: save the best-so-far solution
1150 if ( fasp_dvec_isnan(u) ) {
1151 // If the solution is NAN, restore the best solution
1152 absres = BIGREAL;
1153 goto RESTORE_BESTSOL;
1154 }
1155
1156 if ( absres < absres_best - maxdiff) {
1157 absres_best = absres;
1158 iter_best = iter;
1159 fasp_darray_cp(m,u->val,u_best);
1160 }
1161
1162 // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1163 normuinf = fasp_blas_darray_norminf(m, u->val);
1164 if (normuinf <= sol_inf_tol) {
1165 if ( PrtLvl > PRINT_MIN ) ITS_ZEROSOL;
1166 iter = ERROR_SOLVER_SOLSTAG;
1167 break;
1168 }
1169
1170 // Check II: if staggenated, try to restart
1171 normuu = fasp_blas_darray_norm2(m,p1);
1172 normuu = ABS(alpha)*(normuu/normu2);
1173
1174 if ( normuu < maxdiff ) {
1175
1176 if ( stag < MaxStag ) {
1177 if ( PrtLvl >= PRINT_MORE ) {
1178 ITS_DIFFRES(normuu,relres);
1179 ITS_RESTART;
1180 }
1181 }
1182
1183 fasp_darray_cp(m,b->val,r);
1184 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1185
1186 // compute residuals
1187 switch (StopType) {
1188 case STOP_REL_RES:
1189 temp2 = fasp_blas_darray_dotprod(m,r,r);
1190 absres = sqrt(temp2);
1191 relres = absres/normr0;
1192 break;
1193 case STOP_REL_PRECRES:
1194 if (pc == NULL)
1195 fasp_darray_cp(m,r,t);
1196 else
1197 pc->fct(r,t,pc->data);
1198 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
1199 absres = sqrt(temp2);
1200 relres = absres/normr0;
1201 break;
1202 case STOP_MOD_REL_RES:
1203 temp2 = fasp_blas_darray_dotprod(m,r,r);
1204 absres = sqrt(temp2);
1205 relres = absres/normu2;
1206 break;
1207 }
1208
1209 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1210
1211 if ( relres < tol )
1212 break;
1213 else {
1214 if ( stag >= MaxStag ) {
1215 if ( PrtLvl > PRINT_MIN ) ITS_STAGGED;
1216 iter = ERROR_SOLVER_STAG;
1217 break;
1218 }
1219 fasp_darray_set(m,p0,0.0);
1220 ++stag;
1221 ++restart_step;
1222
1223 // p1 = B(r)
1224 if ( pc != NULL )
1225 pc->fct(r,p1,pc->data); /* Apply preconditioner */
1226 else
1227 fasp_darray_cp(m,r,p1); /* No preconditioner */
1228
1229 // tp = A*p1
1230 fasp_blas_dstr_mxv(A,p1,tp);
1231
1232 // tz = B(tp)
1233 if ( pc != NULL )
1234 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1235 else
1236 fasp_darray_cp(m,tp,tz); /* No preconditioner */
1237
1238 // p1 = p1/normp
1239 normp = fasp_blas_darray_dotprod(m,tz,tp);
1240 normp = sqrt(normp);
1241 fasp_darray_cp(m,p1,t);
1242
1243 // t0 = A*p0=0
1244 fasp_darray_set(m,t0,0.0);
1245 fasp_darray_cp(m,t0,z0);
1246 fasp_darray_cp(m,t0,t1);
1247 fasp_darray_cp(m,t0,z1);
1248 fasp_darray_cp(m,t0,p1);
1249
1250 fasp_blas_darray_axpy(m,1/normp,t,p1);
1251
1252 // t1 = tp/normp, z1 = tz/normp
1253 fasp_blas_darray_axpy(m,1/normp,tp,t1);
1254 fasp_blas_darray_axpy(m,1/normp,tz,z1);
1255 }
1256 }
1257
1258 // Check III: prevent false convergence
1259 if ( relres < tol ) {
1260
1261 if ( PrtLvl >= PRINT_MORE ) ITS_COMPRES(relres);
1262
1263 // compute residual r = b - Ax again
1264 fasp_darray_cp(m,b->val,r);
1265 fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1266
1267 // compute residuals
1268 switch (StopType) {
1269 case STOP_REL_RES:
1270 temp2 = fasp_blas_darray_dotprod(m,r,r);
1271 absres = sqrt(temp2);
1272 relres = absres/normr0;
1273 break;
1274 case STOP_REL_PRECRES:
1275 if (pc == NULL)
1276 fasp_darray_cp(m,r,t);
1277 else
1278 pc->fct(r,t,pc->data);
1279 temp2 = ABS(fasp_blas_darray_dotprod(m,r,t));
1280 absres = sqrt(temp2);
1281 relres = absres/normr0;
1282 break;
1283 case STOP_MOD_REL_RES:
1284 temp2 = fasp_blas_darray_dotprod(m,r,r);
1285 absres = sqrt(temp2);
1286 relres = absres/normu2;
1287 break;
1288 }
1289
1290 if ( PrtLvl >= PRINT_MORE ) ITS_REALRES(relres);
1291
1292 // check convergence
1293 if ( relres < tol ) break;
1294
1295 if ( more_step >= MaxRestartStep ) {
1296 if ( PrtLvl > PRINT_MIN ) ITS_ZEROTOL;
1297 iter = ERROR_SOLVER_TOLSMALL;
1298 break;
1299 }
1300
1301 // prepare for restarting method
1302 fasp_darray_set(m,p0,0.0);
1303 ++more_step;
1304 ++restart_step;
1305
1306 // p1 = B(r)
1307 if ( pc != NULL )
1308 pc->fct(r,p1,pc->data); /* Apply preconditioner */
1309 else
1310 fasp_darray_cp(m,r,p1); /* No preconditioner */
1311
1312 // tp = A*p1
1313 fasp_blas_dstr_mxv(A,p1,tp);
1314
1315 // tz = B(tp)
1316 if ( pc != NULL )
1317 pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1318 else
1319 fasp_darray_cp(m,tp,tz); /* No preconditioner */
1320
1321 // p1 = p1/normp
1322 normp = fasp_blas_darray_dotprod(m,tz,tp);
1323 normp = sqrt(normp);
1324 fasp_darray_cp(m,p1,t);
1325
1326 // t0 = A*p0 = 0
1327 fasp_darray_set(m,t0,0.0);
1328 fasp_darray_cp(m,t0,z0);
1329 fasp_darray_cp(m,t0,t1);
1330 fasp_darray_cp(m,t0,z1);
1331 fasp_darray_cp(m,t0,p1);
1332
1333 fasp_blas_darray_axpy(m,1/normp,t,p1);
1334
1335 // t1=tp/normp,z1=tz/normp
1336 fasp_blas_darray_axpy(m,1/normp,tp,t1);
1337 fasp_blas_darray_axpy(m,1/normp,tz,z1);
1338
1339 } // end of convergence check
1340
1341 // update relative residual here
1342 absres0 = absres;
1343
1344 } // end of the main loop
1345
1346RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1347 if ( iter != iter_best ) {
1348
1349 // compute best residual
1350 fasp_darray_cp(m,b->val,r);
1351 fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
1352
1353 switch ( StopType ) {
1354 case STOP_REL_RES:
1355 absres_best = fasp_blas_darray_norm2(m,r);
1356 break;
1357 case STOP_REL_PRECRES:
1358 if ( pc != NULL )
1359 pc->fct(r,t,pc->data); /* Apply preconditioner */
1360 else
1361 fasp_darray_cp(m,r,t); /* No preconditioner */
1362 absres_best = sqrt(ABS(fasp_blas_darray_dotprod(m,t,r)));
1363 break;
1364 case STOP_MOD_REL_RES:
1365 absres_best = fasp_blas_darray_norm2(m,r);
1366 break;
1367 }
1368
1369 if ( absres > absres_best + maxdiff || isnan(absres) ) {
1370 if ( PrtLvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1371 fasp_darray_cp(m,u_best,u->val);
1372 relres = absres_best / normr0;
1373 }
1374 }
1375
1376FINISHED: // finish iterative method
1377 if ( PrtLvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1378
1379 // clean up temp memory
1380 fasp_mem_free(work); work = NULL;
1381
1382#if DEBUG_MODE > 0
1383 printf("### DEBUG: [--End--] %s ...\n", __FUNCTION__);
1384#endif
1385
1386 if ( iter > MaxIt )
1387 return ERROR_SOLVER_MAXIT;
1388 else
1389 return iter;
1390}
1391
1392/*---------------------------------*/
1393/*-- End of File --*/
1394/*---------------------------------*/
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_dvec_isnan
SHORT fasp_dvec_isnan(const dvector *u)
Check a dvector whether there is NAN.
Definition: AuxVector.c:39
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_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_spminres
INT fasp_solver_dcsr_spminres(const dCSRmat *A, const 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 with safety net.
Definition: KrySPminres.c:60
fasp_solver_dstr_spminres
INT fasp_solver_dstr_spminres(const dSTRmat *A, const 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 with safety net.
Definition: KrySPminres.c:962
fasp_solver_dblc_spminres
INT fasp_solver_dblc_spminres(const dBLCmat *A, const 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 with safety net.
Definition: KrySPminres.c:511
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
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