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: }