Actual source code: plate2.c

```  1: #include <petscdmda.h>
2: #include <petsctao.h>

4: static  char help[] =
5: "This example demonstrates use of the TAO package to \n\
6: solve an unconstrained minimization problem.  This example is based on a \n\
7: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain, \n\
8: boundary values along the edges of the domain, and a plate represented by \n\
9: lower boundary conditions, the objective is to find the\n\
10: surface with the minimal area that satisfies the boundary conditions.\n\
11: The command line options are:\n\
12:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
14:   -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
15:   -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
16:   -bheight <ht>, where <ht> = height of the plate\n\
17:   -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
18:                for an average of the boundary conditions\n\n";

20: /*T
21:    Concepts: TAO^Solving a bound constrained minimization problem
22:    Routines: TaoCreate();
23:    Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
24:    Routines: TaoSetHessianRoutine();
25:    Routines: TaoSetInitialVector();
26:    Routines: TaoSetVariableBounds();
27:    Routines: TaoSetFromOptions();
28:    Routines: TaoSolve(); TaoView();
29:    Routines: TaoDestroy();
30:    Processors: n
31: T*/

33: /*
34:    User-defined application context - contains data needed by the
35:    application-provided call-back routines, FormFunctionGradient(),
36:    FormHessian().
37: */
38: typedef struct {
39:   /* problem parameters */
40:   PetscReal      bheight;                  /* Height of plate under the surface */
41:   PetscInt       mx, my;                   /* discretization in x, y directions */
42:   PetscInt       bmx,bmy;                  /* Size of plate under the surface */
43:   Vec            Bottom, Top, Left, Right; /* boundary values */

45:   /* Working space */
46:   Vec         localX, localV;           /* ghosted local vector */
47:   DM          dm;                       /* distributed array data structure */
48:   Mat         H;
49: } AppCtx;

51: /* -------- User-defined Routines --------- */

53: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
54: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
55: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
56: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
57: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);

59: /* For testing matrix free submatrices */
60: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
61: PetscErrorCode MyMatMult(Mat,Vec,Vec);

63: int main(int argc, char **argv)
64: {
65:   PetscErrorCode         ierr;                 /* used to check for functions returning nonzeros */
66:   PetscInt               Nx, Ny;               /* number of processors in x- and y- directions */
67:   PetscInt               m, N;                 /* number of local and global elements in vectors */
68:   Vec                    x,xl,xu;               /* solution vector  and bounds*/
69:   PetscBool              flg;                /* A return variable when checking for user options */
70:   Tao                    tao;                  /* Tao solver context */
71:   ISLocalToGlobalMapping isltog;   /* local-to-global mapping object */
72:   Mat                    H_shell;                  /* to test matrix-free submatrices */
73:   AppCtx                 user;                 /* user-defined work context */

75:   /* Initialize PETSc, TAO */
76:   PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;

78:   /* Specify default dimension of the problem */
79:   user.mx = 10; user.my = 10; user.bheight=0.1;

81:   /* Check for any command line arguments that override defaults */
82:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
83:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
84:   PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);

86:   user.bmx = user.mx/2; user.bmy = user.my/2;
87:   PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);
88:   PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);

90:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
91:   PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%g\n",user.mx,user.my,user.bmx,user.bmy,(double)user.bheight);

93:   /* Calculate any derived values from parameters */
94:   N    = user.mx*user.my;

96:   /* Let Petsc determine the dimensions of the local vectors */
97:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

99:   /*
100:      A two dimensional distributed array will help define this problem,
101:      which derives from an elliptic PDE on two dimensional domain.  From
102:      the distributed array, Create the vectors.
103:   */
104:   DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
105:   DMSetFromOptions(user.dm);
106:   DMSetUp(user.dm);
107:   /*
108:      Extract global and local vectors from DM; The local vectors are
109:      used solely as work space for the evaluation of the function,
110:      gradient, and Hessian.  Duplicate for remaining vectors that are
111:      the same types.
112:   */
113:   DMCreateGlobalVector(user.dm,&x); /* Solution */
114:   DMCreateLocalVector(user.dm,&user.localX);
115:   VecDuplicate(user.localX,&user.localV);

117:   VecDuplicate(x,&xl);
118:   VecDuplicate(x,&xu);

120:   /* The TAO code begins here */

122:   /*
123:      Create TAO solver and set desired solution method
124:      The method must either be TAOTRON or TAOBLMVM
125:      If TAOBLMVM is used, then hessian function is not called.
126:   */
127:   TaoCreate(PETSC_COMM_WORLD,&tao);
128:   TaoSetType(tao,TAOBLMVM);

130:   /* Set initial solution guess; */
131:   MSA_BoundaryConditions(&user);
132:   MSA_InitialPoint(&user,x);
133:   TaoSetInitialVector(tao,x);

135:   /* Set routines for function, gradient and hessian evaluation */
136:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);

138:   VecGetLocalSize(x,&m);
139:   MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));
140:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);

142:   DMGetLocalToGlobalMapping(user.dm,&isltog);
143:   MatSetLocalToGlobalMapping(user.H,isltog,isltog);
144:   PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);
145:   if (flg) {
146:       MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
147:       MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);
148:       MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);
149:       TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);
150:   } else {
151:       TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);
152:   }

154:   /* Set Variable bounds */
155:   MSA_Plate(xl,xu,(void*)&user);
156:   TaoSetVariableBounds(tao,xl,xu);

158:   /* Check for any tao command line options */
159:   TaoSetFromOptions(tao);

161:   /* SOLVE THE APPLICATION */
162:   TaoSolve(tao);

164:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

166:   /* Free TAO data structures */
167:   TaoDestroy(&tao);

169:   /* Free PETSc data structures */
170:   VecDestroy(&x);
171:   VecDestroy(&xl);
172:   VecDestroy(&xu);
173:   MatDestroy(&user.H);
174:   VecDestroy(&user.localX);
175:   VecDestroy(&user.localV);
176:   VecDestroy(&user.Bottom);
177:   VecDestroy(&user.Top);
178:   VecDestroy(&user.Left);
179:   VecDestroy(&user.Right);
180:   DMDestroy(&user.dm);
181:   if (flg) {
182:     MatDestroy(&H_shell);
183:   }
184:   PetscFinalize();
185:   return ierr;
186: }

188: /*  FormFunctionGradient - Evaluates f(x) and gradient g(x).

190:     Input Parameters:
191: .   tao     - the Tao context
192: .   X      - input vector
193: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

195:     Output Parameters:
196: .   fcn     - the function value
197: .   G      - vector containing the newly evaluated gradient

199:    Notes:
200:    In this case, we discretize the domain and Create triangles. The
201:    surface of each triangle is planar, whose surface area can be easily
202:    computed.  The total surface area is found by sweeping through the grid
203:    and computing the surface area of the two triangles that have their
204:    right angle at the grid point.  The diagonal line segments on the
205:    grid that define the triangles run from top left to lower right.
206:    The numbering of points starts at the lower left and runs left to
207:    right, then bottom to top.
208: */
209: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
210: {
211:   AppCtx         *user = (AppCtx *) userCtx;
213:   PetscInt       i,j,row;
214:   PetscInt       mx=user->mx, my=user->my;
215:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
216:   PetscReal      ft=0;
217:   PetscReal      zero=0.0;
218:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
219:   PetscReal      rhx=mx+1, rhy=my+1;
220:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
221:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
222:   PetscReal      *g, *x,*left,*right,*bottom,*top;
223:   Vec            localX = user->localX, localG = user->localV;

225:   /* Get local mesh boundaries */
226:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
227:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

229:   /* Scatter ghost points to local vector */
230:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
231:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

233:   /* Initialize vector to zero */
234:   VecSet(localG, zero);

236:   /* Get pointers to vector data */
237:   VecGetArray(localX,&x);
238:   VecGetArray(localG,&g);
239:   VecGetArray(user->Top,&top);
240:   VecGetArray(user->Bottom,&bottom);
241:   VecGetArray(user->Left,&left);
242:   VecGetArray(user->Right,&right);

244:   /* Compute function over the locally owned part of the mesh */
245:   for (j=ys; j<ys+ym; j++) {
246:     for (i=xs; i< xs+xm; i++) {
247:       row=(j-gys)*gxm + (i-gxs);

249:       xc = x[row];
250:       xlt=xrb=xl=xr=xb=xt=xc;

252:       if (i==0) { /* left side */
253:         xl= left[j-ys+1];
254:         xlt = left[j-ys+2];
255:       } else {
256:         xl = x[row-1];
257:       }

259:       if (j==0) { /* bottom side */
260:         xb=bottom[i-xs+1];
261:         xrb = bottom[i-xs+2];
262:       } else {
263:         xb = x[row-gxm];
264:       }

266:       if (i+1 == gxs+gxm) { /* right side */
267:         xr=right[j-ys+1];
268:         xrb = right[j-ys];
269:       } else {
270:         xr = x[row+1];
271:       }

273:       if (j+1==gys+gym) { /* top side */
274:         xt=top[i-xs+1];
275:         xlt = top[i-xs];
276:       }else {
277:         xt = x[row+gxm];
278:       }

280:       if (i>gxs && j+1<gys+gym) {
281:         xlt = x[row-1+gxm];
282:       }
283:       if (j>gys && i+1<gxs+gxm) {
284:         xrb = x[row+1-gxm];
285:       }

287:       d1 = (xc-xl);
288:       d2 = (xc-xr);
289:       d3 = (xc-xt);
290:       d4 = (xc-xb);
291:       d5 = (xr-xrb);
292:       d6 = (xrb-xb);
293:       d7 = (xlt-xl);
294:       d8 = (xt-xlt);

296:       df1dxc = d1*hydhx;
297:       df2dxc = (d1*hydhx + d4*hxdhy);
298:       df3dxc = d3*hxdhy;
299:       df4dxc = (d2*hydhx + d3*hxdhy);
300:       df5dxc = d2*hydhx;
301:       df6dxc = d4*hxdhy;

303:       d1 *= rhx;
304:       d2 *= rhx;
305:       d3 *= rhy;
306:       d4 *= rhy;
307:       d5 *= rhy;
308:       d6 *= rhx;
309:       d7 *= rhy;
310:       d8 *= rhx;

312:       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
313:       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
314:       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
315:       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
316:       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
317:       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);

319:       ft = ft + (f2 + f4);

321:       df1dxc /= f1;
322:       df2dxc /= f2;
323:       df3dxc /= f3;
324:       df4dxc /= f4;
325:       df5dxc /= f5;
326:       df6dxc /= f6;

328:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;

330:     }
331:   }

333:   /* Compute triangular areas along the border of the domain. */
334:   if (xs==0) { /* left side */
335:     for (j=ys; j<ys+ym; j++) {
336:       d3=(left[j-ys+1] - left[j-ys+2])*rhy;
337:       d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
338:       ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
339:     }
340:   }
341:   if (ys==0) { /* bottom side */
342:     for (i=xs; i<xs+xm; i++) {
343:       d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
344:       d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
345:       ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
346:     }
347:   }

349:   if (xs+xm==mx) { /* right side */
350:     for (j=ys; j< ys+ym; j++) {
351:       d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
352:       d4=(right[j-ys]-right[j-ys+1])*rhy;
353:       ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
354:     }
355:   }
356:   if (ys+ym==my) { /* top side */
357:     for (i=xs; i<xs+xm; i++) {
358:       d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
359:       d4=(top[i-xs+1] - top[i-xs])*rhx;
360:       ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
361:     }
362:   }

364:   if (ys==0 && xs==0) {
365:     d1=(left[0]-left[1])*rhy;
366:     d2=(bottom[0]-bottom[1])*rhx;
367:     ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
368:   }
369:   if (ys+ym == my && xs+xm == mx) {
370:     d1=(right[ym+1] - right[ym])*rhy;
371:     d2=(top[xm+1] - top[xm])*rhx;
372:     ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
373:   }

375:   ft=ft*area;
376:   MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);

378:   /* Restore vectors */
379:   VecRestoreArray(localX,&x);
380:   VecRestoreArray(localG,&g);
381:   VecRestoreArray(user->Left,&left);
382:   VecRestoreArray(user->Top,&top);
383:   VecRestoreArray(user->Bottom,&bottom);
384:   VecRestoreArray(user->Right,&right);

386:   /* Scatter values to global vector */
387:   DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);
388:   DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);

390:   PetscLogFlops(70.0*xm*ym);

392:   return 0;
393: }

395: /* ------------------------------------------------------------------- */
396: /*
397:    FormHessian - Evaluates Hessian matrix.

399:    Input Parameters:
400: .  tao  - the Tao context
401: .  x    - input vector
402: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

404:    Output Parameters:
405: .  A    - Hessian matrix
406: .  B    - optionally different preconditioning matrix

408:    Notes:
409:    Due to mesh point reordering with DMs, we must always work
410:    with the local mesh points, and then transform them to the new
411:    global numbering with the local-to-global mapping.  We cannot work
412:    directly with the global numbers for the original uniprocessor mesh!

414:    Two methods are available for imposing this transformation
415:    when setting matrix entries:
416:      (A) MatSetValuesLocal(), using the local ordering (including
417:          ghost points!)
418:          - Do the following two steps once, before calling TaoSolve()
419:            - Use DMGetISLocalToGlobalMapping() to extract the
420:              local-to-global map from the DM
421:            - Associate this map with the matrix by calling
422:              MatSetLocalToGlobalMapping()
423:          - Then set matrix entries using the local ordering
424:            by calling MatSetValuesLocal()
425:      (B) MatSetValues(), using the global ordering
426:          - Use DMGetGlobalIndices() to extract the local-to-global map
427:          - Then apply this map explicitly yourself
428:          - Set matrix entries using the global ordering by calling
429:            MatSetValues()
430:    Option (A) seems cleaner/easier in many cases, and is the procedure
431:    used in this example.
432: */
433: PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
434: {
436:   AppCtx         *user = (AppCtx *) ptr;
437:   PetscInt       i,j,k,row;
438:   PetscInt       mx=user->mx, my=user->my;
439:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
440:   PetscReal      hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
441:   PetscReal      rhx=mx+1, rhy=my+1;
442:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
443:   PetscReal      hl,hr,ht,hb,hc,htl,hbr;
444:   PetscReal      *x,*left,*right,*bottom,*top;
445:   PetscReal      v[7];
446:   Vec            localX = user->localX;
447:   PetscBool      assembled;

449:   /* Set various matrix options */
450:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

452:   /* Initialize matrix entries to zero */
453:   MatAssembled(Hessian,&assembled);
454:   if (assembled) {MatZeroEntries(Hessian);}

456:   /* Get local mesh boundaries */
457:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
458:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

460:   /* Scatter ghost points to local vector */
461:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
462:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

464:   /* Get pointers to vector data */
465:   VecGetArray(localX,&x);
466:   VecGetArray(user->Top,&top);
467:   VecGetArray(user->Bottom,&bottom);
468:   VecGetArray(user->Left,&left);
469:   VecGetArray(user->Right,&right);

471:   /* Compute Hessian over the locally owned part of the mesh */

473:   for (i=xs; i< xs+xm; i++) {

475:     for (j=ys; j<ys+ym; j++) {

477:       row=(j-gys)*gxm + (i-gxs);

479:       xc = x[row];
480:       xlt=xrb=xl=xr=xb=xt=xc;

482:       /* Left side */
483:       if (i==gxs) {
484:         xl= left[j-ys+1];
485:         xlt = left[j-ys+2];
486:       } else {
487:         xl = x[row-1];
488:       }

490:       if (j==gys) {
491:         xb=bottom[i-xs+1];
492:         xrb = bottom[i-xs+2];
493:       } else {
494:         xb = x[row-gxm];
495:       }

497:       if (i+1 == gxs+gxm) {
498:         xr=right[j-ys+1];
499:         xrb = right[j-ys];
500:       } else {
501:         xr = x[row+1];
502:       }

504:       if (j+1==gys+gym) {
505:         xt=top[i-xs+1];
506:         xlt = top[i-xs];
507:       }else {
508:         xt = x[row+gxm];
509:       }

511:       if (i>gxs && j+1<gys+gym) {
512:         xlt = x[row-1+gxm];
513:       }
514:       if (j>gys && i+1<gxs+gxm) {
515:         xrb = x[row+1-gxm];
516:       }

518:       d1 = (xc-xl)*rhx;
519:       d2 = (xc-xr)*rhx;
520:       d3 = (xc-xt)*rhy;
521:       d4 = (xc-xb)*rhy;
522:       d5 = (xrb-xr)*rhy;
523:       d6 = (xrb-xb)*rhx;
524:       d7 = (xlt-xl)*rhy;
525:       d8 = (xlt-xt)*rhx;

527:       f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
528:       f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
529:       f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
530:       f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
531:       f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
532:       f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);

534:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
535:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
536:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
537:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
538:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
539:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
540:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
541:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

543:       hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
544:       htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);

546:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
547:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
548:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
549:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

551:       hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5;  hc*=0.5;

553:       k=0;
554:       if (j>0) {
555:         v[k]=hb; col[k]=row - gxm; k++;
556:       }

558:       if (j>0 && i < mx -1) {
559:         v[k]=hbr; col[k]=row - gxm+1; k++;
560:       }

562:       if (i>0) {
563:         v[k]= hl; col[k]=row - 1; k++;
564:       }

566:       v[k]= hc; col[k]=row; k++;

568:       if (i < mx-1) {
569:         v[k]= hr; col[k]=row+1; k++;
570:       }

572:       if (i>0 && j < my-1) {
573:         v[k]= htl; col[k] = row+gxm-1; k++;
574:       }

576:       if (j < my-1) {
577:         v[k]= ht; col[k] = row+gxm; k++;
578:       }

580:       /*
581:          Set matrix values using local numbering, which was defined
582:          earlier, in the main routine.
583:       */
584:       MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);

586:     }
587:   }

589:   /* Restore vectors */
590:   VecRestoreArray(localX,&x);
591:   VecRestoreArray(user->Left,&left);
592:   VecRestoreArray(user->Top,&top);
593:   VecRestoreArray(user->Bottom,&bottom);
594:   VecRestoreArray(user->Right,&right);

596:   /* Assemble the matrix */
597:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
598:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

600:   PetscLogFlops(199.0*xm*ym);
601:   return 0;
602: }

604: /* ------------------------------------------------------------------- */
605: /*
606:    MSA_BoundaryConditions -  Calculates the boundary conditions for
607:    the region.

609:    Input Parameter:
610: .  user - user-defined application context

612:    Output Parameter:
613: .  user - user-defined application context
614: */
615: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
616: {
617:   int        ierr;
618:   PetscInt   i,j,k,maxits=5,limit=0;
619:   PetscInt   xs,ys,xm,ym,gxs,gys,gxm,gym;
620:   PetscInt   mx=user->mx,my=user->my;
621:   PetscInt   bsize=0, lsize=0, tsize=0, rsize=0;
622:   PetscReal  one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
623:   PetscReal  fnorm,det,hx,hy,xt=0,yt=0;
624:   PetscReal  u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
625:   PetscReal  b=-0.5, t=0.5, l=-0.5, r=0.5;
626:   PetscReal  *boundary;
627:   PetscBool  flg;
628:   Vec        Bottom,Top,Right,Left;

630:   /* Get local mesh boundaries */
631:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
632:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

634:   bsize=xm+2;
635:   lsize=ym+2;
636:   rsize=ym+2;
637:   tsize=xm+2;

639:   VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);
640:   VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);
641:   VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);
642:   VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);

644:   user->Top=Top;
645:   user->Left=Left;
646:   user->Bottom=Bottom;
647:   user->Right=Right;

649:   hx= (r-l)/(mx+1); hy=(t-b)/(my+1);

651:   for (j=0; j<4; j++) {
652:     if (j==0) {
653:       yt=b;
654:       xt=l+hx*xs;
655:       limit=bsize;
656:       VecGetArray(Bottom,&boundary);
657:     } else if (j==1) {
658:       yt=t;
659:       xt=l+hx*xs;
660:       limit=tsize;
661:       VecGetArray(Top,&boundary);
662:     } else if (j==2) {
663:       yt=b+hy*ys;
664:       xt=l;
665:       limit=lsize;
666:       VecGetArray(Left,&boundary);
667:     } else if (j==3) {
668:       yt=b+hy*ys;
669:       xt=r;
670:       limit=rsize;
671:       VecGetArray(Right,&boundary);
672:     }

674:     for (i=0; i<limit; i++) {
675:       u1=xt;
676:       u2=-yt;
677:       for (k=0; k<maxits; k++) {
678:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
679:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
680:         fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
681:         if (fnorm <= tol) break;
682:         njac11=one+u2*u2-u1*u1;
683:         njac12=two*u1*u2;
684:         njac21=-two*u1*u2;
685:         njac22=-one - u1*u1 + u2*u2;
686:         det = njac11*njac22-njac21*njac12;
687:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
688:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
689:       }

691:       boundary[i]=u1*u1-u2*u2;
692:       if (j==0 || j==1) {
693:         xt=xt+hx;
694:       } else if (j==2 || j==3) {
695:         yt=yt+hy;
696:       }
697:     }
698:     if (j==0) {
699:       VecRestoreArray(Bottom,&boundary);
700:     } else if (j==1) {
701:       VecRestoreArray(Top,&boundary);
702:     } else if (j==2) {
703:       VecRestoreArray(Left,&boundary);
704:     } else if (j==3) {
705:       VecRestoreArray(Right,&boundary);
706:     }
707:   }

709:   /* Scale the boundary if desired */

711:   PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
712:   if (flg) {
713:     VecScale(Bottom, scl);
714:   }
715:   PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
716:   if (flg) {
717:     VecScale(Top, scl);
718:   }
719:   PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
720:   if (flg) {
721:     VecScale(Right, scl);
722:   }

724:   PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
725:   if (flg) {
726:     VecScale(Left, scl);
727:   }
728:   return 0;
729: }

731: /* ------------------------------------------------------------------- */
732: /*
733:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

735:    Input Parameter:
736: .  user - user-defined application context

738:    Output Parameter:
739: .  user - user-defined application context
740: */
741: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx)
742: {
743:   AppCtx         *user=(AppCtx *)ctx;
745:   PetscInt       i,j,row;
746:   PetscInt       xs,ys,xm,ym;
747:   PetscInt       mx=user->mx, my=user->my, bmy, bmx;
748:   PetscReal      t1,t2,t3;
749:   PetscReal      *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
750:   PetscBool      cylinder;

752:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
753:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
754:   bmy=user->bmy; bmx=user->bmx;

756:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

758:   VecSet(XL, lb);
759:   VecSet(XU, ub);

761:   VecGetArray(XL,&xl);

763:   PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);
764:   /* Compute the optional lower box */
765:   if (cylinder) {
766:     for (i=xs; i< xs+xm; i++) {
767:       for (j=ys; j<ys+ym; j++) {
768:         row=(j-ys)*xm + (i-xs);
769:         t1=(2.0*i-mx)*bmy;
770:         t2=(2.0*j-my)*bmx;
771:         t3=bmx*bmx*bmy*bmy;
772:         if (t1*t1 + t2*t2 <= t3) {
773:           xl[row] = user->bheight;
774:         }
775:       }
776:     }
777:   } else {
778:     /* Compute the optional lower box */
779:     for (i=xs; i< xs+xm; i++) {
780:       for (j=ys; j<ys+ym; j++) {
781:         row=(j-ys)*xm + (i-xs);
782:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
783:             j>=(my-bmy)/2 && j<my-(my-bmy)/2) {
784:           xl[row] = user->bheight;
785:         }
786:       }
787:     }
788:   }
789:     VecRestoreArray(XL,&xl);

791:   return 0;
792: }

794: /* ------------------------------------------------------------------- */
795: /*
796:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

798:    Input Parameters:
799: .  user - user-defined application context
800: .  X - vector for initial guess

802:    Output Parameters:
803: .  X - newly computed initial guess
804: */
805: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
806: {
808:   PetscInt       start=-1,i,j;
809:   PetscReal      zero=0.0;
810:   PetscBool      flg;

812:   PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
813:   if (flg && start==0) { /* The zero vector is reasonable */
814:     VecSet(X, zero);
815:   } else if (flg && start>0) { /* Try a random start between -0.5 and 0.5 */
816:     PetscRandom rctx;  PetscReal np5=-0.5;

818:     PetscRandomCreate(MPI_COMM_WORLD,&rctx);
819:     for (i=0; i<start; i++) {
820:       VecSetRandom(X, rctx);
821:     }
822:     PetscRandomDestroy(&rctx);
823:     VecShift(X, np5);

825:   } else { /* Take an average of the boundary conditions */

827:     PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
828:     PetscInt mx=user->mx,my=user->my;
829:     PetscReal *x,*left,*right,*bottom,*top;
830:     Vec    localX = user->localX;

832:     /* Get local mesh boundaries */
833:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
834:     DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

836:     /* Get pointers to vector data */
837:     VecGetArray(user->Top,&top);
838:     VecGetArray(user->Bottom,&bottom);
839:     VecGetArray(user->Left,&left);
840:     VecGetArray(user->Right,&right);

842:     VecGetArray(localX,&x);
843:     /* Perform local computations */
844:     for (j=ys; j<ys+ym; j++) {
845:       for (i=xs; i< xs+xm; i++) {
846:         row=(j-gys)*gxm + (i-gxs);
847:         x[row] = ((j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+(i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
848:       }
849:     }

851:     /* Restore vectors */
852:     VecRestoreArray(localX,&x);

854:     VecRestoreArray(user->Left,&left);
855:     VecRestoreArray(user->Top,&top);
856:     VecRestoreArray(user->Bottom,&bottom);
857:     VecRestoreArray(user->Right,&right);

859:     /* Scatter values into global vector */
860:     DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);
861:     DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);

863:   }
864:   return 0;
865: }

867: /* For testing matrix free submatrices */
868: PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
869: {
871:   AppCtx         *user = (AppCtx*)ptr;
873:   FormHessian(tao,x,user->H,user->H,ptr);
874:   return(0);
875: }
876: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
877: {
879:   void           *ptr;
880:   AppCtx         *user;
882:   MatShellGetContext(H_shell,&ptr);
883:   user = (AppCtx*)ptr;
884:   MatMult(user->H,X,Y);
885:   return(0);
886: }

888: /*TEST

890:    build:
891:       requires: !complex

893:    test:
894:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
895:       requires: !single

897:    test:
898:       suffix: 2
899:       nsize: 2
900:       args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
901:       requires: !single

903:    test:
904:       suffix: 3
905:       nsize: 3
906:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
907:       requires: !single

909:    test:
910:       suffix: 4
911:       nsize: 3
912:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
913:       requires: !single

915:    test:
916:       suffix: 5
917:       nsize: 3
918:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
919:       requires: !single

921:    test:
922:       suffix: 6
923:       nsize: 3
924:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
925:       requires: !single

927:    test:
928:       suffix: 7
929:       nsize: 3
930:       args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
931:       requires: !single

933:    test:
934:       suffix: 8
935:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
936:       requires: !single

938:    test:
939:       suffix: 9
940:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
941:       requires: !single

943:    test:
944:       suffix: 10
945:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
946:       requires: !single

948:    test:
949:       suffix: 11
950:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
951:       requires: !single

953:    test:
954:       suffix: 12
955:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
956:       requires: !single

958:    test:
959:       suffix: 13
960:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
961:       requires: !single

963:    test:
964:       suffix: 14
965:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
966:       requires: !single

968:    test:
969:       suffix: 15
970:       args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
971:       requires: !single

973:    test:
974:      suffix: 16
975:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
976:      requires: !single

978:    test:
979:      suffix: 17
980:      args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
981:      requires: !single

983:    test:
984:      suffix: 18
985:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
986:      requires: !single

988:    test:
989:      suffix: 19
990:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
991:      requires: !single

993:    test:
994:      suffix: 20
995:      args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian

997: TEST*/
```