Fast Auxiliary Space Preconditioning  1.8.4 Feb/15/2016
pminres.c
Go to the documentation of this file.
1 
57 #include <math.h>
58 
59 #include "fasp.h"
60 #include "fasp_functs.h"
61 #include "itsolver_util.inl"
62 
63 /*---------------------------------*/
64 /*-- Public Functions --*/
65 /*---------------------------------*/
66 
93  dvector *b,
94  dvector *u,
95  precond *pc,
96  const REAL tol,
97  const INT MaxIt,
98  const SHORT stop_type,
99  const SHORT prtlvl)
100 {
101  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
102  const INT m = b->row;
103  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
104  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
105 
106  // local variables
107  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
108  REAL absres0 = BIGREAL, absres = BIGREAL;
109  REAL normr0 = BIGREAL, relres = BIGREAL;
110  REAL normu2, normuu, normp, infnormu, factor;
111  REAL alpha, alpha0, alpha1, temp2;
112 
113  // allocate temp memory (need 11*m REAL)
114  REAL *work=(REAL *)fasp_mem_calloc(11*m,sizeof(REAL));
115  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m;
116  REAL *t0=z1+m, *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m;
117 
118 #if DEBUG_MODE > 0
119  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
120  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
121 #endif
122 
123  // p0 = 0
124  fasp_array_set(m,p0,0.0);
125 
126  // r = b-A*u
127  fasp_array_cp(m,b->val,r);
128  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
129 
130  // p1 = B(r)
131  if ( pc != NULL )
132  pc->fct(r,p1,pc->data); /* Apply preconditioner */
133  else
134  fasp_array_cp(m,r,p1); /* No preconditioner */
135 
136  // compute initial residuals
137  switch ( stop_type ) {
138  case STOP_REL_RES:
139  absres0 = fasp_blas_array_norm2(m,r);
140  normr0 = MAX(SMALLREAL,absres0);
141  relres = absres0/normr0;
142  break;
143  case STOP_REL_PRECRES:
144  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
145  normr0 = MAX(SMALLREAL,absres0);
146  relres = absres0/normr0;
147  break;
148  case STOP_MOD_REL_RES:
149  absres0 = fasp_blas_array_norm2(m,r);
150  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
151  relres = absres0/normu2;
152  break;
153  default:
154  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
155  goto FINISHED;
156  }
157 
158  // if initial residual is small, no need to iterate!
159  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
160 
161  // output iteration information if needed
162  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
163 
164  // tp = A*p1
165  fasp_blas_dcsr_mxv(A,p1,tp);
166 
167  // tz = B(tp)
168  if ( pc != NULL )
169  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
170  else
171  fasp_array_cp(m,tp,tz); /* No preconditioner */
172 
173  // p1 = p1/normp
174  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
175  normp = sqrt(normp);
176  fasp_array_cp(m,p1,t);
177  fasp_array_set(m,p1,0.0);
178  fasp_blas_array_axpy(m,1/normp,t,p1);
179 
180  // t0 = A*p0 = 0
181  fasp_array_set(m,t0,0.0);
182  fasp_array_cp(m,t0,z0);
183  fasp_array_cp(m,t0,t1);
184  fasp_array_cp(m,t0,z1);
185 
186  // t1 = tp/normp, z1 = tz/normp
187  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
188  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
189 
190  // main MinRes loop
191  while ( iter++ < MaxIt ) {
192 
193  // alpha = <r,z1>
194  alpha=fasp_blas_array_dotprod(m,r,z1);
195 
196  // u = u+alpha*p1
197  fasp_blas_array_axpy(m,alpha,p1,u->val);
198 
199  // r = r-alpha*Ap1
200  fasp_blas_array_axpy(m,-alpha,t1,r);
201 
202  // compute t = A*z1 alpha1 = <z1,t>
203  fasp_blas_dcsr_mxv(A,z1,t);
204  alpha1=fasp_blas_array_dotprod(m,z1,t);
205 
206  // compute t = A*z0 alpha0 = <z1,t>
207  fasp_blas_dcsr_mxv(A,z0,t);
208  alpha0=fasp_blas_array_dotprod(m,z1,t);
209 
210  // p2 = z1-alpha1*p1-alpha0*p0
211  fasp_array_cp(m,z1,p2);
212  fasp_blas_array_axpy(m,-alpha1,p1,p2);
213  fasp_blas_array_axpy(m,-alpha0,p0,p2);
214 
215  // tp = A*p2
216  fasp_blas_dcsr_mxv(A,p2,tp);
217 
218  // tz = B(tp)
219  if ( pc != NULL )
220  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
221  else
222  fasp_array_cp(m,tp,tz); /* No preconditioner */
223 
224  // p2 = p2/normp
225  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
226  normp = sqrt(normp);
227  fasp_array_cp(m,p2,t);
228  fasp_array_set(m,p2,0.0);
229  fasp_blas_array_axpy(m,1/normp,t,p2);
230 
231  // prepare for the next iteration
232  fasp_array_cp(m,p1,p0);
233  fasp_array_cp(m,p2,p1);
234  fasp_array_cp(m,t1,t0);
235  fasp_array_cp(m,z1,z0);
236 
237  // t1=tp/normp,z1=tz/normp
238  fasp_array_set(m,t1,0.0);
239  fasp_array_cp(m,t1,z1);
240  fasp_blas_array_axpy(m,1/normp,tp,t1);
241  fasp_blas_array_axpy(m,1/normp,tz,z1);
242 
243  normu2 = fasp_blas_array_norm2(m,u->val);
244 
245  // compute residuals
246  switch ( stop_type ) {
247  case STOP_REL_RES:
248  temp2 = fasp_blas_array_dotprod(m,r,r);
249  absres = sqrt(temp2);
250  relres = absres/normr0;
251  break;
252  case STOP_REL_PRECRES:
253  if (pc == NULL)
254  fasp_array_cp(m,r,t);
255  else
256  pc->fct(r,t,pc->data);
257  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
258  absres = sqrt(temp2);
259  relres = absres/normr0;
260  break;
261  case STOP_MOD_REL_RES:
262  temp2 = fasp_blas_array_dotprod(m,r,r);
263  absres = sqrt(temp2);
264  relres = absres/normu2;
265  break;
266  }
267 
268  // compute reducation factor of residual ||r||
269  factor = absres/absres0;
270 
271  // output iteration information if needed
272  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
273 
274  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
275  infnormu = fasp_blas_array_norminf(m, u->val);
276  if (infnormu <= sol_inf_tol) {
277  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
278  iter = ERROR_SOLVER_SOLSTAG;
279  break;
280  }
281 
282  // Check II: if staggenated, try to restart
283  normuu = fasp_blas_array_norm2(m,p1);
284  normuu = ABS(alpha)*(normuu/normu2);
285 
286  if ( normuu < maxdiff ) {
287 
288  if ( stag < MaxStag ) {
289  if ( prtlvl >= PRINT_MORE ) {
290  ITS_DIFFRES(normuu,relres);
291  ITS_RESTART;
292  }
293  }
294 
295  fasp_array_cp(m,b->val,r);
296  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
297 
298  // compute residuals
299  switch (stop_type) {
300  case STOP_REL_RES:
301  temp2 = fasp_blas_array_dotprod(m,r,r);
302  absres = sqrt(temp2);
303  relres = absres/normr0;
304  break;
305  case STOP_REL_PRECRES:
306  if (pc == NULL)
307  fasp_array_cp(m,r,t);
308  else
309  pc->fct(r,t,pc->data);
310  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
311  absres = sqrt(temp2);
312  relres = absres/normr0;
313  break;
314  case STOP_MOD_REL_RES:
315  temp2 = fasp_blas_array_dotprod(m,r,r);
316  absres = sqrt(temp2);
317  relres = absres/normu2;
318  break;
319  }
320 
321  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
322 
323  if ( relres < tol )
324  break;
325  else {
326  if ( stag >= MaxStag ) {
327  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
328  iter = ERROR_SOLVER_STAG;
329  break;
330  }
331  fasp_array_set(m,p0,0.0);
332  ++stag;
333  ++restart_step;
334 
335  // p1 = B(r)
336  if ( pc != NULL )
337  pc->fct(r,p1,pc->data); /* Apply preconditioner */
338  else
339  fasp_array_cp(m,r,p1); /* No preconditioner */
340 
341  // tp = A*p1
342  fasp_blas_dcsr_mxv(A,p1,tp);
343 
344  // tz = B(tp)
345  if ( pc != NULL )
346  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
347  else
348  fasp_array_cp(m,tp,tz); /* No preconditioner */
349 
350  // p1 = p1/normp
351  normp = fasp_blas_array_dotprod(m,tz,tp);
352  normp = sqrt(normp);
353  fasp_array_cp(m,p1,t);
354 
355  // t0 = A*p0=0
356  fasp_array_set(m,t0,0.0);
357  fasp_array_cp(m,t0,z0);
358  fasp_array_cp(m,t0,t1);
359  fasp_array_cp(m,t0,z1);
360  fasp_array_cp(m,t0,p1);
361 
362  fasp_blas_array_axpy(m,1/normp,t,p1);
363 
364  // t1 = tp/normp, z1 = tz/normp
365  fasp_blas_array_axpy(m,1/normp,tp,t1);
366  fasp_blas_array_axpy(m,1/normp,tz,z1);
367  }
368  }
369 
370  // Check III: prevent false convergence
371  if ( relres < tol ) {
372 
373  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
374 
375  // compute residual r = b - Ax again
376  fasp_array_cp(m,b->val,r);
377  fasp_blas_dcsr_aAxpy(-1.0,A,u->val,r);
378 
379  // compute residuals
380  switch (stop_type) {
381  case STOP_REL_RES:
382  temp2 = fasp_blas_array_dotprod(m,r,r);
383  absres = sqrt(temp2);
384  relres = absres/normr0;
385  break;
386  case STOP_REL_PRECRES:
387  if (pc == NULL)
388  fasp_array_cp(m,r,t);
389  else
390  pc->fct(r,t,pc->data);
391  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
392  absres = sqrt(temp2);
393  relres = absres/normr0;
394  break;
395  case STOP_MOD_REL_RES:
396  temp2 = fasp_blas_array_dotprod(m,r,r);
397  absres = sqrt(temp2);
398  relres = absres/normu2;
399  break;
400  }
401 
402  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
403 
404  // check convergence
405  if ( relres < tol ) break;
406 
407  if ( more_step >= MaxRestartStep ) {
408  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
409  iter = ERROR_SOLVER_TOLSMALL;
410  break;
411  }
412 
413  // prepare for restarting the method
414  fasp_array_set(m,p0,0.0);
415  ++more_step;
416  ++restart_step;
417 
418  // p1 = B(r)
419  if ( pc != NULL )
420  pc->fct(r,p1,pc->data); /* Apply preconditioner */
421  else
422  fasp_array_cp(m,r,p1); /* No preconditioner */
423 
424  // tp = A*p1
425  fasp_blas_dcsr_mxv(A,p1,tp);
426 
427  // tz = B(tp)
428  if ( pc != NULL )
429  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
430  else
431  fasp_array_cp(m,tp,tz); /* No preconditioner */
432 
433  // p1 = p1/normp
434  normp = fasp_blas_array_dotprod(m,tz,tp);
435  normp = sqrt(normp);
436  fasp_array_cp(m,p1,t);
437 
438  // t0 = A*p0 = 0
439  fasp_array_set(m,t0,0.0);
440  fasp_array_cp(m,t0,z0);
441  fasp_array_cp(m,t0,t1);
442  fasp_array_cp(m,t0,z1);
443  fasp_array_cp(m,t0,p1);
444 
445  fasp_blas_array_axpy(m,1/normp,t,p1);
446 
447  // t1=tp/normp,z1=tz/normp
448  fasp_blas_array_axpy(m,1/normp,tp,t1);
449  fasp_blas_array_axpy(m,1/normp,tz,z1);
450 
451  } // end of convergence check
452 
453  // update relative residual here
454  absres0 = absres;
455 
456  } // end of the main loop
457 
458 FINISHED: // finish the iterative method
459  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
460 
461  // clean up temp memory
462  fasp_mem_free(work);
463 
464 #if DEBUG_MODE > 0
465  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
466 #endif
467 
468  if ( iter > MaxIt )
469  return ERROR_SOLVER_MAXIT;
470  else
471  return iter;
472 }
473 
500  dvector *b,
501  dvector *u,
502  precond *pc,
503  const REAL tol,
504  const INT MaxIt,
505  const SHORT stop_type,
506  const SHORT prtlvl)
507 {
508  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
509  const INT m = b->row;
510  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
511  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
512 
513  // local variables
514  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
515  REAL absres0 = BIGREAL, absres = BIGREAL;
516  REAL normr0 = BIGREAL, relres = BIGREAL;
517  REAL normu2, normuu, normp, infnormu, factor;
518  REAL alpha, alpha0, alpha1, temp2;
519 
520  // allocate temp memory (need 11*m REAL)
521  REAL *work=(REAL *)fasp_mem_calloc(11*m,sizeof(REAL));
522  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m;
523  REAL *t0=z1+m, *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m;
524 
525 #if DEBUG_MODE > 0
526  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
527  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
528 #endif
529 
530  // p0 = 0
531  fasp_array_set(m,p0,0.0);
532 
533  // r = b-A*u
534  fasp_array_cp(m,b->val,r);
535  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
536 
537  // p1 = B(r)
538  if ( pc != NULL )
539  pc->fct(r,p1,pc->data); /* Apply preconditioner */
540  else
541  fasp_array_cp(m,r,p1); /* No preconditioner */
542 
543  // compute initial residuals
544  switch ( stop_type ) {
545  case STOP_REL_RES:
546  absres0 = fasp_blas_array_norm2(m,r);
547  normr0 = MAX(SMALLREAL,absres0);
548  relres = absres0/normr0;
549  break;
550  case STOP_REL_PRECRES:
551  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
552  normr0 = MAX(SMALLREAL,absres0);
553  relres = absres0/normr0;
554  break;
555  case STOP_MOD_REL_RES:
556  absres0 = fasp_blas_array_norm2(m,r);
557  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
558  relres = absres0/normu2;
559  break;
560  default:
561  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
562  goto FINISHED;
563  }
564 
565  // if initial residual is small, no need to iterate!
566  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
567 
568  // output iteration information if needed
569  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
570 
571  // tp = A*p1
572  fasp_blas_bdcsr_mxv(A,p1,tp);
573 
574  // tz = B(tp)
575  if ( pc != NULL )
576  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
577  else
578  fasp_array_cp(m,tp,tz); /* No preconditioner */
579 
580  // p1 = p1/normp
581  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
582  normp = sqrt(normp);
583  fasp_array_cp(m,p1,t);
584  fasp_array_set(m,p1,0.0);
585  fasp_blas_array_axpy(m,1/normp,t,p1);
586 
587  // t0 = A*p0 = 0
588  fasp_array_set(m,t0,0.0);
589  fasp_array_cp(m,t0,z0);
590  fasp_array_cp(m,t0,t1);
591  fasp_array_cp(m,t0,z1);
592 
593  // t1 = tp/normp, z1 = tz/normp
594  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
595  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
596 
597  // main MinRes loop
598  while ( iter++ < MaxIt ) {
599 
600  // alpha = <r,z1>
601  alpha=fasp_blas_array_dotprod(m,r,z1);
602 
603  // u = u+alpha*p1
604  fasp_blas_array_axpy(m,alpha,p1,u->val);
605 
606  // r = r-alpha*Ap1
607  fasp_blas_array_axpy(m,-alpha,t1,r);
608 
609  // compute t = A*z1 alpha1 = <z1,t>
610  fasp_blas_bdcsr_mxv(A,z1,t);
611  alpha1=fasp_blas_array_dotprod(m,z1,t);
612 
613  // compute t = A*z0 alpha0 = <z1,t>
614  fasp_blas_bdcsr_mxv(A,z0,t);
615  alpha0=fasp_blas_array_dotprod(m,z1,t);
616 
617  // p2 = z1-alpha1*p1-alpha0*p0
618  fasp_array_cp(m,z1,p2);
619  fasp_blas_array_axpy(m,-alpha1,p1,p2);
620  fasp_blas_array_axpy(m,-alpha0,p0,p2);
621 
622  // tp = A*p2
623  fasp_blas_bdcsr_mxv(A,p2,tp);
624 
625  // tz = B(tp)
626  if ( pc != NULL )
627  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
628  else
629  fasp_array_cp(m,tp,tz); /* No preconditioner */
630 
631  // p2 = p2/normp
632  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
633  normp = sqrt(normp);
634  fasp_array_cp(m,p2,t);
635  fasp_array_set(m,p2,0.0);
636  fasp_blas_array_axpy(m,1/normp,t,p2);
637 
638  // prepare for the next iteration
639  fasp_array_cp(m,p1,p0);
640  fasp_array_cp(m,p2,p1);
641  fasp_array_cp(m,t1,t0);
642  fasp_array_cp(m,z1,z0);
643 
644  // t1=tp/normp,z1=tz/normp
645  fasp_array_set(m,t1,0.0);
646  fasp_array_cp(m,t1,z1);
647  fasp_blas_array_axpy(m,1/normp,tp,t1);
648  fasp_blas_array_axpy(m,1/normp,tz,z1);
649 
650  normu2 = fasp_blas_array_norm2(m,u->val);
651 
652  // compute residuals
653  switch ( stop_type ) {
654  case STOP_REL_RES:
655  temp2 = fasp_blas_array_dotprod(m,r,r);
656  absres = sqrt(temp2);
657  relres = absres/normr0;
658  break;
659  case STOP_REL_PRECRES:
660  if (pc == NULL)
661  fasp_array_cp(m,r,t);
662  else
663  pc->fct(r,t,pc->data);
664  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
665  absres = sqrt(temp2);
666  relres = absres/normr0;
667  break;
668  case STOP_MOD_REL_RES:
669  temp2 = fasp_blas_array_dotprod(m,r,r);
670  absres = sqrt(temp2);
671  relres = absres/normu2;
672  break;
673  }
674 
675  // compute reducation factor of residual ||r||
676  factor = absres/absres0;
677 
678  // output iteration information if needed
679  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
680 
681  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
682  infnormu = fasp_blas_array_norminf(m, u->val);
683  if (infnormu <= sol_inf_tol) {
684  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
685  iter = ERROR_SOLVER_SOLSTAG;
686  break;
687  }
688 
689  // Check II: if staggenated, try to restart
690  normuu = fasp_blas_array_norm2(m,p1);
691  normuu = ABS(alpha)*(normuu/normu2);
692 
693  if ( normuu < maxdiff ) {
694 
695  if ( stag < MaxStag ) {
696  if ( prtlvl >= PRINT_MORE ) {
697  ITS_DIFFRES(normuu,relres);
698  ITS_RESTART;
699  }
700  }
701 
702  fasp_array_cp(m,b->val,r);
703  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
704 
705  // compute residuals
706  switch (stop_type) {
707  case STOP_REL_RES:
708  temp2 = fasp_blas_array_dotprod(m,r,r);
709  absres = sqrt(temp2);
710  relres = absres/normr0;
711  break;
712  case STOP_REL_PRECRES:
713  if (pc == NULL)
714  fasp_array_cp(m,r,t);
715  else
716  pc->fct(r,t,pc->data);
717  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
718  absres = sqrt(temp2);
719  relres = absres/normr0;
720  break;
721  case STOP_MOD_REL_RES:
722  temp2 = fasp_blas_array_dotprod(m,r,r);
723  absres = sqrt(temp2);
724  relres = absres/normu2;
725  break;
726  }
727 
728  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
729 
730  if ( relres < tol )
731  break;
732  else {
733  if ( stag >= MaxStag ) {
734  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
735  iter = ERROR_SOLVER_STAG;
736  break;
737  }
738  fasp_array_set(m,p0,0.0);
739  ++stag;
740  ++restart_step;
741 
742  // p1 = B(r)
743  if ( pc != NULL )
744  pc->fct(r,p1,pc->data); /* Apply preconditioner */
745  else
746  fasp_array_cp(m,r,p1); /* No preconditioner */
747 
748  // tp = A*p1
749  fasp_blas_bdcsr_mxv(A,p1,tp);
750 
751  // tz = B(tp)
752  if ( pc != NULL )
753  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
754  else
755  fasp_array_cp(m,tp,tz); /* No preconditioner */
756 
757  // p1 = p1/normp
758  normp = fasp_blas_array_dotprod(m,tz,tp);
759  normp = sqrt(normp);
760  fasp_array_cp(m,p1,t);
761 
762  // t0 = A*p0=0
763  fasp_array_set(m,t0,0.0);
764  fasp_array_cp(m,t0,z0);
765  fasp_array_cp(m,t0,t1);
766  fasp_array_cp(m,t0,z1);
767  fasp_array_cp(m,t0,p1);
768 
769  fasp_blas_array_axpy(m,1/normp,t,p1);
770 
771  // t1 = tp/normp, z1 = tz/normp
772  fasp_blas_array_axpy(m,1/normp,tp,t1);
773  fasp_blas_array_axpy(m,1/normp,tz,z1);
774  }
775  }
776 
777  // Check III: prevent false convergence
778  if ( relres < tol ) {
779 
780  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
781 
782  // compute residual r = b - Ax again
783  fasp_array_cp(m,b->val,r);
784  fasp_blas_bdcsr_aAxpy(-1.0,A,u->val,r);
785 
786  // compute residuals
787  switch (stop_type) {
788  case STOP_REL_RES:
789  temp2 = fasp_blas_array_dotprod(m,r,r);
790  absres = sqrt(temp2);
791  relres = absres/normr0;
792  break;
793  case STOP_REL_PRECRES:
794  if (pc == NULL)
795  fasp_array_cp(m,r,t);
796  else
797  pc->fct(r,t,pc->data);
798  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
799  absres = sqrt(temp2);
800  relres = absres/normr0;
801  break;
802  case STOP_MOD_REL_RES:
803  temp2 = fasp_blas_array_dotprod(m,r,r);
804  absres = sqrt(temp2);
805  relres = absres/normu2;
806  break;
807  }
808 
809  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
810 
811  // check convergence
812  if ( relres < tol ) break;
813 
814  if ( more_step >= MaxRestartStep ) {
815  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
816  iter = ERROR_SOLVER_TOLSMALL;
817  break;
818  }
819 
820  // prepare for restarting the method
821  fasp_array_set(m,p0,0.0);
822  ++more_step;
823  ++restart_step;
824 
825  // p1 = B(r)
826  if ( pc != NULL )
827  pc->fct(r,p1,pc->data); /* Apply preconditioner */
828  else
829  fasp_array_cp(m,r,p1); /* No preconditioner */
830 
831  // tp = A*p1
832  fasp_blas_bdcsr_mxv(A,p1,tp);
833 
834  // tz = B(tp)
835  if ( pc != NULL )
836  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
837  else
838  fasp_array_cp(m,tp,tz); /* No preconditioner */
839 
840  // p1 = p1/normp
841  normp = fasp_blas_array_dotprod(m,tz,tp);
842  normp = sqrt(normp);
843  fasp_array_cp(m,p1,t);
844 
845  // t0 = A*p0 = 0
846  fasp_array_set(m,t0,0.0);
847  fasp_array_cp(m,t0,z0);
848  fasp_array_cp(m,t0,t1);
849  fasp_array_cp(m,t0,z1);
850  fasp_array_cp(m,t0,p1);
851 
852  fasp_blas_array_axpy(m,1/normp,t,p1);
853 
854  // t1=tp/normp,z1=tz/normp
855  fasp_blas_array_axpy(m,1/normp,tp,t1);
856  fasp_blas_array_axpy(m,1/normp,tz,z1);
857 
858  } // end of convergence check
859 
860  // update relative residual here
861  absres0 = absres;
862 
863  } // end of the main loop
864 
865 FINISHED: // finish the iterative method
866  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
867 
868  // clean up temp memory
869  fasp_mem_free(work);
870 
871 #if DEBUG_MODE > 0
872  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
873 #endif
874 
875  if ( iter > MaxIt )
876  return ERROR_SOLVER_MAXIT;
877  else
878  return iter;
879 }
880 
903  dvector *b,
904  dvector *u,
905  precond *pc,
906  const REAL tol,
907  const INT MaxIt,
908  const SHORT stop_type,
909  const SHORT prtlvl)
910 {
911  const SHORT MaxStag = MAX_STAG, MaxRestartStep = MAX_RESTART;
912  const INT m = b->row;
913  const REAL maxdiff = tol*STAG_RATIO; // staganation tolerance
914  const REAL sol_inf_tol = SMALLREAL; // infinity norm tolerance
915 
916  // local variables
917  INT iter = 0, stag = 1, more_step = 1, restart_step = 1;
918  REAL absres0 = BIGREAL, absres = BIGREAL;
919  REAL normr0 = BIGREAL, relres = BIGREAL;
920  REAL normu2, normuu, normp, infnormu, factor;
921  REAL alpha, alpha0, alpha1, temp2;
922 
923  // allocate temp memory (need 11*m REAL)
924  REAL *work=(REAL *)fasp_mem_calloc(11*m,sizeof(REAL));
925  REAL *p0=work, *p1=work+m, *p2=p1+m, *z0=p2+m, *z1=z0+m;
926  REAL *t0=z1+m, *t1=t0+m, *t=t1+m, *tp=t+m, *tz=tp+m, *r=tz+m;
927 
928 #if DEBUG_MODE > 0
929  printf("### DEBUG: %s ...... [Start]\n", __FUNCTION__);
930  printf("### DEBUG: maxit = %d, tol = %.4le\n", MaxIt, tol);
931 #endif
932 
933  // p0 = 0
934  fasp_array_set(m,p0,0.0);
935 
936  // r = b-A*u
937  fasp_array_cp(m,b->val,r);
938  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
939 
940  // p1 = B(r)
941  if ( pc != NULL )
942  pc->fct(r,p1,pc->data); /* Apply preconditioner */
943  else
944  fasp_array_cp(m,r,p1); /* No preconditioner */
945 
946  // compute initial residuals
947  switch ( stop_type ) {
948  case STOP_REL_RES:
949  absres0 = fasp_blas_array_norm2(m,r);
950  normr0 = MAX(SMALLREAL,absres0);
951  relres = absres0/normr0;
952  break;
953  case STOP_REL_PRECRES:
954  absres0 = sqrt(fasp_blas_array_dotprod(m,r,p1));
955  normr0 = MAX(SMALLREAL,absres0);
956  relres = absres0/normr0;
957  break;
958  case STOP_MOD_REL_RES:
959  absres0 = fasp_blas_array_norm2(m,r);
960  normu2 = MAX(SMALLREAL,fasp_blas_array_norm2(m,u->val));
961  relres = absres0/normu2;
962  break;
963  default:
964  printf("### ERROR: Unrecognized stopping type for %s!\n", __FUNCTION__);
965  goto FINISHED;
966  }
967 
968  // if initial residual is small, no need to iterate!
969  if ( relres < tol || absres0 < 1e-3*tol ) goto FINISHED;
970 
971  // output iteration information if needed
972  print_itinfo(prtlvl,stop_type,iter,relres,absres0,0.0);
973 
974  // tp = A*p1
975  fasp_blas_dstr_mxv(A,p1,tp);
976 
977  // tz = B(tp)
978  if ( pc != NULL )
979  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
980  else
981  fasp_array_cp(m,tp,tz); /* No preconditioner */
982 
983  // p1 = p1/normp
984  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
985  normp = sqrt(normp);
986  fasp_array_cp(m,p1,t);
987  fasp_array_set(m,p1,0.0);
988  fasp_blas_array_axpy(m,1/normp,t,p1);
989 
990  // t0 = A*p0 = 0
991  fasp_array_set(m,t0,0.0);
992  fasp_array_cp(m,t0,z0);
993  fasp_array_cp(m,t0,t1);
994  fasp_array_cp(m,t0,z1);
995 
996  // t1 = tp/normp, z1 = tz/normp
997  fasp_blas_array_axpy(m,1.0/normp,tp,t1);
998  fasp_blas_array_axpy(m,1.0/normp,tz,z1);
999 
1000  // main MinRes loop
1001  while ( iter++ < MaxIt ) {
1002 
1003  // alpha = <r,z1>
1004  alpha=fasp_blas_array_dotprod(m,r,z1);
1005 
1006  // u = u+alpha*p1
1007  fasp_blas_array_axpy(m,alpha,p1,u->val);
1008 
1009  // r = r-alpha*Ap1
1010  fasp_blas_array_axpy(m,-alpha,t1,r);
1011 
1012  // compute t = A*z1 alpha1 = <z1,t>
1013  fasp_blas_dstr_mxv(A,z1,t);
1014  alpha1=fasp_blas_array_dotprod(m,z1,t);
1015 
1016  // compute t = A*z0 alpha0 = <z1,t>
1017  fasp_blas_dstr_mxv(A,z0,t);
1018  alpha0=fasp_blas_array_dotprod(m,z1,t);
1019 
1020  // p2 = z1-alpha1*p1-alpha0*p0
1021  fasp_array_cp(m,z1,p2);
1022  fasp_blas_array_axpy(m,-alpha1,p1,p2);
1023  fasp_blas_array_axpy(m,-alpha0,p0,p2);
1024 
1025  // tp = A*p2
1026  fasp_blas_dstr_mxv(A,p2,tp);
1027 
1028  // tz = B(tp)
1029  if ( pc != NULL )
1030  pc->fct(tp,tz,pc->data); /* Apply preconditioner */
1031  else
1032  fasp_array_cp(m,tp,tz); /* No preconditioner */
1033 
1034  // p2 = p2/normp
1035  normp = ABS(fasp_blas_array_dotprod(m,tz,tp));
1036  normp = sqrt(normp);
1037  fasp_array_cp(m,p2,t);
1038  fasp_array_set(m,p2,0.0);
1039  fasp_blas_array_axpy(m,1/normp,t,p2);
1040 
1041  // prepare for the next iteration
1042  fasp_array_cp(m,p1,p0);
1043  fasp_array_cp(m,p2,p1);
1044  fasp_array_cp(m,t1,t0);
1045  fasp_array_cp(m,z1,z0);
1046 
1047  // t1=tp/normp,z1=tz/normp
1048  fasp_array_set(m,t1,0.0);
1049  fasp_array_cp(m,t1,z1);
1050  fasp_blas_array_axpy(m,1/normp,tp,t1);
1051  fasp_blas_array_axpy(m,1/normp,tz,z1);
1052 
1053  normu2 = fasp_blas_array_norm2(m,u->val);
1054 
1055  // compute residuals
1056  switch ( stop_type ) {
1057  case STOP_REL_RES:
1058  temp2 = fasp_blas_array_dotprod(m,r,r);
1059  absres = sqrt(temp2);
1060  relres = absres/normr0;
1061  break;
1062  case STOP_REL_PRECRES:
1063  if (pc == NULL)
1064  fasp_array_cp(m,r,t);
1065  else
1066  pc->fct(r,t,pc->data);
1067  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1068  absres = sqrt(temp2);
1069  relres = absres/normr0;
1070  break;
1071  case STOP_MOD_REL_RES:
1072  temp2 = fasp_blas_array_dotprod(m,r,r);
1073  absres = sqrt(temp2);
1074  relres = absres/normu2;
1075  break;
1076  }
1077 
1078  // compute reducation factor of residual ||r||
1079  factor = absres/absres0;
1080 
1081  // output iteration information if needed
1082  print_itinfo(prtlvl,stop_type,iter,relres,absres,factor);
1083 
1084  // Check I: if soultion is close to zero, return ERROR_SOLVER_SOLSTAG
1085  infnormu = fasp_blas_array_norminf(m, u->val);
1086  if (infnormu <= sol_inf_tol) {
1087  if ( prtlvl > PRINT_MIN ) ITS_ZEROSOL;
1088  iter = ERROR_SOLVER_SOLSTAG;
1089  break;
1090  }
1091 
1092  // Check II: if staggenated, try to restart
1093  normuu = fasp_blas_array_norm2(m,p1);
1094  normuu = ABS(alpha)*(normuu/normu2);
1095 
1096  if ( normuu < maxdiff ) {
1097 
1098  if ( stag < MaxStag ) {
1099  if ( prtlvl >= PRINT_MORE ) {
1100  ITS_DIFFRES(normuu,relres);
1101  ITS_RESTART;
1102  }
1103  }
1104 
1105  fasp_array_cp(m,b->val,r);
1106  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1107 
1108  // compute residuals
1109  switch (stop_type) {
1110  case STOP_REL_RES:
1111  temp2 = fasp_blas_array_dotprod(m,r,r);
1112  absres = sqrt(temp2);
1113  relres = absres/normr0;
1114  break;
1115  case STOP_REL_PRECRES:
1116  if (pc == NULL)
1117  fasp_array_cp(m,r,t);
1118  else
1119  pc->fct(r,t,pc->data);
1120  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1121  absres = sqrt(temp2);
1122  relres = absres/normr0;
1123  break;
1124  case STOP_MOD_REL_RES:
1125  temp2 = fasp_blas_array_dotprod(m,r,r);
1126  absres = sqrt(temp2);
1127  relres = absres/normu2;
1128  break;
1129  }
1130 
1131  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1132 
1133  if ( relres < tol )
1134  break;
1135  else {
1136  if ( stag >= MaxStag ) {
1137  if ( prtlvl > PRINT_MIN ) ITS_STAGGED;
1138  iter = ERROR_SOLVER_STAG;
1139  break;
1140  }
1141  fasp_array_set(m,p0,0.0);
1142  ++stag;
1143  ++restart_step;
1144 
1145  // p1 = B(r)
1146  if ( pc != NULL )
1147  pc->fct(r,p1,pc->data); /* Apply preconditioner */
1148  else
1149  fasp_array_cp(m,r,p1); /* No preconditioner */
1150 
1151  // tp = A*p1
1152  fasp_blas_dstr_mxv(A,p1,tp);
1153 
1154  // tz = B(tp)
1155  if ( pc != NULL )
1156  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1157  else
1158  fasp_array_cp(m,tp,tz); /* No preconditioner */
1159 
1160  // p1 = p1/normp
1161  normp = fasp_blas_array_dotprod(m,tz,tp);
1162  normp = sqrt(normp);
1163  fasp_array_cp(m,p1,t);
1164 
1165  // t0 = A*p0=0
1166  fasp_array_set(m,t0,0.0);
1167  fasp_array_cp(m,t0,z0);
1168  fasp_array_cp(m,t0,t1);
1169  fasp_array_cp(m,t0,z1);
1170  fasp_array_cp(m,t0,p1);
1171 
1172  fasp_blas_array_axpy(m,1/normp,t,p1);
1173 
1174  // t1 = tp/normp, z1 = tz/normp
1175  fasp_blas_array_axpy(m,1/normp,tp,t1);
1176  fasp_blas_array_axpy(m,1/normp,tz,z1);
1177  }
1178  }
1179 
1180  // Check III: prevent false convergence
1181  if ( relres < tol ) {
1182 
1183  if ( prtlvl >= PRINT_MORE ) ITS_COMPRES(relres);
1184 
1185  // compute residual r = b - Ax again
1186  fasp_array_cp(m,b->val,r);
1187  fasp_blas_dstr_aAxpy(-1.0,A,u->val,r);
1188 
1189  // compute residuals
1190  switch (stop_type) {
1191  case STOP_REL_RES:
1192  temp2 = fasp_blas_array_dotprod(m,r,r);
1193  absres = sqrt(temp2);
1194  relres = absres/normr0;
1195  break;
1196  case STOP_REL_PRECRES:
1197  if (pc == NULL)
1198  fasp_array_cp(m,r,t);
1199  else
1200  pc->fct(r,t,pc->data);
1201  temp2 = ABS(fasp_blas_array_dotprod(m,r,t));
1202  absres = sqrt(temp2);
1203  relres = absres/normr0;
1204  break;
1205  case STOP_MOD_REL_RES:
1206  temp2 = fasp_blas_array_dotprod(m,r,r);
1207  absres = sqrt(temp2);
1208  relres = absres/normu2;
1209  break;
1210  }
1211 
1212  if ( prtlvl >= PRINT_MORE ) ITS_REALRES(relres);
1213 
1214  // check convergence
1215  if ( relres < tol ) break;
1216 
1217  if ( more_step >= MaxRestartStep ) {
1218  if ( prtlvl > PRINT_MIN ) ITS_ZEROTOL;
1219  iter = ERROR_SOLVER_TOLSMALL;
1220  break;
1221  }
1222 
1223  // prepare for restarting the method
1224  fasp_array_set(m,p0,0.0);
1225  ++more_step;
1226  ++restart_step;
1227 
1228  // p1 = B(r)
1229  if ( pc != NULL )
1230  pc->fct(r,p1,pc->data); /* Apply preconditioner */
1231  else
1232  fasp_array_cp(m,r,p1); /* No preconditioner */
1233 
1234  // tp = A*p1
1235  fasp_blas_dstr_mxv(A,p1,tp);
1236 
1237  // tz = B(tp)
1238  if ( pc != NULL )
1239  pc->fct(tp,tz,pc->data); /* Apply rreconditioner */
1240  else
1241  fasp_array_cp(m,tp,tz); /* No preconditioner */
1242 
1243  // p1 = p1/normp
1244  normp = fasp_blas_array_dotprod(m,tz,tp);
1245  normp = sqrt(normp);
1246  fasp_array_cp(m,p1,t);
1247 
1248  // t0 = A*p0 = 0
1249  fasp_array_set(m,t0,0.0);
1250  fasp_array_cp(m,t0,z0);
1251  fasp_array_cp(m,t0,t1);
1252  fasp_array_cp(m,t0,z1);
1253  fasp_array_cp(m,t0,p1);
1254 
1255  fasp_blas_array_axpy(m,1/normp,t,p1);
1256 
1257  // t1=tp/normp,z1=tz/normp
1258  fasp_blas_array_axpy(m,1/normp,tp,t1);
1259  fasp_blas_array_axpy(m,1/normp,tz,z1);
1260 
1261  } // end of convergence check
1262 
1263  // update relative residual here
1264  absres0 = absres;
1265 
1266  } // end of the main loop
1267 
1268 FINISHED: // finish the iterative method
1269  if ( prtlvl > PRINT_NONE ) ITS_FINAL(iter,MaxIt,relres);
1270 
1271  // clean up temp memory
1272  fasp_mem_free(work);
1273 
1274 #if DEBUG_MODE > 0
1275  printf("### DEBUG: %s ...... [Finish]\n", __FUNCTION__);
1276 #endif
1277 
1278  if ( iter > MaxIt )
1279  return ERROR_SOLVER_MAXIT;
1280  else
1281  return iter;
1282 }
1283 
1284 /*---------------------------------*/
1285 /*-- End of File --*/
1286 /*---------------------------------*/
INT fasp_solver_bdcsr_pminres(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.
Definition: pminres.c:499
#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
INT fasp_solver_dstr_pminres(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.
Definition: pminres.c:902
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
INT fasp_solver_dcsr_pminres(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.
Definition: pminres.c:92
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
#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
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
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
#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