1 2 #include <petsc/private/snesimpl.h> 3 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 *); 7 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; 19 20 PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes); 21 PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *); 22 23 static PetscErrorCode JacMatMultCompare(SNES, Vec, Vec, double); 24 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]; 32 33 PetscFunctionBegin; 34 PetscCall(PetscNew(&neP)); 35 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; 42 43 /* Create work vectors */ 44 PetscCall(VecDuplicateVecs(x, 3, &neP->workv)); 45 w = neP->workv[0]; 46 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)); 52 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")); 59 60 *outneP = neP; 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63 64 PetscErrorCode SNESDiffParameterDestroy_More(void *nePv) 65 { 66 DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv; 67 int err; 68 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 } 77 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; 92 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)); 97 98 w = neP->workv[0]; 99 xp = neP->workv[1]; 100 fvec = neP->workv[2]; 101 fp = neP->fp; 102 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; 110 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; 116 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]; 127 128 /* Construct the difference table */ 129 for (i = 0; i < nf; i++) tab[i][0] = fval[i]; 130 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 } 134 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 } 141 142 /* Call the noise estimator */ 143 PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps)); 144 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)); 153 154 /* Save fnoise and fder2. */ 155 if (*fnoise) fnoise_s = *fnoise; 156 if (fder2) fder2_s = fder2; 157 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: 177 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 } 188 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)); 191 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)); 203 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 } 208 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; 223 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)); 231 232 /* Duplicate work vectors */ 233 PetscCall(VecDuplicate(x, &yy2)); 234 PetscCall(VecDuplicate(x, &yy1)); 235 236 /* Compute true matrix-vector product */ 237 PetscCall(MatMult(J, p, yy1)); 238 PetscCall(VecNorm(yy1, NORM_2, &yy1n)); 239 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 } 249 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)); 256 257 /* Compute matrix-vector product via differencing approximation */ 258 PetscCall(MatMult(Jmf, p, yy2)); 259 PetscCall(VecNorm(yy2, NORM_2, &yy2n)); 260 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 } 270 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 } 280 281 static PetscInt lin_its_total = 0; 282 283 PetscErrorCode SNESNoiseMonitor(SNES snes, PetscInt its, double fnorm, void *dummy) 284 { 285 PetscInt lin_its; 286 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)); 291 292 PetscCall(SNESUnSetMatrixFreeParameter(snes)); 293 PetscFunctionReturn(PETSC_SUCCESS); 294 } 295