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: /*
25: User-defined application context - contains data needed by the
26: application-provided call-back routines, FormFunctionGradient(),
27: FormHessian().
28: */
29: typedef struct {
30: /* problem parameters */
31: PetscReal ecc; /* test problem parameter */
32: PetscReal b; /* A dimension of journal bearing */
33: PetscInt nx,ny; /* discretization in x, y directions */
35: /* Working space */
36: DM dm; /* distributed array data structure */
37: Mat A; /* Quadratic Objective term */
38: Vec B; /* Linear Objective term */
39: } AppCtx;
41: /* User-defined routines */
42: static PetscReal p(PetscReal xi, PetscReal ecc);
43: static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *,Vec,void *);
44: static PetscErrorCode FormHessian(Tao,Vec,Mat, Mat, void *);
45: static PetscErrorCode ComputeB(AppCtx*);
46: static PetscErrorCode Monitor(Tao, void*);
47: static PetscErrorCode ConvergenceTest(Tao, void*);
49: int main(int argc, char **argv)
50: {
51: PetscInt Nx, Ny; /* number of processors in x- and y- directions */
52: PetscInt m; /* number of local elements in vectors */
53: Vec x; /* variables vector */
54: Vec xl,xu; /* bounds vectors */
55: PetscReal d1000 = 1000;
56: PetscBool flg,testgetdiag; /* A return variable when checking for user options */
57: Tao tao; /* Tao solver context */
58: KSP ksp;
59: AppCtx user; /* user-defined work context */
60: PetscReal zero = 0.0; /* lower bound on all variables */
62: /* Initialize PETSC and TAO */
63: PetscInitialize(&argc, &argv,(char *)0,help);
65: /* Set the default values for the problem parameters */
66: user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
67: testgetdiag = PETSC_FALSE;
69: /* Check for any command line arguments that override defaults */
70: PetscOptionsGetInt(NULL,NULL,"-mx",&user.nx,&flg);
71: PetscOptionsGetInt(NULL,NULL,"-my",&user.ny,&flg);
72: PetscOptionsGetReal(NULL,NULL,"-ecc",&user.ecc,&flg);
73: PetscOptionsGetReal(NULL,NULL,"-b",&user.b,&flg);
74: PetscOptionsGetBool(NULL,NULL,"-test_getdiagonal",&testgetdiag,NULL);
76: PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem SHB-----\n");
77: PetscPrintf(PETSC_COMM_WORLD,"mx: %D, my: %D, ecc: %g \n\n",user.nx,user.ny,(double)user.ecc);
79: /* Let Petsc determine the grid division */
80: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
82: /*
83: A two dimensional distributed array will help define this problem,
84: which derives from an elliptic PDE on two dimensional domain. From
85: the distributed array, Create the vectors.
86: */
87: 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);
88: DMSetFromOptions(user.dm);
89: DMSetUp(user.dm);
91: /*
92: Extract global and local vectors from DM; the vector user.B is
93: used solely as work space for the evaluation of the function,
94: gradient, and Hessian. Duplicate for remaining vectors that are
95: the same types.
96: */
97: DMCreateGlobalVector(user.dm,&x); /* Solution */
98: VecDuplicate(x,&user.B); /* Linear objective */
100: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
101: VecGetLocalSize(x,&m);
102: DMCreateMatrix(user.dm,&user.A);
104: if (testgetdiag) {
105: MatSetOperation(user.A,MATOP_GET_DIAGONAL,NULL);
106: }
108: /* User defined function -- compute linear term of quadratic */
109: ComputeB(&user);
111: /* The TAO code begins here */
113: /*
114: Create the optimization solver
115: Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
116: */
117: TaoCreate(PETSC_COMM_WORLD,&tao);
118: TaoSetType(tao,TAOBLMVM);
120: /* Set the initial vector */
121: VecSet(x, zero);
122: TaoSetSolution(tao,x);
124: /* Set the user function, gradient, hessian evaluation routines and data structures */
125: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*) &user);
127: TaoSetHessian(tao,user.A,user.A,FormHessian,(void*)&user);
129: /* Set a routine that defines the bounds */
130: VecDuplicate(x,&xl);
131: VecDuplicate(x,&xu);
132: VecSet(xl, zero);
133: VecSet(xu, d1000);
134: TaoSetVariableBounds(tao,xl,xu);
136: TaoGetKSP(tao,&ksp);
137: if (ksp) {
138: KSPSetType(ksp,KSPCG);
139: }
141: PetscOptionsHasName(NULL,NULL,"-testmonitor",&flg);
142: if (flg) {
143: TaoSetMonitor(tao,Monitor,&user,NULL);
144: }
145: PetscOptionsHasName(NULL,NULL,"-testconvergence",&flg);
146: if (flg) {
147: TaoSetConvergenceTest(tao,ConvergenceTest,&user);
148: }
150: /* Check for any tao command line options */
151: TaoSetFromOptions(tao);
153: /* Solve the bound constrained problem */
154: TaoSolve(tao);
156: /* Free PETSc data structures */
157: VecDestroy(&x);
158: VecDestroy(&xl);
159: VecDestroy(&xu);
160: MatDestroy(&user.A);
161: VecDestroy(&user.B);
163: /* Free TAO data structures */
164: TaoDestroy(&tao);
165: DMDestroy(&user.dm);
166: PetscFinalize();
167: return 0;
168: }
170: static PetscReal p(PetscReal xi, PetscReal ecc)
171: {
172: PetscReal t=1.0+ecc*PetscCosScalar(xi);
173: return (t*t*t);
174: }
176: PetscErrorCode ComputeB(AppCtx* user)
177: {
178: PetscInt i,j,k;
179: PetscInt nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
180: PetscReal two=2.0, pi=4.0*atan(1.0);
181: PetscReal hx,hy,ehxhy;
182: PetscReal temp,*b;
183: PetscReal ecc=user->ecc;
185: nx=user->nx;
186: ny=user->ny;
187: hx=two*pi/(nx+1.0);
188: hy=two*user->b/(ny+1.0);
189: ehxhy = ecc*hx*hy;
191: /*
192: Get local grid boundaries
193: */
194: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
195: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
197: /* Compute the linear term in the objective function */
198: VecGetArray(user->B,&b);
199: for (i=xs; i<xs+xm; i++) {
200: temp=PetscSinScalar((i+1)*hx);
201: for (j=ys; j<ys+ym; j++) {
202: k=xm*(j-ys)+(i-xs);
203: b[k]= - ehxhy*temp;
204: }
205: }
206: VecRestoreArray(user->B,&b);
207: PetscLogFlops(5.0*xm*ym+3.0*xm);
208: return 0;
209: }
211: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *ptr)
212: {
213: AppCtx* user=(AppCtx*)ptr;
214: PetscInt i,j,k,kk;
215: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
216: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
217: PetscReal hx,hy,hxhy,hxhx,hyhy;
218: PetscReal xi,v[5];
219: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
220: PetscReal vmiddle, vup, vdown, vleft, vright;
221: PetscReal tt,f1,f2;
222: PetscReal *x,*g,zero=0.0;
223: Vec localX;
225: nx=user->nx;
226: ny=user->ny;
227: hx=two*pi/(nx+1.0);
228: hy=two*user->b/(ny+1.0);
229: hxhy=hx*hy;
230: hxhx=one/(hx*hx);
231: hyhy=one/(hy*hy);
233: DMGetLocalVector(user->dm,&localX);
235: DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
236: DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);
238: VecSet(G, zero);
239: /*
240: Get local grid boundaries
241: */
242: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
243: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
245: VecGetArray(localX,&x);
246: VecGetArray(G,&g);
248: for (i=xs; i< xs+xm; i++) {
249: xi=(i+1)*hx;
250: trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
251: trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
252: trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
253: trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
254: trule5=trule1; /* L(i,j-1) */
255: trule6=trule2; /* U(i,j+1) */
257: vdown=-(trule5+trule2)*hyhy;
258: vleft=-hxhx*(trule2+trule4);
259: vright= -hxhx*(trule1+trule3);
260: vup=-hyhy*(trule1+trule6);
261: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
263: for (j=ys; j<ys+ym; j++) {
265: row=(j-gys)*gxm + (i-gxs);
266: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
268: k=0;
269: if (j>gys) {
270: v[k]=vdown; col[k]=row - gxm; k++;
271: }
273: if (i>gxs) {
274: v[k]= vleft; col[k]=row - 1; k++;
275: }
277: v[k]= vmiddle; col[k]=row; k++;
279: if (i+1 < gxs+gxm) {
280: v[k]= vright; col[k]=row+1; k++;
281: }
283: if (j+1 <gys+gym) {
284: v[k]= vup; col[k] = row+gxm; k++;
285: }
286: tt=0;
287: for (kk=0;kk<k;kk++) {
288: tt+=v[kk]*x[col[kk]];
289: }
290: row=(j-ys)*xm + (i-xs);
291: g[row]=tt;
293: }
295: }
297: VecRestoreArray(localX,&x);
298: VecRestoreArray(G,&g);
300: DMRestoreLocalVector(user->dm,&localX);
302: VecDot(X,G,&f1);
303: VecDot(user->B,X,&f2);
304: VecAXPY(G, one, user->B);
305: *fcn = f1/2.0 + f2;
307: PetscLogFlops((91 + 10.0*ym) * xm);
308: return 0;
310: }
312: /*
313: FormHessian computes the quadratic term in the quadratic objective function
314: Notice that the objective function in this problem is quadratic (therefore a constant
315: hessian). If using a nonquadratic solver, then you might want to reconsider this function
316: */
317: PetscErrorCode FormHessian(Tao tao,Vec X,Mat hes, Mat Hpre, void *ptr)
318: {
319: AppCtx* user=(AppCtx*)ptr;
320: PetscInt i,j,k;
321: PetscInt col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
322: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
323: PetscReal hx,hy,hxhy,hxhx,hyhy;
324: PetscReal xi,v[5];
325: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
326: PetscReal vmiddle, vup, vdown, vleft, vright;
327: PetscBool assembled;
329: nx=user->nx;
330: ny=user->ny;
331: hx=two*pi/(nx+1.0);
332: hy=two*user->b/(ny+1.0);
333: hxhy=hx*hy;
334: hxhx=one/(hx*hx);
335: hyhy=one/(hy*hy);
337: /*
338: Get local grid boundaries
339: */
340: DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
341: DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
342: MatAssembled(hes,&assembled);
343: if (assembled) MatZeroEntries(hes);
345: for (i=xs; i< xs+xm; i++) {
346: xi=(i+1)*hx;
347: trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
348: trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
349: trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
350: trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
351: trule5=trule1; /* L(i,j-1) */
352: trule6=trule2; /* U(i,j+1) */
354: vdown=-(trule5+trule2)*hyhy;
355: vleft=-hxhx*(trule2+trule4);
356: vright= -hxhx*(trule1+trule3);
357: vup=-hyhy*(trule1+trule6);
358: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
359: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
361: for (j=ys; j<ys+ym; j++) {
362: row=(j-gys)*gxm + (i-gxs);
364: k=0;
365: if (j>gys) {
366: v[k]=vdown; col[k]=row - gxm; k++;
367: }
369: if (i>gxs) {
370: v[k]= vleft; col[k]=row - 1; k++;
371: }
373: v[k]= vmiddle; col[k]=row; k++;
375: if (i+1 < gxs+gxm) {
376: v[k]= vright; col[k]=row+1; k++;
377: }
379: if (j+1 <gys+gym) {
380: v[k]= vup; col[k] = row+gxm; k++;
381: }
382: MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES);
384: }
386: }
388: /*
389: Assemble matrix, using the 2-step process:
390: MatAssemblyBegin(), MatAssemblyEnd().
391: By placing code between these two statements, computations can be
392: done while messages are in transition.
393: */
394: MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
395: MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);
397: /*
398: Tell the matrix we will never add a new nonzero location to the
399: matrix. If we do it will generate an error.
400: */
401: MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
402: MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);
404: PetscLogFlops(9.0*xm*ym+49.0*xm);
405: return 0;
406: }
408: PetscErrorCode Monitor(Tao tao, void *ctx)
409: {
410: PetscInt its;
411: PetscReal f,gnorm,cnorm,xdiff;
412: TaoConvergedReason reason;
414: TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
415: if (!(its%5)) {
416: PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%g\n",its,(double)f);
417: }
418: return 0;
419: }
421: PetscErrorCode ConvergenceTest(Tao tao, void *ctx)
422: {
423: PetscInt its;
424: PetscReal f,gnorm,cnorm,xdiff;
425: TaoConvergedReason reason;
427: TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
428: if (its == 100) {
429: TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS);
430: }
431: return 0;
433: }
435: /*TEST
437: build:
438: requires: !complex
440: test:
441: args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gatol 1.e-5
442: requires: !single
444: test:
445: suffix: 2
446: nsize: 2
447: args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gatol 1.e-5
448: requires: !single
450: test:
451: suffix: 3
452: nsize: 2
453: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
454: requires: !single
456: test:
457: suffix: 4
458: nsize: 2
459: args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
460: output_file: output/jbearing2_3.out
461: requires: !single
463: test:
464: suffix: 5
465: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
466: requires: !single
468: test:
469: suffix: 6
470: args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
471: requires: !single
473: test:
474: suffix: 7
475: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5
476: requires: !single
478: test:
479: suffix: 8
480: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5
481: requires: !single
483: test:
484: suffix: 9
485: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5
486: requires: !single
488: test:
489: suffix: 10
490: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
491: requires: !single
493: test:
494: suffix: 11
495: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
496: requires: !single
498: test:
499: suffix: 12
500: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
501: requires: !single
503: test:
504: suffix: 13
505: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls
506: requires: !single
508: test:
509: suffix: 14
510: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type blmvm
511: requires: !single
513: test:
514: suffix: 15
515: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
516: requires: !single
518: test:
519: suffix: 16
520: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
521: requires: !single
523: test:
524: suffix: 17
525: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type scalar -tao_view
526: requires: !single
528: test:
529: suffix: 18
530: args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type none -tao_view
531: requires: !single
533: test:
534: suffix: 19
535: args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
536: requires: !single
538: test:
539: suffix: 20
540: args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
541: requires: !single
543: test:
544: suffix: 21
545: args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
546: requires: !single
547: TEST*/