Actual source code: plate2.c

tao-2.1-p0 2012-07-24
  1: #include "petscdm.h"   /*I "petscdm.h" I*/
  2: #include "taosolver.h" /*I "taosolver.h" I*/

  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: TaoInitialize(); TaoFinalize();
 23:    Routines: TaoCreate();
 24:    Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
 25:    Routines: TaoSetHessianRoutine();
 26:    Routines: TaoSetInitialVector();
 27:    Routines: TaoSetVariableBounds();
 28:    Routines: TaoSetFromOptions();
 29:    Routines: TaoSolve(); TaoView();
 30:    Routines: TaoGetTerminationReason(); TaoDestroy(); 
 31:    Processors: n
 32: T*/ 


 35: /* 
 36:    User-defined application context - contains data needed by the 
 37:    application-provided call-back routines, FormFunctionGradient(),
 38:    FormHessian().
 39: */
 40: typedef struct {
 41:   /* problem parameters */
 42:   PetscReal      bheight;                  /* Height of plate under the surface */
 43:   PetscInt       mx, my;                   /* discretization in x, y directions */
 44:   PetscInt       bmx,bmy;                  /* Size of plate under the surface */
 45:   Vec         Bottom, Top, Left, Right; /* boundary values */
 46:   
 47:   /* Working space */
 48:   Vec         localX, localV;           /* ghosted local vector */
 49:   DM          dm;                       /* distributed array data structure */
 50:   Mat         H;
 51: } AppCtx;

 53: /* -------- User-defined Routines --------- */

 55: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
 56: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
 57: static PetscErrorCode MSA_Plate(Vec,Vec,void*);
 58: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal*,Vec,void*);
 59: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);

 61: /* For testing matrix free submatrices */
 62: PetscErrorCode MatrixFreeHessian(TaoSolver,Vec,Mat*, Mat*,MatStructure*,void*);
 63: PetscErrorCode MyMatMult(Mat,Vec,Vec);

 67: int main( int argc, char **argv )
 68: {
 70:   PetscInt   Nx, Ny;               /* number of processors in x- and y- directions */
 71:   PetscInt   m, N;                 /* number of local and global elements in vectors */
 72:   Vec        x,xl,xu;               /* solution vector  and bounds*/
 73:   PetscBool   flg;                /* A return variable when checking for user options */
 74:   TaoSolver  tao;                  /* TaoSolver solver context */
 75:   ISLocalToGlobalMapping isltog;   /* local-to-global mapping object */
 76:   TaoSolverTerminationReason reason;
 77:   Mat         H_shell;                  /* to test matrix-free submatrices */
 78:   AppCtx     user;                 /* user-defined work context */

 80:   /* Initialize PETSc, TAO */
 81:   PetscInitialize( &argc, &argv,(char *)0,help );
 82:   TaoInitialize( &argc, &argv,(char *)0,help );

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

 87:   /* Check for any command line arguments that override defaults */
 88:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); 
 89:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); 
 90:   PetscOptionsGetReal(PETSC_NULL,"-bheight",&user.bheight,&flg); 

 92:   user.bmx = user.mx/2; user.bmy = user.my/2;
 93:   PetscOptionsGetInt(PETSC_NULL,"-bmx",&user.bmx,&flg); 
 94:   PetscOptionsGetInt(PETSC_NULL,"-bmy",&user.bmy,&flg); 

 96:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");
 97:   PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%G\n",
 98:               user.mx,user.my,user.bmx,user.bmy,user.bheight);

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

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

106:   /*
107:      A two dimensional distributed array will help define this problem,
108:      which derives from an elliptic PDE on two dimensional domain.  From
109:      the distributed array, Create the vectors.
110:   */
111:   DMDACreate2d(MPI_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
112:                       DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,
113:                       PETSC_NULL,PETSC_NULL,&user.dm); 

115:   /*
116:      Extract global and local vectors from DM; The local vectors are
117:      used solely as work space for the evaluation of the function, 
118:      gradient, and Hessian.  Duplicate for remaining vectors that are 
119:      the same types.
120:   */
121:   DMCreateGlobalVector(user.dm,&x);  /* Solution */
122:   DMCreateLocalVector(user.dm,&user.localX); 
123:   VecDuplicate(user.localX,&user.localV); 

