Actual source code: minsurf2.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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*/



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


 49: /* -------- User-defined Routines --------- */

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

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

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

 74:   /* Specify dimension of the problem */
 75:   user.mx = 10; user.my = 10;

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

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


 85:   /* Let PETSc determine the vector distribution */
 86:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 88:   /* Create distributed array (DM) to manage parallel grid and vectors  */
 89:   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);
 90:   DMSetFromOptions(user.dm);
 91:   DMSetUp(user.dm);

 93:   /* Create TAO solver and set desired solution method.*/
 94:   TaoCreate(PETSC_COMM_WORLD,&tao);
 95:   TaoSetType(tao,TAOCG);

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

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

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

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


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

143:   /*
144:      If my_monitor option is in command line, then use the user-provided
145:      monitoring function
146:   */
147:   PetscOptionsHasName(NULL,NULL,"-my_monitor",&viewmat);
148:   if (viewmat){
149:     TaoSetMonitor(tao,My_Monitor,NULL,NULL);
150:   }

152:   /* Check for any tao command line options */
153:   TaoSetFromOptions(tao);

155:   /* SOLVE THE APPLICATION */
156:   TaoSolve(tao);

158:   TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);

160:   /* Free TAO data structures */
161:   TaoDestroy(&tao);

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

178: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G,void *userCtx)
179: {
181:   PetscReal      fcn;

184:   FormFunctionGradient(tao,X,&fcn,G,userCtx);
185:   return(0);
186: }

188: /* -------------------------------------------------------------------- */
189: /*  FormFunctionGradient - Evaluates the function and corresponding gradient.

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

196:     Output Parameters:
197: .   fcn     - the newly evaluated function
198: .   GG       - vector containing the newly evaluated gradient
199: */
200: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *userCtx){

202:   AppCtx         *user = (AppCtx *) userCtx;
204:   PetscInt       i,j;
205:   PetscInt       mx=user->mx, my=user->my;
206:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
207:   PetscReal      ft=0.0;
208:   PetscReal      hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
209:   PetscReal      rhx=mx+1, rhy=my+1;
210:   PetscReal      f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
211:   PetscReal      df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
212:   PetscReal      **g, **x;
213:   Vec            localX;

216:   /* Get local mesh boundaries */
217:   DMGetLocalVector(user->dm,&localX);
218:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
219:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

221:   /* Scatter ghost points to local vector */
222:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
223:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

225:   /* Get pointers to vector data */
226:   DMDAVecGetArray(user->dm,localX,(void**)&x);
227:   DMDAVecGetArray(user->dm,G,(void**)&g);

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

233:       xc = x[j][i];
234:       xlt=xrb=xl=xr=xb=xt=xc;

236:       if (i==0){ /* left side */
237:         xl= user->left[j-ys+1];
238:         xlt = user->left[j-ys+2];
239:       } else {
240:         xl = x[j][i-1];
241:       }

243:       if (j==0){ /* bottom side */
244:         xb=user->bottom[i-xs+1];
245:         xrb =user->bottom[i-xs+2];
246:       } else {
247:         xb = x[j-1][i];
248:       }

250:       if (i+1 == gxs+gxm){ /* right side */
251:         xr=user->right[j-ys+1];
252:         xrb = user->right[j-ys];
253:       } else {
254:         xr = x[j][i+1];
255:       }

257:       if (j+1==gys+gym){ /* top side */
258:         xt=user->top[i-xs+1];
259:         xlt = user->top[i-xs];
260:       }else {
261:         xt = x[j+1][i];
262:       }

264:       if (i>gxs && j+1<gys+gym){
265:         xlt = x[j+1][i-1];
266:       }
267:       if (j>gys && i+1<gxs+gxm){
268:         xrb = x[j-1][i+1];
269:       }

271:       d1 = (xc-xl);
272:       d2 = (xc-xr);
273:       d3 = (xc-xt);
274:       d4 = (xc-xb);
275:       d5 = (xr-xrb);
276:       d6 = (xrb-xb);
277:       d7 = (xlt-xl);
278:       d8 = (xt-xlt);

280:       df1dxc = d1*hydhx;
281:       df2dxc = ( d1*hydhx + d4*hxdhy );
282:       df3dxc = d3*hxdhy;
283:       df4dxc = ( d2*hydhx + d3*hxdhy );
284:       df5dxc = d2*hydhx;
285:       df6dxc = d4*hxdhy;

287:       d1 *= rhx;
288:       d2 *= rhx;
289:       d3 *= rhy;
290:       d4 *= rhy;
291:       d5 *= rhy;
292:       d6 *= rhx;
293:       d7 *= rhy;
294:       d8 *= rhx;

296:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
297:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
298:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
299:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
300:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
301:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);

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

