Actual source code: minsurf1.c
1: #include <petsctao.h>
3: static char help[] =
4: "This example demonstrates use of the TAO package to\n\
5: solve an unconstrained system of equations. This example is based on a\n\
6: problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\
7: boundary values along the edges of the domain, the objective is to find the\n\
8: surface with the minimal area that satisfies the boundary conditions.\n\
9: This application solves this problem using complimentarity -- We are actually\n\
10: solving the system (grad f)_i >= 0, if x_i == l_i \n\
11: (grad f)_i = 0, if l_i < x_i < u_i \n\
12: (grad f)_i <= 0, if x_i == u_i \n\
13: where f is the function to be minimized. \n\
14: \n\
15: The command line options are:\n\
16: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
17: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
18: -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
20: /*T
21: Concepts: TAO^Solving a complementarity problem
22: Routines: TaoCreate(); TaoDestroy();
24: Processors: 1
25: T*/
27: /*
28: User-defined application context - contains data needed by the
29: application-provided call-back routines, FormFunctionGradient(),
30: FormHessian().
31: */
32: typedef struct {
33: PetscInt mx, my;
34: PetscReal *bottom, *top, *left, *right;
35: } AppCtx;
37: /* -------- User-defined Routines --------- */
39: static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
40: static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
41: PetscErrorCode FormConstraints(Tao, Vec, Vec, void *);
42: PetscErrorCode FormJacobian(Tao, Vec, Mat, Mat, void *);
44: int main(int argc, char **argv)
45: {
47: Vec x; /* solution vector */
48: Vec c; /* Constraints function vector */
49: Vec xl,xu; /* Bounds on the variables */
50: PetscBool flg; /* A return variable when checking for user options */
51: Tao tao; /* TAO solver context */
52: Mat J; /* Jacobian matrix */
53: PetscInt N; /* Number of elements in vector */
54: PetscScalar lb = PETSC_NINFINITY; /* lower bound constant */
55: PetscScalar ub = PETSC_INFINITY; /* upper bound constant */
56: AppCtx user; /* user-defined work context */
58: /* Initialize PETSc, TAO */
59: PetscInitialize(&argc, &argv, (char *)0, help);if (ierr) return ierr;
61: /* Specify default dimension of the problem */
62: user.mx = 4; user.my = 4;
64: /* Check for any command line arguments that override defaults */
65: PetscOptionsGetInt(NULL,NULL, "-mx", &user.mx, &flg);
66: PetscOptionsGetInt(NULL,NULL, "-my", &user.my, &flg);
68: /* Calculate any derived values from parameters */
69: N = user.mx*user.my;
71: PetscPrintf(PETSC_COMM_SELF,"\n---- Minimum Surface Area Problem -----\n");
72: PetscPrintf(PETSC_COMM_SELF,"mx:%D, my:%D\n", user.mx,user.my);
74: /* Create appropriate vectors and matrices */
75: VecCreateSeq(MPI_COMM_SELF, N, &x);
76: VecDuplicate(x, &c);
77: MatCreateSeqAIJ(MPI_COMM_SELF, N, N, 7, NULL, &J);
79: /* The TAO code begins here */
81: /* Create TAO solver and set desired solution method */
82: TaoCreate(PETSC_COMM_SELF,&tao);
83: TaoSetType(tao,TAOSSILS);
85: /* Set data structure */
86: TaoSetInitialVector(tao, x);
88: /* Set routines for constraints function and Jacobian evaluation */
89: TaoSetConstraintsRoutine(tao, c, FormConstraints, (void *)&user);
90: TaoSetJacobianRoutine(tao, J, J, FormJacobian, (void *)&user);
92: /* Set the variable bounds */
93: MSA_BoundaryConditions(&user);
95: /* Set initial solution guess */
96: MSA_InitialPoint(&user, x);
98: /* Set Bounds on variables */
99: VecDuplicate(x, &xl);
100: VecDuplicate(x, &xu);
101: VecSet(xl, lb);
102: VecSet(xu, ub);
103: TaoSetVariableBounds(tao,xl,xu);
105: /* Check for any tao command line options */
106: TaoSetFromOptions(tao);
108: /* Solve the application */
109: TaoSolve(tao);
111: /* Free Tao data structures */
112: TaoDestroy(&tao);
114: /* Free PETSc data structures */
115: VecDestroy(&x);
116: VecDestroy(&xl);
117: VecDestroy(&xu);
118: VecDestroy(&c);
119: MatDestroy(&J);
121: /* Free user-created data structures */
122: PetscFree(user.bottom);
123: PetscFree(user.top);
124: PetscFree(user.left);
125: PetscFree(user.right);
127: PetscFinalize();
128: return ierr;
129: }
131: /* -------------------------------------------------------------------- */
133: /* FormConstraints - Evaluates gradient of f.
135: Input Parameters:
136: . tao - the TAO_APPLICATION context
137: . X - input vector
138: . ptr - optional user-defined context, as set by TaoSetConstraintsRoutine()
140: Output Parameters:
141: . G - vector containing the newly evaluated gradient
142: */
143: PetscErrorCode FormConstraints(Tao tao, Vec X, Vec G, void *ptr)
144: {
145: AppCtx *user = (AppCtx *) ptr;
147: PetscInt i,j,row;
148: PetscInt mx=user->mx, my=user->my;
149: PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
150: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
151: PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
152: PetscScalar zero=0.0;
153: PetscScalar *g, *x;
156: /* Initialize vector to zero */
157: VecSet(G, zero);
159: /* Get pointers to vector data */
160: VecGetArray(X, &x);
161: VecGetArray(G, &g);
163: /* Compute function over the locally owned part of the mesh */
164: for (j=0; j<my; j++) {
165: for (i=0; i< mx; i++) {
166: row= j*mx + i;
168: xc = x[row];
169: xlt=xrb=xl=xr=xb=xt=xc;
171: if (i==0) { /* left side */
172: xl= user->left[j+1];
173: xlt = user->left[j+2];
174: } else {
175: xl = x[row-1];
176: }
178: if (j==0) { /* bottom side */
179: xb=user->bottom[i+1];
180: xrb = user->bottom[i+2];
181: } else {
182: xb = x[row-mx];
183: }
185: if (i+1 == mx) { /* right side */
186: xr=user->right[j+1];
187: xrb = user->right[j];
188: } else {
189: xr = x[row+1];
190: }
192: if (j+1==0+my) { /* top side */
193: xt=user->top[i+1];
194: xlt = user->top[i];
195: }else {
196: xt = x[row+mx];
197: }
199: if (i>0 && j+1<my) {
200: xlt = x[row-1+mx];
201: }
202: if (j>0 && i+1<mx) {
203: xrb = x[row+1-mx];
204: }
206: d1 = (xc-xl);
207: d2 = (xc-xr);
208: d3 = (xc-xt);
209: d4 = (xc-xb);
210: d5 = (xr-xrb);
211: d6 = (xrb-xb);
212: d7 = (xlt-xl);
213: d8 = (xt-xlt);
215: df1dxc = d1*hydhx;
216: df2dxc = (d1*hydhx + d4*hxdhy);
217: df3dxc = d3*hxdhy;
218: df4dxc = (d2*hydhx + d3*hxdhy);
219: df5dxc = d2*hydhx;
220: df6dxc = d4*hxdhy;
222: d1 /= hx;
223: d2 /= hx;
224: d3 /= hy;
225: d4 /= hy;
226: d5 /= hy;
227: d6 /= hx;
228: d7 /= hy;
229: d8 /= hx;
231: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
232: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
233: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
234: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
235: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
236: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
238: df1dxc /= f1;
239: df2dxc /= f2;
240: df3dxc /= f3;
241: df4dxc /= f4;
242: df5dxc /= f5;
243: df6dxc /= f6;
245: g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc)/2.0;
246: }
247: }
249: /* Restore vectors */
250: VecRestoreArray(X, &x);
251: VecRestoreArray(G, &g);
252: PetscLogFlops(67*mx*my);
253: return(0);
254: }
256: /* ------------------------------------------------------------------- */
257: /*
258: FormJacobian - Evaluates Jacobian matrix.
260: Input Parameters:
261: . tao - the TAO_APPLICATION context
262: . X - input vector
263: . ptr - optional user-defined context, as set by TaoSetJacobian()
265: Output Parameters:
266: . tH - Jacobian matrix
268: */
269: PetscErrorCode FormJacobian(Tao tao, Vec X, Mat H, Mat tHPre, void *ptr)
270: {
271: AppCtx *user = (AppCtx *) ptr;
272: PetscErrorCode ierr;
273: PetscInt i,j,k,row;
274: PetscInt mx=user->mx, my=user->my;
275: PetscInt col[7];
276: PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
277: PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
278: PetscReal hl,hr,ht,hb,hc,htl,hbr;
279: const PetscScalar *x;
280: PetscScalar v[7];
281: PetscBool assembled;
283: /* Set various matrix options */
285: MatSetOption(H,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);
286: MatAssembled(H,&assembled);
287: if (assembled) {MatZeroEntries(H);}
289: /* Get pointers to vector data */
290: VecGetArrayRead(X, &x);
292: /* Compute Jacobian over the locally owned part of the mesh */
293: for (i=0; i< mx; i++) {
294: for (j=0; j<my; j++) {
295: row= j*mx + i;
297: xc = x[row];
298: xlt=xrb=xl=xr=xb=xt=xc;
300: /* Left side */
301: if (i==0) {
302: xl = user->left[j+1];
303: xlt = user->left[j+2];
304: } else {
305: xl = x[row-1];
306: }
308: if (j==0) {
309: xb = user->bottom[i+1];
310: xrb = user->bottom[i+2];
311: } else {
312: xb = x[row-mx];
313: }
315: if (i+1 == mx) {
316: xr = user->right[j+1];
317: xrb = user->right[j];
318: } else {
319: xr = x[row+1];
320: }
322: if (j+1==my) {
323: xt = user->top[i+1];
324: xlt = user->top[i];
325: }else {
326: xt = x[row+mx];
327: }
329: if (i>0 && j+1<my) {
330: xlt = x[row-1+mx];
331: }
332: if (j>0 && i+1<mx) {
333: xrb = x[row+1-mx];
334: }
336: d1 = (xc-xl)/hx;
337: d2 = (xc-xr)/hx;
338: d3 = (xc-xt)/hy;
339: d4 = (xc-xb)/hy;
340: d5 = (xrb-xr)/hy;
341: d6 = (xrb-xb)/hx;
342: d7 = (xlt-xl)/hy;
343: d8 = (xlt-xt)/hx;
345: f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
346: f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
347: f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
348: f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
349: f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
350: f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
352: hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+(-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
353: hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+(-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
354: ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+(-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
355: hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+(-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
357: hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
358: htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
360: hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
361: (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
363: hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0;
365: k=0;
366: if (j>0) {
367: v[k]=hb; col[k]=row - mx; k++;
368: }
370: if (j>0 && i < mx -1) {
371: v[k]=hbr; col[k]=row - mx+1; k++;
372: }
374: if (i>0) {
375: v[k]= hl; col[k]=row - 1; k++;
376: }
378: v[k]= hc; col[k]=row; k++;
380: if (i < mx-1) {
381: v[k]= hr; col[k]=row+1; k++;
382: }
384: if (i>0 && j < my-1) {
385: v[k]= htl; col[k] = row+mx-1; k++;
386: }
388: if (j < my-1) {
389: v[k]= ht; col[k] = row+mx; k++;
390: }
392: /*
393: Set matrix values using local numbering, which was defined
394: earlier, in the main routine.
395: */
396: MatSetValues(H,1,&row,k,col,v,INSERT_VALUES);
397: }
398: }
400: /* Restore vectors */
401: VecRestoreArrayRead(X,&x);
403: /* Assemble the matrix */
404: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
405: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
406: PetscLogFlops(199*mx*my);
407: return(0);
408: }
410: /* ------------------------------------------------------------------- */
411: /*
412: MSA_BoundaryConditions - Calculates the boundary conditions for
413: the region.
415: Input Parameter:
416: . user - user-defined application context
418: Output Parameter:
419: . user - user-defined application context
420: */
421: static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
422: {
423: PetscErrorCode ierr;
424: PetscInt i,j,k,limit=0,maxits=5;
425: PetscInt mx=user->mx,my=user->my;
426: PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
427: PetscReal one=1.0, two=2.0, three=3.0, tol=1e-10;
428: PetscReal fnorm,det,hx,hy,xt=0,yt=0;
429: PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
430: PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
431: PetscReal *boundary;
434: bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2;
436: PetscMalloc1(bsize, &user->bottom);
437: PetscMalloc1(tsize, &user->top);
438: PetscMalloc1(lsize, &user->left);
439: PetscMalloc1(rsize, &user->right);
441: hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
443: for (j=0; j<4; j++) {
444: if (j==0) {
445: yt=b;
446: xt=l;
447: limit=bsize;
448: boundary=user->bottom;
449: } else if (j==1) {
450: yt=t;
451: xt=l;
452: limit=tsize;
453: boundary=user->top;
454: } else if (j==2) {
455: yt=b;
456: xt=l;
457: limit=lsize;
458: boundary=user->left;
459: } else { /* if (j==3) */
460: yt=b;
461: xt=r;
462: limit=rsize;
463: boundary=user->right;
464: }
466: for (i=0; i<limit; i++) {
467: u1=xt;
468: u2=-yt;
469: for (k=0; k<maxits; k++) {
470: nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
471: nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
472: fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
473: if (fnorm <= tol) break;
474: njac11=one+u2*u2-u1*u1;
475: njac12=two*u1*u2;
476: njac21=-two*u1*u2;
477: njac22=-one - u1*u1 + u2*u2;
478: det = njac11*njac22-njac21*njac12;
479: u1 = u1-(njac22*nf1-njac12*nf2)/det;
480: u2 = u2-(njac11*nf2-njac21*nf1)/det;
481: }
483: boundary[i]=u1*u1-u2*u2;
484: if (j==0 || j==1) {
485: xt=xt+hx;
486: } else { /* if (j==2 || j==3) */
487: yt=yt+hy;
488: }
489: }
490: }
491: return(0);
492: }
494: /* ------------------------------------------------------------------- */
495: /*
496: MSA_InitialPoint - Calculates the initial guess in one of three ways.
498: Input Parameters:
499: . user - user-defined application context
500: . X - vector for initial guess
502: Output Parameters:
503: . X - newly computed initial guess
504: */
505: static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
506: {
508: PetscInt start=-1,i,j;
509: PetscScalar zero=0.0;
510: PetscBool flg;
513: PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);
515: if (flg && start==0) { /* The zero vector is reasonable */
516: VecSet(X, zero);
517: } else { /* Take an average of the boundary conditions */
518: PetscInt row;
519: PetscInt mx=user->mx,my=user->my;
520: PetscScalar *x;
522: /* Get pointers to vector data */
523: VecGetArray(X,&x);
525: /* Perform local computations */
526: for (j=0; j<my; j++) {
527: for (i=0; i< mx; i++) {
528: row=(j)*mx + (i);
529: x[row] = (((j+1)*user->bottom[i+1]+(my-j+1)*user->top[i+1])/(my+2)+ ((i+1)*user->left[j+1]+(mx-i+1)*user->right[j+1])/(mx+2))/2.0;
530: }
531: }
533: /* Restore vectors */
534: VecRestoreArray(X,&x);
535: }
536: return(0);
537: }
539: /*TEST
541: build:
542: requires: !complex
544: test:
545: args: -tao_monitor -tao_view -tao_type ssils -tao_gttol 1.e-5
546: requires: !single
548: test:
549: suffix: 2
550: args: -tao_monitor -tao_view -tao_type ssfls -tao_gttol 1.e-5
552: TEST*/