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