305:       df1dxc /= f1;
306:       df2dxc /= f2;
307:       df3dxc /= f3;
308:       df4dxc /= f4;
309:       df5dxc /= f5;
310:       df6dxc /= f6;

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

314:     }
315:   }

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

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

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

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

362:   /* Restore vectors */
363:   DMDAVecRestoreArray(user->dm,localX,(void**)&x);
364:   DMDAVecRestoreArray(user->dm,G,(void**)&g);

366:   /* Scatter values to global vector */
367:   DMRestoreLocalVector(user->dm,&localX);
368:   PetscLogFlops(67*xm*ym);
369:   return(0);
370: }

372: /* ------------------------------------------------------------------- */
373: /*
374:    FormHessian - Evaluates Hessian matrix.

376:    Input Parameters:
377: .  tao  - the Tao context
378: .  x    - input vector
379: .  ptr  - optional user-defined context, as set by TaoSetHessianRoutine()

381:    Output Parameters:
382: .  H    - Hessian matrix
383: .  Hpre - optionally different preconditioning matrix
384: .  flg  - flag indicating matrix structure

386: */
387: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H, Mat Hpre, void *ptr)
388: {
390:   AppCtx         *user = (AppCtx *) ptr;

393:   /* Evaluate the Hessian entries*/
394:   QuadraticH(user,X,H);
395:   return(0);
396: }

398: /* ------------------------------------------------------------------- */
399: /*
400:    QuadraticH - Evaluates Hessian matrix.

402:    Input Parameters:
403: .  user - user-defined context, as set by TaoSetHessian()
404: .  X    - input vector

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

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

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

430:   /* Scatter ghost points to local vector */
431:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
432:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

437:   /* Initialize matrix entries to zero */
438:   MatAssembled(Hessian,&assembled);
439:   if (assembled){MatZeroEntries(Hessian);}


442:   /* Set various matrix options */
443:   MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);

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

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

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

451:       xc = x[j][i];
452:       xlt=xrb=xl=xr=xb=xt=xc;

454:       /* Left side */
455:       if (i==0){
456:         xl  = user->left[j-ys+1];
457:         xlt = user->left[j-ys+2];
458:       } else {
459:         xl  = x[j][i-1];
460:       }

462:       if (j==0){
463:         xb  = user->bottom[i-xs+1];
464:         xrb = user->bottom[i-xs+2];
465:       } else {
466:         xb  = x[j-1][i];
467:       }

469:       if (i+1 == mx){
470:         xr  = user->right[j-ys+1];
471:         xrb = user->right[j-ys];
472:       } else {
473:         xr  = x[j][i+1];
474:       }

476:       if (j+1==my){
477:         xt  = user->top[i-xs+1];
478:         xlt = user->top[i-xs];
479:       }else {
480:         xt  = x[j+1][i];
481:       }

483:       if (i>0 && j+1<my){
484:         xlt = x[j+1][i-1];
485:       }
486:       if (j>0 && i+1<mx){
487:         xrb = x[j-1][i+1];
488:       }


