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:   PetscFunctionBegin;
 34:   PetscCall(PetscNew(&neP));

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

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

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

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

 60:   *outneP = neP;
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

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

 69:   PetscFunctionBegin;
 70:   /* Destroy work vectors and close output file */
 71:   PetscCall(VecDestroyVecs(3, &neP->workv));
 72:   err = fclose(neP->fp);
 73:   PetscCheck(!err, PETSC_COMM_SELF, PETSC_ERR_SYS, "fclose() failed on file");
 74:   PetscCall(PetscFree(neP));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

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

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

 98:   w    = neP->workv[0];
 99:   xp   = neP->workv[1];
100:   fvec = neP->workv[2];
101:   fp   = neP->fp;

103:   /* Initialize parameters */
104:   hl       = zero;
105:   hu       = zero;
106:   h        = neP->h_first_try;
107:   fnoise_s = zero;
108:   fder2_s  = zero;
109:   fcount   = neP->function_count;

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

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

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

131:     for (j = 0; j < nf - 1; j++) {
132:       for (i = 0; i < nf - j - 1; i++) tab[i][j + 1] = tab[i + 1][j] - tab[i][j];
133:     }

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

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

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

154:     /* Save fnoise and fder2. */
155:     if (*fnoise) fnoise_s = *fnoise;
156:     if (fder2) fder2_s = fder2;

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

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

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

192:   /* For now, compute h **each** MV Mult!! */
193:   /*
194:   PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg));
195:   if (!flg) {
196:     Mat mat;
197:     PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL));
198:     PetscCall(MatSNESMFMoreSetParameters(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt));
199:   }
200:   */
201:   fcount = neP->function_count - fcount;
202:   PetscCall(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));

204:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL));
205:   if (noise_test) PetscCall(JacMatMultCompare(snes, x, p, *hopt));
206:   PetscFunctionReturn(PETSC_SUCCESS);
207: }

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

224:   PetscFunctionBegin;
225:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
226:   /* Compute function and analytic Jacobian at x */
227:   PetscCall(SNESGetJacobian(snes, &Jmf, &J, NULL, NULL));
228:   PetscCall(SNESComputeJacobian(snes, x, Jmf, J));
229:   PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
230:   PetscCall(SNESComputeFunction(snes, x, f));

232:   /* Duplicate work vectors */
233:   PetscCall(VecDuplicate(x, &yy2));
234:   PetscCall(VecDuplicate(x, &yy1));

236:   /* Compute true matrix-vector product */
237:   PetscCall(MatMult(J, p, yy1));
238:   PetscCall(VecNorm(yy1, NORM_2, &yy1n));

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

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

257:     /* Compute matrix-vector product via differencing approximation */
258:     PetscCall(MatMult(Jmf, p, yy2));
259:     PetscCall(VecNorm(yy2, NORM_2, &yy2n));

261:     /* View product vector if desired */
262:     if (printv) {
263:       PetscCall(PetscSNPrintf(filename, PETSC_STATIC_ARRAY_LENGTH(filename), "y2.%d.out", (int)i));
264:       PetscCall(PetscViewerASCIIOpen(comm, filename, &view2));
265:       PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON));
266:       PetscCall(VecView(yy2, view2));
267:       PetscCall(PetscViewerPopFormat(view2));
268:       PetscCall(PetscViewerDestroy(&view2));
269:     }

271:     /* Compute relative error */
272:     PetscCall(VecAXPY(yy2, alpha, yy1));
273:     PetscCall(VecNorm(yy2, NORM_2, &enorm));
274:     enorm = enorm / yy1n;
275:     PetscCall(PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm));
276:     h *= 10.0;
277:   }
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

281: static PetscInt lin_its_total = 0;

283: PetscErrorCode SNESNoiseMonitor(SNES snes, PetscInt its, double fnorm, void *dummy)
284: {
285:   PetscInt lin_its;

287:   PetscFunctionBegin;
288:   PetscCall(SNESGetLinearSolveIterations(snes, &lin_its));
289:   lin_its_total += lin_its;
290:   PetscCall(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));

292:   PetscCall(SNESUnSetMatrixFreeParameter(snes));
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }