Actual source code: minsurf2.c

  1: /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */

  3: /*
  4:   Include "petsctao.h" so we can use TAO solvers.
  5:   petscdmda.h for distributed array
  6: */
  7: #include <petsctao.h>
  8: #include <petscdmda.h>

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

 22: /*T
 23:    Concepts: TAO^Solving an unconstrained minimization problem
 24:    Routines: TaoCreate(); TaoSetType();
 25:    Routines: TaoSetInitialVector();
 26:    Routines: TaoSetObjectiveAndGradientRoutine();
 27:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 28:    Routines: TaoSetMonitor();
 29:    Routines: TaoSolve(); TaoView();
 30:    Routines: TaoDestroy();
 31:    Processors: n
 32: T*/

 34: /*
 35:    User-defined application context - contains data needed by the
 36:    application-provided call-back routines, FormFunctionGradient()
 37:    and FormHessian().
 38: */
 39: typedef struct {
 40:   PetscInt    mx, my;                 /* discretization in x, y directions */
 41:   PetscReal   *bottom, *top, *left, *right;             /* boundary values */
 42:   DM          dm;                      /* distributed array data structure */
 43:   Mat         H;                       /* Hessian */
 44: } AppCtx;

 46: /* -------- User-defined Routines --------- */

 48: static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
 49: static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
 50: PetscErrorCode QuadraticH(AppCtx*,Vec,Mat);
 51: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void*);
 52: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
 53: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
 54: PetscErrorCode My_Monitor(Tao, void *);

 56: int main(int argc, char **argv)
 57: {
 58:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 59:   PetscInt           Nx, Ny;              /* number of processors in x- and y- directions */
 60:   Vec                x;                   /* solution, gradient vectors */
 61:   PetscBool          flg, viewmat;        /* flags */
 62:   PetscBool          fddefault, fdcoloring;   /* flags */
 63:   Tao                tao;                 /* TAO solver context */
 64:   AppCtx             user;                /* user-defined work context */
 65:   ISColoring         iscoloring;
 66:   MatFDColoring      matfdcoloring;

 68:   /* Initialize TAO */
 69:   PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;

 71:   /* Specify dimension of the problem */
 72:   user.mx = 10; user.my = 10;

 74:   /* Check for any command line arguments that override defaults */
 75:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);
 76:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);

 78:   PetscPrintf(MPI_COMM_WORLD,"\n---- Minimum Surface Area Problem -----\n");
 79:   PetscPrintf(MPI_COMM_WORLD,"mx: %D     my: %D   \n\n",user.mx,user.my);

 81:   /* Let PETSc determine the vector distribution */
 82:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 84:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 85:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx, user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
 86:   DMSetFromOptions(user.dm);
 87:   DMSetUp(user.dm);

 89:   /* Create TAO solver and set desired solution method.*/
 90:   TaoCreate(PETSC_COMM_WORLD,&tao);
 91:   TaoSetType(tao,TAOCG);

 93:   /*
 94:      Extract global vector from DA for the vector of variables --  PETSC routine
 95:      Compute the initial solution                              --  application specific, see below
 96:      Set this vector for use by TAO                            --  TAO routine
 97:   */
 98:   DMCreateGlobalVector(user.dm,&x);
 99:   MSA_BoundaryConditions(&user);
100:   MSA_InitialPoint(&user,x);
101:   TaoSetInitialVector(tao,x);

103:   /*
104:      Initialize the Application context for use in function evaluations  --  application specific, see below.
105:      Set routines for function and gradient evaluation
106:   */
107:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

109:   /*
110:      Given the command line arguments, calculate the hessian with either the user-
111:      provided function FormHessian, or the default finite-difference driven Hessian
112:      functions
113:   */
114:   PetscOptionsHasName(NULL,NULL,"-tao_fddefault",&fddefault);
115:   PetscOptionsHasName(NULL,NULL,"-tao_fdcoloring",&fdcoloring);

