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