125:   VecDuplicate(x,&xl); 
126:   VecDuplicate(x,&xu); 



130:   /* The TAO code begins here */

132:   /* 
133:      Create TAO solver and set desired solution method 
134:      The method must either be 'tao_tron' or 'tao_blmvm'
135:      If blmvm is used, then hessian function is not called.
136:   */
137:   TaoCreate(PETSC_COMM_WORLD,&tao); 
138:   TaoSetType(tao,"tao_blmvm"); 

140:   /* Set initial solution guess; */
141:   MSA_BoundaryConditions(&user); 
142:   MSA_InitialPoint(&user,x); 
143:   TaoSetInitialVector(tao,x); 
144:   
145:   /* Set routines for function, gradient and hessian evaluation */
146:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user); 
147:   

149:   VecGetLocalSize(x,&m); 
150:   MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL,
151:                          3,PETSC_NULL,&(user.H)); 
152:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE); 
153:   
154:   DMGetLocalToGlobalMapping(user.dm,&isltog); 
155:   MatSetLocalToGlobalMapping(user.H,isltog,isltog); 
156:   PetscOptionsHasName(PETSC_NULL,"-matrixfree",&flg); 
157:   if (flg) {
158:       MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);
159:       MatShellSetOperation(H_shell,MATOP_MULT,(void(*)())MyMatMult); 
160:       MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE); 
161:       TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user); 
162:   } else {
163:       TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user); 
164:   }
165:   

167:   /* Set Variable bounds */
168:   MSA_Plate(xl,xu,(void*)&user); 
169:   TaoSetVariableBounds(tao,xl,xu); 

171:   /* Check for any tao command line options */
172:   TaoSetFromOptions(tao); 

174:   /* SOLVE THE APPLICATION */
175:   TaoSolve(tao);  

177:   /* Get ierrrmation on termination */
178:   TaoGetTerminationReason(tao,&reason); 
179:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); 
180:   if (reason <= 0){
181:     PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
182:   }
183:   /* Free TAO data structures */
184:   TaoDestroy(&tao); 

186:   /* Free PETSc data structures */
187:   VecDestroy(&x); 
188:   VecDestroy(&xl); 
189:   VecDestroy(&xu); 
190:   MatDestroy(&user.H); 
191:   VecDestroy(&user.localX);  
192:   VecDestroy(&user.localV); 
193:   VecDestroy(&user.Bottom); 
194:   VecDestroy(&user.Top); 
195:   VecDestroy(&user.Left); 
196:   VecDestroy(&user.Right); 
197:   DMDestroy(&user.dm); 
198:   if (flg) {
199:     MatDestroy(&H_shell); 
200:   }
201:   /* Finalize TAO and PETSc */
202:   TaoFinalize();
203:   PetscFinalize();
204:   
205:   return 0;
206: }

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

212:     Input Parameters:
213: .   tao     - the TaoSolver context
214: .   X      - input vector
215: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
216:     
217:     Output Parameters:
218: .   fcn     - the function value
219: .   G      - vector containing the newly evaluated gradient

221:    Notes:
222:    In this case, we discretize the domain and Create triangles. The 
223:    surface of each triangle is planar, whose surface area can be easily 
224:    computed.  The total surface area is found by sweeping through the grid 
225:    and computing the surface area of the two triangles that have their 
226:    right angle at the grid point.  The diagonal line segments on the
227:    grid that define the triangles run from top left to lower right.
228:    The numbering of points starts at the lower left and runs left to
229:    right, then bottom to top.
230: */
231: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
232: {
233:   AppCtx * user = (AppCtx *) userCtx;
234:   PetscErrorCode    ierr;
235:   PetscInt i,j,row;
236:   PetscInt mx=user->mx, my=user->my;
237:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
238:   PetscReal ft=0;
239:   PetscReal zero=0.0;
240:   PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
241:   PetscReal rhx=mx+1, rhy=my+1;
242:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
243:   PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
244:   PetscReal *g, *x,*left,*right,*bottom,*top;
245:   Vec    localX = user->localX, localG = user->localV;

247:   /* Get local mesh boundaries */
248:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
249:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

251:   /* Scatter ghost points to local vector */
252:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
253:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

255:   /* Initialize vector to zero */
256:   VecSet(localG, zero); 

258:   /* Get pointers to vector data */
259:   VecGetArray(localX,&x); 
260:   VecGetArray(localG,&g); 
261:   VecGetArray(user->Top,&top); 
262:   VecGetArray(user->Bottom,&bottom); 
263:   VecGetArray(user->Left,&left); 
264:   VecGetArray(user->Right,&right); 

266:   /* Compute function over the locally owned part of the mesh */
267:   for (j=ys; j<ys+ym; j++){
268:     for (i=xs; i< xs+xm; i++){
269:       row=(j-gys)*gxm + (i-gxs);
270:       
271:       xc = x[row];
272:       xlt=xrb=xl=xr=xb=xt=xc;
273:       
274:       if (i==0){ /* left side */
275:         xl= left[j-ys+1];
276:         xlt = left[j-ys+2];
277:       } else {
278:         xl = x[row-1];
279:       }

281:       if (j==0){ /* bottom side */
282:         xb=bottom[i-xs+1];
283:         xrb = bottom[i-xs+2];
284:       } else {
285:         xb = x[row-gxm];
286:       }
287:       
288:       if (i+1 == gxs+gxm){ /* right side */
289:         xr=right[j-ys+1];
290:         xrb = right[j-ys];
291:       } else {
292:         xr = x[row+1];
293:       }

295:       if (j+1==gys+gym){ /* top side */
296:         xt=top[i-xs+1];
297:         xlt = top[i-xs];
298:       }else {
299:         xt = x[row+gxm];
300:       }

302:       if (i>gxs && j+1<gys+gym){
303:         xlt = x[row-1+gxm];
304:       }
305:       if (j>gys && i+1<gxs+gxm){
306:         xrb = x[row+1-gxm];
307:       }

309:       d1 = (xc-xl);
310:       d2 = (xc-xr);
311:       d3 = (xc-xt);
312:       d4 = (xc-xb);
313:       d5 = (xr-xrb);
314:       d6 = (xrb-xb);
315:       d7 = (xlt-xl);
316:       d8 = (xt-xlt);
317:       
318:       df1dxc = d1*hydhx;
319:       df2dxc = ( d1*hydhx + d4*hxdhy );
320:       df3dxc = d3*hxdhy;
321:       df4dxc = ( d2*hydhx + d3*hxdhy );
322:       df5dxc = d2*hydhx;
323:       df6dxc = d4*hxdhy;

325:       d1 *= rhx;
326:       d2 *= rhx;
327:       d3 *= rhy;
328:       d4 *= rhy;
329:       d5 *= rhy;
330:       d6 *= rhx;
331:       d7 *= rhy;
332:       d8 *= rhx;

334:       f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
335:       f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
336:       f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
337:       f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
338:       f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
339:       f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);
340:       
341:       ft = ft + (f2 + f4);

343:       df1dxc /= f1;
344:       df2dxc /= f2;
345:       df3dxc /= f3;
346:       df4dxc /= f4;
347:       df5dxc /= f5;
348:       df6dxc /= f6;

350:       g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc ) * 0.5;
351:       
352:     }
353:   }


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=(left[j-ys+1] - left[j-ys+2])*rhy;
360:       d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
361:       ft = ft+PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
362:     }
363:   }
364:   if (ys==0){ /* bottom side */
365:     for (i=xs; i<xs+xm; i++){
366:       d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
367:       d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
368:       ft = ft+PetscSqrtScalar( 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+1-gys)*gxm-1]-right[j-ys+1])*rhx;
375:       d4=(right[j-ys]-right[j-ys+1])*rhy;
376:       ft = ft+PetscSqrtScalar( 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[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
382:       d4=(top[i-xs+1] - top[i-xs])*rhx;
383:       ft = ft+PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
384:     }
385:   }

387:   if (ys==0 && xs==0){
388:     d1=(left[0]-left[1])*rhy;
389:     d2=(bottom[0]-bottom[1])*rhx;
390:     ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
391:   }
392:   if (ys+ym == my && xs+xm == mx){
393:     d1=(right[ym+1] - right[ym])*rhy;
394:     d2=(top[xm+1] - top[xm])*rhx;
395:     ft +=PetscSqrtScalar( 1.0 + d1*d1 + d2*d2);
396:   }

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

401:   
402:   /* Restore vectors */
403:   VecRestoreArray(localX,&x); 
404:   VecRestoreArray(localG,&g); 
405:   VecRestoreArray(user->Left,&left); 
406:   VecRestoreArray(user->Top,&top); 
407:   VecRestoreArray(user->Bottom,&bottom); 
408:   VecRestoreArray(user->Right,&right); 

410:   /* Scatter values to global vector */
411:   DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G); 
412:   DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G); 

414:   PetscLogFlops(70*xm*ym); 

416:   return 0;
417: }

419: /* ------------------------------------------------------------------- */
422: /*
423:    FormHessian - Evaluates Hessian matrix.

425:    Input Parameters:
426: .  tao  - the TaoSolver context
427: .  x    - input vector
428: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

430:    Output Parameters:
431: .  A    - Hessian matrix
432: .  B    - optionally different preconditioning matrix
433: .  flag - flag indicating matrix structure

435:    Notes:
436:    Due to mesh point reordering with DMs, we must always work
437:    with the local mesh points, and then transform them to the new
438:    global numbering with the local-to-global mapping.  We cannot work
439:    directly with the global numbers for the original uniprocessor mesh!  

441:    Two methods are available for imposing this transformation
442:    when setting matrix entries:
443:      (A) MatSetValuesLocal(), using the local ordering (including
444:          ghost points!)
445:          - Do the following two steps once, before calling TaoSolve()
446:            - Use DMGetISLocalToGlobalMapping() to extract the
447:              local-to-global map from the DM
448:            - Associate this map with the matrix by calling
449:              MatSetLocalToGlobalMapping() 
450:          - Then set matrix entries using the local ordering
451:            by calling MatSetValuesLocal()
452:      (B) MatSetValues(), using the global ordering 
453:          - Use DMGetGlobalIndices() to extract the local-to-global map
454:          - Then apply this map explicitly yourself
455:          - Set matrix entries using the global ordering by calling
456:            MatSetValues()
457:    Option (A) seems cleaner/easier in many cases, and is the procedure
458:    used in this example.
459: */
460: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *Hptr, Mat *Hpc, MatStructure *flag, void *ptr)
461: { 
462:   PetscErrorCode    ierr;
463:   AppCtx *user = (AppCtx *) ptr;
464:   Mat Hessian = *Hpc;
465:   PetscInt   i,j,k,row;
466:   PetscInt   mx=user->mx, my=user->my;
467:   PetscInt   xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
468:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
469:   PetscReal rhx=mx+1, rhy=my+1;
470:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
471:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
472:   PetscReal *x,*left,*right,*bottom,*top;
473:   PetscReal v[7];
474:   Vec    localX = user->localX;
475:   PetscBool assembled;


478:   /* Set various matrix options */
479:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE); 

481:   /* Initialize matrix entries to zero */
482:   MatAssembled(Hessian,&assembled); 
483:   if (assembled){MatZeroEntries(Hessian);  }

485:   /* Get local mesh boundaries */
486:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
487:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

489:   /* Scatter ghost points to local vector */
490:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
491:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

493:   /* Get pointers to vector data */
494:   VecGetArray(localX,&x); 
495:   VecGetArray(user->Top,&top); 
496:   VecGetArray(user->Bottom,&bottom); 
497:   VecGetArray(user->Left,&left); 
498:   VecGetArray(user->Right,&right); 

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

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

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

506:       row=(j-gys)*gxm + (i-gxs);
507:       
508:       xc = x[row]; 
509:       xlt=xrb=xl=xr=xb=xt=xc;

511:       /* Left side */
512:       if (i==gxs){
513:         xl= left[j-ys+1];
514:         xlt = left[j-ys+2];
515:       } else {
516:         xl = x[row-1];
517:       }
518:       
519:       if (j==gys){
520:         xb=bottom[i-xs+1];
521:         xrb = bottom[i-xs+2];
522:       } else {
523:         xb = x[row-gxm];
524:       }
525:       
526:       if (i+1 == gxs+gxm){
527:         xr=right[j-ys+1];
528:         xrb = right[j-ys];
529:       } else {
530:         xr = x[row+1];
531:       }

533:       if (j+1==gys+gym){
534:         xt=top[i-xs+1];
535:         xlt = top[i-xs];
536:       }else {
537:         xt = x[row+gxm];
538:       }

540:       if (i>gxs && j+1<gys+gym){
541:         xlt = x[row-1+gxm];
542:       }
543:       if (j>gys && i+1<gxs+gxm){
544:         xrb = x[row+1-gxm];
545:       }


548:       d1 = (xc-xl)*rhx;
549:       d2 = (xc-xr)*rhx;
550:       d3 = (xc-xt)*rhy;
551:       d4 = (xc-xb)*rhy;
552:       d5 = (xrb-xr)*rhy;
553:       d6 = (xrb-xb)*rhx;
554:       d7 = (xlt-xl)*rhy;
555:       d8 = (xlt-xt)*rhx;
556:       
557:       f1 = PetscSqrtScalar( 1.0 + d1*d1 + d7*d7);
558:       f2 = PetscSqrtScalar( 1.0 + d1*d1 + d4*d4);
559:       f3 = PetscSqrtScalar( 1.0 + d3*d3 + d8*d8);
560:       f4 = PetscSqrtScalar( 1.0 + d3*d3 + d2*d2);
561:       f5 = PetscSqrtScalar( 1.0 + d2*d2 + d5*d5);
562:       f6 = PetscSqrtScalar( 1.0 + d4*d4 + d6*d6);


565:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
566:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
567:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
568:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
569:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
570:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
571:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
572:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

577:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
578:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
579:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
580:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

584:       k=0;
585:       if (j>0){ 
586:         v[k]=hb; col[k]=row - gxm; k++;
587:       }
588:       
589:       if (j>0 && i < mx -1){
590:         v[k]=hbr; col[k]=row - gxm+1; k++;
591:       }
592:       
593:       if (i>0){
594:         v[k]= hl; col[k]=row - 1; k++;
595:       }
596:       
597:       v[k]= hc; col[k]=row; k++;
598:       
599:       if (i < mx-1 ){
600:         v[k]= hr; col[k]=row+1; k++;
601:       }
602:       
603:       if (i>0 && j < my-1 ){
604:         v[k]= htl; col[k] = row+gxm-1; k++;
605:       }
606:       
607:       if (j < my-1 ){
608:         v[k]= ht; col[k] = row+gxm; k++;
609:       }
610:       
611:       /* 
612:          Set matrix values using local numbering, which was defined
613:          earlier, in the main routine.
614:       */
615:       MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES); 
616:       
617:       
618:     }
619:   }
620:   
621:   /* Restore vectors */
622:   VecRestoreArray(localX,&x); 
623:   VecRestoreArray(user->Left,&left); 
624:   VecRestoreArray(user->Top,&top); 
625:   VecRestoreArray(user->Bottom,&bottom); 
626:   VecRestoreArray(user->Right,&right); 

628:   /* Assemble the matrix */
629:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY); 
630:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY); 

632:   PetscLogFlops(199*xm*ym); 
633:   return 0;
634: }

636: /* ------------------------------------------------------------------- */
639: /* 
640:    MSA_BoundaryConditions -  Calculates the boundary conditions for
641:    the region.

643:    Input Parameter:
644: .  user - user-defined application context

646:    Output Parameter:
647: .  user - user-defined application context
648: */
649: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
650: {
651:   int        ierr;
652:   PetscInt   i,j,k,maxits=5,limit=0;
653:   PetscInt   xs,ys,xm,ym,gxs,gys,gxm,gym;
654:   PetscInt   mx=user->mx,my=user->my;
655:   PetscInt   bsize=0, lsize=0, tsize=0, rsize=0;
656:   PetscReal     one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
657:   PetscReal     fnorm,det,hx,hy,xt=0,yt=0;
658:   PetscReal     u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
659:   PetscReal     b=-0.5, t=0.5, l=-0.5, r=0.5;
660:   PetscReal     *boundary;
661:   PetscBool   flg;
662:   Vec        Bottom,Top,Right,Left;

664:   /* Get local mesh boundaries */
665:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
666:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 


669:   bsize=xm+2;
670:   lsize=ym+2;
671:   rsize=ym+2;
672:   tsize=xm+2;

674:   VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom); 
675:   VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top); 
676:   VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left); 
677:   VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right); 

