Actual source code: minsurf2.c
1: /*$Id$*/
3: /* Program usage: mpirun -np <proc> minsurf2 [-help] [all TAO options] */
5: /*
6: Include "tao.h" so we can use TAO solvers.
7: petscda.h for distributed array
8: */
9: #include "petscda.h"
10: #include "tao.h"
12: static char help[] =
13: "This example demonstrates use of the TAO package to \n\
14: solve an unconstrained minimization problem. This example is based on a \n\
15: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
16: boundary values along the edges of the domain, the objective is to find the\n\
17: surface with the minimal area that satisfies the boundary conditions.\n\
18: The command line options are:\n\
19: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
20: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
21: -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
22: for an average of the boundary conditions\n\n";
24: /*T
25: Concepts: TAO - Solving an unconstrained minimization problem
26: Routines: TaoInitialize(); TaoFinalize();
27: Routines: TaoCreate(); TaoDestroy();
28: Routines: TaoApplicationCreate(); TaoAppDestroy();
29: Routines: TaoAppSetInitialSolutionVec();
30: Routines: TaoAppSetObjectiveAndGradientRoutine();
31: Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
32: Routines: TaoSetOptions();
33: Routines: TaoAppGetKSP(); TaoSolveApplication();
34: Routines: TaoAppSetMonitor(); TaoView();
35: Routines: TaoAppGetSolutionVec();
36: Processors: 1
37: T*/
39: /*
40: User-defined application context - contains data needed by the
41: application-provided call-back routines, FormFunctionGradient()
42: and FormHessian().
43: */
44: typedef struct {
45: PetscInt mx, my; /* discretization in x, y directions */
46: double *bottom, *top, *left, *right; /* boundary values */
47: DA da; /* distributed array data structure */
48: Mat H; /* Hessian */
49: ISColoring iscoloring;
50: } AppCtx;
53: /* -------- User-defined Routines --------- */
55: static int MSA_BoundaryConditions(AppCtx*);
56: static int MSA_InitialPoint(AppCtx*,Vec);
57: int QuadraticH(AppCtx*,Vec,Mat);
58: int FormFunctionGradient(TAO_APPLICATION,Vec,double *,Vec,void*);
59: int FormGradient(TAO_APPLICATION,Vec,Vec,void*);
60: int FormHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure *,void*);
61: int My_Monitor(TAO_APPLICATION, void *);
65: int main( int argc, char **argv )
66: {
67: int info; /* used to check for functions returning nonzeros */
68: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
69: PetscInt iter; /* iteration information */
70: double ff,gnorm;
71: Vec x; /* solution, gradient vectors */
72: PetscTruth flg, viewmat; /* flags */
73: PetscTruth fddefault, fdcoloring; /* flags */
74: KSP ksp; /* Krylov subspace method */
75: TaoMethod method = "tao_nls"; /* minimization method */
76: TaoTerminateReason reason;
77: TAO_SOLVER tao; /* TAO_SOLVER solver context */
78: TAO_APPLICATION minsurfapp; /* The PETSc application */
79: AppCtx user; /* user-defined work context */
81: /* Initialize TAO */
82: PetscInitialize( &argc, &argv,(char *)0,help );
83: TaoInitialize( &argc, &argv,(char *)0,help );
85: /* Specify dimension of the problem */
86: user.mx = 10; user.my = 10;
88: /* Check for any command line arguments that override defaults */
89: info = PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); CHKERRQ(info);
90: info = PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); CHKERRQ(info);
92: PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
93: PetscPrintf(MPI_COMM_WORLD,"mx: %d my: %d \n\n",user.mx,user.my);
96: /* Let PETSc determine the vector distribution */
97: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
99: /* Create distributed array (DA) to manage parallel grid and vectors */
100: info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,user.mx,
101: user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info);
102:
104: /* Create TAO solver and set desired solution method. Create an TAO application structure */
105: info = TaoCreate(PETSC_COMM_WORLD,method,&tao); CHKERRQ(info);
106: info = TaoApplicationCreate(PETSC_COMM_WORLD,&minsurfapp); CHKERRQ(info);
108: /*
109: Extract global vector from DA for the vector of variables -- PETSC routine
110: Compute the initial solution -- application specific, see below
111: Set this vector for use by TAO -- TAO routine
112: */
113: info = DACreateGlobalVector(user.da,&x); CHKERRQ(info);
114: info = MSA_BoundaryConditions(&user); CHKERRQ(info);
115: info = MSA_InitialPoint(&user,x); CHKERRQ(info);
116: info = TaoAppSetInitialSolutionVec(minsurfapp,x); CHKERRQ(info);
118: /*
119: Initialize the Application context for use in function evaluations -- application specific, see below.
120: Set routines for function and gradient evaluation
121: */
122: info = TaoAppSetObjectiveAndGradientRoutine(minsurfapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);
124: /*
125: Given the command line arguments, calculate the hessian with either the user-
126: provided function FormHessian, or the default finite-difference driven Hessian
127: functions
128: */
129: info = PetscOptionsHasName(PETSC_NULL,"-tao_fddefault",&fddefault);CHKERRQ(info);
130: info = PetscOptionsHasName(PETSC_NULL,"-tao_fdcoloring",&fdcoloring);CHKERRQ(info);
133: /*
134: Create a matrix data structure to store the Hessian and set
135: the Hessian evalution routine.
136: Set the matrix structure to be used for Hessian evalutions
137: */
138: info = DAGetMatrix(user.da,MATAIJ,&user.H);CHKERRQ(info);
139: info = MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(info);
141: info = TaoAppSetHessianMat(minsurfapp,user.H,user.H); CHKERRQ(info);
143: if (fdcoloring) {
144: info = DAGetColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&(user.iscoloring));
145: CHKERRQ(info);
146: info = TaoAppSetColoring(minsurfapp, user.iscoloring); CHKERRQ(info);
147: info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessianColor,(void *)&user); CHKERRQ(info);
148: } else if (fddefault){
149: info = TaoAppSetHessianRoutine(minsurfapp,TaoAppDefaultComputeHessian,(void *)&user); CHKERRQ(info);
150: } else {
151: info = TaoAppSetHessianRoutine(minsurfapp,FormHessian,(void *)&user); CHKERRQ(info);
152: }
155: /*
156: If my_monitor option is in command line, then use the user-provided
157: monitoring function
158: */
159: info = PetscOptionsHasName(PETSC_NULL,"-my_monitor",&viewmat); CHKERRQ(info);
160: if (viewmat){
161: info = TaoAppSetMonitor(minsurfapp,My_Monitor,TAO_NULL); CHKERRQ(info);
162: }
164: /* Check for any tao command line options */
165: info = TaoSetOptions(minsurfapp,tao); CHKERRQ(info);
167: /* Limit the number of iterations in the KSP linear solver */
168: info = TaoAppGetKSP(minsurfapp,&ksp); CHKERRQ(info);
169: if (ksp) { /* Modify the PETSc KSP structure */
170: info = KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,user.mx*user.my);
171: CHKERRQ(info);
172: }
174: /* SOLVE THE APPLICATION */
175: info = TaoSolveApplication(minsurfapp,tao); CHKERRQ(info);
177: /* Get information on termination */
178: info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
179: if (reason <= 0 ){
180: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
181: PetscPrintf(MPI_COMM_WORLD," Iterations: %d, Function Value: %4.2e, Residual: %4.2e \n",iter,ff,gnorm);
182: }
184: /*
185: To View TAO solver information use
186: info = TaoView(tao); CHKERRQ(info);
187: */
189: /* Free TAO data structures */
190: info = TaoDestroy(tao); CHKERRQ(info);
191: info = TaoAppDestroy(minsurfapp); CHKERRQ(info);
193: /* Free PETSc data structures */
194: info = VecDestroy(x); CHKERRQ(info);
195: info = MatDestroy(user.H); CHKERRQ(info);
196: if (fdcoloring) {
197: info = ISColoringDestroy(user.iscoloring);CHKERRQ(info);
198: }
199: info = PetscFree(user.bottom); CHKERRQ(info);
200: info = PetscFree(user.top); CHKERRQ(info);
201: info = PetscFree(user.left); CHKERRQ(info);
202: info = PetscFree(user.right); CHKERRQ(info);
203: info = DADestroy(user.da); CHKERRQ(info);
205: /* Finalize TAO */
206: TaoFinalize();
207: PetscFinalize();
208:
209: return 0;
210: }
214: int FormGradient(TAO_APPLICATION taoapp, Vec X, Vec G,void *userCtx){
215: int info;
216: double fcn;
217: TaoFunctionBegin;
218: info = FormFunctionGradient(taoapp,X,&fcn,G,userCtx);CHKERRQ(info);
219: TaoFunctionReturn(0);
220: }
222: /* -------------------------------------------------------------------- */
225: /* FormFunctionGradient - Evaluates the function and corresponding gradient.
227: Input Parameters:
228: . taoapp - the TAO_APPLICATION context
229: . XX - input vector
230: . userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
231:
232: Output Parameters:
233: . fcn - the newly evaluated function
234: . GG - vector containing the newly evaluated gradient
235: */
236: int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *userCtx){
238: AppCtx * user = (AppCtx *) userCtx;
239: int info;
240: PetscInt i,j;
241: PetscInt mx=user->mx, my=user->my;
242: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
243: double ft=0;
244: double hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
245: double rhx=mx+1, rhy=my+1;
246: double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
247: double df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
248: PetscScalar **g, **x;
249: Vec localX;
251: /* Get local mesh boundaries */
252: info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);
254: info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
255: info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);
257: /* Scatter ghost points to local vector */
258: info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
259: info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
261: /* Get pointers to vector data */
262: info = DAVecGetArray(user->da,localX,(void**)&x);
263: info = DAVecGetArray(user->da,G,(void**)&g);
265: /* Compute function and gradient over the locally owned part of the mesh */
266: for (j=ys; j<ys+ym; j++){
267: for (i=xs; i< xs+xm; i++){
268:
269: xc = x[j][i];
270: xlt=xrb=xl=xr=xb=xt=xc;
271:
272: if (i==0){ /* left side */
273: xl= user->left[j-ys+1];
274: xlt = user->left[j-ys+2];
275: } else {
276: xl = x[j][i-1];
277: }
279: if (j==0){ /* bottom side */
280: xb=user->bottom[i-xs+1];
281: xrb =user->bottom[i-xs+2];
282: } else {
283: xb = x[j-1][i];
284: }
285:
286: if (i+1 == gxs+gxm){ /* right side */
287: xr=user->right[j-ys+1];
288: xrb = user->right[j-ys];
289: } else {
290: xr = x[j][i+1];
291: }
293: if (j+1==gys+gym){ /* top side */
294: xt=user->top[i-xs+1];
295: xlt = user->top[i-xs];
296: }else {
297: xt = x[j+1][i];
298: }
300: if (i>gxs && j+1<gys+gym){
301: xlt = x[j+1][i-1];
302: }
303: if (j>gys && i+1<gxs+gxm){
304: xrb = x[j-1][i+1];
305: }
307: d1 = (xc-xl);
308: d2 = (xc-xr);
309: d3 = (xc-xt);
310: d4 = (xc-xb);
311: d5 = (xr-xrb);
312: d6 = (xrb-xb);
313: d7 = (xlt-xl);
314: d8 = (xt-xlt);
315:
316: df1dxc = d1*hydhx;
317: df2dxc = ( d1*hydhx + d4*hxdhy );
318: df3dxc = d3*hxdhy;
319: df4dxc = ( d2*hydhx + d3*hxdhy );
320: df5dxc = d2*hydhx;
321: df6dxc = d4*hxdhy;
323: d1 *= rhx;
324: d2 *= rhx;
325: d3 *= rhy;
326: d4 *= rhy;
327: d5 *= rhy;
328: d6 *= rhx;
329: d7 *= rhy;
330: d8 *= rhx;
332: f1 = sqrt( 1.0 + d1*d1 + d7*d7);
333: f2 = sqrt( 1.0 + d1*d1 + d4*d4);
334: f3 = sqrt( 1.0 + d3*d3 + d8*d8);
335: f4 = sqrt( 1.0 + d3*d3 + d2*d2);
336: f5 = sqrt( 1.0 + d2*d2 + d5*d5);
337: f6 = sqrt( 1.0 + d4*d4 + d6*d6);
338:
339: f2 = sqrt( 1.0 + d1*d1 + d4*d4);
340: f4 = sqrt( 1.0 + d3*d3 + d2*d2);
342: ft = ft + (f2 + f4);
344: df1dxc /= f1;
345: df2dxc /= f2;
346: df3dxc /= f3;
347: df4dxc /= f4;
348: df5dxc /= f5;
349: df6dxc /= f6;
351: g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
352:
353: }
354: }
356: /* Compute triangular areas along the border of the domain. */
357: if (xs==0){ /* left side */
358: for (j=ys; j<ys+ym; j++){
359: d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
360: d2=(user->left[j-ys+1] - x[j][0]) *rhx;
361: ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
362: }
363: }
364: if (ys==0){ /* bottom side */
365: for (i=xs; i<xs+xm; i++){
366: d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
367: d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
368: ft = ft+sqrt( 1.0 + d3*d3 + d2*d2);
369: }
370: }
372: if (xs+xm==mx){ /* right side */
373: for (j=ys; j< ys+ym; j++){
374: d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
375: d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
376: ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
377: }
378: }
379: if (ys+ym==my){ /* top side */
380: for (i=xs; i<xs+xm; i++){
381: d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
382: d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
383: ft = ft+sqrt( 1.0 + d1*d1 + d4*d4);
384: }
385: }
387: if (ys==0 && xs==0){
388: d1=(user->left[0]-user->left[1])/hy;
389: d2=(user->bottom[0]-user->bottom[1])*rhx;
390: ft +=sqrt( 1.0 + d1*d1 + d2*d2);
391: }
392: if (ys+ym == my && xs+xm == mx){
393: d1=(user->right[ym+1] - user->right[ym])*rhy;
394: d2=(user->top[xm+1] - user->top[xm])*rhx;
395: ft +=sqrt( 1.0 + d1*d1 + d2*d2);
396: }
398: ft=ft*area;
399: info = MPI_Allreduce(&ft,fcn,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);CHKERRQ(info);
401: /* Restore vectors */
402: info = DAVecRestoreArray(user->da,localX,(void**)&x);
403: info = DAVecRestoreArray(user->da,G,(void**)&g);
405: /* Scatter values to global vector */
406: info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);
408: info = PetscLogFlops(67*xm*ym); CHKERRQ(info);
410: return 0;
411: }
413: /* ------------------------------------------------------------------- */
416: /*
417: FormHessian - Evaluates Hessian matrix.
419: Input Parameters:
420: . taoapp - the TAO_APPLICATION context
421: . x - input vector
422: . ptr - optional user-defined context, as set by TaoSetHessian()
424: Output Parameters:
425: . H - Hessian matrix
426: . Hpre - optionally different preconditioning matrix
427: . flg - flag indicating matrix structure
429: */
430: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
431: {
432: int info;
433: AppCtx *user = (AppCtx *) ptr;
435: /* Evaluate the Hessian entries*/
436: info = QuadraticH(user,X,*H); CHKERRQ(info);
439: /* Indicate that this matrix has the same sparsity pattern during
440: successive iterations; setting this flag can save significant work
441: in computing the preconditioner for some methods. */
442: *flg=SAME_NONZERO_PATTERN;
444: return 0;
445: }
447: /* ------------------------------------------------------------------- */
450: /*
451: QuadraticH - Evaluates Hessian matrix.
453: Input Parameters:
454: . user - user-defined context, as set by TaoSetHessian()
455: . X - input vector
457: Output Parameter:
458: . H - Hessian matrix
459: */
460: int QuadraticH(AppCtx *user, Vec X, Mat Hessian)
461: {
462: int info;
463: PetscInt i,j,k;
464: PetscInt mx=user->mx, my=user->my;
465: PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
466: double hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
467: double f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
468: double hl,hr,ht,hb,hc,htl,hbr;
469: PetscScalar **x, v[7];
470: MatStencil col[7],row;
471: Vec localX;
472: PetscTruth assembled;
474: /* Get local mesh boundaries */
475: info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);
477: info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
478: info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);
480: /* Scatter ghost points to local vector */
481: info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
482: info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
484: /* Get pointers to vector data */
485: info = DAVecGetArray(user->da,localX,(void**)&x);
487: /* Initialize matrix entries to zero */
488: info = MatAssembled(Hessian,&assembled); CHKERRQ(info);
489: if (assembled){info = MatZeroEntries(Hessian); CHKERRQ(info);}
492: /* Set various matrix options */
493: info = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE); CHKERRQ(info);
495: /* Compute Hessian over the locally owned part of the mesh */
497: for (j=ys; j<ys+ym; j++){
498:
499: for (i=xs; i< xs+xm; i++){
501: xc = x[j][i];
502: xlt=xrb=xl=xr=xb=xt=xc;
504: /* Left side */
505: if (i==0){
506: xl = user->left[j-ys+1];
507: xlt = user->left[j-ys+2];
508: } else {
509: xl = x[j][i-1];
510: }
511:
512: if (j==0){
513: xb = user->bottom[i-xs+1];
514: xrb = user->bottom[i-xs+2];
515: } else {
516: xb = x[j-1][i];
517: }
518:
519: if (i+1 == mx){
520: xr = user->right[j-ys+1];
521: xrb = user->right[j-ys];
522: } else {
523: xr = x[j][i+1];
524: }
526: if (j+1==my){
527: xt = user->top[i-xs+1];
528: xlt = user->top[i-xs];
529: }else {
530: xt = x[j+1][i];
531: }
533: if (i>0 && j+1<my){
534: xlt = x[j+1][i-1];
535: }
536: if (j>0 && i+1<mx){
537: xrb = x[j-1][i+1];
538: }
541: d1 = (xc-xl)/hx;
542: d2 = (xc-xr)/hx;
543: d3 = (xc-xt)/hy;
544: d4 = (xc-xb)/hy;
545: d5 = (xrb-xr)/hy;
546: d6 = (xrb-xb)/hx;
547: d7 = (xlt-xl)/hy;
548: d8 = (xlt-xt)/hx;
549:
550: f1 = sqrt( 1.0 + d1*d1 + d7*d7);
551: f2 = sqrt( 1.0 + d1*d1 + d4*d4);
552: f3 = sqrt( 1.0 + d3*d3 + d8*d8);
553: f4 = sqrt( 1.0 + d3*d3 + d2*d2);
554: f5 = sqrt( 1.0 + d2*d2 + d5*d5);
555: f6 = sqrt( 1.0 + d4*d4 + d6*d6);
558: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
559: (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
560: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
561: (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
562: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
563: (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
564: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
565: (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
567: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
568: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
570: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
571: hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
572: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
573: (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
575: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
577: row.j = j; row.i = i;
578: k=0;
579: if (j>0){
580: v[k]=hb;
581: col[k].j = j - 1; col[k].i = i;
582: k++;
583: }
584:
585: if (j>0 && i < mx -1){
586: v[k]=hbr;
587: col[k].j = j - 1; col[k].i = i+1;
588: k++;
589: }
590:
591: if (i>0){
592: v[k]= hl;
593: col[k].j = j; col[k].i = i-1;
594: k++;
595: }
596:
597: v[k]= hc;
598: col[k].j = j; col[k].i = i;
599: k++;
600:
601: if (i < mx-1 ){
602: v[k]= hr;
603: col[k].j = j; col[k].i = i+1;
604: k++;
605: }
606:
607: if (i>0 && j < my-1 ){
608: v[k]= htl;
609: col[k].j = j+1; col[k].i = i-1;
610: k++;
611: }
612:
613: if (j < my-1 ){
614: v[k]= ht;
615: col[k].j = j+1; col[k].i = i;
616: k++;
617: }
618:
619: /*
620: Set matrix values using local numbering, which was defined
621: earlier, in the main routine.
622: */
623: info = MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
624: CHKERRQ(info);
625:
626: }
627: }
628:
629: /* Restore vectors */
630: info = DAVecRestoreArray(user->da,localX,(void**)&x);
632: info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);
634: /* Assemble the matrix */
635: info = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
636: info = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
638: info = PetscLogFlops(199*xm*ym); CHKERRQ(info);
639: return 0;
640: }
642: /* ------------------------------------------------------------------- */
645: /*
646: MSA_BoundaryConditions - Calculates the boundary conditions for
647: the region.
649: Input Parameter:
650: . user - user-defined application context
652: Output Parameter:
653: . user - user-defined application context
654: */
655: static int MSA_BoundaryConditions(AppCtx * user)
656: {
657: int info;
658: PetscInt i,j,k,limit=0,maxits=5;
659: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
660: PetscInt mx=user->mx,my=user->my;
661: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
662: double one=1.0, two=2.0, three=3.0, tol=1e-10;
663: double fnorm,det,hx,hy,xt=0,yt=0;
664: double u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
665: double b=-0.5, t=0.5, l=-0.5, r=0.5;
666: double *boundary;
667: PetscTruth flg;
669: /* Get local mesh boundaries */
670: info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
671: info = DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); CHKERRQ(info);
673: bsize=xm+2;
674: lsize=ym+2;
675: rsize=ym+2;
676: tsize=xm+2;
678: info = PetscMalloc(bsize*sizeof(double),&user->bottom); CHKERRQ(info);
679: info = PetscMalloc(tsize*sizeof(double),&user->top); CHKERRQ(info);
680: info = PetscMalloc(lsize*sizeof(double),&user->left); CHKERRQ(info);
681: info = PetscMalloc(rsize*sizeof(double),&user->right); CHKERRQ(info);
683: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
685: for (j=0; j<4; j++){
686: if (j==0){
687: yt=b;
688: xt=l+hx*xs;
689: limit=bsize;
690: boundary=user->bottom;
691: } else if (j==1){
692: yt=t;
693: xt=l+hx*xs;
694: limit=tsize;
695: boundary=user->top;
696: } else if (j==2){
697: yt=b+hy*ys;
698: xt=l;
699: limit=lsize;
700: boundary=user->left;
701: } else { //if (j==3)
702: yt=b+hy*ys;
703: xt=r;
704: limit=rsize;
705: boundary=user->right;
706: }
708: for (i=0; i<limit; i++){
709: u1=xt;
710: u2=-yt;
711: for (k=0; k<maxits; k++){
712: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
713: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
714: fnorm=sqrt(nf1*nf1+nf2*nf2);
715: if (fnorm <= tol) break;
716: njac11=one+u2*u2-u1*u1;
717: njac12=two*u1*u2;
718: njac21=-two*u1*u2;
719: njac22=-one - u1*u1 + u2*u2;
720: det = njac11*njac22-njac21*njac12;
721: u1 = u1-(njac22*nf1-njac12*nf2)/det;
722: u2 = u2-(njac11*nf2-njac21*nf1)/det;
723: }
725: boundary[i]=u1*u1-u2*u2;
726: if (j==0 || j==1) {
727: xt=xt+hx;
728: } else { // if (j==2 || j==3)
729: yt=yt+hy;
730: }
731:
732: }
734: }
736: /* Scale the boundary if desired */
737: if (1==1){
738: PetscReal scl = 1.0;
740: info = PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg);
741: CHKERRQ(info);
742: if (flg){
743: for (i=0;i<bsize;i++) user->bottom[i]*=scl;
744: }
746: info = PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg);
747: CHKERRQ(info);
748: if (flg){
749: for (i=0;i<tsize;i++) user->top[i]*=scl;
750: }
752: info = PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg);
753: CHKERRQ(info);
754: if (flg){
755: for (i=0;i<rsize;i++) user->right[i]*=scl;
756: }
758: info = PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg);
759: CHKERRQ(info);
760: if (flg){
761: for (i=0;i<lsize;i++) user->left[i]*=scl;
762: }
763: }
764:
765: return 0;
766: }
768: /* ------------------------------------------------------------------- */
771: /*
772: MSA_InitialPoint - Calculates the initial guess in one of three ways.
774: Input Parameters:
775: . user - user-defined application context
776: . X - vector for initial guess
778: Output Parameters:
779: . X - newly computed initial guess
780: */
781: static int MSA_InitialPoint(AppCtx * user, Vec X)
782: {
783: int info;
784: PetscInt start2=-1,i,j;
785: PetscReal start1=0;
786: PetscTruth flg1,flg2;
788: info = PetscOptionsGetReal(PETSC_NULL,"-start",&start1,&flg1); CHKERRQ(info);
789: info = PetscOptionsGetInt(PETSC_NULL,"-random",&start2,&flg2); CHKERRQ(info);
791: if (flg1){ /* The zero vector is reasonable */
792:
793: info = VecSet(X, start1); CHKERRQ(info);
795: } else if (flg2 && start2>0){ /* Try a random start between -0.5 and 0.5 */
797: PetscRandom rctx; PetscScalar np5=-0.5;
799: info = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
800: CHKERRQ(info);
801: for (i=0; i<start2; i++){
802: info = VecSetRandom(X, rctx); CHKERRQ(info);
803: }
804: info = PetscRandomDestroy(rctx); CHKERRQ(info);
805: info = VecShift(X, np5); CHKERRQ(info);
807: } else { /* Take an average of the boundary conditions */
809: PetscInt xs,xm,ys,ym;
810: PetscInt mx=user->mx,my=user->my;
811: PetscScalar **x;
812:
813: /* Get local mesh boundaries */
814: info = DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); CHKERRQ(info);
815:
816: /* Get pointers to vector data */
817: info = DAVecGetArray(user->da,X,(void**)&x);
819: /* Perform local computations */
820: for (j=ys; j<ys+ym; j++){
821: for (i=xs; i< xs+xm; i++){
822: x[j][i] = ( ((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+
823: ((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
824: }
825: }
826:
827: /* Restore vectors */
828: info = DAVecRestoreArray(user->da,X,(void**)&x); CHKERRQ(info);
830: info = PetscLogFlops(9*xm*ym); CHKERRQ(info);
831:
832: }
833: return 0;
834: }
836: /*-----------------------------------------------------------------------*/
839: int My_Monitor(TAO_APPLICATION minsurfapp, void *ctx){
840: int info;
841: Vec X;
843: info = TaoAppGetSolutionVec(minsurfapp,&X); CHKERRQ(info);
844: info = VecView(X,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(info);
845: return 0;
846: }