Actual source code: jbearing2.c

  1: /*
  2:   Include "petsctao.h" so we can use TAO solvers
  3:   Include "petscdmda.h" so that we can use distributed arrays (DMs) for managing
  4:   Include "petscksp.h" so we can set KSP type
  5:   the parallel mesh.
  6: */

  8: #include <petsctao.h>
  9: #include <petscdmda.h>

 11: static char help[] = "This example demonstrates use of the TAO package to \n\
 12: solve a bound constrained minimization problem.  This example is based on \n\
 13: the problem DPJB from the MINPACK-2 test suite.  This pressure journal \n\
 14: bearing problem is an example of elliptic variational problem defined over \n\
 15: a two dimensional rectangle.  By discretizing the domain into triangular \n\
 16: elements, the pressure surrounding the journal bearing is defined as the \n\
 17: minimum of a quadratic function whose variables are bounded below by zero.\n\
 18: The command line options are:\n\
 19:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 20:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 21:  \n";

 23: /*
 24:    User-defined application context - contains data needed by the
 25:    application-provided call-back routines, FormFunctionGradient(),
 26:    FormHessian().
 27: */
 28: typedef struct {
 29:   /* problem parameters */
 30:   PetscReal ecc;    /* test problem parameter */
 31:   PetscReal b;      /* A dimension of journal bearing */
 32:   PetscInt  nx, ny; /* discretization in x, y directions */

 34:   /* Working space */
 35:   DM  dm; /* distributed array data structure */
 36:   Mat A;  /* Quadratic Objective term */
 37:   Vec B;  /* Linear Objective term */
 38: } AppCtx;

 40: /* User-defined routines */
 41: static PetscReal      p(PetscReal xi, PetscReal ecc);
 42: static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 43: static PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 44: static PetscErrorCode ComputeB(AppCtx *);
 45: static PetscErrorCode Monitor(Tao, void *);
 46: static PetscErrorCode ConvergenceTest(Tao, void *);

 48: int main(int argc, char **argv)
 49: {
 50:   PetscInt  Nx, Ny; /* number of processors in x- and y- directions */
 51:   PetscInt  m;      /* number of local elements in vectors */
 52:   Vec       x;      /* variables vector */
 53:   Vec       xl, xu; /* bounds vectors */
 54:   PetscReal d1000 = 1000;
 55:   PetscBool flg, testgetdiag; /* A return variable when checking for user options */
 56:   Tao       tao;              /* Tao solver context */
 57:   KSP       ksp;
 58:   AppCtx    user;       /* user-defined work context */
 59:   PetscReal zero = 0.0; /* lower bound on all variables */

 61:   /* Initialize PETSC and TAO */
 63:   PetscInitialize(&argc, &argv, (char *)0, help);

 65:   /* Set the default values for the problem parameters */
 66:   user.nx     = 50;
 67:   user.ny     = 50;
 68:   user.ecc    = 0.1;
 69:   user.b      = 10.0;
 70:   testgetdiag = PETSC_FALSE;

 72:   /* Check for any command line arguments that override defaults */
 73:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.nx, &flg);
 74:   PetscOptionsGetInt(NULL, NULL, "-my", &user.ny, &flg);
 75:   PetscOptionsGetReal(NULL, NULL, "-ecc", &user.ecc, &flg);
 76:   PetscOptionsGetReal(NULL, NULL, "-b", &user.b, &flg);
 77:   PetscOptionsGetBool(NULL, NULL, "-test_getdiagonal", &testgetdiag, NULL);

 79:   PetscPrintf(PETSC_COMM_WORLD, "\n---- Journal Bearing Problem SHB-----\n");
 80:   PetscPrintf(PETSC_COMM_WORLD, "mx: %" PetscInt_FMT ",  my: %" PetscInt_FMT ",  ecc: %g \n\n", user.nx, user.ny, (double)user.ecc);

 82:   /* Let Petsc determine the grid division */
 83:   Nx = PETSC_DECIDE;
 84:   Ny = PETSC_DECIDE;

 86:   /*
 87:      A two dimensional distributed array will help define this problem,
 88:      which derives from an elliptic PDE on two dimensional domain.  From
 89:      the distributed array, Create the vectors.
 90:   */
 91:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.nx, user.ny, Nx, Ny, 1, 1, NULL, NULL, &user.dm);
 92:   DMSetFromOptions(user.dm);
 93:   DMSetUp(user.dm);

 95:   /*
 96:      Extract global and local vectors from DM; the vector user.B is
 97:      used solely as work space for the evaluation of the function,
 98:      gradient, and Hessian.  Duplicate for remaining vectors that are
 99:      the same types.
100:   */
101:   DMCreateGlobalVector(user.dm, &x); /* Solution */
102:   VecDuplicate(x, &user.B);          /* Linear objective */