117:   /*
118:      Create a matrix data structure to store the Hessian and set
119:      the Hessian evalution routine.
120:      Set the matrix structure to be used for Hessian evalutions
121:   */
122:   DMCreateMatrix(user.dm,&user.H);
123:   MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);

125:   if (fdcoloring) {
126:     DMCreateColoring(user.dm,IS_COLORING_GLOBAL,&iscoloring);
127:     MatFDColoringCreate(user.H,iscoloring,&matfdcoloring);
128:     ISColoringDestroy(&iscoloring);
129:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormGradient,(void*)&user);
130:     MatFDColoringSetFromOptions(matfdcoloring);
131:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessianColor,(void *)matfdcoloring);
132:   } else if (fddefault) {
133:     TaoSetHessianRoutine(tao,user.H,user.H,TaoDefaultComputeHessian,(void *)NULL);
134:   } else {
135:     TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
136:   }

138:   /*
139:      If my_monitor option is in command line, then use the user-provided
140:      monitoring function
141:   */
142:   PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
143:   if (viewmat) {
144:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
145:   }

147:   /* Check for any tao command line options */
148:   TaoSetFromOptions(tao);

150:   /* SOLVE THE APPLICATION */
151:   TaoSolve(tao);

153:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

155:   /* Free TAO data structures */
156:   TaoDestroy(&tao);

158:   /* Free PETSc data structures */
159:   VecDestroy(&x);
160:   MatDestroy(&user.H);
161:   if (fdcoloring) {
162:     MatFDColoringDestroy(&matfdcoloring);
163:   }
164:   PetscFree(user.bottom);
165:   PetscFree(user.top);
166:   PetscFree(user.left);
167:   PetscFree(user.right);
168:   DMDestroy(&user.dm);
169:   PetscFinalize();
170:   return ierr;
171: }

173: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
174: {
176:   PetscReal      fcn;

179:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
180:   return(0);
181: }

183: /* -------------------------------------------------------------------- */
184: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

186:     Input Parameters:
187: .   tao     - the Tao context
188: .   XX      - input vector
189: .   userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

191:     Output Parameters:
192: .   fcn     - the newly evaluated function
193: .   GG       - vector containing the newly evaluated gradient
194: */
195: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx)
196: {
197:   AppCtx         *user = (AppCtx *) userCtx;
199:   PetscInt       i,j;
200:   PetscInt       mx=user->mx, my=user->my;
201:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
202:   PetscReal      ft=0.0;
203:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
204:   PetscReal      rhx=mx+1, rhy=my+1;
205:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
206:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
207:   PetscReal      **g, **x;
208:   Vec            localX;

211:   /* Get local mesh boundaries */
212:   DMGetLocalVector(user->dm,&localX);
213:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
214:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

216:   /* Scatter ghost points to local vector */
217:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
218:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

220:   /* Get pointers to vector data */
221:   DMDAVecGetArray(user->dm,localX,(void**)&x);
222:   DMDAVecGetArray(user->dm,G,(void**)&g);

224:   /* Compute function and gradient over the locally owned part of the mesh */
225:   for (j=ys; j<ys+ym; j++) {
226:     for (i=xs; i< xs+xm; i++) {

228:       xc = x[j][i];
229:       xlt=xrb=xl=xr=xb=xt=xc;

231:       if (i==0) { /* left side */
232:         xl= user->left[j-ys+1];
233:         xlt = user->left[j-ys+2];
234:       } else {
235:         xl = x[j][i-1];
236:       }

238:       if (j==0) { /* bottom side */
239:         xb=user->bottom[i-xs+1];
240:         xrb =user->bottom[i-xs+2];
241:       } else {
242:         xb = x[j-1][i];
243:       }

245:       if (i+1 == gxs+gxm) { /* right side */
246:         xr=user->right[j-ys+1];
247:         xrb = user->right[j-ys];
248:       } else {
249:         xr = x[j][i+1];
250:       }

252:       if (j+1==gys+gym) { /* top side */
253:         xt=user->top[i-xs+1];
254:         xlt = user->top[i-xs];
255:       }else {
256:         xt = x[j+1][i];
257:       }

259:       if (i>gxs && j+1<gys+gym) {
260:         xlt = x[j+1][i-1];
261:       }
262:       if (j>gys && i+1<gxs+gxm) {
263:         xrb = x[j-1][i+1];
264:       }

266:       d1 = (xc-xl);
267:       d2 = (xc-xr);
268:       d3 = (xc-xt);
269:       d4 = (xc-xb);
270:       d5 = (xr-xrb);
271:       d6 = (xrb-xb);
272:       d7 = (xlt-xl);
273:       d8 = (xt-xlt);

275:       df1dxc = d1*hydhx;
276:       df2dxc = (d1*hydhx + d4*hxdhy);
277:       df3dxc = d3*hxdhy;
278:       df4dxc = (d2*hydhx + d3*hxdhy);
279:       df5dxc = d2*hydhx;
280:       df6dxc = d4*hxdhy;

282:       d1 *= rhx;
283:       d2 *= rhx;
284:       d3 *= rhy;
285:       d4 *= rhy;
286:       d5 *= rhy;
287:       d6 *= rhx;
288:       d7 *= rhy;
289:       d8 *= rhx;

291:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
292:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
293:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
294:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
295:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
296:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

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

300:       df1dxc /= f1;
301:       df2dxc /= f2;
302:       df3dxc /= f3;
303:       df4dxc /= f4;
304:       df5dxc /= f5;
305:       df6dxc /= f6;

307:       g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;

309:     }
310:   }

312:   /* Compute triangular areas along the border of the domain. */
313:   if (xs==0) { /* left side */
314:     for (j=ys; j<ys+ym; j++) {
315:       d3=(user->left[j-ys+1] - user->left[j-ys+2])*rhy;
316:       d2=(user->left[j-ys+1] - x[j][0]) *rhx;
317:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
318:     }
319:   }
320:   if (ys==0) { /* bottom side */
321:     for (i=xs; i<xs+xm; i++) {
322:       d2=(user->bottom[i+1-xs]-user->bottom[i-xs+2])*rhx;
323:       d3=(user->bottom[i-xs+1]-x[0][i])*rhy;
324:       ft = ft+PetscSqrtReal(1.0 + d3*d3 + d2*d2);
325:     }
326:   }

328:   if (xs+xm==mx) { /* right side */
329:     for (j=ys; j< ys+ym; j++) {
330:       d1=(x[j][mx-1] - user->right[j-ys+1])*rhx;
331:       d4=(user->right[j-ys]-user->right[j-ys+1])*rhy;
332:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
333:     }
334:   }
335:   if (ys+ym==my) { /* top side */
336:     for (i=xs; i<xs+xm; i++) {
337:       d1=(x[my-1][i] - user->top[i-xs+1])*rhy;
338:       d4=(user->top[i-xs+1] - user->top[i-xs])*rhx;
339:       ft = ft+PetscSqrtReal(1.0 + d1*d1 + d4*d4);
340:     }
341:   }

343:   if (ys==0 && xs==0) {
344:     d1=(user->left[0]-user->left[1])/hy;
345:     d2=(user->bottom[0]-user->bottom[1])*rhx;
346:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
347:   }
348:   if (ys+ym == my && xs+xm == mx) {
349:     d1=(user->right[ym+1] - user->right[ym])*rhy;
350:     d2=(user->top[xm+1] - user->top[xm])*rhx;
351:     ft +=PetscSqrtReal(1.0 + d1*d1 + d2*d2);
352:   }

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

357:   /* Restore vectors */
358:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
359:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

361:   /* Scatter values to global vector */
362:   DMRestoreLocalVector(user->dm,&localX);
363:   PetscLogFlops(67.0*xm*ym);
364:   return(0);
365: }

