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 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 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(PetscNewLog(snes, &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(0); 61 } 62 63 PetscErrorCode SNESDiffParameterDestroy_More(void *nePv) { 64 DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv; 65 int err; 66 67 PetscFunctionBegin; 68 /* Destroy work vectors and close output file */ 69 PetscCall(VecDestroyVecs(3, &neP->workv)); 70 err = fclose(neP->fp); 71 PetscCheck(!err, PETSC_COMM_SELF, PETSC_ERR_SYS, "fclose() failed on file"); 72 PetscCall(PetscFree(neP)); 73 PetscFunctionReturn(0); 74 } 75 76 PetscErrorCode SNESDiffParameterCompute_More(SNES snes, void *nePv, Vec x, Vec p, double *fnoise, double *hopt) { 77 DIFFPAR_MORE *neP = (DIFFPAR_MORE *)nePv; 78 Vec w, xp, fvec; /* work vectors to use in computing h */ 79 double zero = 0.0, hl, hu, h, fnoise_s, fder2_s; 80 PetscScalar alpha; 81 PetscScalar fval[7], tab[7][7], eps[7], f = -1; 82 double rerrf = -1., fder2; 83 PetscInt iter, k, i, j, info; 84 PetscInt nf = 7; /* number of function evaluations */ 85 PetscInt fcount; 86 MPI_Comm comm; 87 FILE *fp; 88 PetscBool noise_test = PETSC_FALSE; 89 90 PetscFunctionBegin; 91 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 92 /* Call to SNESSetUp() just to set data structures in SNES context */ 93 if (!snes->setupcalled) PetscCall(SNESSetUp(snes)); 94 95 w = neP->workv[0]; 96 xp = neP->workv[1]; 97 fvec = neP->workv[2]; 98 fp = neP->fp; 99 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; 107 108 /* We have 5 tries to attempt to compute a good hopt value */ 109 PetscCall(SNESGetIterationNumber(snes, &i)); 110 PetscCall(PetscFPrintf(comm, fp, "\n ------- SNES iteration %" PetscInt_FMT " ---------\n", i)); 111 for (iter = 0; iter < 5; iter++) { 112 neP->h_first_try = h; 113 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 PetscCall(VecWAXPY(xp, alpha, p, x)); 119 PetscCall(SNESComputeFunction(snes, xp, fvec)); 120 neP->function_count++; 121 PetscCall(VecDot(fvec, w, &fval[k])); 122 } 123 f = fval[(nf + 1) / 2 - 1]; 124 125 /* Construct the difference table */ 126 for (i = 0; i < nf; i++) tab[i][0] = fval[i]; 127 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 } 131 132 /* Print the difference table */ 133 PetscCall(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++) PetscCall(PetscFPrintf(comm, fp, " %10.2e ", tab[i][j])); 136 PetscCall(PetscFPrintf(comm, fp, "\n")); 137 } 138 139 /* Call the noise estimator */ 140 PetscCall(SNESNoise_dnest_(&nf, fval, &h, fnoise, &fder2, hopt, &info, eps)); 141 142 /* Output statements */ 143 rerrf = *fnoise / PetscAbsScalar(f); 144 if (info == 1) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected")); 145 if (info == 2) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too small")); 146 if (info == 3) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise not detected; h is too large")); 147 if (info == 4) PetscCall(PetscFPrintf(comm, fp, "%s\n", "Noise detected, but unreliable hopt")); 148 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])); 149 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)); 150 151 /* Save fnoise and fder2. */ 152 if (*fnoise) fnoise_s = *fnoise; 153 if (fder2) fder2_s = fder2; 154 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: 174 175 if (*fnoise < neP->fnoise_min) { 176 PetscCall(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 PetscCall(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 } 185 186 PetscCall(PetscFPrintf(comm, fp, "Errors in derivative:\n")); 187 PetscCall(PetscFPrintf(comm, fp, "f = %g, fnoise = %g, fder2 = %g, hopt = %g\n", (double)f, (double)*fnoise, (double)fder2, (double)*hopt)); 188 189 /* For now, compute h **each** MV Mult!! */ 190 /* 191 PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg)); 192 if (!flg) { 193 Mat mat; 194 PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL)); 195 PetscCall(MatSNESMFMoreSetParameters(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt)); 196 } 197 */ 198 fcount = neP->function_count - fcount; 199 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)); 200 201 PetscCall(PetscOptionsGetBool(NULL, NULL, "-noise_test", &noise_test, NULL)); 202 if (noise_test) PetscCall(JacMatMultCompare(snes, x, p, *hopt)); 203 PetscFunctionReturn(0); 204 } 205 206 PetscErrorCode JacMatMultCompare(SNES snes, Vec x, Vec p, double hopt) { 207 Vec yy1, yy2; /* work vectors */ 208 PetscViewer view2; /* viewer */ 209 Mat J; /* analytic Jacobian (set as preconditioner matrix) */ 210 Mat Jmf; /* matrix-free Jacobian (set as true system matrix) */ 211 double h; /* differencing parameter */ 212 Vec f; 213 PetscScalar alpha; 214 PetscReal yy1n, yy2n, enorm; 215 PetscInt i; 216 PetscBool printv = PETSC_FALSE; 217 char filename[32]; 218 MPI_Comm comm; 219 220 PetscFunctionBegin; 221 PetscCall(PetscObjectGetComm((PetscObject)snes, &comm)); 222 /* Compute function and analytic Jacobian at x */ 223 PetscCall(SNESGetJacobian(snes, &Jmf, &J, NULL, NULL)); 224 PetscCall(SNESComputeJacobian(snes, x, Jmf, J)); 225 PetscCall(SNESGetFunction(snes, &f, NULL, NULL)); 226 PetscCall(SNESComputeFunction(snes, x, f)); 227 228 /* Duplicate work vectors */ 229 PetscCall(VecDuplicate(x, &yy2)); 230 PetscCall(VecDuplicate(x, &yy1)); 231 232 /* Compute true matrix-vector product */ 233 PetscCall(MatMult(J, p, yy1)); 234 PetscCall(VecNorm(yy1, NORM_2, &yy1n)); 235 236 /* View product vector if desired */ 237 PetscCall(PetscOptionsGetBool(NULL, NULL, "-print_vecs", &printv, NULL)); 238 if (printv) { 239 PetscCall(PetscViewerASCIIOpen(comm, "y1.out", &view2)); 240 PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON)); 241 PetscCall(VecView(yy1, view2)); 242 PetscCall(PetscViewerPopFormat(view2)); 243 PetscCall(PetscViewerDestroy(&view2)); 244 } 245 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 PetscCall(MatSNESMFMoreSetParameters(Jmf, PETSC_DEFAULT, PETSC_DEFAULT, h)); 252 253 /* Compute matrix-vector product via differencing approximation */ 254 PetscCall(MatMult(Jmf, p, yy2)); 255 PetscCall(VecNorm(yy2, NORM_2, &yy2n)); 256 257 /* View product vector if desired */ 258 if (printv) { 259 sprintf(filename, "y2.%d.out", (int)i); 260 PetscCall(PetscViewerASCIIOpen(comm, filename, &view2)); 261 PetscCall(PetscViewerPushFormat(view2, PETSC_VIEWER_ASCII_COMMON)); 262 PetscCall(VecView(yy2, view2)); 263 PetscCall(PetscViewerPopFormat(view2)); 264 PetscCall(PetscViewerDestroy(&view2)); 265 } 266 267 /* Compute relative error */ 268 PetscCall(VecAXPY(yy2, alpha, yy1)); 269 PetscCall(VecNorm(yy2, NORM_2, &enorm)); 270 enorm = enorm / yy1n; 271 PetscCall(PetscFPrintf(comm, stdout, "h = %g: relative error = %g\n", (double)h, (double)enorm)); 272 h *= 10.0; 273 } 274 PetscFunctionReturn(0); 275 } 276 277 static PetscInt lin_its_total = 0; 278 279 PetscErrorCode SNESNoiseMonitor(SNES snes, PetscInt its, double fnorm, void *dummy) { 280 PetscInt lin_its; 281 282 PetscFunctionBegin; 283 PetscCall(SNESGetLinearSolveIterations(snes, &lin_its)); 284 lin_its_total += lin_its; 285 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)); 286 287 PetscCall(SNESUnSetMatrixFreeParameter(snes)); 288 PetscFunctionReturn(0); 289 } 290