491:       d1 = (xc-xl)/hx;
492:       d2 = (xc-xr)/hx;
493:       d3 = (xc-xt)/hy;
494:       d4 = (xc-xb)/hy;
495:       d5 = (xrb-xr)/hy;
496:       d6 = (xrb-xb)/hx;
497:       d7 = (xlt-xl)/hy;
498:       d8 = (xlt-xt)/hx;

500:       f1 = PetscSqrtReal( 1.0 + d1*d1 + d7*d7);
501:       f2 = PetscSqrtReal( 1.0 + d1*d1 + d4*d4);
502:       f3 = PetscSqrtReal( 1.0 + d3*d3 + d8*d8);
503:       f4 = PetscSqrtReal( 1.0 + d3*d3 + d2*d2);
504:       f5 = PetscSqrtReal( 1.0 + d2*d2 + d5*d5);
505:       f6 = PetscSqrtReal( 1.0 + d4*d4 + d6*d6);


508:       hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
509:         (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
510:       hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
511:         (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
512:       ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
513:         (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
514:       hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
515:         (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);

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

520:       hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
521:         hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
522:         (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
523:         (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);

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

527:       row.j = j; row.i = i;
528:       k=0;
529:       if (j>0){
530:         v[k]=hb;
531:         col[k].j = j - 1; col[k].i = i;
532:         k++;
533:       }

535:       if (j>0 && i < mx -1){
536:         v[k]=hbr;
537:         col[k].j = j - 1; col[k].i = i+1;
538:         k++;
539:       }

541:       if (i>0){
542:         v[k]= hl;
543:         col[k].j = j; col[k].i = i-1;
544:         k++;
545:       }

547:       v[k]= hc;
548:       col[k].j = j; col[k].i = i;
549:       k++;

551:       if (i < mx-1 ){
552:         v[k]= hr;
553:         col[k].j = j; col[k].i = i+1;
554:         k++;
555:       }

557:       if (i>0 && j < my-1 ){
558:         v[k]= htl;
559:         col[k].j = j+1; col[k].i = i-1;
560:         k++;
561:       }

563:       if (j < my-1 ){
564:         v[k]= ht;
565:         col[k].j = j+1; col[k].i = i;
566:         k++;
567:       }

569:       /*
570:          Set matrix values using local numbering, which was defined
571:          earlier, in the main routine.
572:       */
573:       MatSetValuesStencil(Hessian,1,&row,k,col,v,INSERT_VALUES);
574:     }
575:   }

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

580:   MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);
581:   MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);

583:   PetscLogFlops(199*xm*ym);
584:   return(0);
585: }

587: /* ------------------------------------------------------------------- */
588: /*
589:    MSA_BoundaryConditions -  Calculates the boundary conditions for
590:    the region.

592:    Input Parameter:
593: .  user - user-defined application context

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

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

617:   bsize=xm+2;
618:   lsize=ym+2;
619:   rsize=ym+2;
620:   tsize=xm+2;

622:   PetscMalloc1(bsize,&user->bottom);
623:   PetscMalloc1(tsize,&user->top);
624:   PetscMalloc1(lsize,&user->left);
625:   PetscMalloc1(rsize,&user->right);

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

629:   for (j=0; j<4; j++){
630:     if (j==0){
631:       yt=b;
632:       xt=l+hx*xs;
633:       limit=bsize;
634:       boundary=user->bottom;
635:     } else if (j==1){
636:       yt=t;
637:       xt=l+hx*xs;
638:       limit=tsize;
639:       boundary=user->top;
640:     } else if (j==2){
641:       yt=b+hy*ys;
642:       xt=l;
643:       limit=lsize;
644:       boundary=user->left;
645:     } else { /* if (j==3) */
646:       yt=b+hy*ys;
647:       xt=r;
648:       limit=rsize;
649:       boundary=user->right;
650:     }

