Actual source code: eptorsion1.c
1: /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */
3: /* ----------------------------------------------------------------------
5: Elastic-plastic torsion problem.
7: The elastic plastic torsion problem arises from the determination
8: of the stress field on an infinitely long cylindrical bar, which is
9: equivalent to the solution of the following problem:
11: min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
13: where C is the torsion angle per unit length.
15: The multiprocessor version of this code is eptorsion2.c.
17: ---------------------------------------------------------------------- */
19: /*
20: Include "petsctao.h" so that we can use TAO solvers. Note that this
21: file automatically includes files for lower-level support, such as those
22: provided by the PETSc library:
23: petsc.h - base PETSc routines petscvec.h - vectors
24: petscsys.h - system routines petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: */
29: #include <petsctao.h>
31: static char help[] = "Demonstrates use of the TAO package to solve \n\
32: unconstrained minimization problems on a single processor. This example \n\
33: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
34: test suite.\n\
35: The command line options are:\n\
36: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
37: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
38: -par <param>, where <param> = angle of twist per unit length\n\n";
40: /*
41: User-defined application context - contains data needed by the
42: application-provided call-back routines, FormFunction(),
43: FormGradient(), and FormHessian().
44: */
46: typedef struct {
47: PetscReal param; /* nonlinearity parameter */
48: PetscInt mx, my; /* discretization in x- and y-directions */
49: PetscInt ndim; /* problem dimension */
50: Vec s, y, xvec; /* work space for computing Hessian */
51: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
52: } AppCtx;
54: /* -------- User-defined Routines --------- */
56: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
57: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
58: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
59: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
60: PetscErrorCode HessianProductMat(Mat, Vec, Vec);
61: PetscErrorCode HessianProduct(void *, Vec, Vec);
62: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
63: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
65: PetscErrorCode main(int argc, char **argv)
66: {
67: PetscInt mx = 10; /* discretization in x-direction */
68: PetscInt my = 10; /* discretization in y-direction */
69: Vec x; /* solution, gradient vectors */
70: PetscBool flg; /* A return value when checking for use options */
71: Tao tao; /* Tao solver context */
72: Mat H; /* Hessian matrix */
73: AppCtx user; /* application context */
74: PetscMPIInt size; /* number of processes */
75: PetscReal one = 1.0;
77: PetscBool test_lmvm = PETSC_FALSE;
78: KSP ksp;
79: PC pc;
80: Mat M;
81: Vec in, out, out2;
82: PetscReal mult_solve_dist;
84: /* Initialize TAO,PETSc */
86: PetscInitialize(&argc, &argv, (char *)0, help);
87: MPI_Comm_size(MPI_COMM_WORLD, &size);
90: /* Specify default parameters for the problem, check for command-line overrides */
91: user.param = 5.0;
92: PetscOptionsGetInt(NULL, NULL, "-my", &my, &flg);
93: PetscOptionsGetInt(NULL, NULL, "-mx", &mx, &flg);
94: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg);
95: PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg);
97: PetscPrintf(PETSC_COMM_SELF, "\n---- Elastic-Plastic Torsion Problem -----\n");
98: PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", mx, my);
99: user.ndim = mx * my;
100: user.mx = mx;
101: user.my = my;
102: user.hx = one / (mx + 1);
103: user.hy = one / (my + 1);
105: /* Allocate vectors */
106: VecCreateSeq(PETSC_COMM_SELF, user.ndim, &user.y);
107: VecDuplicate(user.y, &user.s);
108: VecDuplicate(user.y, &x);
110: /* Create TAO solver and set desired solution method */
111: TaoCreate(PETSC_COMM_SELF, &tao);
112: TaoSetType(tao, TAOLMVM);
114: /* Set solution vector with an initial guess */
115: FormInitialGuess(&user, x);
116: TaoSetSolution(tao, x);
118: /* Set routine for function and gradient evaluation */
119: TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
121: /* From command line options, determine if using matrix-free hessian */
122: PetscOptionsHasName(NULL, NULL, "-my_tao_mf", &flg);
123: if (flg) {
124: MatCreateShell(PETSC_COMM_SELF, user.ndim, user.ndim, user.ndim, user.ndim, (void *)&user, &H);
125: MatShellSetOperation(H, MATOP_MULT, (void (*)(void))HessianProductMat);
126: MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);
128: TaoSetHessian(tao, H, H, MatrixFreeHessian, (void *)&user);
129: } else {
130: MatCreateSeqAIJ(PETSC_COMM_SELF, user.ndim, user.ndim, 5, NULL, &H);
131: MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);
132: TaoSetHessian(tao, H, H, FormHessian, (void *)&user);
133: }
135: /* Test the LMVM matrix */
136: if (test_lmvm) {
137: PetscOptionsSetValue(NULL, "-tao_type", "bntr");
138: PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
139: }
141: /* Check for any TAO command line options */
142: TaoSetFromOptions(tao);
144: /* SOLVE THE APPLICATION */
145: TaoSolve(tao);
147: /* Test the LMVM matrix */
148: if (test_lmvm) {
149: TaoGetKSP(tao, &ksp);
150: KSPGetPC(ksp, &pc);
151: PCLMVMGetMatLMVM(pc, &M);
152: VecDuplicate(x, &in);
153: VecDuplicate(x, &out);
154: VecDuplicate(x, &out2);
155: VecSet(in, 5.0);
156: MatMult(M, in, out);
157: MatSolve(M, out, out2);
158: VecAXPY(out2, -1.0, in);
159: VecNorm(out2, NORM_2, &mult_solve_dist);
160: PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", (double)mult_solve_dist);
161: VecDestroy(&in);
162: VecDestroy(&out);
163: VecDestroy(&out2);
164: }
166: TaoDestroy(&tao);
167: VecDestroy(&user.s);
168: VecDestroy(&user.y);
169: VecDestroy(&x);
170: MatDestroy(&H);
172: PetscFinalize();
173: return 0;
174: }
176: /* ------------------------------------------------------------------- */
177: /*
178: FormInitialGuess - Computes an initial approximation to the solution.
180: Input Parameters:
181: . user - user-defined application context
182: . X - vector
184: Output Parameters:
185: . X - vector
186: */
187: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
188: {
189: PetscReal hx = user->hx, hy = user->hy, temp;
190: PetscReal val;
191: PetscInt i, j, k, nx = user->mx, ny = user->my;
193: /* Compute initial guess */
195: for (j = 0; j < ny; j++) {
196: temp = PetscMin(j + 1, ny - j) * hy;
197: for (i = 0; i < nx; i++) {
198: k = nx * j + i;
199: val = PetscMin((PetscMin(i + 1, nx - i)) * hx, temp);
200: VecSetValues(X, 1, &k, &val, ADD_VALUES);
201: }
202: }
203: VecAssemblyBegin(X);
204: VecAssemblyEnd(X);
205: return 0;
206: }
208: /* ------------------------------------------------------------------- */
209: /*
210: FormFunctionGradient - Evaluates the function and corresponding gradient.
212: Input Parameters:
213: tao - the Tao context
214: X - the input vector
215: ptr - optional user-defined context, as set by TaoSetFunction()
217: Output Parameters:
218: f - the newly evaluated function
219: G - the newly evaluated gradient
220: */
221: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
222: {
224: FormFunction(tao, X, f, ptr);
225: FormGradient(tao, X, G, ptr);
226: return 0;
227: }
229: /* ------------------------------------------------------------------- */
230: /*
231: FormFunction - Evaluates the function, f(X).
233: Input Parameters:
234: . tao - the Tao context
235: . X - the input vector
236: . ptr - optional user-defined context, as set by TaoSetFunction()
238: Output Parameters:
239: . f - the newly evaluated function
240: */
241: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
242: {
243: AppCtx *user = (AppCtx *)ptr;
244: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
245: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
246: PetscReal v, cdiv3 = user->param / three;
247: const PetscScalar *x;
248: PetscInt nx = user->mx, ny = user->my, i, j, k;
251: /* Get pointer to vector data */
252: VecGetArrayRead(X, &x);
254: /* Compute function contributions over the lower triangular elements */
255: for (j = -1; j < ny; j++) {
256: for (i = -1; i < nx; i++) {
257: k = nx * j + i;
258: v = zero;
259: vr = zero;
260: vt = zero;
261: if (i >= 0 && j >= 0) v = x[k];
262: if (i < nx - 1 && j > -1) vr = x[k + 1];
263: if (i > -1 && j < ny - 1) vt = x[k + nx];
264: dvdx = (vr - v) / hx;
265: dvdy = (vt - v) / hy;
266: fquad += dvdx * dvdx + dvdy * dvdy;
267: flin -= cdiv3 * (v + vr + vt);
268: }
269: }
271: /* Compute function contributions over the upper triangular elements */
272: for (j = 0; j <= ny; j++) {
273: for (i = 0; i <= nx; i++) {
274: k = nx * j + i;
275: vb = zero;
276: vl = zero;
277: v = zero;
278: if (i < nx && j > 0) vb = x[k - nx];
279: if (i > 0 && j < ny) vl = x[k - 1];
280: if (i < nx && j < ny) v = x[k];
281: dvdx = (v - vl) / hx;
282: dvdy = (v - vb) / hy;
283: fquad = fquad + dvdx * dvdx + dvdy * dvdy;
284: flin = flin - cdiv3 * (vb + vl + v);
285: }
286: }
288: /* Restore vector */
289: VecRestoreArrayRead(X, &x);
291: /* Scale the function */
292: area = p5 * hx * hy;
293: *f = area * (p5 * fquad + flin);
295: PetscLogFlops(24.0 * nx * ny);
296: return 0;
297: }
299: /* ------------------------------------------------------------------- */
300: /*
301: FormGradient - Evaluates the gradient, G(X).
303: Input Parameters:
304: . tao - the Tao context
305: . X - input vector
306: . ptr - optional user-defined context
308: Output Parameters:
309: . G - vector containing the newly evaluated gradient
310: */
311: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
312: {
313: AppCtx *user = (AppCtx *)ptr;
314: PetscReal zero = 0.0, p5 = 0.5, three = 3.0, area, val;
315: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
316: PetscReal hx = user->hx, hy = user->hy;
317: PetscReal vb, vl, vr, vt, dvdx, dvdy;
318: PetscReal v, cdiv3 = user->param / three;
319: const PetscScalar *x;
322: /* Initialize gradient to zero */
323: VecSet(G, zero);
325: /* Get pointer to vector data */
326: VecGetArrayRead(X, &x);
328: /* Compute gradient contributions over the lower triangular elements */
329: for (j = -1; j < ny; j++) {
330: for (i = -1; i < nx; i++) {
331: k = nx * j + i;
332: v = zero;
333: vr = zero;
334: vt = zero;
335: if (i >= 0 && j >= 0) v = x[k];
336: if (i < nx - 1 && j > -1) vr = x[k + 1];
337: if (i > -1 && j < ny - 1) vt = x[k + nx];
338: dvdx = (vr - v) / hx;
339: dvdy = (vt - v) / hy;
340: if (i != -1 && j != -1) {
341: ind = k;
342: val = -dvdx / hx - dvdy / hy - cdiv3;
343: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
344: }
345: if (i != nx - 1 && j != -1) {
346: ind = k + 1;
347: val = dvdx / hx - cdiv3;
348: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
349: }
350: if (i != -1 && j != ny - 1) {
351: ind = k + nx;
352: val = dvdy / hy - cdiv3;
353: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
354: }
355: }
356: }
358: /* Compute gradient contributions over the upper triangular elements */
359: for (j = 0; j <= ny; j++) {
360: for (i = 0; i <= nx; i++) {
361: k = nx * j + i;
362: vb = zero;
363: vl = zero;
364: v = zero;
365: if (i < nx && j > 0) vb = x[k - nx];
366: if (i > 0 && j < ny) vl = x[k - 1];
367: if (i < nx && j < ny) v = x[k];
368: dvdx = (v - vl) / hx;
369: dvdy = (v - vb) / hy;
370: if (i != nx && j != 0) {
371: ind = k - nx;
372: val = -dvdy / hy - cdiv3;
373: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
374: }
375: if (i != 0 && j != ny) {
376: ind = k - 1;
377: val = -dvdx / hx - cdiv3;
378: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
379: }
380: if (i != nx && j != ny) {
381: ind = k;
382: val = dvdx / hx + dvdy / hy - cdiv3;
383: VecSetValues(G, 1, &ind, &val, ADD_VALUES);
384: }
385: }
386: }
387: VecRestoreArrayRead(X, &x);
389: /* Assemble gradient vector */
390: VecAssemblyBegin(G);
391: VecAssemblyEnd(G);
393: /* Scale the gradient */
394: area = p5 * hx * hy;
395: VecScale(G, area);
396: PetscLogFlops(24.0 * nx * ny);
397: return 0;
398: }
400: /* ------------------------------------------------------------------- */
401: /*
402: FormHessian - Forms the Hessian matrix.
404: Input Parameters:
405: . tao - the Tao context
406: . X - the input vector
407: . ptr - optional user-defined context, as set by TaoSetHessian()
409: Output Parameters:
410: . H - Hessian matrix
411: . PrecH - optionally different preconditioning Hessian
412: . flag - flag indicating matrix structure
414: Notes:
415: This routine is intended simply as an example of the interface
416: to a Hessian evaluation routine. Since this example compute the
417: Hessian a column at a time, it is not particularly efficient and
418: is not recommended.
419: */
420: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
421: {
422: AppCtx *user = (AppCtx *)ptr;
423: PetscInt i, j, ndim = user->ndim;
424: PetscReal *y, zero = 0.0, one = 1.0;
425: PetscBool assembled;
428: user->xvec = X;
430: /* Initialize Hessian entries and work vector to zero */
431: MatAssembled(H, &assembled);
432: if (assembled) MatZeroEntries(H);
434: VecSet(user->s, zero);
436: /* Loop over matrix columns to compute entries of the Hessian */
437: for (j = 0; j < ndim; j++) {
438: VecSetValues(user->s, 1, &j, &one, INSERT_VALUES);
439: VecAssemblyBegin(user->s);
440: VecAssemblyEnd(user->s);
442: HessianProduct(ptr, user->s, user->y);
444: VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES);
445: VecAssemblyBegin(user->s);
446: VecAssemblyEnd(user->s);
448: VecGetArray(user->y, &y);
449: for (i = 0; i < ndim; i++) {
450: if (y[i] != zero) MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES);
451: }
452: VecRestoreArray(user->y, &y);
453: }
454: MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
455: MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
456: return 0;
457: }
459: /* ------------------------------------------------------------------- */
460: /*
461: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
462: products.
464: Input Parameters:
465: . tao - the Tao context
466: . X - the input vector
467: . ptr - optional user-defined context, as set by TaoSetHessian()
469: Output Parameters:
470: . H - Hessian matrix
471: . PrecH - optionally different preconditioning Hessian
472: . flag - flag indicating matrix structure
473: */
474: PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
475: {
476: AppCtx *user = (AppCtx *)ptr;
478: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
480: user->xvec = X;
481: return 0;
482: }
484: /* ------------------------------------------------------------------- */
485: /*
486: HessianProductMat - Computes the matrix-vector product
487: y = mat*svec.
489: Input Parameters:
490: . mat - input matrix
491: . svec - input vector
493: Output Parameters:
494: . y - solution vector
495: */
496: PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
497: {
498: void *ptr;
501: MatShellGetContext(mat, &ptr);
502: HessianProduct(ptr, svec, y);
503: return 0;
504: }
506: /* ------------------------------------------------------------------- */
507: /*
508: Hessian Product - Computes the matrix-vector product:
509: y = f''(x)*svec.
511: Input Parameters:
512: . ptr - pointer to the user-defined context
513: . svec - input vector
515: Output Parameters:
516: . y - product vector
517: */
518: PetscErrorCode HessianProduct(void *ptr, Vec svec, Vec y)
519: {
520: AppCtx *user = (AppCtx *)ptr;
521: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
522: const PetscScalar *x, *s;
523: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
524: PetscInt nx, ny, i, j, k, ind;
527: nx = user->mx;
528: ny = user->my;
529: hx = user->hx;
530: hy = user->hy;
531: hxhx = one / (hx * hx);
532: hyhy = one / (hy * hy);
534: /* Get pointers to vector data */
535: VecGetArrayRead(user->xvec, &x);
536: VecGetArrayRead(svec, &s);
538: /* Initialize product vector to zero */
539: VecSet(y, zero);
541: /* Compute f''(x)*s over the lower triangular elements */
542: for (j = -1; j < ny; j++) {
543: for (i = -1; i < nx; i++) {
544: k = nx * j + i;
545: v = zero;
546: vr = zero;
547: vt = zero;
548: if (i != -1 && j != -1) v = s[k];
549: if (i != nx - 1 && j != -1) {
550: vr = s[k + 1];
551: ind = k + 1;
552: val = hxhx * (vr - v);
553: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
554: }
555: if (i != -1 && j != ny - 1) {
556: vt = s[k + nx];
557: ind = k + nx;
558: val = hyhy * (vt - v);
559: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
560: }
561: if (i != -1 && j != -1) {
562: ind = k;
563: val = hxhx * (v - vr) + hyhy * (v - vt);
564: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
565: }
566: }
567: }
569: /* Compute f''(x)*s over the upper triangular elements */
570: for (j = 0; j <= ny; j++) {
571: for (i = 0; i <= nx; i++) {
572: k = nx * j + i;
573: v = zero;
574: vl = zero;
575: vb = zero;
576: if (i != nx && j != ny) v = s[k];
577: if (i != nx && j != 0) {
578: vb = s[k - nx];
579: ind = k - nx;
580: val = hyhy * (vb - v);
581: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
582: }
583: if (i != 0 && j != ny) {
584: vl = s[k - 1];
585: ind = k - 1;
586: val = hxhx * (vl - v);
587: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
588: }
589: if (i != nx && j != ny) {
590: ind = k;
591: val = hxhx * (v - vl) + hyhy * (v - vb);
592: VecSetValues(y, 1, &ind, &val, ADD_VALUES);
593: }
594: }
595: }
596: /* Restore vector data */
597: VecRestoreArrayRead(svec, &s);
598: VecRestoreArrayRead(user->xvec, &x);
600: /* Assemble vector */
601: VecAssemblyBegin(y);
602: VecAssemblyEnd(y);
604: /* Scale resulting vector by area */
605: area = p5 * hx * hy;
606: VecScale(y, area);
607: PetscLogFlops(18.0 * nx * ny);
608: return 0;
609: }
611: /*TEST
613: build:
614: requires: !complex
616: test:
617: suffix: 1
618: args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4
620: test:
621: suffix: 2
622: args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4
624: test:
625: suffix: 3
626: args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
628: test:
629: suffix: 4
630: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
632: test:
633: suffix: 5
634: args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
636: test:
637: suffix: 6
638: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
640: TEST*/