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 */
62: PetscFunctionBeginUser;
63: PetscCall(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: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.nx, &flg));
74: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.ny, &flg));
75: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ecc", &user.ecc, &flg));
76: PetscCall(PetscOptionsGetReal(NULL, NULL, "-b", &user.b, &flg));
77: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_getdiagonal", &testgetdiag, NULL));
79: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Journal Bearing Problem SHB-----\n"));
80: PetscCall(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: PetscCall(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: PetscCall(DMSetFromOptions(user.dm));
93: PetscCall(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: PetscCall(DMCreateGlobalVector(user.dm, &x)); /* Solution */
102: PetscCall(VecDuplicate(x, &user.B)); /* Linear objective */
104: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
105: PetscCall(VecGetLocalSize(x, &m));
106: PetscCall(DMCreateMatrix(user.dm, &user.A));
108: if (testgetdiag) PetscCall(MatSetOperation(user.A, MATOP_GET_DIAGONAL, NULL));
110: /* User defined function -- compute linear term of quadratic */
111: PetscCall(ComputeB(&user));
113: /* The TAO code begins here */
115: /*
116: Create the optimization solver
117: Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
118: */
119: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
120: PetscCall(TaoSetType(tao, TAOBLMVM));
122: /* Set the initial vector */
123: PetscCall(VecSet(x, zero));
124: PetscCall(TaoSetSolution(tao, x));
126: /* Set the user function, gradient, hessian evaluation routines and data structures */
127: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
129: PetscCall(TaoSetHessian(tao, user.A, user.A, FormHessian, (void *)&user));
131: /* Set a routine that defines the bounds */
132: PetscCall(VecDuplicate(x, &xl));
133: PetscCall(VecDuplicate(x, &xu));
134: PetscCall(VecSet(xl, zero));
135: PetscCall(VecSet(xu, d1000));
136: PetscCall(TaoSetVariableBounds(tao, xl, xu));
138: PetscCall(TaoGetKSP(tao, &ksp));
139: if (ksp) PetscCall(KSPSetType(ksp, KSPCG));
141: PetscCall(PetscOptionsHasName(NULL, NULL, "-testmonitor", &flg));
142: if (flg) PetscCall(TaoSetMonitor(tao, Monitor, &user, NULL));
143: PetscCall(PetscOptionsHasName(NULL, NULL, "-testconvergence", &flg));
144: if (flg) PetscCall(TaoSetConvergenceTest(tao, ConvergenceTest, &user));
146: /* Check for any tao command line options */
147: PetscCall(TaoSetFromOptions(tao));
149: /* Solve the bound constrained problem */
150: PetscCall(TaoSolve(tao));
152: /* Free PETSc data structures */
153: PetscCall(VecDestroy(&x));
154: PetscCall(VecDestroy(&xl));
155: PetscCall(VecDestroy(&xu));
156: PetscCall(MatDestroy(&user.A));
157: PetscCall(VecDestroy(&user.B));
159: /* Free TAO data structures */
160: PetscCall(TaoDestroy(&tao));
161: PetscCall(DMDestroy(&user.dm));
162: PetscCall(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: PetscFunctionBegin;
182: nx = user->nx;
183: ny = user->ny;
184: hx = two * pi / (nx + 1.0);
185: hy = two * user->b / (ny + 1.0);
186: ehxhy = ecc * hx * hy;
188: /*
189: Get local grid boundaries
190: */
191: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
192: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
194: /* Compute the linear term in the objective function */
195: PetscCall(VecGetArray(user->B, &b));
196: for (i = xs; i < xs + xm; i++) {
197: temp = PetscSinScalar((i + 1) * hx);
198: for (j = ys; j < ys + ym; j++) {
199: k = xm * (j - ys) + (i - xs);
200: b[k] = -ehxhy * temp;
201: }
202: }
203: PetscCall(VecRestoreArray(user->B, &b));
204: PetscCall(PetscLogFlops(5.0 * xm * ym + 3.0 * xm));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr)
209: {
210: AppCtx *user = (AppCtx *)ptr;
211: PetscInt i, j, k, kk;
212: PetscInt col[5], row, nx, ny, xs, xm, gxs, gxm, ys, ym, gys, gym;
213: PetscReal one = 1.0, two = 2.0, six = 6.0, pi = 4.0 * atan(1.0);
214: PetscReal hx, hy, hxhy, hxhx, hyhy;
215: PetscReal xi, v[5];
216: PetscReal ecc = user->ecc, trule1, trule2, trule3, trule4, trule5, trule6;
217: PetscReal vmiddle, vup, vdown, vleft, vright;
218: PetscReal tt, f1, f2;
219: PetscReal *x, *g, zero = 0.0;
220: Vec localX;
222: PetscFunctionBegin;
223: nx = user->nx;
224: ny = user->ny;
225: hx = two * pi / (nx + 1.0);
226: hy = two * user->b / (ny + 1.0);
227: hxhy = hx * hy;
228: hxhx = one / (hx * hx);
229: hyhy = one / (hy * hy);
231: PetscCall(DMGetLocalVector(user->dm, &localX));
233: PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
234: PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
236: PetscCall(VecSet(G, zero));
237: /*
238: Get local grid boundaries
239: */
240: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
241: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
243: PetscCall(VecGetArray(localX, &x));
244: PetscCall(VecGetArray(G, &g));
246: for (i = xs; i < xs + xm; i++) {
247: xi = (i + 1) * hx;
248: trule1 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi, ecc)) / six; /* L(i,j) */
249: trule2 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi, ecc)) / six; /* U(i,j) */
250: trule3 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi + hx, ecc)) / six; /* U(i+1,j) */
251: trule4 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi - hx, ecc)) / six; /* L(i-1,j) */
252: trule5 = trule1; /* L(i,j-1) */
253: trule6 = trule2; /* U(i,j+1) */
255: vdown = -(trule5 + trule2) * hyhy;
256: vleft = -hxhx * (trule2 + trule4);
257: vright = -hxhx * (trule1 + trule3);
258: vup = -hyhy * (trule1 + trule6);
259: vmiddle = (hxhx) * (trule1 + trule2 + trule3 + trule4) + hyhy * (trule1 + trule2 + trule5 + trule6);
261: for (j = ys; j < ys + ym; j++) {
262: row = (j - gys) * gxm + (i - gxs);
263: v[0] = 0;
264: v[1] = 0;
265: v[2] = 0;
266: v[3] = 0;
267: v[4] = 0;
269: k = 0;
270: if (j > gys) {
271: v[k] = vdown;
272: col[k] = row - gxm;
273: k++;
274: }
276: if (i > gxs) {
277: v[k] = vleft;
278: col[k] = row - 1;
279: k++;
280: }
282: v[k] = vmiddle;
283: col[k] = row;
284: k++;
286: if (i + 1 < gxs + gxm) {
287: v[k] = vright;
288: col[k] = row + 1;
289: k++;
290: }
292: if (j + 1 < gys + gym) {
293: v[k] = vup;
294: col[k] = row + gxm;
295: k++;
296: }
297: tt = 0;
298: for (kk = 0; kk < k; kk++) tt += v[kk] * x[col[kk]];
299: row = (j - ys) * xm + (i - xs);
300: g[row] = tt;
301: }
302: }
304: PetscCall(VecRestoreArray(localX, &x));
305: PetscCall(VecRestoreArray(G, &g));
307: PetscCall(DMRestoreLocalVector(user->dm, &localX));
309: PetscCall(VecDot(X, G, &f1));
310: PetscCall(VecDot(user->B, X, &f2));
311: PetscCall(VecAXPY(G, one, user->B));
312: *fcn = f1 / 2.0 + f2;
314: PetscCall(PetscLogFlops((91 + 10.0 * ym) * xm));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*
319: FormHessian computes the quadratic term in the quadratic objective function
320: Notice that the objective function in this problem is quadratic (therefore a constant
321: hessian). If using a nonquadratic solver, then you might want to reconsider this function
322: */
323: PetscErrorCode FormHessian(Tao tao, Vec X, Mat hes, Mat Hpre, void *ptr)
324: {
325: AppCtx *user = (AppCtx *)ptr;
326: PetscInt i, j, k;
327: PetscInt col[5], row, nx, ny, xs, xm, gxs, gxm, ys, ym, gys, gym;
328: PetscReal one = 1.0, two = 2.0, six = 6.0, pi = 4.0 * atan(1.0);
329: PetscReal hx, hy, hxhy, hxhx, hyhy;
330: PetscReal xi, v[5];
331: PetscReal ecc = user->ecc, trule1, trule2, trule3, trule4, trule5, trule6;
332: PetscReal vmiddle, vup, vdown, vleft, vright;
333: PetscBool assembled;
335: PetscFunctionBegin;
336: nx = user->nx;
337: ny = user->ny;
338: hx = two * pi / (nx + 1.0);
339: hy = two * user->b / (ny + 1.0);
340: hxhy = hx * hy;
341: hxhx = one / (hx * hx);
342: hyhy = one / (hy * hy);
344: /*
345: Get local grid boundaries
346: */
347: PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
348: PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
349: PetscCall(MatAssembled(hes, &assembled));
350: if (assembled) PetscCall(MatZeroEntries(hes));
352: for (i = xs; i < xs + xm; i++) {
353: xi = (i + 1) * hx;
354: trule1 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi, ecc)) / six; /* L(i,j) */
355: trule2 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi, ecc)) / six; /* U(i,j) */
356: trule3 = hxhy * (p(xi, ecc) + p(xi + hx, ecc) + p(xi + hx, ecc)) / six; /* U(i+1,j) */
357: trule4 = hxhy * (p(xi, ecc) + p(xi - hx, ecc) + p(xi - hx, ecc)) / six; /* L(i-1,j) */
358: trule5 = trule1; /* L(i,j-1) */
359: trule6 = trule2; /* U(i,j+1) */
361: vdown = -(trule5 + trule2) * hyhy;
362: vleft = -hxhx * (trule2 + trule4);
363: vright = -hxhx * (trule1 + trule3);
364: vup = -hyhy * (trule1 + trule6);
365: vmiddle = (hxhx) * (trule1 + trule2 + trule3 + trule4) + hyhy * (trule1 + trule2 + trule5 + trule6);
366: v[0] = 0;
367: v[1] = 0;
368: v[2] = 0;
369: v[3] = 0;
370: v[4] = 0;
372: for (j = ys; j < ys + ym; j++) {
373: row = (j - gys) * gxm + (i - gxs);
375: k = 0;
376: if (j > gys) {
377: v[k] = vdown;
378: col[k] = row - gxm;
379: k++;
380: }
382: if (i > gxs) {
383: v[k] = vleft;
384: col[k] = row - 1;
385: k++;
386: }
388: v[k] = vmiddle;
389: col[k] = row;
390: k++;
392: if (i + 1 < gxs + gxm) {
393: v[k] = vright;
394: col[k] = row + 1;
395: k++;
396: }
398: if (j + 1 < gys + gym) {
399: v[k] = vup;
400: col[k] = row + gxm;
401: k++;
402: }
403: PetscCall(MatSetValuesLocal(hes, 1, &row, k, col, v, INSERT_VALUES));
404: }
405: }
407: /*
408: Assemble matrix, using the 2-step process:
409: MatAssemblyBegin(), MatAssemblyEnd().
410: By placing code between these two statements, computations can be
411: done while messages are in transition.
412: */
413: PetscCall(MatAssemblyBegin(hes, MAT_FINAL_ASSEMBLY));
414: PetscCall(MatAssemblyEnd(hes, MAT_FINAL_ASSEMBLY));
416: /*
417: Tell the matrix we will never add a new nonzero location to the
418: matrix. If we do it will generate an error.
419: */
420: PetscCall(MatSetOption(hes, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
421: PetscCall(MatSetOption(hes, MAT_SYMMETRIC, PETSC_TRUE));
423: PetscCall(PetscLogFlops(9.0 * xm * ym + 49.0 * xm));
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }
427: PetscErrorCode Monitor(Tao tao, void *ctx)
428: {
429: PetscInt its;
430: PetscReal f, gnorm, cnorm, xdiff;
431: TaoConvergedReason reason;
433: PetscFunctionBegin;
434: PetscCall(TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason));
435: if (!(its % 5)) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iteration=%" PetscInt_FMT "\tf=%g\n", its, (double)f));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: PetscErrorCode ConvergenceTest(Tao tao, void *ctx)
440: {
441: PetscInt its;
442: PetscReal f, gnorm, cnorm, xdiff;
443: TaoConvergedReason reason;
445: PetscFunctionBegin;
446: PetscCall(TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason));
447: if (its == 100) PetscCall(TaoSetConvergedReason(tao, TAO_DIVERGED_MAXITS));
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*TEST
453: build:
454: requires: !complex
456: test:
457: args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gatol 1.e-5
458: requires: !single
460: test:
461: suffix: 2
462: nsize: 2
463: args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gatol 1.e-5
464: requires: !single
466: test:
467: suffix: 3
468: nsize: 2
469: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
470: requires: !single
472: test:
473: suffix: 4
474: nsize: 2
475: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
476: output_file: output/jbearing2_3.out
477: requires: !single
479: test:
480: suffix: 5
481: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
482: requires: !single
484: test:
485: suffix: 6
486: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
487: requires: !single
489: test:
490: suffix: 7
491: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5
492: requires: !single
494: test:
495: suffix: 8
496: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5
497: requires: !single
499: test:
500: suffix: 9
501: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5
502: requires: !single
504: test:
505: suffix: 10
506: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
507: requires: !single
509: test:
510: suffix: 11
511: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
512: requires: !single
514: test:
515: suffix: 12
516: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
517: requires: !single
519: test:
520: suffix: 13
521: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls
522: requires: !single
524: test:
525: suffix: 14
526: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type blmvm
527: requires: !single
529: test:
530: suffix: 15
531: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
532: requires: !single
534: test:
535: suffix: 16
536: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
537: requires: !single
539: test:
540: suffix: 17
541: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type scalar -tao_view
542: requires: !single
544: test:
545: suffix: 18
546: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type none -tao_view
547: requires: !single
549: test:
550: suffix: 19
551: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
552: requires: !single
554: test:
555: suffix: 20
556: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
557: requires: !single
559: test:
560: suffix: 21
561: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
562: requires: !single
563: TEST*/