104:   /*  Create matrix user.A to store quadratic, Create a local ordering scheme. */
105:   VecGetLocalSize(x, &m);
106:   DMCreateMatrix(user.dm, &user.A);

108:   if (testgetdiag) MatSetOperation(user.A, MATOP_GET_DIAGONAL, NULL);

110:   /* User defined function -- compute linear term of quadratic */
111:   ComputeB(&user);

113:   /* The TAO code begins here */

115:   /*
116:      Create the optimization solver
117:      Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
118:   */
119:   TaoCreate(PETSC_COMM_WORLD, &tao);
120:   TaoSetType(tao, TAOBLMVM);

122:   /* Set the initial vector */
123:   VecSet(x, zero);
124:   TaoSetSolution(tao, x);

126:   /* Set the user function, gradient, hessian evaluation routines and data structures */
127:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);

129:   TaoSetHessian(tao, user.A, user.A, FormHessian, (void *)&user);

131:   /* Set a routine that defines the bounds */
132:   VecDuplicate(x, &xl);
133:   VecDuplicate(x, &xu);
134:   VecSet(xl, zero);
135:   VecSet(xu, d1000);
136:   TaoSetVariableBounds(tao, xl, xu);

138:   TaoGetKSP(tao, &ksp);
139:   if (ksp) KSPSetType(ksp, KSPCG);

141:   PetscOptionsHasName(NULL, NULL, "-testmonitor", &flg);
142:   if (flg) TaoSetMonitor(tao, Monitor, &user, NULL);
143:   PetscOptionsHasName(NULL, NULL, "-testconvergence", &flg);
144:   if (flg) TaoSetConvergenceTest(tao, ConvergenceTest, &user);

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

149:   /* Solve the bound constrained problem */
150:   TaoSolve(tao);

152:   /* Free PETSc data structures */
153:   VecDestroy(&x);
154:   VecDestroy(&xl);
155:   VecDestroy(&xu);
156:   MatDestroy(&user.A);
157:   VecDestroy(&user.B);

159:   /* Free TAO data structures */
160:   TaoDestroy(&tao);
161:   DMDestroy(&user.dm);
162:   PetscFinalize();
163:   return 0;
164: }

166: static PetscReal p(PetscReal xi, PetscReal ecc)
167: {
168:   PetscReal t = 1.0 + ecc * PetscCosScalar(xi);
169:   return (t * t * t);
170: }

