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; /* minimim 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 SNESDefaultMatrixFreeSetParameters2(Mat, double, double, double); 21 PETSC_INTERN PetscErrorCode SNESUnSetMatrixFreeParameter(SNES snes); 22 PETSC_INTERN PetscErrorCode SNESNoise_dnest_(PetscInt *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscScalar *, PetscInt *, PetscScalar *); 23 24 static PetscErrorCode JacMatMultCompare(SNES, Vec, Vec, double); 25 26 PetscErrorCode SNESDiffParameterCreate_More(SNES snes, Vec x, void **outneP) { 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(PetscNewLog(snes, &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(0); 62 } 63 64 PetscErrorCode SNESDiffParameterDestroy_More(void *nePv) { 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(0); 75 } 76 77 PetscErrorCode SNESDiffParameterCompute_More(SNES snes, void *nePv, Vec x, Vec p, double *fnoise, double *hopt) { 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; 90 91 PetscFunctionBegin; 92 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 93 /* Call to SNESSetUp() just to set data structures in SNES context */ 94 if (!snes->setupcalled) PetscCall(SNESSetUp(snes)); 95 96 w = neP->workv[0]; 97 xp = neP->workv[1]; 98 fvec = neP->workv[2]; 99 fp = neP->fp; 100 101 /* Initialize parameters */ 102 hl = zero; 103 hu = zero; 104 h = neP->h_first_try; 105 fnoise_s = zero; 106 fder2_s = zero; 107 fcount = neP->function_count; 108 109 /* We have 5 tries to attempt to compute a good hopt value */ 110 PetscCall(SNESGetIterationNumber(snes, &i)); 111 PetscCall(PetscFPrintf(comm, fp, "\n ------- SNES iteration %" PetscInt_FMT " ---------\n", i)); 112 for (iter = 0; iter < 5; iter++) { 113 neP->h_first_try = h; 114 115 /* Compute the nf function values needed to estimate the noise from 116 the difference table */ 117 for (k = 0; k < nf; k++) { 118 alpha = h * (k + 1 - (nf + 1) / 2); 119 PetscCall(VecWAXPY(xp, alpha, p, x)); 120 PetscCall(SNESComputeFunction(snes, xp, fvec)); 121 neP->function_count++; 122 PetscCall(VecDot(fvec, w, &fval[k])); 123 } 124 f = fval[(nf + 1) / 2 - 1]; 125 126 /* Construct the difference table */ 127 for (i = 0; i < nf; i++) tab[i][0] = fval[i]; 128 129 for (j = 0; j < nf - 1; j++) { 130 for (i = 0; i < nf - j - 1; i++) { tab[i][j + 1] = tab[i + 1][j] - tab[i][j]; } 131 } 132 133 /* Print the difference table */ 134 PetscCall(PetscFPrintf(comm, fp, "Difference Table: iter = %" PetscInt_FMT "\n", iter)); 135 for (i = 0; i < nf; i++) { 136 for (j = 0; j < nf - i; j++) { PetscCall(PetscFPrintf(comm, fp, " %10.2e ", tab[i][j])); } 137 PetscCall(PetscFPrintf(comm, fp, "\n")); 138 } 139 140 /* Call the noise estimator */ 141 PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps)); 142 143 /* Output statements */ 144 rerrf = *fnoise / PetscAbsScalar(f); 145 if (info == 1) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected")); 146 if (info == 2) { PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too small")); } 147 if (info == 3) { PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too large")); } 148 if (info == 4) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected, but unreliable hopt")); 149 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])); 150 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)); 151 152 /* Save fnoise and fder2. */ 153 if (*fnoise) fnoise_s = *fnoise; 154 if (fder2) fder2_s = fder2; 155 156 /* Check for noise detection. */ 157 if (fnoise_s && fder2_s) { 158 *fnoise = fnoise_s; 159 fder2 = fder2_s; 160 *hopt = 1.68 * sqrt(*fnoise / PetscAbsScalar(fder2)); 161 goto theend; 162 } else { 163 /* Update hl and hu, and determine new h */ 164 if (info == 2 || info == 4) { 165 hl = h; 166 if (hu == zero) h = 100 * h; 167 else h = PetscMin(100 * h, 0.1 * hu); 168 } else if (info == 3) { 169 hu = h; 170 h = PetscMax(1.0e-3, sqrt(hl / hu)) * hu; 171 } 172 } 173 } 174 theend: 175 176 if (*fnoise < neP->fnoise_min) { 177 PetscCall(PetscFPrintf(comm, fp, "Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n", (double)*fnoise, (double)neP->fnoise_min)); 178 *fnoise = neP->fnoise_min; 179 neP->fnoise_resets++; 180 } 181 if (*hopt < neP->hopt_min) { 182 PetscCall(PetscFPrintf(comm, fp, "Resetting hopt: hopt1 = %g, hopt_min = %g\n", (double)*hopt, (double)neP->hopt_min)); 183 *hopt = neP->hopt_min; 184 neP->hopt_resets++; 185 } 186 187 PetscCall(PetscFPrintf(comm, fp, "Errors in derivative:\n")); 188 PetscCall(PetscFPrintf(comm, fp, "f = %g, fnoise = %g, fder2 = %g, hopt = %g\n", (double)f, (double)*fnoise, (double)fder2, (double)*hopt)); 189 190 /* For now, compute h **each** MV Mult!! */ 191 /* 192 PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg)); 193 if (!flg) { 194 Mat mat; 195 PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL)); 196 PetscCall(SNESDefaultMatrixFreeSetParameters2(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt)); 197 } 198 */ 199 fcount = neP->function_count - fcount; 200 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)); 201 202 PetscCall(PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL)); 203 if (noise_test) PetscCall(JacMatMultCompare(snes, x, p, *hopt)); 204 PetscFunctionReturn(0); 205 } 206 207 PetscErrorCode JacMatMultCompare(SNES snes, Vec x, Vec p, double hopt) { 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; 220 221 PetscFunctionBegin; 222 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 223 /* Compute function and analytic Jacobian at x */ 224 PetscCall(SNESGetJacobian(snes, &Jmf, &J, NULL, NULL)); 225 PetscCall(SNESComputeJacobian(snes, x, Jmf, J)); 226 PetscCall(SNESGetFunction(snes, &f, NULL, NULL)); 227 PetscCall(SNESComputeFunction(snes, x, f)); 228 229 /* Duplicate work vectors */ 230 PetscCall(VecDuplicate(x, &yy2)); 231 PetscCall(VecDuplicate(x, &yy1)); 232 233 /* Compute true matrix-vector product */ 234 PetscCall(MatMult(J, p, yy1)); 235 PetscCall(VecNorm(yy1, NORM_2, &yy1n)); 236 237 /* View product vector if desired */ 238 PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_vecs", &printv, NULL)); 239 if (printv) { 240 PetscCall(PetscViewerASCIIOpen(comm, "y1.out", &view2)); 241 PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON)); 242 PetscCall(VecView(yy1, view2)); 243 PetscCall(PetscViewerPopFormat(view2)); 244 PetscCall(PetscViewerDestroy(&view2)); 245 } 246 247 /* Test Jacobian-vector product computation */ 248 alpha = -1.0; 249 h = 0.01 * hopt; 250 for (i = 0; i < 5; i++) { 251 /* Set differencing parameter for matrix-free multiplication */ 252 PetscCall(SNESDefaultMatrixFreeSetParameters2(Jmf, PETSC_DEFAULT, PETSC_DEFAULT, h)); 253 254 /* Compute matrix-vector product via differencing approximation */ 255 PetscCall(MatMult(Jmf, p, yy2)); 256 PetscCall(VecNorm(yy2, NORM_2, &yy2n)); 257 258 /* View product vector if desired */ 259 if (printv) { 260 sprintf(filename, "y2.%d.out", (int)i); 261 PetscCall(PetscViewerASCIIOpen(comm, filename, &view2)); 262 PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON)); 263 PetscCall(VecView(yy2, view2)); 264 PetscCall(PetscViewerPopFormat(view2)); 265 PetscCall(PetscViewerDestroy(&view2)); 266 } 267 268 /* Compute relative error */ 269 PetscCall(VecAXPY(yy2, alpha, yy1)); 270 PetscCall(VecNorm(yy2, NORM_2, &enorm)); 271 enorm = enorm / yy1n; 272 PetscCall(PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm)); 273 h *= 10.0; 274 } 275 PetscFunctionReturn(0); 276 } 277 278 static PetscInt lin_its_total = 0; 279 280 PetscErrorCode SNESNoiseMonitor(SNES snes, PetscInt its, double fnorm, void *dummy) { 281 PetscInt lin_its; 282 283 PetscFunctionBegin; 284 PetscCall(SNESGetLinearSolveIterations(snes, &lin_its)); 285 lin_its_total += lin_its; 286 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)); 287 288 PetscCall(SNESUnSetMatrixFreeParameter(snes)); 289 PetscFunctionReturn(0); 290 } 291