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