367: /* ------------------------------------------------------------------- */
368: /*
369:    FormHessian - Evaluates Hessian matrix.

371:    Input Parameters:
372: .  tao  - the Tao context
373: .  x    - input vector
374: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

376:    Output Parameters:
377: .  H    - Hessian matrix
378: .  Hpre - optionally different preconditioning matrix
379: .  flg  - flag indicating matrix structure

381: */
382: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
383: {
385:   AppCtx         *user = (AppCtx *) ptr;

388:   /* Evaluate the Hessian entries*/
389:   QuadraticH(user,X,H);
390:   return(0);
391: }

393: /* ------------------------------------------------------------------- */
394: /*
395:    QuadraticH - Evaluates Hessian matrix.

397:    Input Parameters:
398: .  user - user-defined context, as set by TaoSetHessian()
399: .  X    - input vector

401:    Output Parameter:
402: .  H    - Hessian matrix
403: */
404: PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
405: {
406:   PetscErrorCode    ierr;
407:   PetscInt i,j,k;
408:   PetscInt mx=user->mx, my=user->my;
409:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
410:   PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
411:   PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
412:   PetscReal hl,hr,ht,hb,hc,htl,hbr;
413:   PetscReal **x, v[7];
414:   MatStencil col[7],row;
415:   Vec    localX;
416:   PetscBool assembled;

419:   /* Get local mesh boundaries */
420:   DMGetLocalVector(user->dm,&localX);

422:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
423:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

425:   /* Scatter ghost points to local vector */
426:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
427:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

429:   /* Get pointers to vector data */
430:   DMDAVecGetArray(user->dm,localX,(void**)&x);

432:   /* Initialize matrix entries to zero */
433:   MatAssembled(Hessian,&assembled);
434:   if (assembled) {MatZeroEntries(Hessian);}

436:   /* Set various matrix options */
437:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

445:       xc = x[j][i];
446:       xlt=xrb=xl=xr=xb=xt=xc;

448:       /* Left side */
449:       if (i==0) {
450:         xl  = user->left[j-ys+1];
451:         xlt = user->left[j-ys+2];
452:       } else {
453:         xl  = x[j][i-1];
454:       }

456:       if (j==0) {
457:         xb  = user->bottom[i-xs+1];
458:         xrb = user->bottom[i-xs+2];
459:       } else {
460:         xb  = x[j-1][i];
461:       }

463:       if (i+1 == mx) {
464:         xr  = user->right[j-ys+1];
465:         xrb = user->right[j-ys];
466:       } else {
467:         xr  = x[j][i+1];
468:       }

470:       if (j+1==my) {
471:         xt  = user->top[i-xs+1];
472:         xlt = user->top[i-xs];
473:       }else {
474:         xt  = x[j+1][i];
475:       }

477:       if (i>0 && j+1<my) {
478:         xlt = x[j+1][i-1];
479:       }
480:       if (j>0 && i+1<mx) {
481:         xrb = x[j-1][i+1];
482:       }

484:       d1 = (xc-xl)/hx;
485:       d2 = (xc-xr)/hx;
486:       d3 = (xc-xt)/hy;
487:       d4 = (xc-xb)/hy;
488:       d5 = (xrb-xr)/hy;
489:       d6 = (xrb-xb)/hx;
490:       d7 = (xlt-xl)/hy;
491:       d8 = (xlt-xt)/hx;

493:       f1 = PetscSqrtReal(1.0 + d1*d1 + d7*d7);
494:       f2 = PetscSqrtReal(1.0 + d1*d1 + d4*d4);
495:       f3 = PetscSqrtReal(1.0 + d3*d3 + d8*d8);
496:       f4 = PetscSqrtReal(1.0 + d3*d3 + d2*d2);
497:       f5 = PetscSqrtReal(1.0 + d2*d2 + d5*d5);
498:       f6 = PetscSqrtReal(1.0 + d4*d4 + d6*d6);

500:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
501:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
502:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
503:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
504:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
505:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
506:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
507:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

512:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
513:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
514:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
515:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

517:       hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0;  hc/=2.0;

519:       row.j = j; row.i = i;
520:       k=0;
521:       if (j>0) {
522:         v[k]=hb;
523:         col[k].j = j - 1; col[k].i = i;
524:         k++;
525:       }

527:       if (j>0 && i < mx -1) {
528:         v[k]=hbr;
529:         col[k].j = j - 1; col[k].i = i+1;
530:         k++;
531:       }

533:       if (i>0) {
534:         v[k]= hl;
535:         col[k].j = j; col[k].i = i-1;
536:         k++;
537:       }

539:       v[k]= hc;
540:       col[k].j = j; col[k].i = i;
541:       k++;

543:       if (i < mx-1) {
544:         v[k]= hr;
545:         col[k].j = j; col[k].i = i+1;
546:         k++;
547:       }

549:       if (i>0 && j < my-1) {
550:         v[k]= htl;
551:         col[k].j = j+1; col[k].i = i-1;
552:         k++;
553:       }

555:       if (j < my-1) {
556:         v[k]= ht;
557:         col[k].j = j+1; col[k].i = i;
558:         k++;
559:       }

561:       /*
562:          Set matrix values using local numbering, which was defined
563:          earlier, in the main routine.
564:       */
565:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
566:     }
567:   }

569:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
570:   DMRestoreLocalVector(user->dm,&localX);

572:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
573:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

575:   PetscLogFlops(199.0*xm*ym);
576:   return(0);
577: }

579: /* ------------------------------------------------------------------- */
580: /*
581:    MSA_BoundaryConditions -  Calculates the boundary conditions for
582:    the region.

584:    Input Parameter:
585: .  user - user-defined application context

587:    Output Parameter:
588: .  user - user-defined application context
589: */
590: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
591: {
593:   PetscInt       i,j,k,limit=0,maxits=5;
594:   PetscInt       xs,ys,xm,ym,gxs,gys,gxm,gym;
595:   PetscInt       mx=user->mx,my=user->my;
596:   PetscInt       bsize=0, lsize=0, tsize=0, rsize=0;
597:   PetscReal      one=1.0, two=2.0, three=3.0, tol=1e-10;
598:   PetscReal      fnorm,det,hx,hy,xt=0,yt=0;
599:   PetscReal      u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
600:   PetscReal      b=-0.5, t=0.5, l=-0.5, r=0.5;
601:   PetscReal      *boundary;
602:   PetscBool      flg;

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

609:   bsize=xm+2;
610:   lsize=ym+2;
611:   rsize=ym+2;
612:   tsize=xm+2;

614:   PetscMalloc1(bsize,&user->bottom);
615:   PetscMalloc1(tsize,&user->top);
616:   PetscMalloc1(lsize,&user->left);
617:   PetscMalloc1(rsize,&user->right);

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

621:   for (j=0; j<4; j++) {
622:     if (j==0) {
623:       yt=b;
624:       xt=l+hx*xs;
625:       limit=bsize;
626:       boundary=user->bottom;
627:     } else if (j==1) {
628:       yt=t;
629:       xt=l+hx*xs;
630:       limit=tsize;
631:       boundary=user->top;
632:     } else if (j==2) {
633:       yt=b+hy*ys;
634:       xt=l;
635:       limit=lsize;
636:       boundary=user->left;
637:     } else { /* if (j==3) */
638:       yt=b+hy*ys;
639:       xt=r;
640:       limit=rsize;
641:       boundary=user->right;
642:     }

644:     for (i=0; i<limit; i++) {
645:       u1=xt;
646:       u2=-yt;
647:       for (k=0; k<maxits; k++) {
648:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
649:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
650:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
651:         if (fnorm <= tol) break;
652:         njac11=one+u2*u2-u1*u1;
653:         njac12=two*u1*u2;
654:         njac21=-two*u1*u2;
655:         njac22=-one - u1*u1 + u2*u2;
656:         det = njac11*njac22-njac21*njac12;
657:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
658:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
659:       }

661:       boundary[i]=u1*u1-u2*u2;
662:       if (j==0 || j==1) {
663:         xt=xt+hx;
664:       } else { /*  if (j==2 || j==3) */
665:         yt=yt+hy;
666:       }

668:     }

670:   }

672:   /* Scale the boundary if desired */
673:   if (1==1) {
674:     PetscReal scl = 1.0;

676:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
677:     if (flg) {
678:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
679:     }

681:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
682:     if (flg) {
683:       for (i=0;i<tsize;i++) user->top[i]*=scl;
684:     }

686:     PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
687:     if (flg) {
688:       for (i=0;i<rsize;i++) user->right[i]*=scl;
689:     }

691:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
692:     if (flg) {
693:       for (i=0;i<lsize;i++) user->left[i]*=scl;
694:     }
695:   }
696:   return(0);
697: }