652:     for (i=0; i<limit; i++){
653:       u1=xt;
654:       u2=-yt;
655:       for (k=0; k<maxits; k++){
656:         nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
657:         nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
658:         fnorm=PetscSqrtReal(nf1*nf1+nf2*nf2);
659:         if (fnorm <= tol) break;
660:         njac11=one+u2*u2-u1*u1;
661:         njac12=two*u1*u2;
662:         njac21=-two*u1*u2;
663:         njac22=-one - u1*u1 + u2*u2;
664:         det = njac11*njac22-njac21*njac12;
665:         u1 = u1-(njac22*nf1-njac12*nf2)/det;
666:         u2 = u2-(njac11*nf2-njac21*nf1)/det;
667:       }

669:       boundary[i]=u1*u1-u2*u2;
670:       if (j==0 || j==1) {
671:         xt=xt+hx;
672:       } else { /*  if (j==2 || j==3) */
673:         yt=yt+hy;
674:       }

676:     }

678:   }

680:   /* Scale the boundary if desired */
681:   if (1==1){
682:     PetscReal scl = 1.0;

684:     PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);
685:     if (flg){
686:       for (i=0;i<bsize;i++) user->bottom[i]*=scl;
687:     }

689:     PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);
690:     if (flg){
691:       for (i=0;i<tsize;i++) user->top[i]*=scl;
692:     }

694:     PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);
695:     if (flg){
696:       for (i=0;i<rsize;i++) user->right[i]*=scl;
697:     }

699:     PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);
700:     if (flg){
701:       for (i=0;i<lsize;i++) user->left[i]*=scl;
702:     }
703:   }
704:   return(0);
705: }

707: /* ------------------------------------------------------------------- */
708: /*
709:    MSA_InitialPoint - Calculates the initial guess in one of three ways.

711:    Input Parameters:
712: .  user - user-defined application context
713: .  X - vector for initial guess

715:    Output Parameters:
716: .  X - newly computed initial guess
717: */
718: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
719: {
721:   PetscInt       start2=-1,i,j;
722:   PetscReal      start1=0;
723:   PetscBool      flg1,flg2;

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

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

731:     VecSet(X, start1);

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

736:     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
737:     VecSetRandom(X, rctx);
738:     PetscRandomDestroy(&rctx);
739:     VecShift(X, np5);

741:   } else { /* Take an average of the boundary conditions */
742:     PetscInt  xs,xm,ys,ym;
743:     PetscInt  mx=user->mx,my=user->my;
744:     PetscReal **x;

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

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

752:     /* Perform local computations */
753:     for (j=ys; j<ys+ym; j++){
754:       for (i=xs; i< xs+xm; i++){
755:         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;
756:       }
757:     }
758:     DMDAVecRestoreArray(user->dm,X,(void**)&x);
759:     PetscLogFlops(9*xm*ym);
760:   }
761:   return(0);
762: }

764: /*-----------------------------------------------------------------------*/
765: PetscErrorCode My_Monitor(Tao tao, void *ctx)
766: {
768:   Vec            X;

771:   TaoGetSolutionVector(tao,&X);
772:   VecView(X,PETSC_VIEWER_STDOUT_WORLD);
773:   return(0);
774: }


777: /*TEST

779:    build:
780:       requires: !complex

782:    test:
783:       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gttol 1.e-2
784:       requires: !single

786:    test:
787:       suffix: 2
788:       nsize: 2
789:       args: -tao_smonitor -tao_type nls -ksp_max_it 15 -tao_gttol 1.e-2
790:       filter: grep -v "nls ksp"
791:       requires: !single

793:    test:
794:       suffix: 3
795:       nsize: 3
796:       args: -tao_smonitor -tao_type cg -tao_cg_type fr -mx 10 -my 10 -tao_gttol 1.e-2
797:       requires: !single

799:    test:
800:       suffix: 5
801:       nsize: 2
802:       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gttol 1.e-2
803:       requires: !single

805: TEST*/