System Preprocessors
kspmonitor.c
Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include <stdio.h>
00003 #include <math.h>
00004 #include "syspro.h"
00005 #include "sysprolinear.h"
00006 
00007 #define ITER_STAGNATION -21
00008 #define ITER_DIVERGENCE -22
00009 
00010 #undef __FUNCT__
00011 #define __FUNCT__ "estimate_completion_from_hist"
00012 static PetscErrorCode estimate_completion_from_hist
00013 (double *hist,int n,double rtol,int tracing,int *est)
00014 {
00015   PetscReal t; int i,j,m;
00016   PetscFunctionBegin;
00017     
00018   *est = 0; /* default case: no estimate returned */
00019 
00020   /*
00021    * first detect stagnation
00022    */
00023   {
00024     double l_avg=0., rlo,rhi; int p;
00025     if (n<10) p=n; else p=10;
00026     for (i=n-p; i<n; i++) l_avg += hist[i];
00027     l_avg /= p; rlo = (1-1.e-3)*l_avg; rhi = (1+1.e-3)*l_avg;
00028     for (i=n-p; i<n; i++)
00029       if (hist[i]<rlo || hist[i]>rhi) goto go_on;
00030     *est = ITER_STAGNATION;
00031     if (tracing)
00032       printf(".. method seems to stagnate (at %d:%e over %d its)\n",
00033              n,l_avg,p);
00034     PetscFunctionReturn(0);
00035   }
00036   go_on:
00037   
00038   /* detect highest error value */
00039   t = 0; m = 0;
00040   for (i=0; i<n; i++) {
00041     if (hist[i]>t) {t=hist[i]; m=i;}
00042   }
00043   /*
00044    * if highest error was in last iteration; we call it quits; otherwise:
00045    */
00046   if (m==n-1) {
00047     *est = ITER_DIVERGENCE;
00048     if (tracing) printf(".. method seems to diverge\n");
00049     return 0;
00050   } else {
00051     if (tracing) 
00052       printf(".. error decrease %d--%d: %e--%e\n",m,n-1,t,hist[n-1]);
00053     /* let's tighten the interval */
00054     for (i=m+1; i<(m+n-1)/2; i++) {
00055       for (j=i+1; j<n; j++)
00056         if (hist[j]>hist[i]) goto out;
00057       m = i; t = hist[m];}
00058   out:
00059     if (tracing) 
00060       printf(".. error decrease %d--%d: %e--%e\n",m,n-1,t,hist[n-1]);
00061     t = pow(hist[n-1]/t,1./(n-m)); /* reduction per iteration */
00062     if (tracing) 
00063       printf(".. reduction estimated as %e\n",t);
00064     if (t>.999) {
00065       *est = ITER_STAGNATION;
00066       if (tracing) 
00067         printf(".. method seems to stagnate (reduction too low)\n");
00068     } else {
00069       *est = log(rtol)/log(t)+m; /* use access function into info? */
00070       if (*est<1.1*n) *est=1.1*n;
00071     }
00072   }
00073 
00074   PetscFunctionReturn(0);
00075 }
00076 
00077 extern int gmrescycleid;
00078 #undef __FUNCT__
00079 #define __FUNCT__ "MonitorAdjustMaxit"
00080 /*! This routine analyzes the convergence history, and if the iterative
00081   method is still making progress, extends the maximum number of iterations.
00082  */
00083 PetscErrorCode MonitorAdjustMaxit(KSP ksp,int it,PetscReal cg_err,void *data)
00084 {
00085   PC pc; Mat A; MPI_Comm comm; PetscReal rtol,*h;
00086   int M,nh,est,i,maxit; PetscErrorCode ierr;
00087 
00088   PetscFunctionBegin;
00089   ierr = PetscObjectGetComm((PetscObject)ksp,&comm); CHKERRQ(ierr);
00090 
00091   ierr = KSPGetTolerances
00092     (ksp,&rtol,PETSC_NULL,PETSC_NULL,&maxit); CHKERRQ(ierr);
00093   if (it<maxit) PetscFunctionReturn(0);
00094 
00095   PetscPrintf(comm,".. end of rope @ it=%d\n",it);
00096   ierr = KSPGetResidualHistory(ksp,&h,&nh); CHKERRQ(ierr);
00097 
00098   /*
00099    * if this is gmres, make sure we do whole cycles;
00100    * for instance, stagnation at the start of a cycle is not fatal
00101    */
00102   {
00103     PetscBool flg; int cycle;
00104     ierr = PetscObjectComposedDataGetInt
00105       ((PetscObject)ksp,gmrescycleid,cycle,flg); CHKERRQ(ierr);
00106     if (flg && cycle>0 && it%cycle != 0) {
00107       est = cycle*( 1+(int)(it/cycle) );
00108       goto adjust;
00109     }
00110   }
00111 
00112   ierr = KSPGetPC(ksp,&pc); CHKERRQ(ierr);
00113   ierr = PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00114   ierr = MatGetSize(A,&M,PETSC_NULL); CHKERRQ(ierr);
00115 
00116   if (nh<10)
00117     est = 20; /* this is a kludge: someone is mucking with maxit */
00118   else {
00119     ierr = estimate_completion_from_hist
00120       (h,nh,rtol,1,&est); CHKERRQ(ierr);
00121     switch (est) {
00122     case ITER_STAGNATION :
00123       PetscPrintf(comm,".. method stagnates\n");
00124       PetscFunctionReturn(0); break;
00125     case ITER_DIVERGENCE :
00126       PetscPrintf(comm,".. method diverges\n");
00127       PetscFunctionReturn(0); break;
00128     }
00129   }
00130 
00131   if (est<=it) {
00132     /*PetscPrintf(comm,".. this should not happen: %d\n",est);*/
00133     est = 1.5 * it;
00134   }
00135  adjust:
00136   if (est>M) {
00137     PetscPrintf(comm,".. this will take too long: %d\n",est);
00138   } else if (est>it) {
00139     PetscReal *new_h;
00140     PetscPrintf(comm,".. extending to %d\n",est);
00141     ierr = KSPSetTolerances
00142       (ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,est); CHKERRQ(ierr);
00143     ierr = PetscMalloc(est*sizeof(PetscReal),&new_h); CHKERRQ(ierr);
00144     for (i=0; i<it; i++) new_h[i] = h[i];
00145     ierr = PetscFree(h); CHKERRQ(ierr);
00146     ierr = KSPSetResidualHistory(ksp,new_h,est,PETSC_TRUE); CHKERRQ(ierr);
00147   }
00148   PetscFunctionReturn(0);
00149 }
00150