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