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