1 #include <petsc/private/snesimpl.h>
2
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 *);
6
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;
18
19 PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes);
20 PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *);
21
22 static PetscErrorCode JacMatMultCompare(SNES, Vec, Vec, double);
23
SNESDiffParameterCreate_More(SNES snes,Vec x,void ** outneP)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];
31
32 PetscFunctionBegin;
33 PetscCall(PetscNew(&neP));
34
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;
41
42 /* Create work vectors */
43 PetscCall(VecDuplicateVecs(x, 3, &neP->workv));
44 w = neP->workv[0];
45
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));
51
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"));
58
59 *outneP = neP;
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
SNESDiffParameterDestroy_More(void * nePv)63 PetscErrorCode SNESDiffParameterDestroy_More(void *nePv)
64 {
65 DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv;
66 int err;
67
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 }
76
SNESDiffParameterCompute_More(SNES snes,void * nePv,Vec x,Vec p,double * fnoise,double * hopt)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;
91
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));
96
97 w = neP->workv[0];
98 xp = neP->workv[1];
99 fvec = neP->workv[2];
100 fp = neP->fp;
101
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;
109
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;
115
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];
126
127 /* Construct the difference table */
128 for (i = 0; i < nf; i++) tab[i][0] = fval[i];
129
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 }
133
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 }
140
141 /* Call the noise estimator */
142 PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps));
143
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));
152
153 /* Save fnoise and fder2. */
154 if (*fnoise) fnoise_s = *fnoise;
155 if (fder2) fder2_s = fder2;
156
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:
176
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 }
187
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));
190
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));
202
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 }
207
JacMatMultCompare(SNES snes,Vec x,Vec p,double hopt)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 matrix from which to compute the preconditioner) */
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;
222
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));
230
231 /* Duplicate work vectors */
232 PetscCall(VecDuplicate(x, &yy2));
233 PetscCall(VecDuplicate(x, &yy1));
234
235 /* Compute true matrix-vector product */
236 PetscCall(MatMult(J, p, yy1));
237 PetscCall(VecNorm(yy1, NORM_2, &yy1n));
238
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 }
248
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));
255
256 /* Compute matrix-vector product via differencing approximation */
257 PetscCall(MatMult(Jmf, p, yy2));
258 PetscCall(VecNorm(yy2, NORM_2, &yy2n));
259
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 }
269
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 }
279