System Preprocessors
|
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