1 #include <petscdmda.h>
2 #include <petsctao.h>
3
4 static char help[] =
5 "This example demonstrates use of the TAO package to \n\
6 solve an unconstrained minimization problem. This example is based on a \n\
7 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain, \n\
8 boundary values along the edges of the domain, and a plate represented by \n\
9 lower boundary conditions, the objective is to find the\n\
10 surface with the minimal area that satisfies the boundary conditions.\n\
11 The command line options are:\n\
12 -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13 -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
14 -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction\n\
15 -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction\n\
16 -bheight <ht>, where <ht> = height of the plate\n\
17 -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
18 for an average of the boundary conditions\n\n";
19
20 /*T
21 Concepts: TAO^Solving a bound constrained minimization problem
22 Routines: TaoCreate();
23 Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
24 Routines: TaoSetHessianRoutine();
25 Routines: TaoSetInitialVector();
26 Routines: TaoSetVariableBounds();
27 Routines: TaoSetFromOptions();
28 Routines: TaoSolve(); TaoView();
29 Routines: TaoDestroy();
30 Processors: n
31 T*/
32
33
34
35
36 /*
37 User-defined application context - contains data needed by the
38 application-provided call-back routines, FormFunctionGradient(),
39 FormHessian().
40 */
41 typedef struct {
42 /* problem parameters */
43 PetscReal bheight; /* Height of plate under the surface */
44 PetscInt mx, my; /* discretization in x, y directions */
45 PetscInt bmx,bmy; /* Size of plate under the surface */
46 Vec Bottom, Top, Left, Right; /* boundary values */
47
48 /* Working space */
49 Vec localX, localV; /* ghosted local vector */
50 DM dm; /* distributed array data structure */
51 Mat H;
52 } AppCtx;
53
54 /* -------- User-defined Routines --------- */
55
56 static PetscErrorCode MSA_BoundaryConditions(AppCtx*);
57 static PetscErrorCode MSA_InitialPoint(AppCtx*,Vec);
58 static PetscErrorCode MSA_Plate(Vec,Vec,void*);
59 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
60 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
61
62 /* For testing matrix free submatrices */
63 PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat, Mat,void*);
64 PetscErrorCode MyMatMult(Mat,Vec,Vec);
65
main(int argc,char ** argv)66 int main(int argc, char **argv)
67 {
68 PetscErrorCode ierr; /* used to check for functions returning nonzeros */
69 PetscInt Nx, Ny; /* number of processors in x- and y- directions */
70 PetscInt m, N; /* number of local and global elements in vectors */
71 Vec x,xl,xu; /* solution vector and bounds*/
72 PetscBool flg; /* A return variable when checking for user options */
73 Tao tao; /* Tao solver context */
74 ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
75 Mat H_shell; /* to test matrix-free submatrices */
76 AppCtx user; /* user-defined work context */
77
78 /* Initialize PETSc, TAO */
79 ierr = PetscInitialize(&argc, &argv,(char *)0,help);if (ierr) return ierr;
80
81 /* Specify default dimension of the problem */
82 user.mx = 10; user.my = 10; user.bheight=0.1;
83
84 /* Check for any command line arguments that override defaults */
85 ierr = PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);CHKERRQ(ierr);
86 ierr = PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);CHKERRQ(ierr);
87 ierr = PetscOptionsGetReal(NULL,NULL,"-bheight",&user.bheight,&flg);CHKERRQ(ierr);
88
89 user.bmx = user.mx/2; user.bmy = user.my/2;
90 ierr = PetscOptionsGetInt(NULL,NULL,"-bmx",&user.bmx,&flg);CHKERRQ(ierr);
91 ierr = PetscOptionsGetInt(NULL,NULL,"-bmy",&user.bmy,&flg);CHKERRQ(ierr);
92
93 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n---- Minimum Surface Area With Plate Problem -----\n");CHKERRQ(ierr);
94 ierr = PetscPrintf(PETSC_COMM_WORLD,"mx:%D, my:%D, bmx:%D, bmy:%D, height:%g\n",user.mx,user.my,user.bmx,user.bmy,(double)user.bheight);CHKERRQ(ierr);
95
96 /* Calculate any derived values from parameters */
97 N = user.mx*user.my;
98
99 /* Let Petsc determine the dimensions of the local vectors */
100 Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
101
102 /*
103 A two dimensional distributed array will help define this problem,
104 which derives from an elliptic PDE on two dimensional domain. From
105 the distributed array, Create the vectors.
106 */
107 ierr = DMDACreate2d(MPI_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);CHKERRQ(ierr);
108 ierr = DMSetFromOptions(user.dm);CHKERRQ(ierr);
109 ierr = DMSetUp(user.dm);CHKERRQ(ierr);
110 /*
111 Extract global and local vectors from DM; The local vectors are
112 used solely as work space for the evaluation of the function,
113 gradient, and Hessian. Duplicate for remaining vectors that are
114 the same types.
115 */
116 ierr = DMCreateGlobalVector(user.dm,&x);CHKERRQ(ierr); /* Solution */
117 ierr = DMCreateLocalVector(user.dm,&user.localX);CHKERRQ(ierr);
118 ierr = VecDuplicate(user.localX,&user.localV);CHKERRQ(ierr);
119
120 ierr = VecDuplicate(x,&xl);CHKERRQ(ierr);
121 ierr = VecDuplicate(x,&xu);CHKERRQ(ierr);
122
123 /* The TAO code begins here */
124
125 /*
126 Create TAO solver and set desired solution method
127 The method must either be TAOTRON or TAOBLMVM
128 If TAOBLMVM is used, then hessian function is not called.
129 */
130 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
131 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr);
132
133 /* Set initial solution guess; */
134 ierr = MSA_BoundaryConditions(&user);CHKERRQ(ierr);
135 ierr = MSA_InitialPoint(&user,x);CHKERRQ(ierr);
136 ierr = TaoSetInitialVector(tao,x);CHKERRQ(ierr);
137
138 /* Set routines for function, gradient and hessian evaluation */
139 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);CHKERRQ(ierr);
140
141 ierr = VecGetLocalSize(x,&m);CHKERRQ(ierr);
142 ierr = MatCreateAIJ(MPI_COMM_WORLD,m,m,N,N,7,NULL,3,NULL,&(user.H));CHKERRQ(ierr);
143 ierr = MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
144
145 ierr = DMGetLocalToGlobalMapping(user.dm,&isltog);CHKERRQ(ierr);
146 ierr = MatSetLocalToGlobalMapping(user.H,isltog,isltog);CHKERRQ(ierr);
147 ierr = PetscOptionsHasName(NULL,NULL,"-matrixfree",&flg);CHKERRQ(ierr);
148 if (flg) {
149 ierr = MatCreateShell(PETSC_COMM_WORLD,m,m,N,N,(void*)&user,&H_shell);CHKERRQ(ierr);
150 ierr = MatShellSetOperation(H_shell,MATOP_MULT,(void(*)(void))MyMatMult);CHKERRQ(ierr);
151 ierr = MatSetOption(H_shell,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
152 ierr = TaoSetHessianRoutine(tao,H_shell,H_shell,MatrixFreeHessian,(void*)&user);CHKERRQ(ierr);
153 } else {
154 ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void*)&user);CHKERRQ(ierr);
155 }
156
157 /* Set Variable bounds */
158 ierr = MSA_Plate(xl,xu,(void*)&user);CHKERRQ(ierr);
159 ierr = TaoSetVariableBounds(tao,xl,xu);CHKERRQ(ierr);
160
161 /* Check for any tao command line options */
162 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
163
164 /* SOLVE THE APPLICATION */
165 ierr = TaoSolve(tao);CHKERRQ(ierr);
166
167 ierr = TaoView(tao,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
168
169 /* Free TAO data structures */
170 ierr = TaoDestroy(&tao);CHKERRQ(ierr);
171
172 /* Free PETSc data structures */
173 ierr = VecDestroy(&x);CHKERRQ(ierr);
174 ierr = VecDestroy(&xl);CHKERRQ(ierr);
175 ierr = VecDestroy(&xu);CHKERRQ(ierr);
176 ierr = MatDestroy(&user.H);CHKERRQ(ierr);
177 ierr = VecDestroy(&user.localX);CHKERRQ(ierr);
178 ierr = VecDestroy(&user.localV);CHKERRQ(ierr);
179 ierr = VecDestroy(&user.Bottom);CHKERRQ(ierr);
180 ierr = VecDestroy(&user.Top);CHKERRQ(ierr);
181 ierr = VecDestroy(&user.Left);CHKERRQ(ierr);
182 ierr = VecDestroy(&user.Right);CHKERRQ(ierr);
183 ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
184 if (flg) {
185 ierr = MatDestroy(&H_shell);CHKERRQ(ierr);
186 }
187 ierr = PetscFinalize();
188 return ierr;
189 }
190
191 /* FormFunctionGradient - Evaluates f(x) and gradient g(x).
192
193 Input Parameters:
194 . tao - the Tao context
195 . X - input vector
196 . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
197
198 Output Parameters:
199 . fcn - the function value
200 . G - vector containing the newly evaluated gradient
201
202 Notes:
203 In this case, we discretize the domain and Create triangles. The
204 surface of each triangle is planar, whose surface area can be easily
205 computed. The total surface area is found by sweeping through the grid
206 and computing the surface area of the two triangles that have their
207 right angle at the grid point. The diagonal line segments on the
208 grid that define the triangles run from top left to lower right.
209 The numbering of points starts at the lower left and runs left to
210 right, then bottom to top.
211 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * fcn,Vec G,void * userCtx)212 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G,void *userCtx)
213 {
214 AppCtx *user = (AppCtx *) userCtx;
215 PetscErrorCode ierr;
216 PetscInt i,j,row;
217 PetscInt mx=user->mx, my=user->my;
218 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
219 PetscReal ft=0;
220 PetscReal zero=0.0;
221 PetscReal hx=1.0/(mx+1),hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy, area=0.5*hx*hy;
222 PetscReal rhx=mx+1, rhy=my+1;
223 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
224 PetscReal df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc;
225 PetscReal *g, *x,*left,*right,*bottom,*top;
226 Vec localX = user->localX, localG = user->localV;
227
228 /* Get local mesh boundaries */
229 ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
230 ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);
231
232 /* Scatter ghost points to local vector */
233 ierr = DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
234 ierr = DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
235
236 /* Initialize vector to zero */
237 ierr = VecSet(localG, zero);CHKERRQ(ierr);
238
239 /* Get pointers to vector data */
240 ierr = VecGetArray(localX,&x);CHKERRQ(ierr);
241 ierr = VecGetArray(localG,&g);CHKERRQ(ierr);
242 ierr = VecGetArray(user->Top,&top);CHKERRQ(ierr);
243 ierr = VecGetArray(user->Bottom,&bottom);CHKERRQ(ierr);
244 ierr = VecGetArray(user->Left,&left);CHKERRQ(ierr);
245 ierr = VecGetArray(user->Right,&right);CHKERRQ(ierr);
246
247 /* Compute function over the locally owned part of the mesh */
248 for (j=ys; j<ys+ym; j++){
249 for (i=xs; i< xs+xm; i++){
250 row=(j-gys)*gxm + (i-gxs);
251
252 xc = x[row];
253 xlt=xrb=xl=xr=xb=xt=xc;
254
255 if (i==0){ /* left side */
256 xl= left[j-ys+1];
257 xlt = left[j-ys+2];
258 } else {
259 xl = x[row-1];
260 }
261
262 if (j==0){ /* bottom side */
263 xb=bottom[i-xs+1];
264 xrb = bottom[i-xs+2];
265 } else {
266 xb = x[row-gxm];
267 }
268
269 if (i+1 == gxs+gxm){ /* right side */
270 xr=right[j-ys+1];
271 xrb = right[j-ys];
272 } else {
273 xr = x[row+1];
274 }
275
276 if (j+1==gys+gym){ /* top side */
277 xt=top[i-xs+1];
278 xlt = top[i-xs];
279 }else {
280 xt = x[row+gxm];
281 }
282
283 if (i>gxs && j+1<gys+gym){
284 xlt = x[row-1+gxm];
285 }
286 if (j>gys && i+1<gxs+gxm){
287 xrb = x[row+1-gxm];
288 }
289
290 d1 = (xc-xl);
291 d2 = (xc-xr);
292 d3 = (xc-xt);
293 d4 = (xc-xb);
294 d5 = (xr-xrb);
295 d6 = (xrb-xb);
296 d7 = (xlt-xl);
297 d8 = (xt-xlt);
298
299 df1dxc = d1*hydhx;
300 df2dxc = (d1*hydhx + d4*hxdhy);
301 df3dxc = d3*hxdhy;
302 df4dxc = (d2*hydhx + d3*hxdhy);
303 df5dxc = d2*hydhx;
304 df6dxc = d4*hxdhy;
305
306 d1 *= rhx;
307 d2 *= rhx;
308 d3 *= rhy;
309 d4 *= rhy;
310 d5 *= rhy;
311 d6 *= rhx;
312 d7 *= rhy;
313 d8 *= rhx;
314
315 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
316 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
317 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
318 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
319 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
320 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
321
322 ft = ft + (f2 + f4);
323
324 df1dxc /= f1;
325 df2dxc /= f2;
326 df3dxc /= f3;
327 df4dxc /= f4;
328 df5dxc /= f5;
329 df6dxc /= f6;
330
331 g[row] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc) * 0.5;
332
333 }
334 }
335
336
337 /* Compute triangular areas along the border of the domain. */
338 if (xs==0){ /* left side */
339 for (j=ys; j<ys+ym; j++){
340 d3=(left[j-ys+1] - left[j-ys+2])*rhy;
341 d2=(left[j-ys+1] - x[(j-gys)*gxm])*rhx;
342 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
343 }
344 }
345 if (ys==0){ /* bottom side */
346 for (i=xs; i<xs+xm; i++){
347 d2=(bottom[i+1-xs]-bottom[i-xs+2])*rhx;
348 d3=(bottom[i-xs+1]-x[i-gxs])*rhy;
349 ft = ft+PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
350 }
351 }
352
353 if (xs+xm==mx){ /* right side */
354 for (j=ys; j< ys+ym; j++){
355 d1=(x[(j+1-gys)*gxm-1]-right[j-ys+1])*rhx;
356 d4=(right[j-ys]-right[j-ys+1])*rhy;
357 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
358 }
359 }
360 if (ys+ym==my){ /* top side */
361 for (i=xs; i<xs+xm; i++){
362 d1=(x[(gym-1)*gxm + i-gxs] - top[i-xs+1])*rhy;
363 d4=(top[i-xs+1] - top[i-xs])*rhx;
364 ft = ft+PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
365 }
366 }
367
368 if (ys==0 && xs==0){
369 d1=(left[0]-left[1])*rhy;
370 d2=(bottom[0]-bottom[1])*rhx;
371 ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
372 }
373 if (ys+ym == my && xs+xm == mx){
374 d1=(right[ym+1] - right[ym])*rhy;
375 d2=(top[xm+1] - top[xm])*rhx;
376 ft +=PetscSqrtScalar(1.0 + d1*d1 + d2*d2);
377 }
378
379 ft=ft*area;
380 ierr = MPI_Allreduce(&ft,fcn,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);CHKERRQ(ierr);
381
382
383 /* Restore vectors */
384 ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);
385 ierr = VecRestoreArray(localG,&g);CHKERRQ(ierr);
386 ierr = VecRestoreArray(user->Left,&left);CHKERRQ(ierr);
387 ierr = VecRestoreArray(user->Top,&top);CHKERRQ(ierr);
388 ierr = VecRestoreArray(user->Bottom,&bottom);CHKERRQ(ierr);
389 ierr = VecRestoreArray(user->Right,&right);CHKERRQ(ierr);
390
391 /* Scatter values to global vector */
392 ierr = DMLocalToGlobalBegin(user->dm,localG,INSERT_VALUES,G);CHKERRQ(ierr);
393 ierr = DMLocalToGlobalEnd(user->dm,localG,INSERT_VALUES,G);CHKERRQ(ierr);
394
395 ierr = PetscLogFlops(70.0*xm*ym);CHKERRQ(ierr);
396
397 return 0;
398 }
399
400 /* ------------------------------------------------------------------- */
401 /*
402 FormHessian - Evaluates Hessian matrix.
403
404 Input Parameters:
405 . tao - the Tao context
406 . x - input vector
407 . ptr - optional user-defined context, as set by TaoSetHessianRoutine()
408
409 Output Parameters:
410 . A - Hessian matrix
411 . B - optionally different preconditioning matrix
412
413 Notes:
414 Due to mesh point reordering with DMs, we must always work
415 with the local mesh points, and then transform them to the new
416 global numbering with the local-to-global mapping. We cannot work
417 directly with the global numbers for the original uniprocessor mesh!
418
419 Two methods are available for imposing this transformation
420 when setting matrix entries:
421 (A) MatSetValuesLocal(), using the local ordering (including
422 ghost points!)
423 - Do the following two steps once, before calling TaoSolve()
424 - Use DMGetISLocalToGlobalMapping() to extract the
425 local-to-global map from the DM
426 - Associate this map with the matrix by calling
427 MatSetLocalToGlobalMapping()
428 - Then set matrix entries using the local ordering
429 by calling MatSetValuesLocal()
430 (B) MatSetValues(), using the global ordering
431 - Use DMGetGlobalIndices() to extract the local-to-global map
432 - Then apply this map explicitly yourself
433 - Set matrix entries using the global ordering by calling
434 MatSetValues()
435 Option (A) seems cleaner/easier in many cases, and is the procedure
436 used in this example.
437 */
FormHessian(Tao tao,Vec X,Mat Hptr,Mat Hessian,void * ptr)438 PetscErrorCode FormHessian(Tao tao,Vec X,Mat Hptr, Mat Hessian, void *ptr)
439 {
440 PetscErrorCode ierr;
441 AppCtx *user = (AppCtx *) ptr;
442 PetscInt i,j,k,row;
443 PetscInt mx=user->mx, my=user->my;
444 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym,col[7];
445 PetscReal hx=1.0/(mx+1), hy=1.0/(my+1), hydhx=hy/hx, hxdhy=hx/hy;
446 PetscReal rhx=mx+1, rhy=my+1;
447 PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb;
448 PetscReal hl,hr,ht,hb,hc,htl,hbr;
449 PetscReal *x,*left,*right,*bottom,*top;
450 PetscReal v[7];
451 Vec localX = user->localX;
452 PetscBool assembled;
453
454
455 /* Set various matrix options */
456 ierr = MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
457
458 /* Initialize matrix entries to zero */
459 ierr = MatAssembled(Hessian,&assembled);CHKERRQ(ierr);
460 if (assembled){ierr = MatZeroEntries(Hessian);CHKERRQ(ierr);}
461
462 /* Get local mesh boundaries */
463 ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
464 ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);
465
466 /* Scatter ghost points to local vector */
467 ierr = DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
468 ierr = DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);CHKERRQ(ierr);
469
470 /* Get pointers to vector data */
471 ierr = VecGetArray(localX,&x);CHKERRQ(ierr);
472 ierr = VecGetArray(user->Top,&top);CHKERRQ(ierr);
473 ierr = VecGetArray(user->Bottom,&bottom);CHKERRQ(ierr);
474 ierr = VecGetArray(user->Left,&left);CHKERRQ(ierr);
475 ierr = VecGetArray(user->Right,&right);CHKERRQ(ierr);
476
477 /* Compute Hessian over the locally owned part of the mesh */
478
479 for (i=xs; i< xs+xm; i++){
480
481 for (j=ys; j<ys+ym; j++){
482
483 row=(j-gys)*gxm + (i-gxs);
484
485 xc = x[row];
486 xlt=xrb=xl=xr=xb=xt=xc;
487
488 /* Left side */
489 if (i==gxs){
490 xl= left[j-ys+1];
491 xlt = left[j-ys+2];
492 } else {
493 xl = x[row-1];
494 }
495
496 if (j==gys){
497 xb=bottom[i-xs+1];
498 xrb = bottom[i-xs+2];
499 } else {
500 xb = x[row-gxm];
501 }
502
503 if (i+1 == gxs+gxm){
504 xr=right[j-ys+1];
505 xrb = right[j-ys];
506 } else {
507 xr = x[row+1];
508 }
509
510 if (j+1==gys+gym){
511 xt=top[i-xs+1];
512 xlt = top[i-xs];
513 }else {
514 xt = x[row+gxm];
515 }
516
517 if (i>gxs && j+1<gys+gym){
518 xlt = x[row-1+gxm];
519 }
520 if (j>gys && i+1<gxs+gxm){
521 xrb = x[row+1-gxm];
522 }
523
524
525 d1 = (xc-xl)*rhx;
526 d2 = (xc-xr)*rhx;
527 d3 = (xc-xt)*rhy;
528 d4 = (xc-xb)*rhy;
529 d5 = (xrb-xr)*rhy;
530 d6 = (xrb-xb)*rhx;
531 d7 = (xlt-xl)*rhy;
532 d8 = (xlt-xt)*rhx;
533
534 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7);
535 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4);
536 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8);
537 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2);
538 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5);
539 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6);
540
541
542 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+
543 (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2);
544 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+
545 (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4);
546 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+
547 (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4);
548 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+
549 (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2);
550
551 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6);
552 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3);
553
554 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) +
555 hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) +
556 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2*d1*d4)/(f2*f2*f2) +
557 (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4);
558
559 hl*=0.5; hr*=0.5; ht*=0.5; hb*=0.5; hbr*=0.5; htl*=0.5; hc*=0.5;
560
561 k=0;
562 if (j>0){
563 v[k]=hb; col[k]=row - gxm; k++;
564 }
565
566 if (j>0 && i < mx -1){
567 v[k]=hbr; col[k]=row - gxm+1; k++;
568 }
569
570 if (i>0){
571 v[k]= hl; col[k]=row - 1; k++;
572 }
573
574 v[k]= hc; col[k]=row; k++;
575
576 if (i < mx-1){
577 v[k]= hr; col[k]=row+1; k++;
578 }
579
580 if (i>0 && j < my-1){
581 v[k]= htl; col[k] = row+gxm-1; k++;
582 }
583
584 if (j < my-1){
585 v[k]= ht; col[k] = row+gxm; k++;
586 }
587
588 /*
589 Set matrix values using local numbering, which was defined
590 earlier, in the main routine.
591 */
592 ierr = MatSetValuesLocal(Hessian,1,&row,k,col,v,INSERT_VALUES);CHKERRQ(ierr);
593
594 }
595 }
596
597 /* Restore vectors */
598 ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);
599 ierr = VecRestoreArray(user->Left,&left);CHKERRQ(ierr);
600 ierr = VecRestoreArray(user->Top,&top);CHKERRQ(ierr);
601 ierr = VecRestoreArray(user->Bottom,&bottom);CHKERRQ(ierr);
602 ierr = VecRestoreArray(user->Right,&right);CHKERRQ(ierr);
603
604 /* Assemble the matrix */
605 ierr = MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
606 ierr = MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
607
608 ierr = PetscLogFlops(199.0*xm*ym);CHKERRQ(ierr);
609 return 0;
610 }
611
612 /* ------------------------------------------------------------------- */
613 /*
614 MSA_BoundaryConditions - Calculates the boundary conditions for
615 the region.
616
617 Input Parameter:
618 . user - user-defined application context
619
620 Output Parameter:
621 . user - user-defined application context
622 */
MSA_BoundaryConditions(AppCtx * user)623 static PetscErrorCode MSA_BoundaryConditions(AppCtx * user)
624 {
625 int ierr;
626 PetscInt i,j,k,maxits=5,limit=0;
627 PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym;
628 PetscInt mx=user->mx,my=user->my;
629 PetscInt bsize=0, lsize=0, tsize=0, rsize=0;
630 PetscReal one=1.0, two=2.0, three=3.0, scl=1.0, tol=1e-10;
631 PetscReal fnorm,det,hx,hy,xt=0,yt=0;
632 PetscReal u1,u2,nf1,nf2,njac11,njac12,njac21,njac22;
633 PetscReal b=-0.5, t=0.5, l=-0.5, r=0.5;
634 PetscReal *boundary;
635 PetscBool flg;
636 Vec Bottom,Top,Right,Left;
637
638 /* Get local mesh boundaries */
639 ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
640 ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);
641
642
643 bsize=xm+2;
644 lsize=ym+2;
645 rsize=ym+2;
646 tsize=xm+2;
647
648 ierr = VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,&Bottom);CHKERRQ(ierr);
649 ierr = VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,&Top);CHKERRQ(ierr);
650 ierr = VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,&Left);CHKERRQ(ierr);
651 ierr = VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,&Right);CHKERRQ(ierr);
652
653 user->Top=Top;
654 user->Left=Left;
655 user->Bottom=Bottom;
656 user->Right=Right;
657
658 hx= (r-l)/(mx+1); hy=(t-b)/(my+1);
659
660 for (j=0; j<4; j++){
661 if (j==0){
662 yt=b;
663 xt=l+hx*xs;
664 limit=bsize;
665 VecGetArray(Bottom,&boundary);
666 } else if (j==1){
667 yt=t;
668 xt=l+hx*xs;
669 limit=tsize;
670 VecGetArray(Top,&boundary);
671 } else if (j==2){
672 yt=b+hy*ys;
673 xt=l;
674 limit=lsize;
675 VecGetArray(Left,&boundary);
676 } else if (j==3){
677 yt=b+hy*ys;
678 xt=r;
679 limit=rsize;
680 VecGetArray(Right,&boundary);
681 }
682
683 for (i=0; i<limit; i++){
684 u1=xt;
685 u2=-yt;
686 for (k=0; k<maxits; k++){
687 nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt;
688 nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt;
689 fnorm=PetscSqrtScalar(nf1*nf1+nf2*nf2);
690 if (fnorm <= tol) break;
691 njac11=one+u2*u2-u1*u1;
692 njac12=two*u1*u2;
693 njac21=-two*u1*u2;
694 njac22=-one - u1*u1 + u2*u2;
695 det = njac11*njac22-njac21*njac12;
696 u1 = u1-(njac22*nf1-njac12*nf2)/det;
697 u2 = u2-(njac11*nf2-njac21*nf1)/det;
698 }
699
700 boundary[i]=u1*u1-u2*u2;
701 if (j==0 || j==1) {
702 xt=xt+hx;
703 } else if (j==2 || j==3){
704 yt=yt+hy;
705 }
706 }
707 if (j==0){
708 ierr = VecRestoreArray(Bottom,&boundary);CHKERRQ(ierr);
709 } else if (j==1){
710 ierr = VecRestoreArray(Top,&boundary);CHKERRQ(ierr);
711 } else if (j==2){
712 ierr = VecRestoreArray(Left,&boundary);CHKERRQ(ierr);
713 } else if (j==3){
714 ierr = VecRestoreArray(Right,&boundary);CHKERRQ(ierr);
715 }
716 }
717
718 /* Scale the boundary if desired */
719
720 ierr = PetscOptionsGetReal(NULL,NULL,"-bottom",&scl,&flg);CHKERRQ(ierr);
721 if (flg){
722 ierr = VecScale(Bottom, scl);CHKERRQ(ierr);
723 }
724 ierr = PetscOptionsGetReal(NULL,NULL,"-top",&scl,&flg);CHKERRQ(ierr);
725 if (flg){
726 ierr = VecScale(Top, scl);CHKERRQ(ierr);
727 }
728 ierr = PetscOptionsGetReal(NULL,NULL,"-right",&scl,&flg);CHKERRQ(ierr);
729 if (flg){
730 ierr = VecScale(Right, scl);CHKERRQ(ierr);
731 }
732
733 ierr = PetscOptionsGetReal(NULL,NULL,"-left",&scl,&flg);CHKERRQ(ierr);
734 if (flg){
735 ierr = VecScale(Left, scl);CHKERRQ(ierr);
736 }
737 return 0;
738 }
739
740
741 /* ------------------------------------------------------------------- */
742 /*
743 MSA_Plate - Calculates an obstacle for surface to stretch over.
744
745 Input Parameter:
746 . user - user-defined application context
747
748 Output Parameter:
749 . user - user-defined application context
750 */
MSA_Plate(Vec XL,Vec XU,void * ctx)751 static PetscErrorCode MSA_Plate(Vec XL,Vec XU,void *ctx){
752
753 AppCtx *user=(AppCtx *)ctx;
754 PetscErrorCode ierr;
755 PetscInt i,j,row;
756 PetscInt xs,ys,xm,ym;
757 PetscInt mx=user->mx, my=user->my, bmy, bmx;
758 PetscReal t1,t2,t3;
759 PetscReal *xl, lb=PETSC_NINFINITY, ub=PETSC_INFINITY;
760 PetscBool cylinder;
761
762 user->bmy = PetscMax(0,user->bmy);user->bmy = PetscMin(my,user->bmy);
763 user->bmx = PetscMax(0,user->bmx);user->bmx = PetscMin(mx,user->bmx);
764 bmy=user->bmy; bmx=user->bmx;
765
766 ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
767
768 ierr = VecSet(XL, lb);CHKERRQ(ierr);
769 ierr = VecSet(XU, ub);CHKERRQ(ierr);
770
771 ierr = VecGetArray(XL,&xl);CHKERRQ(ierr);
772
773 ierr = PetscOptionsHasName(NULL,NULL,"-cylinder",&cylinder);CHKERRQ(ierr);
774 /* Compute the optional lower box */
775 if (cylinder){
776 for (i=xs; i< xs+xm; i++){
777 for (j=ys; j<ys+ym; j++){
778 row=(j-ys)*xm + (i-xs);
779 t1=(2.0*i-mx)*bmy;
780 t2=(2.0*j-my)*bmx;
781 t3=bmx*bmx*bmy*bmy;
782 if (t1*t1 + t2*t2 <= t3){
783 xl[row] = user->bheight;
784 }
785 }
786 }
787 } else {
788 /* Compute the optional lower box */
789 for (i=xs; i< xs+xm; i++){
790 for (j=ys; j<ys+ym; j++){
791 row=(j-ys)*xm + (i-xs);
792 if (i>=(mx-bmx)/2 && i<mx-(mx-bmx)/2 &&
793 j>=(my-bmy)/2 && j<my-(my-bmy)/2){
794 xl[row] = user->bheight;
795 }
796 }
797 }
798 }
799 ierr = VecRestoreArray(XL,&xl);CHKERRQ(ierr);
800
801 return 0;
802 }
803
804
805 /* ------------------------------------------------------------------- */
806 /*
807 MSA_InitialPoint - Calculates the initial guess in one of three ways.
808
809 Input Parameters:
810 . user - user-defined application context
811 . X - vector for initial guess
812
813 Output Parameters:
814 . X - newly computed initial guess
815 */
MSA_InitialPoint(AppCtx * user,Vec X)816 static PetscErrorCode MSA_InitialPoint(AppCtx * user, Vec X)
817 {
818 PetscErrorCode ierr;
819 PetscInt start=-1,i,j;
820 PetscReal zero=0.0;
821 PetscBool flg;
822
823 ierr = PetscOptionsGetInt(NULL,NULL,"-start",&start,&flg);CHKERRQ(ierr);
824 if (flg && start==0){ /* The zero vector is reasonable */
825 ierr = VecSet(X, zero);CHKERRQ(ierr);
826 } else if (flg && start>0){ /* Try a random start between -0.5 and 0.5 */
827 PetscRandom rctx; PetscReal np5=-0.5;
828
829 ierr = PetscRandomCreate(MPI_COMM_WORLD,&rctx);CHKERRQ(ierr);
830 for (i=0; i<start; i++){
831 ierr = VecSetRandom(X, rctx);CHKERRQ(ierr);
832 }
833 ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
834 ierr = VecShift(X, np5);CHKERRQ(ierr);
835
836 } else { /* Take an average of the boundary conditions */
837
838 PetscInt row,xs,xm,gxs,gxm,ys,ym,gys,gym;
839 PetscInt mx=user->mx,my=user->my;
840 PetscReal *x,*left,*right,*bottom,*top;
841 Vec localX = user->localX;
842
843 /* Get local mesh boundaries */
844 ierr = DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
845 ierr = DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);
846
847 /* Get pointers to vector data */
848 ierr = VecGetArray(user->Top,&top);CHKERRQ(ierr);
849 ierr = VecGetArray(user->Bottom,&bottom);CHKERRQ(ierr);
850 ierr = VecGetArray(user->Left,&left);CHKERRQ(ierr);
851 ierr = VecGetArray(user->Right,&right);CHKERRQ(ierr);
852
853 ierr = VecGetArray(localX,&x);CHKERRQ(ierr);
854 /* Perform local computations */
855 for (j=ys; j<ys+ym; j++){
856 for (i=xs; i< xs+xm; i++){
857 row=(j-gys)*gxm + (i-gxs);
858 x[row] = ((j+1)*bottom[i-xs+1]/my + (my-j+1)*top[i-xs+1]/(my+2)+(i+1)*left[j-ys+1]/mx + (mx-i+1)*right[j-ys+1]/(mx+2))/2.0;
859 }
860 }
861
862 /* Restore vectors */
863 ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr);
864
865 ierr = VecRestoreArray(user->Left,&left);CHKERRQ(ierr);
866 ierr = VecRestoreArray(user->Top,&top);CHKERRQ(ierr);
867 ierr = VecRestoreArray(user->Bottom,&bottom);CHKERRQ(ierr);
868 ierr = VecRestoreArray(user->Right,&right);CHKERRQ(ierr);
869
870 /* Scatter values into global vector */
871 ierr = DMLocalToGlobalBegin(user->dm,localX,INSERT_VALUES,X);CHKERRQ(ierr);
872 ierr = DMLocalToGlobalEnd(user->dm,localX,INSERT_VALUES,X);CHKERRQ(ierr);
873
874 }
875 return 0;
876 }
877
878 /* For testing matrix free submatrices */
MatrixFreeHessian(Tao tao,Vec x,Mat H,Mat Hpre,void * ptr)879 PetscErrorCode MatrixFreeHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ptr)
880 {
881 PetscErrorCode ierr;
882 AppCtx *user = (AppCtx*)ptr;
883 PetscFunctionBegin;
884 ierr = FormHessian(tao,x,user->H,user->H,ptr);CHKERRQ(ierr);
885 PetscFunctionReturn(0);
886 }
MyMatMult(Mat H_shell,Vec X,Vec Y)887 PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
888 {
889 PetscErrorCode ierr;
890 void *ptr;
891 AppCtx *user;
892 PetscFunctionBegin;
893 ierr = MatShellGetContext(H_shell,&ptr);CHKERRQ(ierr);
894 user = (AppCtx*)ptr;
895 ierr = MatMult(user->H,X,Y);CHKERRQ(ierr);
896 PetscFunctionReturn(0);
897 }
898
899
900 /*TEST
901
902 build:
903 requires: !complex
904
905 test:
906 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type tron -tao_gatol 1.e-5
907 requires: !single
908
909 test:
910 suffix: 2
911 nsize: 2
912 args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_type blmvm -tao_gatol 1.e-4
913 requires: !single
914
915 test:
916 suffix: 3
917 nsize: 3
918 args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_type tron -tao_gatol 1.e-5
919 requires: !single
920
921 test:
922 suffix: 4
923 nsize: 3
924 args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type mask -tao_type tron -tao_gatol 1.e-5
925 requires: !single
926
927 test:
928 suffix: 5
929 nsize: 3
930 args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -pc_type none -tao_type tron -tao_gatol 1.e-5
931 requires: !single
932
933 test:
934 suffix: 6
935 nsize: 3
936 args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -matrixfree -tao_type blmvm -tao_gatol 1.e-4
937 requires: !single
938
939 test:
940 suffix: 7
941 nsize: 3
942 args: -tao_smonitor -mx 8 -my 12 -bmx 4 -bmy 10 -bheight 0.1 -tao_subset_type matrixfree -pc_type none -tao_type gpcg -tao_gatol 1.e-5
943 requires: !single
944
945 test:
946 suffix: 8
947 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
948 requires: !single
949
950 test:
951 suffix: 9
952 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bncg -tao_gatol 1e-4
953 requires: !single
954
955 test:
956 suffix: 10
957 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5
958 requires: !single
959
960 test:
961 suffix: 11
962 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5
963 requires: !single
964
965 test:
966 suffix: 12
967 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5
968 requires: !single
969
970 test:
971 suffix: 13
972 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
973 requires: !single
974
975 test:
976 suffix: 14
977 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
978 requires: !single
979
980 test:
981 suffix: 15
982 args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
983 requires: !single
984
985 test:
986 suffix: 16
987 args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnls
988 requires: !single
989
990 test:
991 suffix: 17
992 args: -tao_smonitor -mx 8 -my 8 -bmx 2 -bmy 5 -bheight 0.3 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
993 requires: !single
994
995 TEST*/
996