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[]=
12: "This example demonstrates use of the TAO package to \n\
13: solve a bound constrained minimization problem. This example is based on \n\
14: the problem DPJB from the MINPACK-2 test suite. This pressure journal \n\
15: bearing problem is an example of elliptic variational problem defined over \n\
16: a two dimensional rectangle. By discretizing the domain into triangular \n\
17: elements, the pressure surrounding the journal bearing is defined as the \n\
18: minimum of a quadratic function whose variables are bounded below by zero.\n\
19: The command line options are:\n\
20: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
21: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
22: \n";
24: /*T
25: Concepts: TAO^Solving a bound constrained minimization problem
26: Routines: TaoCreate();
27: Routines: TaoSetType(); TaoSetObjectiveAndGradient();
28: Routines: TaoSetHessian();
29: Routines: TaoSetVariableBounds();
30: Routines: TaoSetMonitor(); TaoSetConvergenceTest();
31: Routines: TaoSetSolution();
32: Routines: TaoSetFromOptions();
33: Routines: TaoSolve();
34: Routines: TaoDestroy();
35: Processors: n
36: T*/
38: /*
39: User-defined application context - contains data needed by the
40: application-provided call-back routines, FormFunctionGradient(),
41: FormHessian().
42: */
43: typedef struct {
44: /* problem parameters */
45: PetscReal ecc; /* test problem parameter */
46: PetscReal b; /* A dimension of journal bearing */
47: PetscInt nx,ny; /* discretization in x, y directions */
49: /* Working space */
50: DM dm; /* distributed array data structure */
51: Mat A; /* Quadratic Objective term */
52: Vec B; /* Linear Objective term */
53: } AppCtx;
55: /* User-defined routines */
56: static PetscReal p(PetscReal xi, PetscReal ecc);
57: static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *,Vec,void *);
58: static PetscErrorCode FormHessian(Tao,Vec,Mat, Mat, void *);
59: static PetscErrorCode ComputeB(AppCtx*);
60: static PetscErrorCode Monitor(Tao, void*);
61: static PetscErrorCode ConvergenceTest(Tao, void*);
63: int main(int argc, char **argv)
64: {
65: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
66: PetscInt m; /* number of local elements in vectors */
67: Vec x; /* variables vector */
68: Vec xl,xu; /* bounds vectors */
69: PetscReal d1000 = 1000;
70: PetscBool flg,testgetdiag; /* A return variable when checking for user options */
71: Tao tao; /* Tao solver context */
72: KSP ksp;
73: AppCtx user; /* user-defined work context */
74: PetscReal zero = 0.0; /* lower bound on all variables */
76: /* Initialize PETSC and TAO */
77: PetscInitialize(&argc, &argv,(char *)0,help);
79: /* Set the default values for the problem parameters */
80: user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
81: testgetdiag = PETSC_FALSE;
83: /* Check for any command line arguments that override defaults */
84: PetscOptionsGetInt(NULL,NULL,"-mx",&user.nx,&flg);
85: PetscOptionsGetInt(NULL,NULL,"-my",&user.ny,&flg);
86: PetscOptionsGetReal(NULL,NULL,"-ecc",&user.ecc,&flg);
87: PetscOptionsGetReal(NULL,NULL,"-b",&user.b,&flg);
88: PetscOptionsGetBool(NULL,NULL,"-test_getdiagonal",&testgetdiag,NULL);
90: PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem SHB-----\n");
91: PetscPrintf(PETSC_COMM_WORLD,"mx: %D, my: %D, ecc: %g \n\n",user.nx,user.ny,(double)user.ecc);
93: /* Let Petsc determine the grid division */
94: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
96: /*
97: A two dimensional distributed array will help define this problem,
98: which derives from an elliptic PDE on two dimensional domain. From
99: the distributed array, Create the vectors.
100: */
101: 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);
102: DMSetFromOptions(user.dm);
103: DMSetUp(user.dm);
105: /*
106: Extract global and local vectors from DM; the vector user.B is
107: used solely as work space for the evaluation of the function,
108: gradient, and Hessian. Duplicate for remaining vectors that are
109: the same types.
110: */
111: DMCreateGlobalVector(user.dm,&x); /* Solution */
112: VecDuplicate(x,&user.B); /* Linear objective */
114: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
115: VecGetLocalSize(x,&m);
116: DMCreateMatrix(user.dm,&user.A);
118: if (testgetdiag) {
119: MatSetOperation(user.A,MATOP_GET_DIAGONAL,NULL);
120: }
122: /* User defined function -- compute linear term of quadratic */
123: ComputeB(&user);
125: /* The TAO code begins here */
127: /*
128: Create the optimization solver
129: Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
130: */
131: TaoCreate(PETSC_COMM_WORLD,&tao);
132: TaoSetType(tao,TAOBLMVM);
134: /* Set the initial vector */
135: VecSet(x, zero);
136: TaoSetSolution(tao,x);
138: /* Set the user function, gradient, hessian evaluation routines and data structures */
139: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*) &user);
141: TaoSetHessian(tao,user.A,user.A,FormHessian,(void*)&user);
143: /* Set a routine that defines the bounds */
144: VecDuplicate(x,&xl);
145: VecDuplicate(x,&xu);
146: VecSet(xl, zero);
147: VecSet(xu, d1000);
148: TaoSetVariableBounds(tao,xl,xu);
150: TaoGetKSP(tao,&ksp);
151: if (ksp) {
152: KSPSetType(ksp,KSPCG);
153: }
155: PetscOptionsHasName(NULL,NULL,"-testmonitor",&flg);
156: if (flg) {
157: TaoSetMonitor(tao,Monitor,&user,NULL);
158: }
159: PetscOptionsHasName(NULL,NULL,"-testconvergence",&flg);
160: if (flg) {
161: TaoSetConvergenceTest(tao,ConvergenceTest,&user);
162: }
164: /* Check for any tao command line options */
165: TaoSetFromOptions(tao);
167: /* Solve the bound constrained problem */
168: TaoSolve(tao);
170: /* Free PETSc data structures */
171: VecDestroy(&x);
172: VecDestroy(&xl);
173: VecDestroy(&xu);
174: MatDestroy(&user.A);
175: VecDestroy(&user.B);
177: /* Free TAO data structures */
178: TaoDestroy(&tao);
179: DMDestroy(&user.dm);
180: PetscFinalize();
181: return 0;
182: }
184: static PetscReal p(PetscReal xi, PetscReal ecc)
185: {
186: PetscReal t=1.0+ecc*PetscCosScalar(xi);
187: return (t*t*t);
188: }
190: PetscErrorCode ComputeB(AppCtx* user)
191: {
192: PetscInt i,j,k;
193: PetscInt nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
194: PetscReal two=2.0, pi=4.0*atan(1.0);
195: PetscReal hx,hy,ehxhy;
196: PetscReal temp,*b;
197: PetscReal ecc=user->ecc;
199: nx=user->nx;
200: ny=user->ny;
201: hx=two*pi/(nx+1.0);
202: hy=two*user->b/(ny+1.0);
203: ehxhy = ecc*hx*hy;
205: /*
206: Get local grid boundaries
207: */
208: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
209: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
211: /* Compute the linear term in the objective function */
212: VecGetArray(user->B,&b);
213: for (i=xs; i<xs+xm; i++) {
214: temp=PetscSinScalar((i+1)*hx);
215: for (j=ys; j<ys+ym; j++) {
216: k=xm*(j-ys)+(i-xs);
217: b[k]= - ehxhy*temp;
218: }
219: }
220: VecRestoreArray(user->B,&b);
221: PetscLogFlops(5.0*xm*ym+3.0*xm);
222: return 0;
223: }
225: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *ptr)
226: {
227: AppCtx* user=(AppCtx*)ptr;
228: PetscInt i,j,k,kk;
229: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
230: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
231: PetscReal hx,hy,hxhy,hxhx,hyhy;
232: PetscReal xi,v[5];
233: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
234: PetscReal vmiddle, vup, vdown, vleft, vright;
235: PetscReal tt,f1,f2;
236: PetscReal *x,*g,zero=0.0;
237: Vec localX;
239: nx=user->nx;
240: ny=user->ny;
241: hx=two*pi/(nx+1.0);
242: hy=two*user->b/(ny+1.0);
243: hxhy=hx*hy;
244: hxhx=one/(hx*hx);
245: hyhy=one/(hy*hy);
247: DMGetLocalVector(user->dm,&localX);
249: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
250: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
252: VecSet(G, zero);
253: /*
254: Get local grid boundaries
255: */
256: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
257: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
259: VecGetArray(localX,&x);
260: VecGetArray(G,&g);
262: for (i=xs; i< xs+xm; i++) {
263: xi=(i+1)*hx;
264: trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
265: trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
266: trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
267: trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
268: trule5=trule1; /* L(i,j-1) */
269: trule6=trule2; /* U(i,j+1) */
271: vdown=-(trule5+trule2)*hyhy;
272: vleft=-hxhx*(trule2+trule4);
273: vright= -hxhx*(trule1+trule3);
274: vup=-hyhy*(trule1+trule6);
275: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
277: for (j=ys; j<ys+ym; j++) {
279: row=(j-gys)*gxm + (i-gxs);
280: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
282: k=0;
283: if (j>gys) {
284: v[k]=vdown; col[k]=row - gxm; k++;
285: }
287: if (i>gxs) {
288: v[k]= vleft; col[k]=row - 1; k++;
289: }
291: v[k]= vmiddle; col[k]=row; k++;
293: if (i+1 < gxs+gxm) {
294: v[k]= vright; col[k]=row+1; k++;
295: }
297: if (j+1 <gys+gym) {
298: v[k]= vup; col[k] = row+gxm; k++;
299: }
300: tt=0;
301: for (kk=0;kk<k;kk++) {
302: tt+=v[kk]*x[col[kk]];
303: }
304: row=(j-ys)*xm + (i-xs);
305: g[row]=tt;
307: }
309: }
311: VecRestoreArray(localX,&x);
312: VecRestoreArray(G,&g);
314: DMRestoreLocalVector(user->dm,&localX);
316: VecDot(X,G,&f1);
317: VecDot(user->B,X,&f2);
318: VecAXPY(G, one, user->B);
319: *fcn = f1/2.0 + f2;
321: PetscLogFlops((91 + 10.0*ym) * xm);
322: return 0;
324: }
326: /*
327: FormHessian computes the quadratic term in the quadratic objective function
328: Notice that the objective function in this problem is quadratic (therefore a constant
329: hessian). If using a nonquadratic solver, then you might want to reconsider this function
330: */
331: PetscErrorCode FormHessian(Tao tao,Vec X,Mat hes, Mat Hpre, void *ptr)
332: {
333: AppCtx* user=(AppCtx*)ptr;
334: PetscInt i,j,k;
335: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
336: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
337: PetscReal hx,hy,hxhy,hxhx,hyhy;
338: PetscReal xi,v[5];
339: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
340: PetscReal vmiddle, vup, vdown, vleft, vright;
341: PetscBool assembled;
343: nx=user->nx;
344: ny=user->ny;
345: hx=two*pi/(nx+1.0);
346: hy=two*user->b/(ny+1.0);
347: hxhy=hx*hy;
348: hxhx=one/(hx*hx);
349: hyhy=one/(hy*hy);
351: /*
352: Get local grid boundaries
353: */
354: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
355: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
356: MatAssembled(hes,&assembled);
357: if (assembled) MatZeroEntries(hes);
359: for (i=xs; i< xs+xm; i++) {
360: xi=(i+1)*hx;
361: trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
362: trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
363: trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
364: trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
365: trule5=trule1; /* L(i,j-1) */
366: trule6=trule2; /* U(i,j+1) */
368: vdown=-(trule5+trule2)*hyhy;
369: vleft=-hxhx*(trule2+trule4);
370: vright= -hxhx*(trule1+trule3);
371: vup=-hyhy*(trule1+trule6);
372: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
373: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
375: for (j=ys; j<ys+ym; j++) {
376: row=(j-gys)*gxm + (i-gxs);
378: k=0;
379: if (j>gys) {
380: v[k]=vdown; col[k]=row - gxm; k++;
381: }
383: if (i>gxs) {
384: v[k]= vleft; col[k]=row - 1; k++;
385: }
387: v[k]= vmiddle; col[k]=row; k++;
389: if (i+1 < gxs+gxm) {
390: v[k]= vright; col[k]=row+1; k++;
391: }
393: if (j+1 <gys+gym) {
394: v[k]= vup; col[k] = row+gxm; k++;
395: }
396: MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES);
398: }
400: }
402: /*
403: Assemble matrix, using the 2-step process:
404: MatAssemblyBegin(), MatAssemblyEnd().
405: By placing code between these two statements, computations can be
406: done while messages are in transition.
407: */
408: MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
409: MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);
411: /*
412: Tell the matrix we will never add a new nonzero location to the
413: matrix. If we do it will generate an error.
414: */
415: MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
416: MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);
418: PetscLogFlops(9.0*xm*ym+49.0*xm);
419: return 0;
420: }
422: PetscErrorCode Monitor(Tao tao, void *ctx)
423: {
424: PetscInt its;
425: PetscReal f,gnorm,cnorm,xdiff;
426: TaoConvergedReason reason;
428: TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
429: if (!(its%5)) {
430: PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%g\n",its,(double)f);
431: }
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) {
443: TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS);
444: }
445: return 0;
447: }
449: /*TEST
451: build:
452: requires: !complex
454: test:
455: args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gatol 1.e-5
456: requires: !single
458: test:
459: suffix: 2
460: nsize: 2
461: args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gatol 1.e-5
462: requires: !single
464: test:
465: suffix: 3
466: nsize: 2
467: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
468: requires: !single
470: test:
471: suffix: 4
472: nsize: 2
473: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
474: output_file: output/jbearing2_3.out
475: requires: !single
477: test:
478: suffix: 5
479: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
480: requires: !single
482: test:
483: suffix: 6
484: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
485: requires: !single
487: test:
488: suffix: 7
489: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5
490: requires: !single
492: test:
493: suffix: 8
494: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5
495: requires: !single
497: test:
498: suffix: 9
499: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5
500: requires: !single
502: test:
503: suffix: 10
504: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
505: requires: !single
507: test:
508: suffix: 11
509: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
510: requires: !single
512: test:
513: suffix: 12
514: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
515: requires: !single
517: test:
518: suffix: 13
519: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls
520: requires: !single
522: test:
523: suffix: 14
524: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type blmvm
525: requires: !single
527: test:
528: suffix: 15
529: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
530: requires: !single
532: test:
533: suffix: 16
534: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
535: requires: !single
537: test:
538: suffix: 17
539: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type scalar -tao_view
540: requires: !single
542: test:
543: suffix: 18
544: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type none -tao_view
545: requires: !single
547: test:
548: suffix: 19
549: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
550: requires: !single
552: test:
553: suffix: 20
554: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
555: requires: !single
557: test:
558: suffix: 21
559: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
560: requires: !single
561: TEST*/