172: PetscErrorCode ComputeB(AppCtx *user)
173: {
174:   PetscInt  i, j, k;
175:   PetscInt  nx, ny, xs, xm, gxs, gxm, ys, ym, gys, gym;
176:   PetscReal two = 2.0, pi = 4.0 * atan(1.0);
177:   PetscReal hx, hy, ehxhy;
178:   PetscReal temp, *b;
179:   PetscReal ecc = user->ecc;

181:   nx    = user->nx;
182:   ny    = user->ny;
183:   hx    = two * pi / (nx + 1.0);
184:   hy    = two * user->b / (ny + 1.0);
185:   ehxhy = ecc * hx * hy;

187:   /*
188:      Get local grid boundaries
189:   */
190:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
191:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

193:   /* Compute the linear term in the objective function */
194:   VecGetArray(user->B, &b);
195:   for (i = xs; i < xs + xm; i++) {
196:     temp = PetscSinScalar((i + 1) * hx);
197:     for (j = ys; j < ys + ym; j++) {
198:       k    = xm * (j - ys) + (i - xs);
199:       b[k] = -ehxhy * temp;
200:     }
201:   }
202:   VecRestoreArray(user->B, &b);
203:   PetscLogFlops(5.0 * xm * ym + 3.0 * xm);
204:   return 0;
205: }

207: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr)
208: {
209:   AppCtx    *user = (AppCtx *)ptr;
210:   PetscInt   i, j, k, kk;
211:   PetscInt   col[5], row, nx, ny, xs, xm, gxs, gxm, ys, ym, gys, gym;
212:   PetscReal  one = 1.0, two = 2.0, six = 6.0, pi = 4.0 * atan(1.0);
213:   PetscReal  hx, hy, hxhy, hxhx, hyhy;
214:   PetscReal  xi, v[5];
215:   PetscReal  ecc = user->ecc, trule1, trule2, trule3, trule4, trule5, trule6;
216:   PetscReal  vmiddle, vup, vdown, vleft, vright;
217:   PetscReal  tt, f1, f2;
218:   PetscReal *x, *g, zero = 0.0;
219:   Vec        localX;

221:   nx   = user->nx;
222:   ny   = user->ny;
223:   hx   = two * pi / (nx + 1.0);
224:   hy   = two * user->b / (ny + 1.0);
225:   hxhy = hx * hy;
226:   hxhx = one / (hx * hx);
227:   hyhy = one / (hy * hy);

229:   DMGetLocalVector(user->dm, &localX);

231:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
232:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);

234:   VecSet(G, zero);
235:   /*
236:     Get local grid boundaries
237:   */
238:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
239:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);

241:   VecGetArray(localX, &x);
242:   VecGetArray(G, &g);

244:   for (i = xs; i < xs + xm; i++) {
245:     xi     = (i + 1) * hx;
246:     trule1 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi, ecc)) / six;      /* L(i,j) */
247:     trule2 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi, ecc)) / six;      /* U(i,j) */
248:     trule3 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi + hx, ecc)) / six; /* U(i+1,j) */
249:     trule4 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi - hx, ecc)) / six; /* L(i-1,j) */
250:     trule5 = trule1;                                                        /* L(i,j-1) */
251:     trule6 = trule2;                                                        /* U(i,j+1) */

253:     vdown   = -(trule5 + trule2) * hyhy;
254:     vleft   = -hxhx * (trule2 + trule4);
255:     vright  = -hxhx * (trule1 + trule3);
256:     vup     = -hyhy * (trule1 + trule6);
257:     vmiddle = (hxhx) * (trule1 + trule2 + trule3 + trule4) + hyhy * (trule1 + trule2 + trule5 + trule6);

259:     for (j = ys; j < ys + ym; j++) {
260:       row  = (j - gys) * gxm + (i - gxs);
261:       v[0] = 0;
262:       v[1] = 0;
263:       v[2] = 0;
264:       v[3] = 0;
265:       v[4] = 0;

267:       k = 0;
268:       if (j > gys) {
269:         v[k]   = vdown;
270:         col[k] = row - gxm;
271:         k++;
272:       }

274:       if (i > gxs) {
275:         v[k]   = vleft;
276:         col[k] = row - 1;
277:         k++;
278:       }

280:       v[k]   = vmiddle;
281:       col[k] = row;
282:       k++;

284:       if (i + 1 < gxs + gxm) {
285:         v[k]   = vright;
286:         col[k] = row + 1;
287:         k++;
288:       }

290:       if (j + 1 < gys + gym) {
291:         v[k]   = vup;
292:         col[k] = row + gxm;
293:         k++;
294:       }
295:       tt = 0;
296:       for (kk = 0; kk < k; kk++) tt += v[kk] * x[col[kk]];
297:       row    = (j - ys) * xm + (i - xs);
298:       g[row] = tt;
299:     }
300:   }

