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