679:   user->Top=Top;
680:   user->Left=Left;
681:   user->Bottom=Bottom;
682:   user->Right=Right;

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

686:   for (j=0; j<4; j++){
687:     if (j==0){
688:       yt=b;
689:       xt=l+hx*xs;
690:       limit=bsize;
691:       VecGetArray(Bottom,&boundary);
692:     } else if (j==1){
693:       yt=t;
694:       xt=l+hx*xs;
695:       limit=tsize;
696:       VecGetArray(Top,&boundary);
697:     } else if (j==2){
698:       yt=b+hy*ys;
699:       xt=l;
700:       limit=lsize;
701:       VecGetArray(Left,&boundary);
702:     } else if (j==3){
703:       yt=b+hy*ys;
704:       xt=r;
705:       limit=rsize;
706:       VecGetArray(Right,&boundary);
707:     }

709:     for (i=0; i<limit; i++){
710:       u1=xt;
711:       u2=-yt;
712:       for (k=0; k<maxits; k++){
713:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
714:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
715:         fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
716:         if (fnorm <= tol) break;
717:         njac11=one+u2*u2-u1*u1;
718:         njac12=two*u1*u2;
719:         njac21=-two*u1*u2;
720:         njac22=-one - u1*u1 + u2*u2;
721:         det = njac11*njac22-njac21*njac12;
722:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
723:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
724:       }

726:       boundary[i]=u1*u1-u2*u2;
727:       if (j==0 || j==1) {
728:         xt=xt+hx;
729:       } else if (j==2 || j==3){
730:         yt=yt+hy;
731:       }
732:       
733:     }
734:     
735:     if (j==0){
736:       VecRestoreArray(Bottom,&boundary); 
737:     } else if (j==1){
738:       VecRestoreArray(Top,&boundary); 
739:     } else if (j==2){
740:       VecRestoreArray(Left,&boundary); 
741:     } else if (j==3){
742:       VecRestoreArray(Right,&boundary); 
743:     }

745:   }

747:   /* Scale the boundary if desired */

749:   PetscOptionsGetReal(PETSC_NULL,"-bottom",&scl,&flg); 
750:   
751:   if (flg){
752:     VecScale(Bottom, scl); 
753:   }
754:   
755:   PetscOptionsGetReal(PETSC_NULL,"-top",&scl,&flg); 
756:   
757:   if (flg){
758:     VecScale(Top, scl); 
759:   }
760:   
761:   PetscOptionsGetReal(PETSC_NULL,"-right",&scl,&flg); 
762:   
763:   if (flg){
764:     VecScale(Right, scl); 
765:   }
766:   
767:   PetscOptionsGetReal(PETSC_NULL,"-left",&scl,&flg); 
768:   
769:   if (flg){
770:     VecScale(Left, scl); 
771:   }

773:   return 0;
774: }


