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;
00019
00020
00021
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
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
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
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));
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;
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
00081
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
00100
00101
00102 {
00103 PetscTruth 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;
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
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