302:   VecRestoreArray(localX, &x);
303:   VecRestoreArray(G, &g);

305:   DMRestoreLocalVector(user->dm, &localX);

307:   VecDot(X, G, &f1);
308:   VecDot(user->B, X, &f2);
309:   VecAXPY(G, one, user->B);
310:   *fcn = f1 / 2.0 + f2;

312:   PetscLogFlops((91 + 10.0 * ym) * xm);
313:   return 0;
314: }

316: /*
317:    FormHessian computes the quadratic term in the quadratic objective function
318:    Notice that the objective function in this problem is quadratic (therefore a constant
319:    hessian).  If using a nonquadratic solver, then you might want to reconsider this function
320: */
321: PetscErrorCode FormHessian(Tao tao, Vec X, Mat hes, Mat Hpre, void *ptr)
322: {
323:   AppCtx   *user = (AppCtx *)ptr;
324:   PetscInt  i, j, k;
325:   PetscInt  col[5], row, nx, ny, xs, xm, gxs, gxm, ys, ym, gys, gym;
326:   PetscReal one = 1.0, two = 2.0, six = 6.0, pi = 4.0 * atan(1.0);
327:   PetscReal hx, hy, hxhy, hxhx, hyhy;
328:   PetscReal xi, v[5];
329:   PetscReal ecc = user->ecc, trule1, trule2, trule3, trule4, trule5, trule6;
330:   PetscReal vmiddle, vup, vdown, vleft, vright;
331:   PetscBool assembled;

333:   nx   = user->nx;
334:   ny   = user->ny;
335:   hx   = two * pi / (nx + 1.0);
336:   hy   = two * user->b / (ny + 1.0);
337:   hxhy = hx * hy;
338:   hxhx = one / (hx * hx);
339:   hyhy = one / (hy * hy);

341:   /*
342:     Get local grid boundaries
343:   */
344:   DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL);
345:   DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL);
346:   MatAssembled(hes, &assembled);
347:   if (assembled) MatZeroEntries(hes);

349:   for (i = xs; i < xs + xm; i++) {
350:     xi     = (i + 1) * hx;
351:     trule1 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi, ecc)) / six;      /* L(i,j) */
352:     trule2 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi, ecc)) / six;      /* U(i,j) */
353:     trule3 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi + hx, ecc)) / six; /* U(i+1,j) */
354:     trule4 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi - hx, ecc)) / six; /* L(i-1,j) */
355:     trule5 = trule1;                                                        /* L(i,j-1) */
356:     trule6 = trule2;                                                        /* U(i,j+1) */

358:     vdown   = -(trule5 + trule2) * hyhy;
359:     vleft   = -hxhx * (trule2 + trule4);
360:     vright  = -hxhx * (trule1 + trule3);
361:     vup     = -hyhy * (trule1 + trule6);
362:     vmiddle = (hxhx) * (trule1 + trule2 + trule3 + trule4) + hyhy * (trule1 + trule2 + trule5 + trule6);
363:     v[0]    = 0;
364:     v[1]    = 0;
365:     v[2]    = 0;
366:     v[3]    = 0;
367:     v[4]    = 0;

