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