Actual source code: jbearing2.c
1: /*
2: Include "taosolver.h" so we can use TAO solvers
3: Include "petscdm.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 "petscdm.h"
9: #include "petscksp.h"
10: #include "taosolver.h"
12: static char help[]=
13: "This example demonstrates use of the TAO package to \n\
14: solve a bound constrained minimization problem. This example is based on \n\
15: the problem DPJB from the MINPACK-2 test suite. This pressure journal \n\
16: bearing problem is an example of elliptic variational problem defined over \n\
17: a two dimensional rectangle. By discretizing the domain into triangular \n\
18: elements, the pressure surrounding the journal bearing is defined as the \n\
19: minimum of a quadratic function whose variables are bounded below by zero.\n\
20: The command line options are:\n\
21: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
22: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
23: \n";
25: /*T
26: Concepts: TAO - Solving a bound constrained minimization problem
27: Routines: TaoInitialize(); TaoFinalize();
28: Routines: TaoCreate();
29: Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
30: Routines: TaoSetHessianRoutine();
31: Routines: TaoSetVariableBounds();
32: Routines: TaoSetMonitor(); TaoSetConvergenceTest();
33: Routines: TaoSetInitialVector();
34: Routines: TaoSetFromOptions();
35: Routines: TaoSolve();
36: Routines: TaoGetTerminationReason(); TaoDestroy();
37: Processors: n
38: T*/
40: /*
41: User-defined application context - contains data needed by the
42: application-provided call-back routines, FormFunctionGradient(),
43: FormHessian().
44: */
45: typedef struct {
46: /* problem parameters */
47: PetscReal ecc; /* test problem parameter */
48: PetscReal b; /* A dimension of journal bearing */
49: PetscInt nx,ny; /* discretization in x, y directions */
51: /* Working space */
52: DM dm; /* distributed array data structure */
53: Mat A; /* Quadratic Objective term */
54: Vec B; /* Linear Objective term */
55: } AppCtx;
57: /* User-defined routines */
58: static PetscReal p(PetscReal xi, PetscReal ecc);
59: static PetscErrorCode FormFunctionGradient(TaoSolver, Vec, PetscReal *,Vec,void *);
60: static PetscErrorCode FormHessian(TaoSolver,Vec,Mat *, Mat *, MatStructure *, void *);
61: static PetscErrorCode ComputeB(AppCtx*);
62: static PetscErrorCode Monitor(TaoSolver, void*);
63: static PetscErrorCode ConvergenceTest(TaoSolver, void*);
67: int main( int argc, char **argv )
68: {
69: PetscErrorCode ierr; /* used to check for functions returning nonzeros */
70: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
71: PetscInt m; /* number of local elements in vectors */
72: Vec x; /* variables vector */
73: Vec xl,xu; /* bounds vectors */
74: PetscReal d1000 = 1000;
75: PetscBool flg; /* A return variable when checking for user options */
76: TaoSolver tao; /* TaoSolver solver context */
78: TaoSolverTerminationReason reason;
79: KSP ksp;
80: AppCtx user; /* user-defined work context */
81: PetscReal zero=0.0; /* lower bound on all variables */
84:
85: /* Initialize PETSC and TAO */
86: PetscInitialize( &argc, &argv,(char *)0,help );
87: TaoInitialize( &argc, &argv,(char *)0,help );
89: /* Set the default values for the problem parameters */
90: user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
92: /* Check for any command line arguments that override defaults */
93: PetscOptionsGetInt(PETSC_NULL,"-mx",&user.nx,&flg);
94: PetscOptionsGetInt(PETSC_NULL,"-my",&user.ny,&flg);
95: PetscOptionsGetReal(PETSC_NULL,"-ecc",&user.ecc,&flg);
96: PetscOptionsGetReal(PETSC_NULL,"-b",&user.b,&flg);
99: PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem SHB-----\n");
100: PetscPrintf(PETSC_COMM_WORLD,"mx: %D, my: %D, ecc: %G \n\n",
101: user.nx,user.ny,user.ecc);
103: /* Let Petsc determine the grid division */
104: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
106: /*
107: A two dimensional distributed array will help define this problem,
108: which derives from an elliptic PDE on two dimensional domain. From
109: the distributed array, Create the vectors.
110: */
111: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
112: user.nx,user.ny,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
113: &user.dm);
115: /*
116: Extract global and local vectors from DM; the vector user.B is
117: used solely as work space for the evaluation of the function,
118: gradient, and Hessian. Duplicate for remaining vectors that are
119: the same types.
120: */
121: DMCreateGlobalVector(user.dm,&x); /* Solution */
122: VecDuplicate(x,&user.B); /* Linear objective */
125: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
126: VecGetLocalSize(x,&m);
127: DMGetMatrix(user.dm,MATAIJ,&user.A);
129: /* User defined function -- compute linear term of quadratic */
130: ComputeB(&user);
132: /* The TAO code begins here */
134: /*
135: Create the optimization solver, Petsc application
136: Suitable methods: "tao_gpcg","tao_bqpip","tao_tron","tao_blmvm"
137: */
138: TaoCreate(PETSC_COMM_WORLD,&tao);
139: TaoSetType(tao,"tao_blmvm");
142: /* Set the initial vector */
143: VecSet(x, zero);
144: TaoSetInitialVector(tao,x);
146: /* Set the user function, gradient, hessian evaluation routines and data structures */
147: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);
148:
149:
150: TaoSetHessianRoutine(tao,user.A,user.A,FormHessian,(void*)&user);
152: /* Set a routine that defines the bounds */
153: VecDuplicate(x,&xl);
154: VecDuplicate(x,&xu);
155: VecSet(xl, zero);
156: VecSet(xu, d1000);
157: TaoSetVariableBounds(tao,xl,xu);
159: TaoGetKSP(tao,&ksp);
160: if (ksp) {
161: KSPSetType(ksp,KSPCG);
162: }
164: PetscOptionsHasName(PETSC_NULL,"-testmonitor",&flg);
165: if (flg) {
166: TaoSetMonitor(tao,Monitor,&user,PETSC_NULL);
167: }
168: PetscOptionsHasName(PETSC_NULL,"-testconvergence",&flg);
169: if (flg) {
170: TaoSetConvergenceTest(tao,ConvergenceTest,&user);
171: }
173: /* Check for any tao command line options */
174: TaoSetFromOptions(tao);
176: /* Solve the bound constrained problem */
177: TaoSolve(tao);
179: TaoGetTerminationReason(tao,&reason);
180: if (reason <= 0)
181: PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
184: /* Free PETSc data structures */
185: VecDestroy(&x);
186: VecDestroy(&xl);
187: VecDestroy(&xu);
188: MatDestroy(&user.A);
189: VecDestroy(&user.B);
190: /* Free TAO data structures */
191: TaoDestroy(&tao);
193: DMDestroy(&user.dm);
195: TaoFinalize();
196: PetscFinalize();
198: return 0;
199: }
202: static PetscReal p(PetscReal xi, PetscReal ecc)
203: {
204: PetscReal t=1.0+ecc*PetscCosScalar(xi);
205: return (t*t*t);
206: }
210: PetscErrorCode ComputeB(AppCtx* user)
211: {
213: PetscInt i,j,k;
214: PetscInt nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
215: PetscReal two=2.0, pi=4.0*atan(1.0);
216: PetscReal hx,hy,ehxhy;
217: PetscReal temp,*b;
218: PetscReal ecc=user->ecc;
220: nx=user->nx;
221: ny=user->ny;
222: hx=two*pi/(nx+1.0);
223: hy=two*user->b/(ny+1.0);
224: ehxhy = ecc*hx*hy;
227: /*
228: Get local grid boundaries
229: */
230: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
231: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
232:
234: /* Compute the linear term in the objective function */
235: VecGetArray(user->B,&b);
236: for (i=xs; i<xs+xm; i++){
237: temp=PetscSinScalar((i+1)*hx);
238: for (j=ys; j<ys+ym; j++){
239: k=xm*(j-ys)+(i-xs);
240: b[k]= - ehxhy*temp;
241: }
242: }
243: VecRestoreArray(user->B,&b);
244: PetscLogFlops(5*xm*ym+3*xm);
246: return 0;
247: }
251: PetscErrorCode FormFunctionGradient(TaoSolver tao, Vec X, PetscReal *fcn,Vec G,void *ptr)
252: {
253: AppCtx* user=(AppCtx*)ptr;
255: PetscInt i,j,k,kk;
256: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
257: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
258: PetscReal hx,hy,hxhy,hxhx,hyhy;
259: PetscReal xi,v[5];
260: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
261: PetscReal vmiddle, vup, vdown, vleft, vright;
262: PetscReal tt,f1,f2;
263: PetscReal *x,*g,zero=0.0;
264: Vec localX;
266: nx=user->nx;
267: ny=user->ny;
268: hx=two*pi/(nx+1.0);
269: hy=two*user->b/(ny+1.0);
270: hxhy=hx*hy;
271: hxhx=one/(hx*hx);
272: hyhy=one/(hy*hy);
274: DMGetLocalVector(user->dm,&localX);
276: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
277: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
279: VecSet(G, zero);
280: /*
281: Get local grid boundaries
282: */
283: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
284: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
285:
286: VecGetArray(localX,&x);
287: VecGetArray(G,&g);
289: for (i=xs; i< xs+xm; i++){
290: xi=(i+1)*hx;
291: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
292: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
293: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
294: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
295: trule5=trule1; /* L(i,j-1) */
296: trule6=trule2; /* U(i,j+1) */
298: vdown=-(trule5+trule2)*hyhy;
299: vleft=-hxhx*(trule2+trule4);
300: vright= -hxhx*(trule1+trule3);
301: vup=-hyhy*(trule1+trule6);
302: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
304: for (j=ys; j<ys+ym; j++){
305:
306: row=(j-gys)*gxm + (i-gxs);
307: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
308:
309: k=0;
310: if (j>gys){
311: v[k]=vdown; col[k]=row - gxm; k++;
312: }
313:
314: if (i>gxs){
315: v[k]= vleft; col[k]=row - 1; k++;
316: }
318: v[k]= vmiddle; col[k]=row; k++;
319:
320: if (i+1 < gxs+gxm){
321: v[k]= vright; col[k]=row+1; k++;
322: }
323:
324: if (j+1 <gys+gym){
325: v[k]= vup; col[k] = row+gxm; k++;
326: }
327: tt=0;
328: for (kk=0;kk<k;kk++){
329: tt+=v[kk]*x[col[kk]];
330: }
331: row=(j-ys)*xm + (i-xs);
332: g[row]=tt;
334: }
336: }
338: VecRestoreArray(localX,&x);
339: VecRestoreArray(G,&g);
341: DMRestoreLocalVector(user->dm,&localX);
343: VecDot(X,G,&f1);
344: VecDot(user->B,X,&f2);
345: VecAXPY(G, one, user->B);
346: *fcn = f1/2.0 + f2;
347:
349: PetscLogFlops((91 + 10*ym) * xm);
350: return 0;
352: }
357: /*
358: FormHessian computes the quadratic term in the quadratic objective function
359: Notice that the objective function in this problem is quadratic (therefore a constant
360: hessian). If using a nonquadratic solver, then you might want to reconsider this function
361: */
362: PetscErrorCode FormHessian(TaoSolver tao,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
363: {
364: AppCtx* user=(AppCtx*)ptr;
366: PetscInt i,j,k;
367: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
368: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
369: PetscReal hx,hy,hxhy,hxhx,hyhy;
370: PetscReal xi,v[5];
371: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
372: PetscReal vmiddle, vup, vdown, vleft, vright;
373: Mat hes=*H;
374: PetscBool assembled;
376: nx=user->nx;
377: ny=user->ny;
378: hx=two*pi/(nx+1.0);
379: hy=two*user->b/(ny+1.0);
380: hxhy=hx*hy;
381: hxhx=one/(hx*hx);
382: hyhy=one/(hy*hy);
384: *flg=SAME_NONZERO_PATTERN;
385: /*
386: Get local grid boundaries
387: */
388: DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
389: DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);
390:
391: MatAssembled(hes,&assembled);
392: if (assembled){MatZeroEntries(hes); }
394: for (i=xs; i< xs+xm; i++){
395: xi=(i+1)*hx;
396: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
397: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
398: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
399: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
400: trule5=trule1; /* L(i,j-1) */
401: trule6=trule2; /* U(i,j+1) */
403: vdown=-(trule5+trule2)*hyhy;
404: vleft=-hxhx*(trule2+trule4);
405: vright= -hxhx*(trule1+trule3);
406: vup=-hyhy*(trule1+trule6);
407: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
408: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
410: for (j=ys; j<ys+ym; j++){
411: row=(j-gys)*gxm + (i-gxs);
412:
413: k=0;
414: if (j>gys){
415: v[k]=vdown; col[k]=row - gxm; k++;
416: }
417:
418: if (i>gxs){
419: v[k]= vleft; col[k]=row - 1; k++;
420: }
422: v[k]= vmiddle; col[k]=row; k++;
423:
424: if (i+1 < gxs+gxm){
425: v[k]= vright; col[k]=row+1; k++;
426: }
427:
428: if (j+1 <gys+gym){
429: v[k]= vup; col[k] = row+gxm; k++;
430: }
431: MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES);
432:
433: }
435: }
437: /*
438: Assemble matrix, using the 2-step process:
439: MatAssemblyBegin(), MatAssemblyEnd().
440: By placing code between these two statements, computations can be
441: done while messages are in transition.
442: */
443: MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
444: MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);
446: /*
447: Tell the matrix we will never add a new nonzero location to the
448: matrix. If we do it will generate an error.
449: */
450: MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
451: MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);
453: PetscLogFlops(9*xm*ym+49*xm);
454: MatNorm(hes,NORM_1,&hx);
455: return 0;
456: }
461: PetscErrorCode Monitor(TaoSolver tao, void *ctx)
462: {
464: PetscInt its;
465: PetscReal f,gnorm,cnorm,xdiff;
466: TaoSolverTerminationReason reason;
468: TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
469: if (!(its%5)) {
470: PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%G\n",its,f);
471: }
472: return(0);
473: }
477: PetscErrorCode ConvergenceTest(TaoSolver tao, void *ctx)
478: {
480: PetscInt its;
481: PetscReal f,gnorm,cnorm,xdiff;
482: TaoSolverTerminationReason reason;
484: TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
485: if (its == 100) {
486: TaoSetTerminationReason(tao,TAO_DIVERGED_MAXITS);
487: }
488: return(0);
489:
490: }