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[]=
32: "Demonstrates use of the TAO package to solve \n\
33: unconstrained minimization problems on a single processor. This example \n\
34: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
35: test suite.\n\
36: The command line options are:\n\
37: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
38: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
39: -par <param>, where <param> = angle of twist per unit length\n\n";
41: /*
42: User-defined application context - contains data needed by the
43: application-provided call-back routines, FormFunction(),
44: FormGradient(), and FormHessian().
45: */
47: typedef struct {
48: PetscReal param; /* nonlinearity parameter */
49: PetscInt mx, my; /* discretization in x- and y-directions */
50: PetscInt ndim; /* problem dimension */
51: Vec s, y, xvec; /* work space for computing Hessian */
52: PetscReal hx, hy; /* mesh spacing in x- and y-directions */
53: } AppCtx;
55: /* -------- User-defined Routines --------- */
57: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
58: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
59: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
60: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
61: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
62: PetscErrorCode HessianProduct(void*,Vec,Vec);
63: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
64: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
66: PetscErrorCode main(int argc,char **argv)
67: {
68: PetscInt mx=10; /* discretization in x-direction */
69: PetscInt my=10; /* discretization in y-direction */
70: Vec x; /* solution, gradient vectors */
71: PetscBool flg; /* A return value when checking for use options */
72: Tao tao; /* Tao solver context */
73: Mat H; /* Hessian matrix */
74: AppCtx user; /* application context */
75: PetscMPIInt size; /* number of processes */
76: PetscReal one=1.0;
78: PetscBool test_lmvm = PETSC_FALSE;
79: KSP ksp;
80: PC pc;
81: Mat M;
82: Vec in, out, out2;
83: PetscReal mult_solve_dist;
85: /* 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: %D my: %D \n\n",mx,my);
99: user.ndim = mx * my; user.mx = mx; user.my = my;
100: user.hx = one/(mx+1); user.hy = one/(my+1);
102: /* Allocate vectors */
103: VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
104: VecDuplicate(user.y,&user.s);
105: VecDuplicate(user.y,&x);
107: /* Create TAO solver and set desired solution method */
108: TaoCreate(PETSC_COMM_SELF,&tao);
109: TaoSetType(tao,TAOLMVM);
111: /* Set solution vector with an initial guess */
112: FormInitialGuess(&user,x);
113: TaoSetSolution(tao,x);
115: /* Set routine for function and gradient evaluation */
116: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
118: /* From command line options, determine if using matrix-free hessian */
119: PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg);
120: if (flg) {
121: MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
122: MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
123: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
125: TaoSetHessian(tao,H,H,MatrixFreeHessian,(void *)&user);
126: } else {
127: MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
128: MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
129: TaoSetHessian(tao,H,H,FormHessian,(void *)&user);
130: }
132: /* Test the LMVM matrix */
133: if (test_lmvm) {
134: PetscOptionsSetValue(NULL, "-tao_type", "bntr");
135: PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
136: }
138: /* Check for any TAO command line options */
139: TaoSetFromOptions(tao);
141: /* SOLVE THE APPLICATION */
142: TaoSolve(tao);
144: /* Test the LMVM matrix */
145: if (test_lmvm) {
146: TaoGetKSP(tao, &ksp);
147: KSPGetPC(ksp, &pc);
148: PCLMVMGetMatLMVM(pc, &M);
149: VecDuplicate(x, &in);
150: VecDuplicate(x, &out);
151: VecDuplicate(x, &out2);
152: VecSet(in, 5.0);
153: MatMult(M, in, out);
154: MatSolve(M, out, out2);
155: VecAXPY(out2, -1.0, in);
156: VecNorm(out2, NORM_2, &mult_solve_dist);
157: PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist);
158: VecDestroy(&in);
159: VecDestroy(&out);
160: VecDestroy(&out2);
161: }
163: TaoDestroy(&tao);
164: VecDestroy(&user.s);
165: VecDestroy(&user.y);
166: VecDestroy(&x);
167: MatDestroy(&H);
169: PetscFinalize();
170: return 0;
171: }
173: /* ------------------------------------------------------------------- */
174: /*
175: FormInitialGuess - Computes an initial approximation to the solution.
177: Input Parameters:
178: . user - user-defined application context
179: . X - vector
181: Output Parameters:
182: . X - vector
183: */
184: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
185: {
186: PetscReal hx = user->hx, hy = user->hy, temp;
187: PetscReal val;
188: PetscInt i, j, k, nx = user->mx, ny = user->my;
190: /* Compute initial guess */
192: for (j=0; j<ny; j++) {
193: temp = PetscMin(j+1,ny-j)*hy;
194: for (i=0; i<nx; i++) {
195: k = nx*j + i;
196: val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
197: VecSetValues(X,1,&k,&val,ADD_VALUES);
198: }
199: }
200: VecAssemblyBegin(X);
201: VecAssemblyEnd(X);
202: return 0;
203: }
205: /* ------------------------------------------------------------------- */
206: /*
207: FormFunctionGradient - Evaluates the function and corresponding gradient.
209: Input Parameters:
210: tao - the Tao context
211: X - the input vector
212: ptr - optional user-defined context, as set by TaoSetFunction()
214: Output Parameters:
215: f - the newly evaluated function
216: G - the newly evaluated gradient
217: */
218: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
219: {
221: FormFunction(tao,X,f,ptr);
222: FormGradient(tao,X,G,ptr);
223: return 0;
224: }
226: /* ------------------------------------------------------------------- */
227: /*
228: FormFunction - Evaluates the function, f(X).
230: Input Parameters:
231: . tao - the Tao context
232: . X - the input vector
233: . ptr - optional user-defined context, as set by TaoSetFunction()
235: Output Parameters:
236: . f - the newly evaluated function
237: */
238: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
239: {
240: AppCtx *user = (AppCtx *) ptr;
241: PetscReal hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
242: PetscReal zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
243: PetscReal v, cdiv3 = user->param/three;
244: const PetscScalar *x;
245: PetscInt nx = user->mx, ny = user->my, i, j, k;
248: /* Get pointer to vector data */
249: VecGetArrayRead(X,&x);
251: /* Compute function contributions over the lower triangular elements */
252: for (j=-1; j<ny; j++) {
253: for (i=-1; i<nx; i++) {
254: k = nx*j + i;
255: v = zero;
256: vr = zero;
257: vt = zero;
258: if (i >= 0 && j >= 0) v = x[k];
259: if (i < nx-1 && j > -1) vr = x[k+1];
260: if (i > -1 && j < ny-1) vt = x[k+nx];
261: dvdx = (vr-v)/hx;
262: dvdy = (vt-v)/hy;
263: fquad += dvdx*dvdx + dvdy*dvdy;
264: flin -= cdiv3*(v+vr+vt);
265: }
266: }
268: /* Compute function contributions over the upper triangular elements */
269: for (j=0; j<=ny; j++) {
270: for (i=0; i<=nx; i++) {
271: k = nx*j + i;
272: vb = zero;
273: vl = zero;
274: v = zero;
275: if (i < nx && j > 0) vb = x[k-nx];
276: if (i > 0 && j < ny) vl = x[k-1];
277: if (i < nx && j < ny) v = x[k];
278: dvdx = (v-vl)/hx;
279: dvdy = (v-vb)/hy;
280: fquad = fquad + dvdx*dvdx + dvdy*dvdy;
281: flin = flin - cdiv3*(vb+vl+v);
282: }
283: }
285: /* Restore vector */
286: VecRestoreArrayRead(X,&x);
288: /* Scale the function */
289: area = p5*hx*hy;
290: *f = area*(p5*fquad+flin);
292: PetscLogFlops(24.0*nx*ny);
293: return 0;
294: }
296: /* ------------------------------------------------------------------- */
297: /*
298: FormGradient - Evaluates the gradient, G(X).
300: Input Parameters:
301: . tao - the Tao context
302: . X - input vector
303: . ptr - optional user-defined context
305: Output Parameters:
306: . G - vector containing the newly evaluated gradient
307: */
308: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
309: {
310: AppCtx *user = (AppCtx *) ptr;
311: PetscReal zero=0.0, p5=0.5, three = 3.0, area, val;
312: PetscInt nx = user->mx, ny = user->my, ind, i, j, k;
313: PetscReal hx = user->hx, hy = user->hy;
314: PetscReal vb, vl, vr, vt, dvdx, dvdy;
315: PetscReal v, cdiv3 = user->param/three;
316: const PetscScalar *x;
319: /* Initialize gradient to zero */
320: VecSet(G, zero);
322: /* Get pointer to vector data */
323: VecGetArrayRead(X,&x);
325: /* Compute gradient contributions over the lower triangular elements */
326: for (j=-1; j<ny; j++) {
327: for (i=-1; i<nx; i++) {
328: k = nx*j + i;
329: v = zero;
330: vr = zero;
331: vt = zero;
332: if (i >= 0 && j >= 0) v = x[k];
333: if (i < nx-1 && j > -1) vr = x[k+1];
334: if (i > -1 && j < ny-1) vt = x[k+nx];
335: dvdx = (vr-v)/hx;
336: dvdy = (vt-v)/hy;
337: if (i != -1 && j != -1) {
338: ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
339: VecSetValues(G,1,&ind,&val,ADD_VALUES);
340: }
341: if (i != nx-1 && j != -1) {
342: ind = k+1; val = dvdx/hx - cdiv3;
343: VecSetValues(G,1,&ind,&val,ADD_VALUES);
344: }
345: if (i != -1 && j != ny-1) {
346: ind = k+nx; val = dvdy/hy - cdiv3;
347: VecSetValues(G,1,&ind,&val,ADD_VALUES);
348: }
349: }
350: }
352: /* Compute gradient contributions over the upper triangular elements */
353: for (j=0; j<=ny; j++) {
354: for (i=0; i<=nx; i++) {
355: k = nx*j + i;
356: vb = zero;
357: vl = zero;
358: v = zero;
359: if (i < nx && j > 0) vb = x[k-nx];
360: if (i > 0 && j < ny) vl = x[k-1];
361: if (i < nx && j < ny) v = x[k];
362: dvdx = (v-vl)/hx;
363: dvdy = (v-vb)/hy;
364: if (i != nx && j != 0) {
365: ind = k-nx; val = - dvdy/hy - cdiv3;
366: VecSetValues(G,1,&ind,&val,ADD_VALUES);
367: }
368: if (i != 0 && j != ny) {
369: ind = k-1; val = - dvdx/hx - cdiv3;
370: VecSetValues(G,1,&ind,&val,ADD_VALUES);
371: }
372: if (i != nx && j != ny) {
373: ind = k; val = dvdx/hx + dvdy/hy - cdiv3;
374: VecSetValues(G,1,&ind,&val,ADD_VALUES);
375: }
376: }
377: }
378: VecRestoreArrayRead(X,&x);
380: /* Assemble gradient vector */
381: VecAssemblyBegin(G);
382: VecAssemblyEnd(G);
384: /* Scale the gradient */
385: area = p5*hx*hy;
386: VecScale(G, area);
387: PetscLogFlops(24.0*nx*ny);
388: return 0;
389: }
391: /* ------------------------------------------------------------------- */
392: /*
393: FormHessian - Forms the Hessian matrix.
395: Input Parameters:
396: . tao - the Tao context
397: . X - the input vector
398: . ptr - optional user-defined context, as set by TaoSetHessian()
400: Output Parameters:
401: . H - Hessian matrix
402: . PrecH - optionally different preconditioning Hessian
403: . flag - flag indicating matrix structure
405: Notes:
406: This routine is intended simply as an example of the interface
407: to a Hessian evaluation routine. Since this example compute the
408: Hessian a column at a time, it is not particularly efficient and
409: is not recommended.
410: */
411: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
412: {
413: AppCtx *user = (AppCtx *) ptr;
414: PetscInt i,j, ndim = user->ndim;
415: PetscReal *y, zero = 0.0, one = 1.0;
416: PetscBool assembled;
419: user->xvec = X;
421: /* Initialize Hessian entries and work vector to zero */
422: MatAssembled(H,&assembled);
423: if (assembled)MatZeroEntries(H);
425: VecSet(user->s, zero);
427: /* Loop over matrix columns to compute entries of the Hessian */
428: for (j=0; j<ndim; j++) {
429: VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
430: VecAssemblyBegin(user->s);
431: VecAssemblyEnd(user->s);
433: HessianProduct(ptr,user->s,user->y);
435: VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
436: VecAssemblyBegin(user->s);
437: VecAssemblyEnd(user->s);
439: VecGetArray(user->y,&y);
440: for (i=0; i<ndim; i++) {
441: if (y[i] != zero) {
442: MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
443: }
444: }
445: VecRestoreArray(user->y,&y);
446: }
447: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
448: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
449: return 0;
450: }
452: /* ------------------------------------------------------------------- */
453: /*
454: MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
455: products.
457: Input Parameters:
458: . tao - the Tao context
459: . X - the input vector
460: . ptr - optional user-defined context, as set by TaoSetHessian()
462: Output Parameters:
463: . H - Hessian matrix
464: . PrecH - optionally different preconditioning Hessian
465: . flag - flag indicating matrix structure
466: */
467: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
468: {
469: AppCtx *user = (AppCtx *) ptr;
471: /* Sets location of vector for use in computing matrix-vector products of the form H(X)*y */
473: user->xvec = X;
474: return 0;
475: }
477: /* ------------------------------------------------------------------- */
478: /*
479: HessianProductMat - Computes the matrix-vector product
480: y = mat*svec.
482: Input Parameters:
483: . mat - input matrix
484: . svec - input vector
486: Output Parameters:
487: . y - solution vector
488: */
489: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
490: {
491: void *ptr;
494: MatShellGetContext(mat,&ptr);
495: HessianProduct(ptr,svec,y);
496: return 0;
497: }
499: /* ------------------------------------------------------------------- */
500: /*
501: Hessian Product - Computes the matrix-vector product:
502: y = f''(x)*svec.
504: Input Parameters:
505: . ptr - pointer to the user-defined context
506: . svec - input vector
508: Output Parameters:
509: . y - product vector
510: */
511: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
512: {
513: AppCtx *user = (AppCtx *)ptr;
514: PetscReal p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
515: const PetscScalar *x, *s;
516: PetscReal v, vb, vl, vr, vt, hxhx, hyhy;
517: PetscInt nx, ny, i, j, k, ind;
520: nx = user->mx;
521: ny = user->my;
522: hx = user->hx;
523: hy = user->hy;
524: hxhx = one/(hx*hx);
525: hyhy = one/(hy*hy);
527: /* Get pointers to vector data */
528: VecGetArrayRead(user->xvec,&x);
529: VecGetArrayRead(svec,&s);
531: /* Initialize product vector to zero */
532: VecSet(y, zero);
534: /* Compute f''(x)*s over the lower triangular elements */
535: for (j=-1; j<ny; j++) {
536: for (i=-1; i<nx; i++) {
537: k = nx*j + i;
538: v = zero;
539: vr = zero;
540: vt = zero;
541: if (i != -1 && j != -1) v = s[k];
542: if (i != nx-1 && j != -1) {
543: vr = s[k+1];
544: ind = k+1; val = hxhx*(vr-v);
545: VecSetValues(y,1,&ind,&val,ADD_VALUES);
546: }
547: if (i != -1 && j != ny-1) {
548: vt = s[k+nx];
549: ind = k+nx; val = hyhy*(vt-v);
550: VecSetValues(y,1,&ind,&val,ADD_VALUES);
551: }
552: if (i != -1 && j != -1) {
553: ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
554: VecSetValues(y,1,&ind,&val,ADD_VALUES);
555: }
556: }
557: }
559: /* Compute f''(x)*s over the upper triangular elements */
560: for (j=0; j<=ny; j++) {
561: for (i=0; i<=nx; i++) {
562: k = nx*j + i;
563: v = zero;
564: vl = zero;
565: vb = zero;
566: if (i != nx && j != ny) v = s[k];
567: if (i != nx && j != 0) {
568: vb = s[k-nx];
569: ind = k-nx; val = hyhy*(vb-v);
570: VecSetValues(y,1,&ind,&val,ADD_VALUES);
571: }
572: if (i != 0 && j != ny) {
573: vl = s[k-1];
574: ind = k-1; val = hxhx*(vl-v);
575: VecSetValues(y,1,&ind,&val,ADD_VALUES);
576: }
577: if (i != nx && j != ny) {
578: ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
579: VecSetValues(y,1,&ind,&val,ADD_VALUES);
580: }
581: }
582: }
583: /* Restore vector data */
584: VecRestoreArrayRead(svec,&s);
585: VecRestoreArrayRead(user->xvec,&x);
587: /* Assemble vector */
588: VecAssemblyBegin(y);
589: VecAssemblyEnd(y);
591: /* Scale resulting vector by area */
592: area = p5*hx*hy;
593: VecScale(y, area);
594: PetscLogFlops(18.0*nx*ny);
595: return 0;
596: }
598: /*TEST
600: build:
601: requires: !complex
603: test:
604: suffix: 1
605: args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4
607: test:
608: suffix: 2
609: args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4
611: test:
612: suffix: 3
613: args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
615: test:
616: suffix: 4
617: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
619: test:
620: suffix: 5
621: args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
623: test:
624: suffix: 6
625: args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
627: TEST*/