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;     /* minimum 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 SNESUnSetMatrixFreeParameter(SNES snes);
 21: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);

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

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

 33:   PetscNew(&neP);

 35:   neP->function_count = 0;
 36:   neP->fnoise_min     = 1.0e-20;
 37:   neP->hopt_min       = 1.0e-8;
 38:   neP->h_first_try    = 1.0e-3;
 39:   neP->fnoise_resets  = 0;
 40:   neP->hopt_resets    = 0;

 42:   /* Create work vectors */
 43:   VecDuplicateVecs(x, 3, &neP->workv);
 44:   w = neP->workv[0];

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

 52:   /* Open output file */
 53:   PetscOptionsGetString(((PetscObject)snes)->options, ((PetscObject)snes)->prefix, "-snes_mf_noise_file", noise_file, sizeof(noise_file), &flg);
 54:   if (flg) neP->fp = fopen(noise_file, "w");
 55:   else neP->fp = fopen("noise.out", "w");
 57:   PetscInfo(snes, "Creating Jorge's differencing parameter context\n");

 59:   *outneP = neP;
 60:   return 0;
 61: }

 63: PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
 64: {
 65:   DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
 66:   int           err;

 68:   /* Destroy work vectors and close output file */
 69:   VecDestroyVecs(3, &neP->workv);
 70:   err = fclose(neP->fp);
 72:   PetscFree(neP);
 73:   return 0;
 74: }

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

 91:   PetscObjectGetComm((PetscObject)snes, &comm);
 92:   /* Call to SNESSetUp() just to set data structures in SNES context */
 93:   if (!snes->setupcalled) SNESSetUp(snes);

 95:   w    = neP->workv[0];
 96:   xp   = neP->workv[1];
 97:   fvec = neP->workv[2];
 98:   fp   = neP->fp;

100:   /* Initialize parameters */
101:   hl       = zero;
102:   hu       = zero;
103:   h        = neP->h_first_try;
104:   fnoise_s = zero;
105:   fder2_s  = zero;
106:   fcount   = neP->function_count;

108:   /* We have 5 tries to attempt to compute a good hopt value */
109:   SNESGetIterationNumber(snes, &i);
110:   PetscFPrintf(comm, fp, "\n ------- SNES iteration %" PetscInt_FMT " ---------\n", i);
111:   for (iter = 0; iter < 5; iter++) {
112:     neP->h_first_try = h;

114:     /* Compute the nf function values needed to estimate the noise from
115:        the difference table */
116:     for (k = 0; k < nf; k++) {
117:       alpha = h * (k + 1 - (nf + 1) / 2);
118:       VecWAXPY(xp, alpha, p, x);
119:       SNESComputeFunction(snes, xp, fvec);
120:       neP->function_count++;
121:       VecDot(fvec, w, &fval[k]);
122:     }
123:     f = fval[(nf + 1) / 2 - 1];

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

128:     for (j = 0; j < nf - 1; j++) {
129:       for (i = 0; i < nf - j - 1; i++) tab[i][j + 1] = tab[i + 1][j] - tab[i][j];
130:     }

132:     /* Print the difference table */
133:     PetscFPrintf(comm, fp, "Difference Table: iter = %" PetscInt_FMT "\n", iter);
134:     for (i = 0; i < nf; i++) {
135:       for (j = 0; j < nf - i; j++) PetscFPrintf(comm, fp, " %10.2e ", tab[i][j]);
136:       PetscFPrintf(comm, fp, "\n");
137:     }

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

142:     /* Output statements */
143:     rerrf = *fnoise / PetscAbsScalar(f);
144:     if (info == 1) PetscFPrintf(comm, fp, "%s\n", "Noise detected");
145:     if (info == 2) PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too small");
146:     if (info == 3) PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too large");
147:     if (info == 4) PetscFPrintf(comm, fp, "%s\n", "Noise detected, but unreliable hopt");
148:     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]);
149:     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);

151:     /* Save fnoise and fder2. */
152:     if (*fnoise) fnoise_s = *fnoise;
153:     if (fder2) fder2_s = fder2;

155:     /* Check for noise detection. */
156:     if (fnoise_s && fder2_s) {
157:       *fnoise = fnoise_s;
158:       fder2   = fder2_s;
159:       *hopt   = 1.68 * sqrt(*fnoise / PetscAbsScalar(fder2));
160:       goto theend;
161:     } else {
162:       /* Update hl and hu, and determine new h */
163:       if (info == 2 || info == 4) {
164:         hl = h;
165:         if (hu == zero) h = 100 * h;
166:         else h = PetscMin(100 * h, 0.1 * hu);
167:       } else if (info == 3) {
168:         hu = h;
169:         h  = PetscMax(1.0e-3, sqrt(hl / hu)) * hu;
170:       }
171:     }
172:   }
173: theend:

