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