777: /* ------------------------------------------------------------------- */
780: /* 
781:    MSA_Plate -  Calculates an obstacle for surface to stretch over.

783:    Input Parameter:
784: .  user - user-defined application context

786:    Output Parameter:
787: .  user - user-defined application context
788: */
789: static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){

791:   AppCtx *user=(AppCtx *)ctx;
793:   PetscInt i,j,row;
794:   PetscInt xs,ys,xm,ym;
795:   PetscInt mx=user->mx, my=user->my, bmy, bmx;
796:   PetscReal t1,t2,t3;
797:   PetscReal *xl, lb=TAO_NINFINITY, ub=TAO_INFINITY;
798:   PetscBool cylinder;

800:   user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
801:   user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
802:   bmy=user->bmy, bmx=user->bmx;

804:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 

806:   VecSet(XL, lb); 
807:   VecSet(XU, ub); 

809:   VecGetArray(XL,&xl); 

811:   PetscOptionsHasName(PETSC_NULL,"-cylinder",&cylinder); 
812:   /* Compute the optional lower box */
813:   if (cylinder){
814:     for (i=xs; i< xs+xm; i++){    
815:       for (j=ys; j<ys+ym; j++){
816:         row=(j-ys)*xm + (i-xs);
817:         t1=(2.0*i-mx)*bmy;
818:         t2=(2.0*j-my)*bmx;
819:         t3=bmx*bmx*bmy*bmy;
820:         if ( t1*t1 + t2*t2 <= t3 ){
821:           xl[row] = user->bheight;
822:         }
823:       }
824:     }
825:   } else {
826:     /* Compute the optional lower box */
827:     for (i=xs; i< xs+xm; i++){    
828:       for (j=ys; j<ys+ym; j++){
829:         row=(j-ys)*xm + (i-xs);
830:         if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 && 
831:             j>=(my-bmy)/2 && j<my-(my-bmy)/2 ){
832:           xl[row] = user->bheight;
833:         }
834:       }
835:     }
836:   }
837:     VecRestoreArray(XL,&xl); 

839:   return 0;
840: }


843: /* ------------------------------------------------------------------- */
846: /*
847:    MSA_InitialPoint - Calculates the initial guess in one of three ways. 

849:    Input Parameters:
850: .  user - user-defined application context
851: .  X - vector for initial guess

853:    Output Parameters:
854: .  X - newly computed initial guess
855: */
856: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
857: {
858:   PetscErrorCode      ierr;
859:   PetscInt start=-1,i,j;
860:   PetscReal   zero=0.0;
861:   PetscBool flg;

863:   PetscOptionsGetInt(PETSC_NULL,"-start",&start,&flg); 

865:   if (flg && start==0){ /* The zero vector is reasonable */
866:  
867:     VecSet(X, zero); 

869:   } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */

871:     PetscRandom rctx;  PetscReal np5=-0.5;

873:     PetscRandomCreate(MPI_COMM_WORLD,&rctx); 
874:     
875:     for (i=0; i<start; i++){
876:       VecSetRandom(X, rctx); 
877:     }
878:     PetscRandomDestroy(&rctx); 
879:     VecShift(X, np5);

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

883:     PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
884:     PetscInt mx=user->mx,my=user->my;
885:     PetscReal *x,*left,*right,*bottom,*top;
886:     Vec    localX = user->localX;
887:     
888:     /* Get local mesh boundaries */
889:     DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
890:     DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 
891:     
892:     /* Get pointers to vector data */
893:     VecGetArray(user->Top,&top); 
894:     VecGetArray(user->Bottom,&bottom); 
895:     VecGetArray(user->Left,&left); 
896:     VecGetArray(user->Right,&right); 

898:     VecGetArray(localX,&x); 
899:     /* Perform local computations */    
900:     for (j=ys; j<ys+ym; j++){
901:       for (i=xs; i< xs+xm; i++){
902:         row=(j-gys)*gxm + (i-gxs);
903:         x[row] = ( (j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+
904:                    (i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0; 
905:       }
906:     }
907:     
908:     /* Restore vectors */
909:     VecRestoreArray(localX,&x); 

911:     VecRestoreArray(user->Left,&left); 
912:     VecRestoreArray(user->Top,&top); 
913:     VecRestoreArray(user->Bottom,&bottom); 
914:     VecRestoreArray(user->Right,&right); 
915:     
916:     /* Scatter values into global vector */
917:     DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X); 
918:     DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X); 
919:     
920:   }
921:   return 0;
922: }

924: /* For testing matrix free submatrices */
927: PetscErrorCode MatrixFreeHessian(TaoSolver tao, Vec x, Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
928: {
930:   AppCtx *user = (AppCtx*)ptr;
932:   FormHessian(tao,x,&user->H,&user->H,flg,ptr); 
933:   return(0);
934: }
937: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
938: {
940:   void *ptr;
941:   AppCtx *user;
943:   MatShellGetContext(H_shell,&ptr); 
944:   user = (AppCtx*)ptr;
945:   MatMult(user->H,X,Y); 
946:   return(0);
947: }