175:   if (*fnoise < neP->fnoise_min) {
176:     PetscFPrintf(comm, fp, "Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n", (double)*fnoise, (double)neP->fnoise_min);
177:     *fnoise = neP->fnoise_min;
178:     neP->fnoise_resets++;
179:   }
180:   if (*hopt < neP->hopt_min) {
181:     PetscFPrintf(comm, fp, "Resetting hopt: hopt1 = %g, hopt_min = %g\n", (double)*hopt, (double)neP->hopt_min);
182:     *hopt = neP->hopt_min;
183:     neP->hopt_resets++;
184:   }

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

189:   /* For now, compute h **each** MV Mult!! */
190:   /*
191:   PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg);
192:   if (!flg) {
193:     Mat mat;
194:     SNESGetJacobian(snes,&mat,NULL,NULL);
195:     MatSNESMFMoreSetParameters(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt);
196:   }
197:   */
198:   fcount = neP->function_count - fcount;
199:   PetscInfo(snes, "fct_now = %" PetscInt_FMT ", fct_cum = %" PetscInt_FMT ", rerrf=%g, sqrt(noise)=%g, h_more=%g\n", fcount, neP->function_count, (double)rerrf, (double)PetscSqrtReal(*fnoise), (double)*hopt);

201:   PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL);
202:   if (noise_test) JacMatMultCompare(snes, x, p, *hopt);
203:   return 0;
204: }

206: PetscErrorCode JacMatMultCompare(SNES snes, Vec x, Vec p, double hopt)
207: {
208:   Vec         yy1, yy2; /* work vectors */
209:   PetscViewer view2;    /* viewer */
210:   Mat         J;        /* analytic Jacobian (set as preconditioner matrix) */
211:   Mat         Jmf;      /* matrix-free Jacobian (set as true system matrix) */
212:   double      h;        /* differencing parameter */
213:   Vec         f;
214:   PetscScalar alpha;
215:   PetscReal   yy1n, yy2n, enorm;
216:   PetscInt    i;
217:   PetscBool   printv = PETSC_FALSE;
218:   char        filename[32];
219:   MPI_Comm    comm;

221:   PetscObjectGetComm((PetscObject)snes, &comm);
222:   /* Compute function and analytic Jacobian at x */
223:   SNESGetJacobian(snes, &Jmf, &J, NULL, NULL);
224:   SNESComputeJacobian(snes, x, Jmf, J);
225:   SNESGetFunction(snes, &f, NULL, NULL);
226:   SNESComputeFunction(snes, x, f);

228:   /* Duplicate work vectors */
229:   VecDuplicate(x, &yy2);
230:   VecDuplicate(x, &yy1);

232:   /* Compute true matrix-vector product */
233:   MatMult(J, p, yy1);
234:   VecNorm(yy1, NORM_2, &yy1n);

236:   /* View product vector if desired */
237:   PetscOptionsGetBool(NULL, NULL, "-print_vecs", &printv, NULL);
238:   if (printv) {
239:     PetscViewerASCIIOpen(comm, "y1.out", &view2);
240:     PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON);
241:     VecView(yy1, view2);
242:     PetscViewerPopFormat(view2);
243:     PetscViewerDestroy(&view2);
244:   }

246:   /* Test Jacobian-vector product computation */
247:   alpha = -1.0;
248:   h     = 0.01 * hopt;
249:   for (i = 0; i < 5; i++) {
250:     /* Set differencing parameter for matrix-free multiplication */
251:     MatSNESMFMoreSetParameters(Jmf, PETSC_DEFAULT, PETSC_DEFAULT, h);

253:     /* Compute matrix-vector product via differencing approximation */
254:     MatMult(Jmf, p, yy2);
255:     VecNorm(yy2, NORM_2, &yy2n);

257:     /* View product vector if desired */
258:     if (printv) {
259:       sprintf(filename, "y2.%d.out", (int)i);
260:       PetscViewerASCIIOpen(comm, filename, &view2);
261:       PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON);
262:       VecView(yy2, view2);
263:       PetscViewerPopFormat(view2);
264:       PetscViewerDestroy(&view2);
265:     }

267:     /* Compute relative error */
268:     VecAXPY(yy2, alpha, yy1);
269:     VecNorm(yy2, NORM_2, &enorm);
270:     enorm = enorm / yy1n;
271:     PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm);
272:     h *= 10.0;
273:   }
274:   return 0;
275: }

277: static PetscInt lin_its_total = 0;

279: PetscErrorCode SNESNoiseMonitor(SNES snes, PetscInt its, double fnorm, void *dummy)
280: {
281:   PetscInt lin_its;

283:   SNESGetLinearSolveIterations(snes, &lin_its);
284:   lin_its_total += lin_its;
285:   PetscPrintf(PetscObjectComm((PetscObject)snes), "iter = %" PetscInt_FMT ", SNES Function norm = %g, lin_its = %" PetscInt_FMT ", total_lin_its = %" PetscInt_FMT "\n", its, (double)fnorm, lin_its, lin_its_total);

287:   SNESUnSetMatrixFreeParameter(snes);
288:   return 0;
289: }