Actual source code: snesnoise.c


  2: #include <petsc/private/snesimpl.h>

  4: PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES,Vec,void**);
  5: PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES,void*,Vec,Vec,PetscReal*,PetscReal*);
  6: PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void*);

  8: /* Data used by Jorge's diff parameter computation method */
  9: typedef struct {
 10:   Vec      *workv;           /* work vectors */
 11:   FILE     *fp;              /* output file */
 12:   PetscInt function_count;   /* count of function evaluations for diff param estimation */
 13:   double   fnoise_min;       /* minimim allowable noise */
 14:   double   hopt_min;         /* minimum allowable hopt */
 15:   double   h_first_try;      /* first try for h used in diff parameter estimate */
 16:   PetscInt fnoise_resets;    /* number of times we've reset the noise estimate */
 17:   PetscInt hopt_resets;      /* number of times we've reset the hopt estimate */
 18: } DIFFPAR_MORE;

 20: PETSC_INTERN PetscErrorCode SNESDefaultMatrixFreeSetParameters2(Mat,double,double,double);
 21: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes);
 22: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscScalar*,PetscInt*,PetscScalar*);

 24: static PetscErrorCode JacMatMultCompare(SNES,Vec,Vec,double);

 26: PetscErrorCode SNESDiffParameterCreate_More(SNES snes,Vec x,void **outneP)
 27: {
 28:   DIFFPAR_MORE   *neP;
 29:   Vec            w;
 30:   PetscRandom    rctx;  /* random number generator context */
 32:   PetscBool      flg;
 33:   char           noise_file[PETSC_MAX_PATH_LEN];

 36:   PetscNewLog(snes,&neP);

 38:   neP->function_count = 0;
 39:   neP->fnoise_min     = 1.0e-20;
 40:   neP->hopt_min       = 1.0e-8;
 41:   neP->h_first_try    = 1.0e-3;
 42:   neP->fnoise_resets  = 0;
 43:   neP->hopt_resets    = 0;

 45:   /* Create work vectors */
 46:   VecDuplicateVecs(x,3,&neP->workv);
 47:   w    = neP->workv[0];

 49:   /* Set components of vector w to random numbers */
 50:   PetscRandomCreate(PetscObjectComm((PetscObject)snes),&rctx);
 51:   PetscRandomSetFromOptions(rctx);
 52:   VecSetRandom(w,rctx);
 53:   PetscRandomDestroy(&rctx);

 55:   /* Open output file */
 56:   PetscOptionsGetString(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_noise_file",noise_file,sizeof(noise_file),&flg);
 57:   if (flg) neP->fp = fopen(noise_file,"w");
 58:   else     neP->fp = fopen("noise.out","w");
 59:   if (!neP->fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file");
 60:   PetscInfo(snes,"Creating Jorge's differencing parameter context\n");

 62:   *outneP = neP;
 63:   return(0);
 64: }

 66: PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
 67: {
 68:   DIFFPAR_MORE   *neP = (DIFFPAR_MORE*)nePv;
 70:   int            err;

 73:   /* Destroy work vectors and close output file */
 74:   VecDestroyVecs(3,&neP->workv);
 75:   err  = fclose(neP->fp);
 76:   if (err) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SYS,"fclose() failed on file");
 77:   PetscFree(neP);
 78:   return(0);
 79: }

 81: PetscErrorCode SNESDiffParameterCompute_More(SNES snes,void *nePv,Vec x,Vec p,double *fnoise,double *hopt)
 82: {
 83:   DIFFPAR_MORE   *neP = (DIFFPAR_MORE*)nePv;
 84:   Vec            w, xp, fvec;    /* work vectors to use in computing h */
 85:   double         zero = 0.0, hl, hu, h, fnoise_s, fder2_s;
 86:   PetscScalar    alpha;
 87:   PetscScalar    fval[7], tab[7][7], eps[7], f = -1;
 88:   double         rerrf = -1., fder2;
 90:   PetscInt       iter, k, i, j,  info;
 91:   PetscInt       nf = 7;         /* number of function evaluations */
 92:   PetscInt       fcount;
 93:   MPI_Comm       comm;
 94:   FILE           *fp;
 95:   PetscBool      noise_test = PETSC_FALSE;

 98:   PetscObjectGetComm((PetscObject)snes,&comm);
 99:   /* Call to SNESSetUp() just to set data structures in SNES context */
100:   if (!snes->setupcalled) {SNESSetUp(snes);}

102:   w    = neP->workv[0];
103:   xp   = neP->workv[1];
104:   fvec = neP->workv[2];
105:   fp   = neP->fp;

107:   /* Initialize parameters */
108:   hl       = zero;
109:   hu       = zero;
110:   h        = neP->h_first_try;
111:   fnoise_s = zero;
112:   fder2_s  = zero;
113:   fcount   = neP->function_count;

115:   /* We have 5 tries to attempt to compute a good hopt value */
116:   SNESGetIterationNumber(snes,&i);
117:   PetscFPrintf(comm,fp,"\n ------- SNES iteration %D ---------\n",i);
118:   for (iter=0; iter<5; iter++) {
119:     neP->h_first_try = h;

121:     /* Compute the nf function values needed to estimate the noise from
122:        the difference table */
123:     for (k=0; k<nf; k++) {
124:       alpha = h * (k+1 - (nf+1)/2);
125:       VecWAXPY(xp,alpha,p,x);
126:       SNESComputeFunction(snes,xp,fvec);
127:       neP->function_count++;
128:       VecDot(fvec,w,&fval[k]);
129:     }
130:     f = fval[(nf+1)/2 - 1];

132:     /* Construct the difference table */
133:     for (i=0; i<nf; i++) tab[i][0] = fval[i];

135:     for (j=0; j<nf-1; j++) {
136:       for (i=0; i<nf-j-1; i++) {
137:         tab[i][j+1] = tab[i+1][j] - tab[i][j];
138:       }
139:     }

141:     /* Print the difference table */
142:     PetscFPrintf(comm,fp,"Difference Table: iter = %D\n",iter);
143:     for (i=0; i<nf; i++) {
144:       for (j=0; j<nf-i; j++) {
145:         PetscFPrintf(comm,fp," %10.2e ",tab[i][j]);
146:       }
147:       PetscFPrintf(comm,fp,"\n");
148:     }

150:     /* Call the noise estimator */
151:     SNESNoise_dnest_(&nf,fval,&h,fnoise,&fder2,hopt,&info,eps);

153:     /* Output statements */
154:     rerrf = *fnoise/PetscAbsScalar(f);
155:     if (info == 1) {PetscFPrintf(comm,fp,"%s\n","Noise detected");}
156:     if (info == 2) {PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too small");}
157:     if (info == 3) {PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too large");}
158:     if (info == 4) {PetscFPrintf(comm,fp,"%s\n","Noise detected, but unreliable hopt");}
159:     PetscFPrintf(comm,fp,"Approximate epsfcn %g  %g  %g  %g  %g  %g\n",(double)eps[0],(double)eps[1],(double)eps[2],(double)eps[3],(double)eps[4],(double)eps[5]);
160:     PetscFPrintf(comm,fp,"h = %g, fnoise = %g, fder2 = %g, rerrf = %g, hopt = %g\n\n",(double)h, (double)*fnoise, (double)fder2, (double)rerrf, (double)*hopt);

162:     /* Save fnoise and fder2. */
163:     if (*fnoise) fnoise_s = *fnoise;
164:     if (fder2) fder2_s = fder2;

166:     /* Check for noise detection. */
167:     if (fnoise_s && fder2_s) {
168:       *fnoise = fnoise_s;
169:       fder2   = fder2_s;
170:       *hopt   = 1.68*sqrt(*fnoise/PetscAbsScalar(fder2));
171:       goto theend;
172:     } else {

174:       /* Update hl and hu, and determine new h */
175:       if (info == 2 || info == 4) {
176:         hl = h;
177:         if (hu == zero) h = 100*h;
178:         else            h = PetscMin(100*h,0.1*hu);
179:       } else if (info == 3) {
180:         hu = h;
181:         h  = PetscMax(1.0e-3,sqrt(hl/hu))*hu;
182:       }
183:     }
184:   }
185: theend:

187:   if (*fnoise < neP->fnoise_min) {
188:     PetscFPrintf(comm,fp,"Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n",(double)*fnoise,(double)neP->fnoise_min);
189:     *fnoise = neP->fnoise_min;
190:     neP->fnoise_resets++;
191:   }
192:   if (*hopt < neP->hopt_min) {
193:     PetscFPrintf(comm,fp,"Resetting hopt: hopt1 = %g, hopt_min = %g\n",(double)*hopt,(double)neP->hopt_min);
194:     *hopt = neP->hopt_min;
195:     neP->hopt_resets++;
196:   }

198:   PetscFPrintf(comm,fp,"Errors in derivative:\n");
199:   PetscFPrintf(comm,fp,"f = %g, fnoise = %g, fder2 = %g, hopt = %g\n",(double)f,(double)*fnoise,(double)fder2,(double)*hopt);

201:   /* For now, compute h **each** MV Mult!! */
202:   /*
203:   PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg);
204:   if (!flg) {
205:     Mat mat;
206:     SNESGetJacobian(snes,&mat,NULL,NULL);
207:     SNESDefaultMatrixFreeSetParameters2(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt);
208:   }
209:   */
210:   fcount = neP->function_count - fcount;
211:   PetscInfo5(snes,"fct_now = %D, fct_cum = %D, rerrf=%g, sqrt(noise)=%g, h_more=%g\n",fcount,neP->function_count,(double)rerrf,(double)PetscSqrtReal(*fnoise),(double)*hopt);

213:   PetscOptionsGetBool(NULL,NULL,"-noise_test",&noise_test,NULL);
214:   if (noise_test) {
215:     JacMatMultCompare(snes,x,p,*hopt);
216:   }
217:   return(0);
218: }

220: PetscErrorCode JacMatMultCompare(SNES snes,Vec x,Vec p,double hopt)
221: {
222:   Vec            yy1, yy2; /* work vectors */
223:   PetscViewer    view2;    /* viewer */
224:   Mat            J;        /* analytic Jacobian (set as preconditioner matrix) */
225:   Mat            Jmf;      /* matrix-free Jacobian (set as true system matrix) */
226:   double         h;        /* differencing parameter */
227:   Vec            f;
228:   PetscScalar    alpha;
229:   PetscReal      yy1n,yy2n,enorm;
231:   PetscInt       i;
232:   PetscBool      printv = PETSC_FALSE;
233:   char           filename[32];
234:   MPI_Comm       comm;

237:   PetscObjectGetComm((PetscObject)snes,&comm);
238:   /* Compute function and analytic Jacobian at x */
239:   SNESGetJacobian(snes,&Jmf,&J,NULL,NULL);
240:   SNESComputeJacobian(snes,x,Jmf,J);
241:   SNESGetFunction(snes,&f,NULL,NULL);
242:   SNESComputeFunction(snes,x,f);

244:   /* Duplicate work vectors */
245:   VecDuplicate(x,&yy2);
246:   VecDuplicate(x,&yy1);

248:   /* Compute true matrix-vector product */
249:   MatMult(J,p,yy1);
250:   VecNorm(yy1,NORM_2,&yy1n);

252:   /* View product vector if desired */
253:   PetscOptionsGetBool(NULL,NULL,"-print_vecs",&printv,NULL);
254:   if (printv) {
255:     PetscViewerASCIIOpen(comm,"y1.out",&view2);
256:     PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON);
257:     VecView(yy1,view2);
258:     PetscViewerPopFormat(view2);
259:     PetscViewerDestroy(&view2);
260:   }

262:   /* Test Jacobian-vector product computation */
263:   alpha = -1.0;
264:   h     = 0.01 * hopt;
265:   for (i=0; i<5; i++) {
266:     /* Set differencing parameter for matrix-free multiplication */
267:     SNESDefaultMatrixFreeSetParameters2(Jmf,PETSC_DEFAULT,PETSC_DEFAULT,h);

269:     /* Compute matrix-vector product via differencing approximation */
270:     MatMult(Jmf,p,yy2);
271:     VecNorm(yy2,NORM_2,&yy2n);

273:     /* View product vector if desired */
274:     if (printv) {
275:       sprintf(filename,"y2.%d.out",(int)i);
276:       PetscViewerASCIIOpen(comm,filename,&view2);
277:       PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON);
278:       VecView(yy2,view2);
279:       PetscViewerPopFormat(view2);
280:       PetscViewerDestroy(&view2);
281:     }

283:     /* Compute relative error */
284:     VecAXPY(yy2,alpha,yy1);
285:     VecNorm(yy2,NORM_2,&enorm);
286:     enorm = enorm/yy1n;
287:     PetscFPrintf(comm,stdout,"h = %g: relative error = %g\n",(double)h,(double)enorm);
288:     h    *= 10.0;
289:   }
290:   return(0);
291: }

293: static PetscInt lin_its_total = 0;

295: PetscErrorCode SNESNoiseMonitor(SNES snes,PetscInt its,double fnorm,void *dummy)
296: {
298:   PetscInt       lin_its;

301:   SNESGetLinearSolveIterations(snes,&lin_its);
302:   lin_its_total += lin_its;
303:   PetscPrintf(PetscObjectComm((PetscObject)snes), "iter = %D, SNES Function norm = %g, lin_its = %D, total_lin_its = %D\n",its,(double)fnorm,lin_its,lin_its_total);

305:   SNESUnSetMatrixFreeParameter(snes);
306:   return(0);
307: }