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*/