Actual source code: snesnoise.c
1: #include <petsc/private/snesimpl.h>
3: PETSC_INTERN PetscErrorCode SNESDiffParameterCreate_More(SNES, Vec, void **);
4: PETSC_INTERN PetscErrorCode SNESDiffParameterCompute_More(SNES, void *, Vec, Vec, PetscReal *, PetscReal *);
5: PETSC_INTERN PetscErrorCode SNESDiffParameterDestroy_More(void *);
7: /* Data used by Jorge's diff parameter computation method */
8: typedef struct {
9: Vec *workv; /* work vectors */
10: FILE *fp; /* output file */
11: PetscInt function_count; /* count of function evaluations for diff param estimation */
12: double fnoise_min; /* minimum allowable noise */
13: double hopt_min; /* minimum allowable hopt */
14: double h_first_try; /* first try for h used in diff parameter estimate */
15: PetscInt fnoise_resets; /* number of times we've reset the noise estimate */
16: PetscInt hopt_resets; /* number of times we've reset the hopt estimate */
17: } DIFFPAR_MORE;
19: PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes);
20: PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);
22: static PetscErrorCode JacMatMultCompare(SNES, Vec, Vec, double);
24: PetscErrorCode SNESDiffParameterCreate_More(SNES snes, Vec x, void **outneP)
25: {
26: DIFFPAR_MORE *neP;
27: Vec w;
28: PetscRandom rctx; /* random number generator context */
29: PetscBool flg;
30: char noise_file[PETSC_MAX_PATH_LEN];
32: PetscFunctionBegin;
33: PetscCall(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: PetscCall(VecDuplicateVecs(x, 3, &neP->workv));
44: w = neP->workv[0];
46: /* Set components of vector w to random numbers */
47: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)snes), &rctx));
48: PetscCall(PetscRandomSetFromOptions(rctx));
49: PetscCall(VecSetRandom(w, rctx));
50: PetscCall(PetscRandomDestroy(&rctx));
52: /* Open output file */
53: PetscCall(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");
56: PetscCheck(neP->fp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Cannot open file");
57: PetscCall(PetscInfo(snes, "Creating Jorge's differencing parameter context\n"));
59: *outneP = neP;
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
64: {
65: DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
66: int err;
68: PetscFunctionBegin;
69: /* Destroy work vectors and close output file */
70: PetscCall(VecDestroyVecs(3, &neP->workv));
71: err = fclose(neP->fp);
72: PetscCheck(!err, PETSC_COMM_SELF, PETSC_ERR_SYS, "fclose() failed on file");
73: PetscCall(PetscFree(neP));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: PetscErrorCode SNESDiffParameterCompute_More(SNES snes, void *nePv, Vec x, Vec p, double *fnoise, double *hopt)
78: {
79: DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
80: Vec w, xp, fvec; /* work vectors to use in computing h */
81: double zero = 0.0, hl, hu, h, fnoise_s, fder2_s;
82: PetscScalar alpha;
83: PetscScalar fval[7], tab[7][7], eps[7], f = -1;
84: double rerrf = -1., fder2;
85: PetscInt iter, k, i, j, info;
86: PetscInt nf = 7; /* number of function evaluations */
87: PetscInt fcount;
88: MPI_Comm comm;
89: FILE *fp;
90: PetscBool noise_test = PETSC_FALSE;
92: PetscFunctionBegin;
93: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
94: /* Call to SNESSetUp() just to set data structures in SNES context */
95: if (!snes->setupcalled) PetscCall(SNESSetUp(snes));
97: w = neP->workv[0];
98: xp = neP->workv[1];
99: fvec = neP->workv[2];
100: fp = neP->fp;
102: /* Initialize parameters */
103: hl = zero;
104: hu = zero;
105: h = neP->h_first_try;
106: fnoise_s = zero;
107: fder2_s = zero;
108: fcount = neP->function_count;
110: /* We have 5 tries to attempt to compute a good hopt value */
111: PetscCall(SNESGetIterationNumber(snes, &i));
112: PetscCall(PetscFPrintf(comm, fp, "\n ------- SNES iteration %" PetscInt_FMT " ---------\n", i));
113: for (iter = 0; iter < 5; iter++) {
114: neP->h_first_try = h;
116: /* Compute the nf function values needed to estimate the noise from
117: the difference table */
118: for (k = 0; k < nf; k++) {
119: alpha = h * (k + 1 - (nf + 1) / 2);
120: PetscCall(VecWAXPY(xp, alpha, p, x));
121: PetscCall(SNESComputeFunction(snes, xp, fvec));
122: neP->function_count++;
123: PetscCall(VecDot(fvec, w, &fval[k]));
124: }
125: f = fval[(nf + 1) / 2 - 1];
127: /* Construct the difference table */
128: for (i = 0; i < nf; i++) tab[i][0] = fval[i];
130: for (j = 0; j < nf - 1; j++) {
131: for (i = 0; i < nf - j - 1; i++) tab[i][j + 1] = tab[i + 1][j] - tab[i][j];
132: }
134: /* Print the difference table */
135: PetscCall(PetscFPrintf(comm, fp, "Difference Table: iter = %" PetscInt_FMT "\n", iter));
136: for (i = 0; i < nf; i++) {
137: for (j = 0; j < nf - i; j++) PetscCall(PetscFPrintf(comm, fp, " %10.2e ", tab[i][j]));
138: PetscCall(PetscFPrintf(comm, fp, "\n"));
139: }
141: /* Call the noise estimator */
142: PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps));
144: /* Output statements */
145: rerrf = *fnoise / PetscAbsScalar(f);
146: if (info == 1) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected"));
147: if (info == 2) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too small"));
148: if (info == 3) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too large"));
149: if (info == 4) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected, but unreliable hopt"));
150: 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]));
151: 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));
153: /* Save fnoise and fder2. */
154: if (*fnoise) fnoise_s = *fnoise;
155: if (fder2) fder2_s = fder2;
157: /* Check for noise detection. */
158: if (fnoise_s && fder2_s) {
159: *fnoise = fnoise_s;
160: fder2 = fder2_s;
161: *hopt = 1.68 * sqrt(*fnoise / PetscAbsScalar(fder2));
162: goto theend;
163: } else {
164: /* Update hl and hu, and determine new h */
165: if (info == 2 || info == 4) {
166: hl = h;
167: if (hu == zero) h = 100 * h;
168: else h = PetscMin(100 * h, 0.1 * hu);
169: } else if (info == 3) {
170: hu = h;
171: h = PetscMax(1.0e-3, sqrt(hl / hu)) * hu;
172: }
173: }
174: }
175: theend:
177: if (*fnoise < neP->fnoise_min) {
178: PetscCall(PetscFPrintf(comm, fp, "Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n", (double)*fnoise, (double)neP->fnoise_min));
179: *fnoise = neP->fnoise_min;
180: neP->fnoise_resets++;
181: }
182: if (*hopt < neP->hopt_min) {
183: PetscCall(PetscFPrintf(comm, fp, "Resetting hopt: hopt1 = %g, hopt_min = %g\n", (double)*hopt, (double)neP->hopt_min));
184: *hopt = neP->hopt_min;
185: neP->hopt_resets++;
186: }
188: PetscCall(PetscFPrintf(comm, fp, "Errors in derivative:\n"));
189: PetscCall(PetscFPrintf(comm, fp, "f = %g, fnoise = %g, fder2 = %g, hopt = %g\n", (double)f, (double)*fnoise, (double)fder2, (double)*hopt));
191: /* For now, compute h **each** MV Mult!! */
192: /*
193: PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg));
194: if (!flg) {
195: Mat mat;
196: PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL));
197: PetscCall(MatSNESMFMoreSetParameters(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt));
198: }
199: */
200: fcount = neP->function_count - fcount;
201: 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));
203: PetscCall(PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL));
204: if (noise_test) PetscCall(JacMatMultCompare(snes, x, p, *hopt));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode JacMatMultCompare(SNES snes, Vec x, Vec p, double hopt)
209: {
210: Vec yy1, yy2; /* work vectors */
211: PetscViewer view2; /* viewer */
212: Mat J; /* analytic Jacobian (set as preconditioner matrix) */
213: Mat Jmf; /* matrix-free Jacobian (set as true system matrix) */
214: double h; /* differencing parameter */
215: Vec f;
216: PetscScalar alpha;
217: PetscReal yy1n, yy2n, enorm;
218: PetscInt i;
219: PetscBool printv = PETSC_FALSE;
220: char filename[32];
221: MPI_Comm comm;
223: PetscFunctionBegin;
224: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
225: /* Compute function and analytic Jacobian at x */
226: PetscCall(SNESGetJacobian(snes, &Jmf, &J, NULL, NULL));
227: PetscCall(SNESComputeJacobian(snes, x, Jmf, J));
228: PetscCall(SNESGetFunction(snes, &f, NULL, NULL));
229: PetscCall(SNESComputeFunction(snes, x, f));
231: /* Duplicate work vectors */
232: PetscCall(VecDuplicate(x, &yy2));
233: PetscCall(VecDuplicate(x, &yy1));
235: /* Compute true matrix-vector product */
236: PetscCall(MatMult(J, p, yy1));
237: PetscCall(VecNorm(yy1, NORM_2, &yy1n));
239: /* View product vector if desired */
240: PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_vecs", &printv, NULL));
241: if (printv) {
242: PetscCall(PetscViewerASCIIOpen(comm, "y1.out", &view2));
243: PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON));
244: PetscCall(VecView(yy1, view2));
245: PetscCall(PetscViewerPopFormat(view2));
246: PetscCall(PetscViewerDestroy(&view2));
247: }
249: /* Test Jacobian-vector product computation */
250: alpha = -1.0;
251: h = 0.01 * hopt;
252: for (i = 0; i < 5; i++) {
253: /* Set differencing parameter for matrix-free multiplication */
254: PetscCall(MatSNESMFMoreSetParameters(Jmf, PETSC_DEFAULT, PETSC_DEFAULT, h));
256: /* Compute matrix-vector product via differencing approximation */
257: PetscCall(MatMult(Jmf, p, yy2));
258: PetscCall(VecNorm(yy2, NORM_2, &yy2n));
260: /* View product vector if desired */
261: if (printv) {
262: PetscCall(PetscSNPrintf(filename, PETSC_STATIC_ARRAY_LENGTH(filename), "y2.%d.out", (int)i));
263: PetscCall(PetscViewerASCIIOpen(comm, filename, &view2));
264: PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON));
265: PetscCall(VecView(yy2, view2));
266: PetscCall(PetscViewerPopFormat(view2));
267: PetscCall(PetscViewerDestroy(&view2));
268: }
270: /* Compute relative error */
271: PetscCall(VecAXPY(yy2, alpha, yy1));
272: PetscCall(VecNorm(yy2, NORM_2, &enorm));
273: enorm = enorm / yy1n;
274: PetscCall(PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm));
275: h *= 10.0;
276: }
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }