Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
spminres.c
Go to the documentation of this file.
1 
64 #include <math.h>
65 
66 #include "fasp.h"
67 #include "fasp_functs.h"
68 #include "itsolver_util.inl"
69 
70 /*---------------------------------*/
71 /*-- Public Functions --*/
72 /*---------------------------------*/
73 
96  dvector *b,
97  dvector *u,
98  precond *pc,
99  const REAL tol,
100  const INT MaxIt,
101  const SHORT stop_type,
102  const SHORT prtlvl)
103 {
104  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
105  const INT m = b->row;
106  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
107  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
108 
109  // local variables
110  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
111  REAL absres0 = BIGREAL, absres = BIGREAL;
112  REAL normr0 = BIGREAL, relres = BIGREAL;
113  REAL normu2, normuu, normp, infnormu, factor;
114  REAL alpha, alpha0, alpha1, temp2;
115 
116  INT iter_best = 0; // initial best known iteration
117  REAL absres_best = BIGREAL; // initial best known residual
118 
119  // allocate temp memory (need 12*m REAL)
120  REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
121  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
122  REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
123 
124 #if DEBUG_MODE > 0
125  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
126  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
127 #endif
128 
129  // p0 = 0
130  fasp_array_set(m,p0,0.0);
131 
132  // r = b-A*u
133  fasp_array_cp(m,b->val,r);
134  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
135 
136  // p1 = B(r)
137  if ( pc != NULL )
138  pc->fct(r,p1,pc->data); /* Apply preconditioner */
139  else
140  fasp_array_cp(m,r,p1); /* No preconditioner */
141 
142  // compute initial residuals
143  switch ( stop_type ) {
144  case STOP_REL_RES:
145  absres0 = fasp_blas_array_norm2(m,r);
146  normr0 = MAX(SMALLREAL,absres0);
147  relres = absres0/normr0;
148  break;
149  case STOP_REL_PRECRES:
150  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
151  normr0 = MAX(SMALLREAL,absres0);
152  relres = absres0/normr0;
153  break;
154  case STOP_MOD_REL_RES:
155  absres0 = fasp_blas_array_norm2(m,r);
156  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
157  relres = absres0/normu2;
158  break;
159  default:
160  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
161  goto FINISHED;
162  }
163 
164  // if initial residual is small, no need to iterate!
165  if ( relres < tol ) goto FINISHED;
166 
167  // output iteration information if needed
168  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
169 
170  // tp = A*p1
171  fasp_blas_dcsr_mxv(A,p1,tp);
172 
173  // tz = B(tp)
174  if ( pc != NULL )
175  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
176  else
177  fasp_array_cp(m,tp,tz); /* No preconditioner */
178 
179  // p1 = p1/normp
180  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
181  normp = sqrt(normp);
182  fasp_array_cp(m,p1,t);
183  fasp_array_set(m,p1,0.0);
184  fasp_blas_array_axpy(m,1/normp,t,p1);
185 
186  // t0 = A*p0 = 0
187  fasp_array_set(m,t0,0.0);
188  fasp_array_cp(m,t0,z0);
189  fasp_array_cp(m,t0,t1);
190  fasp_array_cp(m,t0,z1);
191 
192  // t1 = tp/normp, z1 = tz/normp
193  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
194  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
195 
196  // main MinRes loop
197  while ( iter++ < MaxIt ) {
198 
199  // alpha = <r,z1>
200  alpha=fasp_blas_array_dotprod(m,r,z1);
201 
202  // u = u+alpha*p1
203  fasp_blas_array_axpy(m,alpha,p1,u->val);
204 
205  // r = r-alpha*Ap1
206  fasp_blas_array_axpy(m,-alpha,t1,r);
207 
208  // compute t = A*z1 alpha1 = <z1,t>
209  fasp_blas_dcsr_mxv(A,z1,t);
210  alpha1=fasp_blas_array_dotprod(m,z1,t);
211 
212  // compute t = A*z0 alpha0 = <z1,t>
213  fasp_blas_dcsr_mxv(A,z0,t);
214  alpha0=fasp_blas_array_dotprod(m,z1,t);
215 
216  // p2 = z1-alpha1*p1-alpha0*p0
217  fasp_array_cp(m,z1,p2);
218  fasp_blas_array_axpy(m,-alpha1,p1,p2);
219  fasp_blas_array_axpy(m,-alpha0,p0,p2);
220 
221  // tp = A*p2
222  fasp_blas_dcsr_mxv(A,p2,tp);
223 
224  // tz = B(tp)
225  if ( pc != NULL )
226  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
227  else
228  fasp_array_cp(m,tp,tz); /* No preconditioner */
229 
230  // p2 = p2/normp
231  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
232  normp = sqrt(normp);
233  fasp_array_cp(m,p2,t);
234  fasp_array_set(m,p2,0.0);
235  fasp_blas_array_axpy(m,1/normp,t,p2);
236 
237  // prepare for the next iteration
238  fasp_array_cp(m,p1,p0);
239  fasp_array_cp(m,p2,p1);
240  fasp_array_cp(m,t1,t0);
241  fasp_array_cp(m,z1,z0);
242 
243  // t1=tp/normp,z1=tz/normp
244  fasp_array_set(m,t1,0.0);
245  fasp_array_cp(m,t1,z1);
246  fasp_blas_array_axpy(m,1/normp,tp,t1);
247  fasp_blas_array_axpy(m,1/normp,tz,z1);
248 
249  normu2 = fasp_blas_array_norm2(m,u->val);
250 
251  // compute residuals
252  switch ( stop_type ) {
253  case STOP_REL_RES:
254  temp2 = fasp_blas_array_dotprod(m,r,r);
255  absres = sqrt(temp2);
256  relres = absres/normr0;
257  break;
258  case STOP_REL_PRECRES:
259  if (pc == NULL)
260  fasp_array_cp(m,r,t);
261  else
262  pc->fct(r,t,pc->data);
263  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
264  absres = sqrt(temp2);
265  relres = absres/normr0;
266  break;
267  case STOP_MOD_REL_RES:
268  temp2 = fasp_blas_array_dotprod(m,r,r);
269  absres = sqrt(temp2);
270  relres = absres/normu2;
271  break;
272  }
273 
274  // compute reducation factor of residual ||r||
275  factor = absres/absres0;
276 
277  // output iteration information if needed
278  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
279 
280  // safety net check: save the best-so-far solution
281  if ( fasp_dvec_isnan(u) ) {
282  // If the solution is NAN, restrore the best solution
283  absres = BIGREAL;
284  goto RESTORE_BESTSOL;
285  }
286 
287  if ( absres < absres_best - maxdiff) {
288  absres_best = absres;
289  iter_best = iter;
290  fasp_array_cp(m,u->val,u_best);
291  }
292 
293  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
294  infnormu = fasp_blas_array_norminf(m, u->val);
295  if (infnormu <= sol_inf_tol) {
296  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
297  iter = ERROR_SOLVER_SOLSTAG;
298  break;
299  }
300 
301  // Check II: if staggenated, try to restart
302  normuu = fasp_blas_array_norm2(m,p1);
303  normuu = ABS(alpha)*(normuu/normu2);
304 
305  if ( normuu < maxdiff ) {
306 
307  if ( stag < MaxStag ) {
308  if ( prtlvl >= PRINT_MORE ) {
309  ITS_DIFFRES(normuu,relres);
310  ITS_RESTART;
311  }
312  }
313 
314  fasp_array_cp(m,b->val,r);
315  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
316 
317  // compute residuals
318  switch (stop_type) {
319  case STOP_REL_RES:
320  temp2 = fasp_blas_array_dotprod(m,r,r);
321  absres = sqrt(temp2);
322  relres = absres/normr0;
323  break;
324  case STOP_REL_PRECRES:
325  if (pc == NULL)
326  fasp_array_cp(m,r,t);
327  else
328  pc->fct(r,t,pc->data);
329  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
330  absres = sqrt(temp2);
331  relres = absres/normr0;
332  break;
333  case STOP_MOD_REL_RES:
334  temp2 = fasp_blas_array_dotprod(m,r,r);
335  absres = sqrt(temp2);
336  relres = absres/normu2;
337  break;
338  }
339 
340  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
341 
342  if ( relres < tol )
343  break;
344  else {
345  if ( stag >= MaxStag ) {
346  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
347  iter = ERROR_SOLVER_STAG;
348  break;
349  }
350  fasp_array_set(m,p0,0.0);
351  ++stag;
352  ++restart_step;
353 
354  // p1 = B(r)
355  if ( pc != NULL )
356  pc->fct(r,p1,pc->data); /* Apply preconditioner */
357  else
358  fasp_array_cp(m,r,p1); /* No preconditioner */
359 
360  // tp = A*p1
361  fasp_blas_dcsr_mxv(A,p1,tp);
362 
363  // tz = B(tp)
364  if ( pc != NULL )
365  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
366  else
367  fasp_array_cp(m,tp,tz); /* No preconditioner */
368 
369  // p1 = p1/normp
370  normp = fasp_blas_array_dotprod(m,tz,tp);
371  normp = sqrt(normp);
372  fasp_array_cp(m,p1,t);
373 
374  // t0 = A*p0=0
375  fasp_array_set(m,t0,0.0);
376  fasp_array_cp(m,t0,z0);
377  fasp_array_cp(m,t0,t1);
378  fasp_array_cp(m,t0,z1);
379  fasp_array_cp(m,t0,p1);
380 
381  fasp_blas_array_axpy(m,1/normp,t,p1);
382 
383  // t1 = tp/normp, z1 = tz/normp
384  fasp_blas_array_axpy(m,1/normp,tp,t1);
385  fasp_blas_array_axpy(m,1/normp,tz,z1);
386  }
387  }
388 
389  // Check III: prevent false convergence
390  if ( relres < tol ) {
391 
392  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
393 
394  // compute residual r = b - Ax again
395  fasp_array_cp(m,b->val,r);
396  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
397 
398  // compute residuals
399  switch (stop_type) {
400  case STOP_REL_RES:
401  temp2 = fasp_blas_array_dotprod(m,r,r);
402  absres = sqrt(temp2);
403  relres = absres/normr0;
404  break;
405  case STOP_REL_PRECRES:
406  if (pc == NULL)
407  fasp_array_cp(m,r,t);
408  else
409  pc->fct(r,t,pc->data);
410  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
411  absres = sqrt(temp2);
412  relres = absres/normr0;
413  break;
414  case STOP_MOD_REL_RES:
415  temp2 = fasp_blas_array_dotprod(m,r,r);
416  absres = sqrt(temp2);
417  relres = absres/normu2;
418  break;
419  }
420 
421  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
422 
423  // check convergence
424  if ( relres < tol ) break;
425 
426  if ( more_step >= MaxRestartStep ) {
427  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
428  iter = ERROR_SOLVER_TOLSMALL;
429  break;
430  }
431 
432  // prepare for restarting the method
433  fasp_array_set(m,p0,0.0);
434  ++more_step;
435  ++restart_step;
436 
437  // p1 = B(r)
438  if ( pc != NULL )
439  pc->fct(r,p1,pc->data); /* Apply preconditioner */
440  else
441  fasp_array_cp(m,r,p1); /* No preconditioner */
442 
443  // tp = A*p1
444  fasp_blas_dcsr_mxv(A,p1,tp);
445 
446  // tz = B(tp)
447  if ( pc != NULL )
448  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
449  else
450  fasp_array_cp(m,tp,tz); /* No preconditioner */
451 
452  // p1 = p1/normp
453  normp = fasp_blas_array_dotprod(m,tz,tp);
454  normp = sqrt(normp);
455  fasp_array_cp(m,p1,t);
456 
457  // t0 = A*p0 = 0
458  fasp_array_set(m,t0,0.0);
459  fasp_array_cp(m,t0,z0);
460  fasp_array_cp(m,t0,t1);
461  fasp_array_cp(m,t0,z1);
462  fasp_array_cp(m,t0,p1);
463 
464  fasp_blas_array_axpy(m,1/normp,t,p1);
465 
466  // t1=tp/normp,z1=tz/normp
467  fasp_blas_array_axpy(m,1/normp,tp,t1);
468  fasp_blas_array_axpy(m,1/normp,tz,z1);
469 
470  } // end of convergence check
471 
472  // update relative residual here
473  absres0 = absres;
474 
475  } // end of the main loop
476 
477 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
478  if ( iter != iter_best ) {
479 
480  // compute best residual
481  fasp_array_cp(m,b->val,r);
482  fasp_blas_dcsr_aAxpy(-1.0,A,u_best,r);
483 
484  switch ( stop_type ) {
485  case STOP_REL_RES:
486  absres_best = fasp_blas_array_norm2(m,r);
487  break;
488  case STOP_REL_PRECRES:
489  if ( pc != NULL )
490  pc->fct(r,t,pc->data); /* Apply preconditioner */
491  else
492  fasp_array_cp(m,r,t); /* No preconditioner */
493  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,t,r)));
494  break;
495  case STOP_MOD_REL_RES:
496  absres_best = fasp_blas_array_norm2(m,r);
497  break;
498  }
499 
500  if ( absres > absres_best + maxdiff ) {
501  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
502  fasp_array_cp(m,u_best,u->val);
503  relres = absres_best / normr0;
504  }
505  }
506 
507 FINISHED: // finish the iterative method
508  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
509 
510  // clean up temp memory
511  fasp_mem_free(work);
512 
513 #if DEBUG_MODE > 0
514  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
515 #endif
516 
517  if ( iter > MaxIt )
518  return ERROR_SOLVER_MAXIT;
519  else
520  return iter;
521 }
522 
545  dvector *b,
546  dvector *u,
547  precond *pc,
548  const REAL tol,
549  const INT MaxIt,
550  const SHORT stop_type,
551  const SHORT prtlvl)
552 {
553  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
554  const INT m = b->row;
555  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
556  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
557 
558  // local variables
559  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
560  REAL absres0 = BIGREAL, absres = BIGREAL;
561  REAL normr0 = BIGREAL, relres = BIGREAL;
562  REAL normu2, normuu, normp, infnormu, factor;
563  REAL alpha, alpha0, alpha1, temp2;
564 
565  INT iter_best = 0; // initial best known iteration
566  REAL absres_best = BIGREAL; // initial best known residual
567 
568  // allocate temp memory (need 12*m REAL)
569  REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
570  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
571  REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
572 
573 #if DEBUG_MODE > 0
574  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
575  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
576 #endif
577 
578  // p0 = 0
579  fasp_array_set(m,p0,0.0);
580 
581  // r = b-A*u
582  fasp_array_cp(m,b->val,r);
583  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
584 
585  // p1 = B(r)
586  if ( pc != NULL )
587  pc->fct(r,p1,pc->data); /* Apply preconditioner */
588  else
589  fasp_array_cp(m,r,p1); /* No preconditioner */
590 
591  // compute initial residuals
592  switch ( stop_type ) {
593  case STOP_REL_RES:
594  absres0 = fasp_blas_array_norm2(m,r);
595  normr0 = MAX(SMALLREAL,absres0);
596  relres = absres0/normr0;
597  break;
598  case STOP_REL_PRECRES:
599  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
600  normr0 = MAX(SMALLREAL,absres0);
601  relres = absres0/normr0;
602  break;
603  case STOP_MOD_REL_RES:
604  absres0 = fasp_blas_array_norm2(m,r);
605  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
606  relres = absres0/normu2;
607  break;
608  default:
609  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
610  goto FINISHED;
611  }
612 
613  // if initial residual is small, no need to iterate!
614  if ( relres < tol ) goto FINISHED;
615 
616  // output iteration information if needed
617  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
618 
619  // tp = A*p1
620  fasp_blas_bdcsr_mxv(A,p1,tp);
621 
622  // tz = B(tp)
623  if ( pc != NULL )
624  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
625  else
626  fasp_array_cp(m,tp,tz); /* No preconditioner */
627 
628  // p1 = p1/normp
629  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
630  normp = sqrt(normp);
631  fasp_array_cp(m,p1,t);
632  fasp_array_set(m,p1,0.0);
633  fasp_blas_array_axpy(m,1/normp,t,p1);
634 
635  // t0 = A*p0 = 0
636  fasp_array_set(m,t0,0.0);
637  fasp_array_cp(m,t0,z0);
638  fasp_array_cp(m,t0,t1);
639  fasp_array_cp(m,t0,z1);
640 
641  // t1 = tp/normp, z1 = tz/normp
642  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
643  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
644 
645  // main MinRes loop
646  while ( iter++ < MaxIt ) {
647 
648  // alpha = <r,z1>
649  alpha=fasp_blas_array_dotprod(m,r,z1);
650 
651  // u = u+alpha*p1
652  fasp_blas_array_axpy(m,alpha,p1,u->val);
653 
654  // r = r-alpha*Ap1
655  fasp_blas_array_axpy(m,-alpha,t1,r);
656 
657  // compute t = A*z1 alpha1 = <z1,t>
658  fasp_blas_bdcsr_mxv(A,z1,t);
659  alpha1=fasp_blas_array_dotprod(m,z1,t);
660 
661  // compute t = A*z0 alpha0 = <z1,t>
662  fasp_blas_bdcsr_mxv(A,z0,t);
663  alpha0=fasp_blas_array_dotprod(m,z1,t);
664 
665  // p2 = z1-alpha1*p1-alpha0*p0
666  fasp_array_cp(m,z1,p2);
667  fasp_blas_array_axpy(m,-alpha1,p1,p2);
668  fasp_blas_array_axpy(m,-alpha0,p0,p2);
669 
670  // tp = A*p2
671  fasp_blas_bdcsr_mxv(A,p2,tp);
672 
673  // tz = B(tp)
674  if ( pc != NULL )
675  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
676  else
677  fasp_array_cp(m,tp,tz); /* No preconditioner */
678 
679  // p2 = p2/normp
680  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
681  normp = sqrt(normp);
682  fasp_array_cp(m,p2,t);
683  fasp_array_set(m,p2,0.0);
684  fasp_blas_array_axpy(m,1/normp,t,p2);
685 
686  // prepare for the next iteration
687  fasp_array_cp(m,p1,p0);
688  fasp_array_cp(m,p2,p1);
689  fasp_array_cp(m,t1,t0);
690  fasp_array_cp(m,z1,z0);
691 
692  // t1=tp/normp,z1=tz/normp
693  fasp_array_set(m,t1,0.0);
694  fasp_array_cp(m,t1,z1);
695  fasp_blas_array_axpy(m,1/normp,tp,t1);
696  fasp_blas_array_axpy(m,1/normp,tz,z1);
697 
698  normu2 = fasp_blas_array_norm2(m,u->val);
699 
700  // compute residuals
701  switch ( stop_type ) {
702  case STOP_REL_RES:
703  temp2 = fasp_blas_array_dotprod(m,r,r);
704  absres = sqrt(temp2);
705  relres = absres/normr0;
706  break;
707  case STOP_REL_PRECRES:
708  if (pc == NULL)
709  fasp_array_cp(m,r,t);
710  else
711  pc->fct(r,t,pc->data);
712  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
713  absres = sqrt(temp2);
714  relres = absres/normr0;
715  break;
716  case STOP_MOD_REL_RES:
717  temp2 = fasp_blas_array_dotprod(m,r,r);
718  absres = sqrt(temp2);
719  relres = absres/normu2;
720  break;
721  }
722 
723  // compute reducation factor of residual ||r||
724  factor = absres/absres0;
725 
726  // output iteration information if needed
727  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
728 
729  // safety net check: save the best-so-far solution
730  if ( fasp_dvec_isnan(u) ) {
731  // If the solution is NAN, restrore the best solution
732  absres = BIGREAL;
733  goto RESTORE_BESTSOL;
734  }
735 
736  if ( absres < absres_best - maxdiff) {
737  absres_best = absres;
738  iter_best = iter;
739  fasp_array_cp(m,u->val,u_best);
740  }
741 
742  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
743  infnormu = fasp_blas_array_norminf(m, u->val);
744  if (infnormu <= sol_inf_tol) {
745  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
746  iter = ERROR_SOLVER_SOLSTAG;
747  break;
748  }
749 
750  // Check II: if staggenated, try to restart
751  normuu = fasp_blas_array_norm2(m,p1);
752  normuu = ABS(alpha)*(normuu/normu2);
753 
754  if ( normuu < maxdiff ) {
755 
756  if ( stag < MaxStag ) {
757  if ( prtlvl >= PRINT_MORE ) {
758  ITS_DIFFRES(normuu,relres);
759  ITS_RESTART;
760  }
761  }
762 
763  fasp_array_cp(m,b->val,r);
764  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
765 
766  // compute residuals
767  switch (stop_type) {
768  case STOP_REL_RES:
769  temp2 = fasp_blas_array_dotprod(m,r,r);
770  absres = sqrt(temp2);
771  relres = absres/normr0;
772  break;
773  case STOP_REL_PRECRES:
774  if (pc == NULL)
775  fasp_array_cp(m,r,t);
776  else
777  pc->fct(r,t,pc->data);
778  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
779  absres = sqrt(temp2);
780  relres = absres/normr0;
781  break;
782  case STOP_MOD_REL_RES:
783  temp2 = fasp_blas_array_dotprod(m,r,r);
784  absres = sqrt(temp2);
785  relres = absres/normu2;
786  break;
787  }
788 
789  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
790 
791  if ( relres < tol )
792  break;
793  else {
794  if ( stag >= MaxStag ) {
795  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
796  iter = ERROR_SOLVER_STAG;
797  break;
798  }
799  fasp_array_set(m,p0,0.0);
800  ++stag;
801  ++restart_step;
802 
803  // p1 = B(r)
804  if ( pc != NULL )
805  pc->fct(r,p1,pc->data); /* Apply preconditioner */
806  else
807  fasp_array_cp(m,r,p1); /* No preconditioner */
808 
809  // tp = A*p1
810  fasp_blas_bdcsr_mxv(A,p1,tp);
811 
812  // tz = B(tp)
813  if ( pc != NULL )
814  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
815  else
816  fasp_array_cp(m,tp,tz); /* No preconditioner */
817 
818  // p1 = p1/normp
819  normp = fasp_blas_array_dotprod(m,tz,tp);
820  normp = sqrt(normp);
821  fasp_array_cp(m,p1,t);
822 
823  // t0 = A*p0=0
824  fasp_array_set(m,t0,0.0);
825  fasp_array_cp(m,t0,z0);
826  fasp_array_cp(m,t0,t1);
827  fasp_array_cp(m,t0,z1);
828  fasp_array_cp(m,t0,p1);
829 
830  fasp_blas_array_axpy(m,1/normp,t,p1);
831 
832  // t1 = tp/normp, z1 = tz/normp
833  fasp_blas_array_axpy(m,1/normp,tp,t1);
834  fasp_blas_array_axpy(m,1/normp,tz,z1);
835  }
836  }
837 
838  // Check III: prevent false convergence
839  if ( relres < tol ) {
840 
841  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
842 
843  // compute residual r = b - Ax again
844  fasp_array_cp(m,b->val,r);
845  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
846 
847  // compute residuals
848  switch (stop_type) {
849  case STOP_REL_RES:
850  temp2 = fasp_blas_array_dotprod(m,r,r);
851  absres = sqrt(temp2);
852  relres = absres/normr0;
853  break;
854  case STOP_REL_PRECRES:
855  if (pc == NULL)
856  fasp_array_cp(m,r,t);
857  else
858  pc->fct(r,t,pc->data);
859  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
860  absres = sqrt(temp2);
861  relres = absres/normr0;
862  break;
863  case STOP_MOD_REL_RES:
864  temp2 = fasp_blas_array_dotprod(m,r,r);
865  absres = sqrt(temp2);
866  relres = absres/normu2;
867  break;
868  }
869 
870  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
871 
872  // check convergence
873  if ( relres < tol ) break;
874 
875  if ( more_step >= MaxRestartStep ) {
876  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
877  iter = ERROR_SOLVER_TOLSMALL;
878  break;
879  }
880 
881  // prepare for restarting the method
882  fasp_array_set(m,p0,0.0);
883  ++more_step;
884  ++restart_step;
885 
886  // p1 = B(r)
887  if ( pc != NULL )
888  pc->fct(r,p1,pc->data); /* Apply preconditioner */
889  else
890  fasp_array_cp(m,r,p1); /* No preconditioner */
891 
892  // tp = A*p1
893  fasp_blas_bdcsr_mxv(A,p1,tp);
894 
895  // tz = B(tp)
896  if ( pc != NULL )
897  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
898  else
899  fasp_array_cp(m,tp,tz); /* No preconditioner */
900 
901  // p1 = p1/normp
902  normp = fasp_blas_array_dotprod(m,tz,tp);
903  normp = sqrt(normp);
904  fasp_array_cp(m,p1,t);
905 
906  // t0 = A*p0 = 0
907  fasp_array_set(m,t0,0.0);
908  fasp_array_cp(m,t0,z0);
909  fasp_array_cp(m,t0,t1);
910  fasp_array_cp(m,t0,z1);
911  fasp_array_cp(m,t0,p1);
912 
913  fasp_blas_array_axpy(m,1/normp,t,p1);
914 
915  // t1=tp/normp,z1=tz/normp
916  fasp_blas_array_axpy(m,1/normp,tp,t1);
917  fasp_blas_array_axpy(m,1/normp,tz,z1);
918 
919  } // end of convergence check
920 
921  // update relative residual here
922  absres0 = absres;
923 
924  } // end of the main loop
925 
926 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
927  if ( iter != iter_best ) {
928 
929  // compute best residual
930  fasp_array_cp(m,b->val,r);
931  fasp_blas_bdcsr_aAxpy(-1.0,A,u_best,r);
932 
933  switch ( stop_type ) {
934  case STOP_REL_RES:
935  absres_best = fasp_blas_array_norm2(m,r);
936  break;
937  case STOP_REL_PRECRES:
938  if ( pc != NULL )
939  pc->fct(r,t,pc->data); /* Apply preconditioner */
940  else
941  fasp_array_cp(m,r,t); /* No preconditioner */
942  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,t,r)));
943  break;
944  case STOP_MOD_REL_RES:
945  absres_best = fasp_blas_array_norm2(m,r);
946  break;
947  }
948 
949  if ( absres > absres_best + maxdiff ) {
950  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
951  fasp_array_cp(m,u_best,u->val);
952  relres = absres_best / normr0;
953  }
954  }
955 
956 FINISHED: // finish the iterative method
957  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
958 
959  // clean up temp memory
960  fasp_mem_free(work);
961 
962 #if DEBUG_MODE > 0
963  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
964 #endif
965 
966  if ( iter > MaxIt )
967  return ERROR_SOLVER_MAXIT;
968  else
969  return iter;
970 }
971 
994  dvector *b,
995  dvector *u,
996  precond *pc,
997  const REAL tol,
998  const INT MaxIt,
999  const SHORT stop_type,
1000  const SHORT prtlvl)
1001 {
1002  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
1003  const INT m = b->row;
1004  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
1005  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
1006 
1007  // local variables
1008  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
1009  REAL absres0 = BIGREAL, absres = BIGREAL;
1010  REAL normr0 = BIGREAL, relres = BIGREAL;
1011  REAL normu2, normuu, normp, infnormu, factor;
1012  REAL alpha, alpha0, alpha1, temp2;
1013 
1014  INT iter_best = 0; // initial best known iteration
1015  REAL absres_best = BIGREAL; // initial best known residual
1016 
1017  // allocate temp memory (need 12*m REAL)
1018  REAL *work=(REAL *)fasp_mem_calloc(12*m,sizeof(REAL));
1019  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m, *t0=z1+m;
1020  REAL *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m, *u_best = r+m;
1021 
1022 #if DEBUG_MODE > 0
1023  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
1024  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
1025 #endif
1026 
1027  // p0 = 0
1028  fasp_array_set(m,p0,0.0);
1029 
1030  // r = b-A*u
1031  fasp_array_cp(m,b->val,r);
1032  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1033 
1034  // p1 = B(r)
1035  if ( pc != NULL )
1036  pc->fct(r,p1,pc->data); /* Apply preconditioner */
1037  else
1038  fasp_array_cp(m,r,p1); /* No preconditioner */
1039 
1040  // compute initial residuals
1041  switch ( stop_type ) {
1042  case STOP_REL_RES:
1043  absres0 = fasp_blas_array_norm2(m,r);
1044  normr0 = MAX(SMALLREAL,absres0);
1045  relres = absres0/normr0;
1046  break;
1047  case STOP_REL_PRECRES:
1048  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
1049  normr0 = MAX(SMALLREAL,absres0);
1050  relres = absres0/normr0;
1051  break;
1052  case STOP_MOD_REL_RES:
1053  absres0 = fasp_blas_array_norm2(m,r);
1054  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
1055  relres = absres0/normu2;
1056  break;
1057  default:
1058  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
1059  goto FINISHED;
1060  }
1061 
1062  // if initial residual is small, no need to iterate!
1063  if ( relres < tol ) goto FINISHED;
1064 
1065  // output iteration information if needed
1066  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
1067 
1068  // tp = A*p1
1069  fasp_blas_dstr_mxv(A,p1,tp);
1070 
1071  // tz = B(tp)
1072  if ( pc != NULL )
1073  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
1074  else
1075  fasp_array_cp(m,tp,tz); /* No preconditioner */
1076 
1077  // p1 = p1/normp
1078  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
1079  normp = sqrt(normp);
1080  fasp_array_cp(m,p1,t);
1081  fasp_array_set(m,p1,0.0);
1082  fasp_blas_array_axpy(m,1/normp,t,p1);
1083 
1084  // t0 = A*p0 = 0
1085  fasp_array_set(m,t0,0.0);
1086  fasp_array_cp(m,t0,z0);
1087  fasp_array_cp(m,t0,t1);
1088  fasp_array_cp(m,t0,z1);
1089 
1090  // t1 = tp/normp, z1 = tz/normp
1091  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
1092  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
1093 
1094  // main MinRes loop
1095  while ( iter++ < MaxIt ) {
1096 
1097  // alpha = <r,z1>
1098  alpha=fasp_blas_array_dotprod(m,r,z1);
1099 
1100  // u = u+alpha*p1
1101  fasp_blas_array_axpy(m,alpha,p1,u->val);
1102 
1103  // r = r-alpha*Ap1
1104  fasp_blas_array_axpy(m,-alpha,t1,r);
1105 
1106  // compute t = A*z1 alpha1 = <z1,t>
1107  fasp_blas_dstr_mxv(A,z1,t);
1108  alpha1=fasp_blas_array_dotprod(m,z1,t);
1109 
1110  // compute t = A*z0 alpha0 = <z1,t>
1111  fasp_blas_dstr_mxv(A,z0,t);
1112  alpha0=fasp_blas_array_dotprod(m,z1,t);
1113 
1114  // p2 = z1-alpha1*p1-alpha0*p0
1115  fasp_array_cp(m,z1,p2);
1116  fasp_blas_array_axpy(m,-alpha1,p1,p2);
1117  fasp_blas_array_axpy(m,-alpha0,p0,p2);
1118 
1119  // tp = A*p2
1120  fasp_blas_dstr_mxv(A,p2,tp);
1121 
1122  // tz = B(tp)
1123  if ( pc != NULL )
1124  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
1125  else
1126  fasp_array_cp(m,tp,tz); /* No preconditioner */
1127 
1128  // p2 = p2/normp
1129  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
1130  normp = sqrt(normp);
1131  fasp_array_cp(m,p2,t);
1132  fasp_array_set(m,p2,0.0);
1133  fasp_blas_array_axpy(m,1/normp,t,p2);
1134 
1135  // prepare for the next iteration
1136  fasp_array_cp(m,p1,p0);
1137  fasp_array_cp(m,p2,p1);
1138  fasp_array_cp(m,t1,t0);
1139  fasp_array_cp(m,z1,z0);
1140 
1141  // t1=tp/normp,z1=tz/normp
1142  fasp_array_set(m,t1,0.0);
1143  fasp_array_cp(m,t1,z1);
1144  fasp_blas_array_axpy(m,1/normp,tp,t1);
1145  fasp_blas_array_axpy(m,1/normp,tz,z1);
1146 
1147  normu2 = fasp_blas_array_norm2(m,u->val);
1148 
1149  // compute residuals
1150  switch ( stop_type ) {
1151  case STOP_REL_RES:
1152  temp2 = fasp_blas_array_dotprod(m,r,r);
1153  absres = sqrt(temp2);
1154  relres = absres/normr0;
1155  break;
1156  case STOP_REL_PRECRES:
1157  if (pc == NULL)
1158  fasp_array_cp(m,r,t);
1159  else
1160  pc->fct(r,t,pc->data);
1161  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1162  absres = sqrt(temp2);
1163  relres = absres/normr0;
1164  break;
1165  case STOP_MOD_REL_RES:
1166  temp2 = fasp_blas_array_dotprod(m,r,r);
1167  absres = sqrt(temp2);
1168  relres = absres/normu2;
1169  break;
1170  }
1171 
1172  // compute reducation factor of residual ||r||
1173  factor = absres/absres0;
1174 
1175  // output iteration information if needed
1176  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1177 
1178  // safety net check: save the best-so-far solution
1179  if ( fasp_dvec_isnan(u) ) {
1180  // If the solution is NAN, restrore the best solution
1181  absres = BIGREAL;
1182  goto RESTORE_BESTSOL;
1183  }
1184 
1185  if ( absres < absres_best - maxdiff) {
1186  absres_best = absres;
1187  iter_best = iter;
1188  fasp_array_cp(m,u->val,u_best);
1189  }
1190 
1191  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1192  infnormu = fasp_blas_array_norminf(m, u->val);
1193  if (infnormu <= sol_inf_tol) {
1194  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1195  iter = ERROR_SOLVER_SOLSTAG;
1196  break;
1197  }
1198 
1199  // Check II: if staggenated, try to restart
1200  normuu = fasp_blas_array_norm2(m,p1);
1201  normuu = ABS(alpha)*(normuu/normu2);
1202 
1203  if ( normuu < maxdiff ) {
1204 
1205  if ( stag < MaxStag ) {
1206  if ( prtlvl >= PRINT_MORE ) {
1207  ITS_DIFFRES(normuu,relres);
1208  ITS_RESTART;
1209  }
1210  }
1211 
1212  fasp_array_cp(m,b->val,r);
1213  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1214 
1215  // compute residuals
1216  switch (stop_type) {
1217  case STOP_REL_RES:
1218  temp2 = fasp_blas_array_dotprod(m,r,r);
1219  absres = sqrt(temp2);
1220  relres = absres/normr0;
1221  break;
1222  case STOP_REL_PRECRES:
1223  if (pc == NULL)
1224  fasp_array_cp(m,r,t);
1225  else
1226  pc->fct(r,t,pc->data);
1227  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1228  absres = sqrt(temp2);
1229  relres = absres/normr0;
1230  break;
1231  case STOP_MOD_REL_RES:
1232  temp2 = fasp_blas_array_dotprod(m,r,r);
1233  absres = sqrt(temp2);
1234  relres = absres/normu2;
1235  break;
1236  }
1237 
1238  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1239 
1240  if ( relres < tol )
1241  break;
1242  else {
1243  if ( stag >= MaxStag ) {
1244  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1245  iter = ERROR_SOLVER_STAG;
1246  break;
1247  }
1248  fasp_array_set(m,p0,0.0);
1249  ++stag;
1250  ++restart_step;
1251 
1252  // p1 = B(r)
1253  if ( pc != NULL )
1254  pc->fct(r,p1,pc->data); /* Apply preconditioner */
1255  else
1256  fasp_array_cp(m,r,p1); /* No preconditioner */
1257 
1258  // tp = A*p1
1259  fasp_blas_dstr_mxv(A,p1,tp);
1260 
1261  // tz = B(tp)
1262  if ( pc != NULL )
1263  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1264  else
1265  fasp_array_cp(m,tp,tz); /* No preconditioner */
1266 
1267  // p1 = p1/normp
1268  normp = fasp_blas_array_dotprod(m,tz,tp);
1269  normp = sqrt(normp);
1270  fasp_array_cp(m,p1,t);
1271 
1272  // t0 = A*p0=0
1273  fasp_array_set(m,t0,0.0);
1274  fasp_array_cp(m,t0,z0);
1275  fasp_array_cp(m,t0,t1);
1276  fasp_array_cp(m,t0,z1);
1277  fasp_array_cp(m,t0,p1);
1278 
1279  fasp_blas_array_axpy(m,1/normp,t,p1);
1280 
1281  // t1 = tp/normp, z1 = tz/normp
1282  fasp_blas_array_axpy(m,1/normp,tp,t1);
1283  fasp_blas_array_axpy(m,1/normp,tz,z1);
1284  }
1285  }
1286 
1287  // Check III: prevent false convergence
1288  if ( relres < tol ) {
1289 
1290  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
1291 
1292  // compute residual r = b - Ax again
1293  fasp_array_cp(m,b->val,r);
1294  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1295 
1296  // compute residuals
1297  switch (stop_type) {
1298  case STOP_REL_RES:
1299  temp2 = fasp_blas_array_dotprod(m,r,r);
1300  absres = sqrt(temp2);
1301  relres = absres/normr0;
1302  break;
1303  case STOP_REL_PRECRES:
1304  if (pc == NULL)
1305  fasp_array_cp(m,r,t);
1306  else
1307  pc->fct(r,t,pc->data);
1308  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1309  absres = sqrt(temp2);
1310  relres = absres/normr0;
1311  break;
1312  case STOP_MOD_REL_RES:
1313  temp2 = fasp_blas_array_dotprod(m,r,r);
1314  absres = sqrt(temp2);
1315  relres = absres/normu2;
1316  break;
1317  }
1318 
1319  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1320 
1321  // check convergence
1322  if ( relres < tol ) break;
1323 
1324  if ( more_step >= MaxRestartStep ) {
1325  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1326  iter = ERROR_SOLVER_TOLSMALL;
1327  break;
1328  }
1329 
1330  // prepare for restarting the method
1331  fasp_array_set(m,p0,0.0);
1332  ++more_step;
1333  ++restart_step;
1334 
1335  // p1 = B(r)
1336  if ( pc != NULL )
1337  pc->fct(r,p1,pc->data); /* Apply preconditioner */
1338  else
1339  fasp_array_cp(m,r,p1); /* No preconditioner */
1340 
1341  // tp = A*p1
1342  fasp_blas_dstr_mxv(A,p1,tp);
1343 
1344  // tz = B(tp)
1345  if ( pc != NULL )
1346  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1347  else
1348  fasp_array_cp(m,tp,tz); /* No preconditioner */
1349 
1350  // p1 = p1/normp
1351  normp = fasp_blas_array_dotprod(m,tz,tp);
1352  normp = sqrt(normp);
1353  fasp_array_cp(m,p1,t);
1354 
1355  // t0 = A*p0 = 0
1356  fasp_array_set(m,t0,0.0);
1357  fasp_array_cp(m,t0,z0);
1358  fasp_array_cp(m,t0,t1);
1359  fasp_array_cp(m,t0,z1);
1360  fasp_array_cp(m,t0,p1);
1361 
1362  fasp_blas_array_axpy(m,1/normp,t,p1);
1363 
1364  // t1=tp/normp,z1=tz/normp
1365  fasp_blas_array_axpy(m,1/normp,tp,t1);
1366  fasp_blas_array_axpy(m,1/normp,tz,z1);
1367 
1368  } // end of convergence check
1369 
1370  // update relative residual here
1371  absres0 = absres;
1372 
1373  } // end of the main loop
1374 
1375 RESTORE_BESTSOL: // restore the best-so-far solution if necessary
1376  if ( iter != iter_best ) {
1377 
1378  // compute best residual
1379  fasp_array_cp(m,b->val,r);
1380  fasp_blas_dstr_aAxpy(-1.0,A,u_best,r);
1381 
1382  switch ( stop_type ) {
1383  case STOP_REL_RES:
1384  absres_best = fasp_blas_array_norm2(m,r);
1385  break;
1386  case STOP_REL_PRECRES:
1387  if ( pc != NULL )
1388  pc->fct(r,t,pc->data); /* Apply preconditioner */
1389  else
1390  fasp_array_cp(m,r,t); /* No preconditioner */
1391  absres_best = sqrt(ABS(fasp_blas_array_dotprod(m,t,r)));
1392  break;
1393  case STOP_MOD_REL_RES:
1394  absres_best = fasp_blas_array_norm2(m,r);
1395  break;
1396  }
1397 
1398  if ( absres > absres_best + maxdiff ) {
1399  if ( prtlvl > PRINT_NONE ) ITS_RESTORE(iter_best);
1400  fasp_array_cp(m,u_best,u->val);
1401  relres = absres_best / normr0;
1402  }
1403  }
1404 
1405 FINISHED: // finish the iterative method
1406  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1407 
1408  // clean up temp memory
1409  fasp_mem_free(work);
1410 
1411 #if DEBUG_MODE > 0
1412  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1413 #endif
1414 
1415  if ( iter > MaxIt )
1416  return ERROR_SOLVER_MAXIT;
1417  else
1418  return iter;
1419 }
1420 
1421 /*---------------------------------*/
1422 /*-- End of File --*/
1423 /*---------------------------------*/
#define REAL
Definition: fasp.h:67
void fasp_blas_dstr_aAxpy(const REAL alpha, dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_str.c:47
void fasp_blas_bdcsr_mxv(block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_bcsr.c:155
void fasp_blas_bdcsr_aAxpy(const REAL alpha, block_dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_bcsr.c:30
#define STOP_MOD_REL_RES
Definition: fasp_const.h:133
Structure matrix of REAL type.
Definition: fasp.h:304
REAL * val
actual vector entries
Definition: fasp.h:348
void * fasp_mem_calloc(LONGLONG size, INT type)
1M = 1024*1024
Definition: memory.c:60
void(* fct)(REAL *, REAL *, void *)
action for preconditioner, void function pointer
Definition: fasp.h:1005
Vector with n entries of REAL type.
Definition: fasp.h:342
INT fasp_dvec_isnan(dvector *u)
Check a dvector whether there is NAN.
Definition: vec.c:33
#define ERROR_SOLVER_MAXIT
Definition: fasp_const.h:54
#define INT
Definition: fasp.h:64
void print_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: message.c:36
Preconditioner data and action.
Definition: fasp.h:999
void fasp_mem_free(void *mem)
Free up previous allocated memory body.
Definition: memory.c:150
void fasp_array_set(const INT n, REAL *x, const REAL val)
Set initial value for an array to be x=val.
Definition: array.c:48
void * data
data for preconditioner, void pointer
Definition: fasp.h:1002
#define STOP_REL_PRECRES
Definition: fasp_const.h:132
void fasp_blas_dstr_mxv(dSTRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_str.c:117
Main header file for FASP.
#define PRINT_MORE
Definition: fasp_const.h:82
#define PRINT_NONE
Print level for all subroutines – not including DEBUG output.
Definition: fasp_const.h:79
#define ERROR_SOLVER_TOLSMALL
Definition: fasp_const.h:51
#define MAX_RESTART
Definition: fasp_const.h:245
Sparse matrix of REAL type in CSR format.
Definition: fasp.h:148
void fasp_blas_array_axpy(const INT n, const REAL a, REAL *x, REAL *y)
y = a*x + y
Definition: blas_array.c:87
#define ERROR_SOLVER_SOLSTAG
Definition: fasp_const.h:50
#define ERROR_SOLVER_STAG
Definition: fasp_const.h:49
#define STOP_REL_RES
Definition of iterative solver stopping criteria types.
Definition: fasp_const.h:131
REAL fasp_blas_array_norminf(const INT n, const REAL *x)
Linf norm of array x.
Definition: blas_array.c:388
INT row
number of rows
Definition: fasp.h:345
Block REAL CSR matrix format.
Definition: fasp_block.h:84
#define ABS(a)
Definition: fasp.h:74
#define SHORT
FASP integer and floating point numbers.
Definition: fasp.h:63
INT fasp_solver_dcsr_spminres(dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
A preconditioned minimal residual (Minres) method for solving Au=b with safety net.
Definition: spminres.c:95
REAL fasp_blas_array_dotprod(const INT n, const REAL *x, const REAL *y)
Inner product of two arraies (x,y)
Definition: blas_array.c:267
REAL fasp_blas_array_norm2(const INT n, const REAL *x)
L2 norm of array x.
Definition: blas_array.c:347
#define BIGREAL
Some global constants.
Definition: fasp_const.h:237
INT fasp_solver_bdcsr_spminres(block_dCSRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
A preconditioned minimal residual (Minres) method for solving Au=b with safety net.
Definition: spminres.c:544
void fasp_blas_dcsr_aAxpy(const REAL alpha, dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = alpha*A*x + y.
Definition: blas_csr.c:479
void fasp_array_cp(const INT n, REAL *x, REAL *y)
Copy an array to the other y=x.
Definition: array.c:165
#define MAX_STAG
Definition: fasp_const.h:246
#define STAG_RATIO
Definition: fasp_const.h:247
#define MAX(a, b)
Definition of max, min, abs.
Definition: fasp.h:72
#define SMALLREAL
Definition: fasp_const.h:238
INT fasp_solver_dstr_spminres(dSTRmat *A, dvector *b, dvector *u, precond *pc, const REAL tol, const INT MaxIt, const SHORT stop_type, const SHORT prtlvl)
A preconditioned minimal residual (Minres) method for solving Au=b with safety net.
Definition: spminres.c:993
#define PRINT_MIN
Definition: fasp_const.h:80
void fasp_blas_dcsr_mxv(dCSRmat *A, REAL *x, REAL *y)
Matrix-vector multiplication y = A*x.
Definition: blas_csr.c:225