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 { 28 DIFFPAR_MORE *neP; 29 Vec w; 30 PetscRandom rctx; /* random number generator context */ 31 PetscBool flg; 32 char noise_file[PETSC_MAX_PATH_LEN]; 33 34 PetscFunctionBegin; 35 PetscCall(PetscNewLog(snes,&neP)); 36 37 neP->function_count = 0; 38 neP->fnoise_min = 1.0e-20; 39 neP->hopt_min = 1.0e-8; 40 neP->h_first_try = 1.0e-3; 41 neP->fnoise_resets = 0; 42 neP->hopt_resets = 0; 43 44 /* Create work vectors */ 45 PetscCall(VecDuplicateVecs(x,3,&neP->workv)); 46 w = neP->workv[0]; 47 48 /* Set components of vector w to random numbers */ 49 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)snes),&rctx)); 50 PetscCall(PetscRandomSetFromOptions(rctx)); 51 PetscCall(VecSetRandom(w,rctx)); 52 PetscCall(PetscRandomDestroy(&rctx)); 53 54 /* Open output file */ 55 PetscCall(PetscOptionsGetString(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_mf_noise_file",noise_file,sizeof(noise_file),&flg)); 56 if (flg) neP->fp = fopen(noise_file,"w"); 57 else neP->fp = fopen("noise.out","w"); 58 PetscCheck(neP->fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file"); 59 PetscCall(PetscInfo(snes,"Creating Jorge's differencing parameter context\n")); 60 61 *outneP = neP; 62 PetscFunctionReturn(0); 63 } 64 65 PetscErrorCode SNESDiffParameterDestroy_More(void *nePv) 66 { 67 DIFFPAR_MORE *neP = (DIFFPAR_MORE*)nePv; 68 int err; 69 70 PetscFunctionBegin; 71 /* Destroy work vectors and close output file */ 72 PetscCall(VecDestroyVecs(3,&neP->workv)); 73 err = fclose(neP->fp); 74 PetscCheck(!err,PETSC_COMM_SELF,PETSC_ERR_SYS,"fclose() failed on file"); 75 PetscCall(PetscFree(neP)); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode SNESDiffParameterCompute_More(SNES snes,void *nePv,Vec x,Vec p,double *fnoise,double *hopt) 80 { 81 DIFFPAR_MORE *neP = (DIFFPAR_MORE*)nePv; 82 Vec w, xp, fvec; /* work vectors to use in computing h */ 83 double zero = 0.0, hl, hu, h, fnoise_s, fder2_s; 84 PetscScalar alpha; 85 PetscScalar fval[7], tab[7][7], eps[7], f = -1; 86 double rerrf = -1., fder2; 87 PetscInt iter, k, i, j, info; 88 PetscInt nf = 7; /* number of function evaluations */ 89 PetscInt fcount; 90 MPI_Comm comm; 91 FILE *fp; 92 PetscBool noise_test = PETSC_FALSE; 93 94 PetscFunctionBegin; 95 PetscCall(PetscObjectGetComm((PetscObject)snes,&comm)); 96 /* Call to SNESSetUp() just to set data structures in SNES context */ 97 if (!snes->setupcalled) PetscCall(SNESSetUp(snes)); 98 99 w = neP->workv[0]; 100 xp = neP->workv[1]; 101 fvec = neP->workv[2]; 102 fp = neP->fp; 103 104 /* Initialize parameters */ 105 hl = zero; 106 hu = zero; 107 h = neP->h_first_try; 108 fnoise_s = zero; 109 fder2_s = zero; 110 fcount = neP->function_count; 111 112 /* We have 5 tries to attempt to compute a good hopt value */ 113 PetscCall(SNESGetIterationNumber(snes,&i)); 114 PetscCall(PetscFPrintf(comm,fp,"\n ------- SNES iteration %" PetscInt_FMT " ---------\n",i)); 115 for (iter=0; iter<5; iter++) { 116 neP->h_first_try = h; 117 118 /* Compute the nf function values needed to estimate the noise from 119 the difference table */ 120 for (k=0; k<nf; k++) { 121 alpha = h * (k+1 - (nf+1)/2); 122 PetscCall(VecWAXPY(xp,alpha,p,x)); 123 PetscCall(SNESComputeFunction(snes,xp,fvec)); 124 neP->function_count++; 125 PetscCall(VecDot(fvec,w,&fval[k])); 126 } 127 f = fval[(nf+1)/2 - 1]; 128 129 /* Construct the difference table */ 130 for (i=0; i<nf; i++) tab[i][0] = fval[i]; 131 132 for (j=0; j<nf-1; j++) { 133 for (i=0; i<nf-j-1; i++) { 134 tab[i][j+1] = tab[i+1][j] - tab[i][j]; 135 } 136 } 137 138 /* Print the difference table */ 139 PetscCall(PetscFPrintf(comm,fp,"Difference Table: iter = %" PetscInt_FMT "\n",iter)); 140 for (i=0; i<nf; i++) { 141 for (j=0; j<nf-i; j++) { 142 PetscCall(PetscFPrintf(comm,fp," %10.2e ",tab[i][j])); 143 } 144 PetscCall(PetscFPrintf(comm,fp,"\n")); 145 } 146 147 /* Call the noise estimator */ 148 PetscCall(SNESNoise_dnest_(&nf,fval,&h,fnoise,&fder2,hopt,&info,eps)); 149 150 /* Output statements */ 151 rerrf = *fnoise/PetscAbsScalar(f); 152 if (info == 1) PetscCall(PetscFPrintf(comm,fp,"%s\n","Noise detected")); 153 if (info == 2) {PetscCall(PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too small"));} 154 if (info == 3) {PetscCall(PetscFPrintf(comm,fp,"%s\n","Noise not detected; h is too large"));} 155 if (info == 4) PetscCall(PetscFPrintf(comm,fp,"%s\n","Noise detected, but unreliable hopt")); 156 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])); 157 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)); 158 159 /* Save fnoise and fder2. */ 160 if (*fnoise) fnoise_s = *fnoise; 161 if (fder2) fder2_s = fder2; 162 163 /* Check for noise detection. */ 164 if (fnoise_s && fder2_s) { 165 *fnoise = fnoise_s; 166 fder2 = fder2_s; 167 *hopt = 1.68*sqrt(*fnoise/PetscAbsScalar(fder2)); 168 goto theend; 169 } else { 170 171 /* Update hl and hu, and determine new h */ 172 if (info == 2 || info == 4) { 173 hl = h; 174 if (hu == zero) h = 100*h; 175 else h = PetscMin(100*h,0.1*hu); 176 } else if (info == 3) { 177 hu = h; 178 h = PetscMax(1.0e-3,sqrt(hl/hu))*hu; 179 } 180 } 181 } 182 theend: 183 184 if (*fnoise < neP->fnoise_min) { 185 PetscCall(PetscFPrintf(comm,fp,"Resetting fnoise: fnoise1 = %g, fnoise_min = %g\n",(double)*fnoise,(double)neP->fnoise_min)); 186 *fnoise = neP->fnoise_min; 187 neP->fnoise_resets++; 188 } 189 if (*hopt < neP->hopt_min) { 190 PetscCall(PetscFPrintf(comm,fp,"Resetting hopt: hopt1 = %g, hopt_min = %g\n",(double)*hopt,(double)neP->hopt_min)); 191 *hopt = neP->hopt_min; 192 neP->hopt_resets++; 193 } 194 195 PetscCall(PetscFPrintf(comm,fp,"Errors in derivative:\n")); 196 PetscCall(PetscFPrintf(comm,fp,"f = %g, fnoise = %g, fder2 = %g, hopt = %g\n",(double)f,(double)*fnoise,(double)fder2,(double)*hopt)); 197 198 /* For now, compute h **each** MV Mult!! */ 199 /* 200 PetscCall(PetscOptionsHasName(NULL,"-matrix_free_jorge_each_mvp",&flg)); 201 if (!flg) { 202 Mat mat; 203 PetscCall(SNESGetJacobian(snes,&mat,NULL,NULL)); 204 PetscCall(SNESDefaultMatrixFreeSetParameters2(mat,PETSC_DEFAULT,PETSC_DEFAULT,*hopt)); 205 } 206 */ 207 fcount = neP->function_count - fcount; 208 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)); 209 210 PetscCall(PetscOptionsGetBool(NULL,NULL,"-noise_test",&noise_test,NULL)); 211 if (noise_test) { 212 PetscCall(JacMatMultCompare(snes,x,p,*hopt)); 213 } 214 PetscFunctionReturn(0); 215 } 216 217 PetscErrorCode JacMatMultCompare(SNES snes,Vec x,Vec p,double hopt) 218 { 219 Vec yy1, yy2; /* work vectors */ 220 PetscViewer view2; /* viewer */ 221 Mat J; /* analytic Jacobian (set as preconditioner matrix) */ 222 Mat Jmf; /* matrix-free Jacobian (set as true system matrix) */ 223 double h; /* differencing parameter */ 224 Vec f; 225 PetscScalar alpha; 226 PetscReal yy1n,yy2n,enorm; 227 PetscInt i; 228 PetscBool printv = PETSC_FALSE; 229 char filename[32]; 230 MPI_Comm comm; 231 232 PetscFunctionBegin; 233 PetscCall(PetscObjectGetComm((PetscObject)snes,&comm)); 234 /* Compute function and analytic Jacobian at x */ 235 PetscCall(SNESGetJacobian(snes,&Jmf,&J,NULL,NULL)); 236 PetscCall(SNESComputeJacobian(snes,x,Jmf,J)); 237 PetscCall(SNESGetFunction(snes,&f,NULL,NULL)); 238 PetscCall(SNESComputeFunction(snes,x,f)); 239 240 /* Duplicate work vectors */ 241 PetscCall(VecDuplicate(x,&yy2)); 242 PetscCall(VecDuplicate(x,&yy1)); 243 244 /* Compute true matrix-vector product */ 245 PetscCall(MatMult(J,p,yy1)); 246 PetscCall(VecNorm(yy1,NORM_2,&yy1n)); 247 248 /* View product vector if desired */ 249 PetscCall(PetscOptionsGetBool(NULL,NULL,"-print_vecs",&printv,NULL)); 250 if (printv) { 251 PetscCall(PetscViewerASCIIOpen(comm,"y1.out",&view2)); 252 PetscCall(PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON)); 253 PetscCall(VecView(yy1,view2)); 254 PetscCall(PetscViewerPopFormat(view2)); 255 PetscCall(PetscViewerDestroy(&view2)); 256 } 257 258 /* Test Jacobian-vector product computation */ 259 alpha = -1.0; 260 h = 0.01 * hopt; 261 for (i=0; i<5; i++) { 262 /* Set differencing parameter for matrix-free multiplication */ 263 PetscCall(SNESDefaultMatrixFreeSetParameters2(Jmf,PETSC_DEFAULT,PETSC_DEFAULT,h)); 264 265 /* Compute matrix-vector product via differencing approximation */ 266 PetscCall(MatMult(Jmf,p,yy2)); 267 PetscCall(VecNorm(yy2,NORM_2,&yy2n)); 268 269 /* View product vector if desired */ 270 if (printv) { 271 sprintf(filename,"y2.%d.out",(int)i); 272 PetscCall(PetscViewerASCIIOpen(comm,filename,&view2)); 273 PetscCall(PetscViewerPushFormat(view2,PETSC_VIEWER_ASCII_COMMON)); 274 PetscCall(VecView(yy2,view2)); 275 PetscCall(PetscViewerPopFormat(view2)); 276 PetscCall(PetscViewerDestroy(&view2)); 277 } 278 279 /* Compute relative error */ 280 PetscCall(VecAXPY(yy2,alpha,yy1)); 281 PetscCall(VecNorm(yy2,NORM_2,&enorm)); 282 enorm = enorm/yy1n; 283 PetscCall(PetscFPrintf(comm,stdout,"h = %g: relative error = %g\n",(double)h,(double)enorm)); 284 h *= 10.0; 285 } 286 PetscFunctionReturn(0); 287 } 288 289 static PetscInt lin_its_total = 0; 290 291 PetscErrorCode SNESNoiseMonitor(SNES snes,PetscInt its,double fnorm,void *dummy) 292 { 293 PetscInt lin_its; 294 295 PetscFunctionBegin; 296 PetscCall(SNESGetLinearSolveIterations(snes,&lin_its)); 297 lin_its_total += lin_its; 298 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)); 299 300 PetscCall(SNESUnSetMatrixFreeParameter(snes)); 301 PetscFunctionReturn(0); 302 } 303