699: /* ------------------------------------------------------------------- */
700: /*
701:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

703:    Input Parameters:
704: .  user - user-defined application context
705: .  X - vector for initial guess

707:    Output Parameters:
708: .  X - newly computed initial guess
709: */
710: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
711: {
713:   PetscInt       start2=-1,i,j;
714:   PetscReal      start1=0;
715:   PetscBool      flg1,flg2;

718:   PetscOptionsGetReal(NULL,NULL,"-start",&start1,&flg1);
719:   PetscOptionsGetInt(NULL,NULL,"-random",&start2,&flg2);

721:   if (flg1) { /* The zero vector is reasonable */

723:     VecSet(X, start1);

725:   } else if (flg2 && start2>0) { /* Try a random start between -0.5 and 0.5 */
726:     PetscRandom rctx;  PetscReal np5=-0.5;

728:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
729:     VecSetRandom(X, rctx);
730:     PetscRandomDestroy(&rctx);
731:     VecShift(X, np5);

733:   } else { /* Take an average of the boundary conditions */
734:     PetscInt  xs,xm,ys,ym;
735:     PetscInt  mx=user->mx,my=user->my;
736:     PetscReal **x;

738:     /* Get local mesh boundaries */
739:     DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);

741:     /* Get pointers to vector data */
742:     DMDAVecGetArray(user->dm,X,(void**)&x);

744:     /* Perform local computations */
745:     for (j=ys; j<ys+ym; j++) {
746:       for (i=xs; i< xs+xm; i++) {
747:         x[j][i] = (((j+1)*user->bottom[i-xs+1]+(my-j+1)*user->top[i-xs+1])/(my+2)+((i+1)*user->left[j-ys+1]+(mx-i+1)*user->right[j-ys+1])/(mx+2))/2.0;
748:       }
749:     }
750:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
751:     PetscLogFlops(9.0*xm*ym);
752:   }
753:   return(0);
754: }

756: /*-----------------------------------------------------------------------*/
757: PetscErrorCode My_Monitor(Tao tao, void *ctx)
758: {
760:   Vec            X;

763:   TaoGetSolutionVector(tao,&X);
764:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
765:   return(0);
766: }

768: /*TEST

770:    build:
771:       requires: !complex

773:    test:
774:       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
775:       requires: !single

777:    test:
778:       suffix: 2
779:       nsize: 2
780:       args: -tao_smonitor -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
781:       filter: grep -v "nls ksp"
782:       requires: !single

784:    test:
785:       suffix: 3
786:       nsize: 3
787:       args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gatol 1.e-3
788:       requires: !single

790:    test:
791:       suffix: 5
792:       nsize: 2
793:       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
794:       requires: !single

796: TEST*/