369:     for (j = ys; j < ys + ym; j++) {
370:       row = (j - gys) * gxm + (i - gxs);

372:       k = 0;
373:       if (j > gys) {
374:         v[k]   = vdown;
375:         col[k] = row - gxm;
376:         k++;
377:       }

379:       if (i > gxs) {
380:         v[k]   = vleft;
381:         col[k] = row - 1;
382:         k++;
383:       }

385:       v[k]   = vmiddle;
386:       col[k] = row;
387:       k++;

389:       if (i + 1 < gxs + gxm) {
390:         v[k]   = vright;
391:         col[k] = row + 1;
392:         k++;
393:       }

395:       if (j + 1 < gys + gym) {
396:         v[k]   = vup;
397:         col[k] = row + gxm;
398:         k++;
399:       }
400:       MatSetValuesLocal(hes, 1, &row, k, col, v, INSERT_VALUES);
401:     }
402:   }

404:   /*
405:      Assemble matrix, using the 2-step process:
406:      MatAssemblyBegin(), MatAssemblyEnd().
407:      By placing code between these two statements, computations can be
408:      done while messages are in transition.
409:   */
410:   MatAssemblyBegin(hes, MAT_FINAL_ASSEMBLY);
411:   MatAssemblyEnd(hes, MAT_FINAL_ASSEMBLY);

413:   /*
414:     Tell the matrix we will never add a new nonzero location to the
415:     matrix. If we do it will generate an error.
416:   */
417:   MatSetOption(hes, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
418:   MatSetOption(hes, MAT_SYMMETRIC, PETSC_TRUE);

420:   PetscLogFlops(9.0 * xm * ym + 49.0 * xm);
421:   return 0;
422: }

424: PetscErrorCode Monitor(Tao tao, void *ctx)
425: {
426:   PetscInt           its;
427:   PetscReal          f, gnorm, cnorm, xdiff;
428:   TaoConvergedReason reason;

430:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
431:   if (!(its % 5)) PetscPrintf(PETSC_COMM_WORLD, "iteration=%" PetscInt_FMT "\tf=%g\n", its, (double)f);
432:   return 0;
433: }

435: PetscErrorCode ConvergenceTest(Tao tao, void *ctx)
436: {
437:   PetscInt           its;
438:   PetscReal          f, gnorm, cnorm, xdiff;
439:   TaoConvergedReason reason;

441:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
442:   if (its == 100) TaoSetConvergedReason(tao, TAO_DIVERGED_MAXITS);
443:   return 0;
444: }

446: /*TEST

448:    build:
449:       requires: !complex

451:    test:
452:       args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gatol 1.e-5
453:       requires: !single

455:    test:
456:       suffix: 2
457:       nsize: 2
458:       args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gatol 1.e-5
459:       requires: !single

461:    test:
462:       suffix: 3
463:       nsize: 2
464:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
465:       requires: !single

467:    test:
468:       suffix: 4
469:       nsize: 2
470:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
471:       output_file: output/jbearing2_3.out
472:       requires: !single

474:    test:
475:       suffix: 5
476:       args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
477:       requires: !single

479:    test:
480:       suffix: 6
481:       args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
482:       requires: !single

484:    test:
485:       suffix: 7
486:       args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5
487:       requires: !single

489:    test:
490:       suffix: 8
491:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5
492:       requires: !single

494:    test:
495:       suffix: 9
496:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5
497:       requires: !single

499:    test:
500:       suffix: 10
501:       args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
502:       requires: !single

504:    test:
505:       suffix: 11
506:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
507:       requires: !single

509:    test:
510:       suffix: 12
511:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
512:       requires: !single

514:    test:
515:      suffix: 13
516:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls
517:      requires: !single

519:    test:
520:      suffix: 14
521:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type blmvm
522:      requires: !single

524:    test:
525:      suffix: 15
526:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
527:      requires: !single

529:    test:
530:      suffix: 16
531:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
532:      requires: !single

534:    test:
535:      suffix: 17
536:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type scalar -tao_view
537:      requires: !single

539:    test:
540:      suffix: 18
541:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type none -tao_view
542:      requires: !single

544:    test:
545:      suffix: 19
546:      args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
547:      requires: !single

549:    test:
550:       suffix: 20
551:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
552:       requires: !single

554:    test:
555:       suffix: 21
556:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
557:       requires: !single
558: TEST*/