1 /* Program usage: mpiexec -n <proc> minsurf2 [-help] [all TAO options] */
2
3 /*
4 Include "petsctao.h" so we can use TAO solvers.
5 petscdmda.h for distributed array
6 */
7 #include <petsctao.h>
8 #include <petscdmda.h>
9
10 static char help[] = "This example demonstrates use of the TAO package to \n\
11 solve an unconstrained minimization problem. This example is based on a \n\
12 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and \n\
13 boundary values along the edges of the domain, the objective is to find the\n\
14 surface with the minimal area that satisfies the boundary conditions.\n\
15 The command line options are:\n\
16 -da_grid_x <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
17 -da_grid_y <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
18 -start <st>, where <st> =0 for zero vector, <st> >0 for random start, and <st> <0 \n\
19 for an average of the boundary conditions\n\n";
20
21 /*
22 User-defined application context - contains data needed by the
23 application-provided call-back routines, FormFunction(),
24 FormFunctionGradient(), and FormHessian().
25 */
26 typedef struct {
27 PetscInt mx, my; /* discretization in x, y directions */
28 PetscReal *bottom, *top, *left, *right; /* boundary values */
29 DM dm; /* distributed array data structure */
30 Mat H; /* Hessian */
31 } AppCtx;
32
33 /* -------- User-defined Routines --------- */
34
35 static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
36 static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
37 static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
38 static PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
39 static PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
40 static PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
41 static PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
42 static PetscErrorCode My_Monitor(Tao, void *);
43
main(int argc,char ** argv)44 int main(int argc, char **argv)
45 {
46 Vec x; /* solution, gradient vectors */
47 PetscBool viewmat; /* flags */
48 PetscBool fddefault, fdcoloring; /* flags */
49 Tao tao; /* TAO solver context */
50 AppCtx user; /* user-defined work context */
51 ISColoring iscoloring;
52 MatFDColoring matfdcoloring;
53
54 /* Initialize TAO */
55 PetscFunctionBeginUser;
56 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
57
58 /* Create distributed array (DM) to manage parallel grid and vectors */
59 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &user.dm));
60 PetscCall(DMSetFromOptions(user.dm));
61 PetscCall(DMSetUp(user.dm));
62 PetscCall(DMDASetUniformCoordinates(user.dm, -0.5, 0.5, -0.5, 0.5, PETSC_DECIDE, PETSC_DECIDE));
63 PetscCall(DMDAGetInfo(user.dm, PETSC_IGNORE, &user.mx, &user.my, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
64
65 PetscCall(PetscPrintf(MPI_COMM_WORLD, "\n---- Minimum Surface Area Problem -----\n"));
66 PetscCall(PetscPrintf(MPI_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my));
67
68 /* Create TAO solver and set desired solution method.*/
69 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
70 PetscCall(TaoSetType(tao, TAOCG));
71
72 /*
73 Extract global vector from DA for the vector of variables -- PETSc routine
74 Compute the initial solution -- application specific, see below
75 Set this vector for use by TAO -- TAO routine
76 */
77 PetscCall(DMCreateGlobalVector(user.dm, &x));
78 PetscCall(MSA_BoundaryConditions(&user));
79 PetscCall(MSA_InitialPoint(&user, x));
80 PetscCall(TaoSetSolution(tao, x));
81
82 /*
83 Initialize the Application context for use in function evaluations -- application specific, see below.
84 Set routines for function and gradient evaluation
85 */
86 PetscCall(TaoSetObjective(tao, FormFunction, (void *)&user));
87 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
88
89 /*
90 Given the command line arguments, calculate the hessian with either the user-
91 provided function FormHessian, or the default finite-difference driven Hessian
92 functions
93 */
94 PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fddefault", &fddefault));
95 PetscCall(PetscOptionsHasName(NULL, NULL, "-tao_fdcoloring", &fdcoloring));
96
97 /*
98 Create a matrix data structure to store the Hessian and set
99 the Hessian evaluation routine.
100 Set the matrix nonzero structure to be used for Hessian evaluations
101 */
102 PetscCall(DMCreateMatrix(user.dm, &user.H));
103 PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
104
105 if (fdcoloring) {
106 PetscCall(DMCreateColoring(user.dm, IS_COLORING_GLOBAL, &iscoloring));
107 PetscCall(MatFDColoringCreate(user.H, iscoloring, &matfdcoloring));
108 PetscCall(ISColoringDestroy(&iscoloring));
109 PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)FormGradient, (void *)&user));
110 PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
111 PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessianColor, (void *)matfdcoloring));
112 } else if (fddefault) {
113 PetscCall(TaoSetHessian(tao, user.H, user.H, TaoDefaultComputeHessian, (void *)NULL));
114 } else {
115 PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
116 }
117
118 /*
119 If my_monitor option is in command line, then use the user-provided
120 monitoring function
121 */
122 PetscCall(PetscOptionsHasName(NULL, NULL, "-my_monitor", &viewmat));
123 if (viewmat) PetscCall(TaoMonitorSet(tao, My_Monitor, NULL, NULL));
124
125 /* Check for any tao command line options */
126 PetscCall(TaoSetFromOptions(tao));
127
128 /* SOLVE THE APPLICATION */
129 PetscCall(TaoSolve(tao));
130
131 PetscCall(TaoView(tao, PETSC_VIEWER_STDOUT_WORLD));
132
133 /* Free TAO data structures */
134 PetscCall(TaoDestroy(&tao));
135
136 /* Free PETSc data structures */
137 PetscCall(VecDestroy(&x));
138 PetscCall(MatDestroy(&user.H));
139 if (fdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
140 PetscCall(PetscFree(user.bottom));
141 PetscCall(PetscFree(user.top));
142 PetscCall(PetscFree(user.left));
143 PetscCall(PetscFree(user.right));
144 PetscCall(DMDestroy(&user.dm));
145 PetscCall(PetscFinalize());
146 return 0;
147 }
148
149 /* -------------------------------------------------------------------- */
150 /*
151 FormFunction - Evaluates the objective function.
152
153 Input Parameters:
154 . tao - the Tao context
155 . X - input vector
156 . userCtx - optional user-defined context, as set by TaoSetObjective()
157
158 Output Parameters:
159 . fcn - the newly evaluated function
160 */
FormFunction(Tao tao,Vec X,PetscReal * fcn,void * userCtx)161 PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *fcn, void *userCtx)
162 {
163 AppCtx *user = (AppCtx *)userCtx;
164 PetscInt i, j;
165 PetscInt mx = user->mx, my = user->my;
166 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
167 PetscReal ft = 0.0;
168 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), area = 0.5 * hx * hy;
169 PetscReal rhx = mx + 1, rhy = my + 1;
170 PetscReal f2, f4, d1, d2, d3, d4, xc, xl, xr, xt, xb;
171 PetscReal **x;
172 Vec localX;
173
174 PetscFunctionBegin;
175 /* Get local mesh boundaries */
176 PetscCall(DMGetLocalVector(user->dm, &localX));
177 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
178 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
179
180 /* Scatter ghost points to local vector */
181 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
182 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
183
184 /* Get pointers to vector data */
185 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
186
187 /* Compute function and gradient over the locally owned part of the mesh */
188 for (j = ys; j < ys + ym; j++) {
189 for (i = xs; i < xs + xm; i++) {
190 xc = x[j][i];
191
192 if (i == 0) { /* left side */
193 xl = user->left[j - ys + 1];
194 } else {
195 xl = x[j][i - 1];
196 }
197
198 if (j == 0) { /* bottom side */
199 xb = user->bottom[i - xs + 1];
200 } else {
201 xb = x[j - 1][i];
202 }
203
204 if (i + 1 == gxs + gxm) { /* right side */
205 xr = user->right[j - ys + 1];
206 } else {
207 xr = x[j][i + 1];
208 }
209
210 if (j + 1 == gys + gym) { /* top side */
211 xt = user->top[i - xs + 1];
212 } else {
213 xt = x[j + 1][i];
214 }
215
216 d1 = (xc - xl);
217 d2 = (xc - xr);
218 d3 = (xc - xt);
219 d4 = (xc - xb);
220
221 d1 *= rhx;
222 d2 *= rhx;
223 d3 *= rhy;
224 d4 *= rhy;
225
226 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
227 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
228
229 ft = ft + (f2 + f4);
230 }
231 }
232
233 /* Compute triangular areas along the border of the domain. */
234 if (xs == 0) { /* left side */
235 for (j = ys; j < ys + ym; j++) {
236 d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
237 d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
238 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
239 }
240 }
241 if (ys == 0) { /* bottom side */
242 for (i = xs; i < xs + xm; i++) {
243 d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
244 d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
245 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
246 }
247 }
248 if (xs + xm == mx) { /* right side */
249 for (j = ys; j < ys + ym; j++) {
250 d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
251 d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
252 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
253 }
254 }
255 if (ys + ym == my) { /* top side */
256 for (i = xs; i < xs + xm; i++) {
257 d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
258 d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
259 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
260 }
261 }
262 if (ys == 0 && xs == 0) {
263 d1 = (user->left[0] - user->left[1]) / hy;
264 d2 = (user->bottom[0] - user->bottom[1]) * rhx;
265 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
266 }
267 if (ys + ym == my && xs + xm == mx) {
268 d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
269 d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
270 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
271 }
272
273 ft = ft * area;
274 PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
275
276 /* Restore vectors */
277 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
278 PetscCall(DMRestoreLocalVector(user->dm, &localX));
279 PetscFunctionReturn(PETSC_SUCCESS);
280 }
281
282 /* -------------------------------------------------------------------- */
283 /*
284 FormFunctionGradient - Evaluates the function and corresponding gradient.
285
286 Input Parameters:
287 . tao - the Tao context
288 . X - input vector
289 . userCtx - optional user-defined context, as set by TaoSetObjectiveAndGradient()
290
291 Output Parameters:
292 . fcn - the newly evaluated function
293 . G - vector containing the newly evaluated gradient
294 */
FormFunctionGradient(Tao tao,Vec X,PetscReal * fcn,Vec G,void * userCtx)295 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx)
296 {
297 AppCtx *user = (AppCtx *)userCtx;
298 PetscInt i, j;
299 PetscInt mx = user->mx, my = user->my;
300 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
301 PetscReal ft = 0.0;
302 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy;
303 PetscReal rhx = mx + 1, rhy = my + 1;
304 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
305 PetscReal df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
306 PetscReal **g, **x;
307 Vec localX;
308
309 PetscFunctionBegin;
310 /* Get local mesh boundaries */
311 PetscCall(DMGetLocalVector(user->dm, &localX));
312 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
313 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
314
315 /* Scatter ghost points to local vector */
316 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
317 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
318
319 /* Get pointers to vector data */
320 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
321 PetscCall(DMDAVecGetArray(user->dm, G, (void **)&g));
322
323 /* Compute function and gradient over the locally owned part of the mesh */
324 for (j = ys; j < ys + ym; j++) {
325 for (i = xs; i < xs + xm; i++) {
326 xc = x[j][i];
327 xlt = xrb = xl = xr = xb = xt = xc;
328
329 if (i == 0) { /* left side */
330 xl = user->left[j - ys + 1];
331 xlt = user->left[j - ys + 2];
332 } else {
333 xl = x[j][i - 1];
334 }
335
336 if (j == 0) { /* bottom side */
337 xb = user->bottom[i - xs + 1];
338 xrb = user->bottom[i - xs + 2];
339 } else {
340 xb = x[j - 1][i];
341 }
342
343 if (i + 1 == gxs + gxm) { /* right side */
344 xr = user->right[j - ys + 1];
345 xrb = user->right[j - ys];
346 } else {
347 xr = x[j][i + 1];
348 }
349
350 if (j + 1 == gys + gym) { /* top side */
351 xt = user->top[i - xs + 1];
352 xlt = user->top[i - xs];
353 } else {
354 xt = x[j + 1][i];
355 }
356
357 if (i > gxs && j + 1 < gys + gym) xlt = x[j + 1][i - 1];
358 if (j > gys && i + 1 < gxs + gxm) xrb = x[j - 1][i + 1];
359
360 d1 = (xc - xl);
361 d2 = (xc - xr);
362 d3 = (xc - xt);
363 d4 = (xc - xb);
364 d5 = (xr - xrb);
365 d6 = (xrb - xb);
366 d7 = (xlt - xl);
367 d8 = (xt - xlt);
368
369 df1dxc = d1 * hydhx;
370 df2dxc = (d1 * hydhx + d4 * hxdhy);
371 df3dxc = d3 * hxdhy;
372 df4dxc = (d2 * hydhx + d3 * hxdhy);
373 df5dxc = d2 * hydhx;
374 df6dxc = d4 * hxdhy;
375
376 d1 *= rhx;
377 d2 *= rhx;
378 d3 *= rhy;
379 d4 *= rhy;
380 d5 *= rhy;
381 d6 *= rhx;
382 d7 *= rhy;
383 d8 *= rhx;
384
385 f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
386 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
387 f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
388 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
389 f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
390 f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
391
392 ft = ft + (f2 + f4);
393
394 df1dxc /= f1;
395 df2dxc /= f2;
396 df3dxc /= f3;
397 df4dxc /= f4;
398 df5dxc /= f5;
399 df6dxc /= f6;
400
401 g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) * 0.5;
402 }
403 }
404
405 /* Compute triangular areas along the border of the domain. */
406 if (xs == 0) { /* left side */
407 for (j = ys; j < ys + ym; j++) {
408 d3 = (user->left[j - ys + 1] - user->left[j - ys + 2]) * rhy;
409 d2 = (user->left[j - ys + 1] - x[j][0]) * rhx;
410 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
411 }
412 }
413 if (ys == 0) { /* bottom side */
414 for (i = xs; i < xs + xm; i++) {
415 d2 = (user->bottom[i + 1 - xs] - user->bottom[i - xs + 2]) * rhx;
416 d3 = (user->bottom[i - xs + 1] - x[0][i]) * rhy;
417 ft = ft + PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
418 }
419 }
420
421 if (xs + xm == mx) { /* right side */
422 for (j = ys; j < ys + ym; j++) {
423 d1 = (x[j][mx - 1] - user->right[j - ys + 1]) * rhx;
424 d4 = (user->right[j - ys] - user->right[j - ys + 1]) * rhy;
425 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
426 }
427 }
428 if (ys + ym == my) { /* top side */
429 for (i = xs; i < xs + xm; i++) {
430 d1 = (x[my - 1][i] - user->top[i - xs + 1]) * rhy;
431 d4 = (user->top[i - xs + 1] - user->top[i - xs]) * rhx;
432 ft = ft + PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
433 }
434 }
435
436 if (ys == 0 && xs == 0) {
437 d1 = (user->left[0] - user->left[1]) / hy;
438 d2 = (user->bottom[0] - user->bottom[1]) * rhx;
439 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
440 }
441 if (ys + ym == my && xs + xm == mx) {
442 d1 = (user->right[ym + 1] - user->right[ym]) * rhy;
443 d2 = (user->top[xm + 1] - user->top[xm]) * rhx;
444 ft += PetscSqrtReal(1.0 + d1 * d1 + d2 * d2);
445 }
446
447 ft = ft * area;
448 PetscCallMPI(MPIU_Allreduce(&ft, fcn, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD));
449
450 /* Restore vectors */
451 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
452 PetscCall(DMDAVecRestoreArray(user->dm, G, (void **)&g));
453 PetscCall(DMRestoreLocalVector(user->dm, &localX));
454 PetscCall(PetscLogFlops(67.0 * xm * ym));
455 PetscFunctionReturn(PETSC_SUCCESS);
456 }
457
FormGradient(Tao tao,Vec X,Vec G,void * userCtx)458 PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *userCtx)
459 {
460 PetscReal fcn;
461
462 PetscFunctionBegin;
463 PetscCall(FormFunctionGradient(tao, X, &fcn, G, userCtx));
464 PetscFunctionReturn(PETSC_SUCCESS);
465 }
466
467 /* ------------------------------------------------------------------- */
468 /*
469 FormHessian - Evaluates Hessian matrix.
470
471 Input Parameters:
472 . tao - the Tao context
473 . x - input vector
474 . ptr - optional user-defined context, as set by TaoSetHessian()
475
476 Output Parameters:
477 . H - Hessian matrix
478 . Hpre - optionally different matrix used to compute the preconditioner
479
480 */
FormHessian(Tao tao,Vec X,Mat H,Mat Hpre,void * ptr)481 PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
482 {
483 AppCtx *user = (AppCtx *)ptr;
484
485 PetscFunctionBegin;
486 /* Evaluate the Hessian entries*/
487 PetscCall(QuadraticH(user, X, H));
488 PetscFunctionReturn(PETSC_SUCCESS);
489 }
490
491 /* ------------------------------------------------------------------- */
492 /*
493 QuadraticH - Evaluates Hessian matrix.
494
495 Input Parameters:
496 . user - user-defined context, as set by TaoSetHessian()
497 . X - input vector
498
499 Output Parameter:
500 . H - Hessian matrix
501 */
QuadraticH(AppCtx * user,Vec X,Mat Hessian)502 PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian)
503 {
504 PetscInt i, j, k;
505 PetscInt mx = user->mx, my = user->my;
506 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym;
507 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
508 PetscReal f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
509 PetscReal hl, hr, ht, hb, hc, htl, hbr;
510 PetscReal **x, v[7];
511 MatStencil col[7], row;
512 Vec localX;
513 PetscBool assembled;
514
515 PetscFunctionBegin;
516 /* Get local mesh boundaries */
517 PetscCall(DMGetLocalVector(user->dm, &localX));
518
519 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
520 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
521
522 /* Scatter ghost points to local vector */
523 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX));
524 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX));
525
526 /* Get pointers to vector data */
527 PetscCall(DMDAVecGetArray(user->dm, localX, (void **)&x));
528
529 /* Initialize matrix entries to zero */
530 PetscCall(MatAssembled(Hessian, &assembled));
531 if (assembled) PetscCall(MatZeroEntries(Hessian));
532
533 /* Set various matrix options */
534 PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
535
536 /* Compute Hessian over the locally owned part of the mesh */
537 for (j = ys; j < ys + ym; j++) {
538 for (i = xs; i < xs + xm; i++) {
539 xc = x[j][i];
540 xlt = xrb = xl = xr = xb = xt = xc;
541
542 /* Left side */
543 if (i == 0) {
544 xl = user->left[j - ys + 1];
545 xlt = user->left[j - ys + 2];
546 } else {
547 xl = x[j][i - 1];
548 }
549
550 if (j == 0) {
551 xb = user->bottom[i - xs + 1];
552 xrb = user->bottom[i - xs + 2];
553 } else {
554 xb = x[j - 1][i];
555 }
556
557 if (i + 1 == mx) {
558 xr = user->right[j - ys + 1];
559 xrb = user->right[j - ys];
560 } else {
561 xr = x[j][i + 1];
562 }
563
564 if (j + 1 == my) {
565 xt = user->top[i - xs + 1];
566 xlt = user->top[i - xs];
567 } else {
568 xt = x[j + 1][i];
569 }
570
571 if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];
572 if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];
573
574 d1 = (xc - xl) / hx;
575 d2 = (xc - xr) / hx;
576 d3 = (xc - xt) / hy;
577 d4 = (xc - xb) / hy;
578 d5 = (xrb - xr) / hy;
579 d6 = (xrb - xb) / hx;
580 d7 = (xlt - xl) / hy;
581 d8 = (xlt - xt) / hx;
582
583 f1 = PetscSqrtReal(1.0 + d1 * d1 + d7 * d7);
584 f2 = PetscSqrtReal(1.0 + d1 * d1 + d4 * d4);
585 f3 = PetscSqrtReal(1.0 + d3 * d3 + d8 * d8);
586 f4 = PetscSqrtReal(1.0 + d3 * d3 + d2 * d2);
587 f5 = PetscSqrtReal(1.0 + d2 * d2 + d5 * d5);
588 f6 = PetscSqrtReal(1.0 + d4 * d4 + d6 * d6);
589
590 hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
591 hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
592 ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
593 hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
594
595 hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
596 htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
597
598 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);
599
600 hl /= 2.0;
601 hr /= 2.0;
602 ht /= 2.0;
603 hb /= 2.0;
604 hbr /= 2.0;
605 htl /= 2.0;
606 hc /= 2.0;
607
608 row.j = j;
609 row.i = i;
610 k = 0;
611 if (j > 0) {
612 v[k] = hb;
613 col[k].j = j - 1;
614 col[k].i = i;
615 k++;
616 }
617
618 if (j > 0 && i < mx - 1) {
619 v[k] = hbr;
620 col[k].j = j - 1;
621 col[k].i = i + 1;
622 k++;
623 }
624
625 if (i > 0) {
626 v[k] = hl;
627 col[k].j = j;
628 col[k].i = i - 1;
629 k++;
630 }
631
632 v[k] = hc;
633 col[k].j = j;
634 col[k].i = i;
635 k++;
636
637 if (i < mx - 1) {
638 v[k] = hr;
639 col[k].j = j;
640 col[k].i = i + 1;
641 k++;
642 }
643
644 if (i > 0 && j < my - 1) {
645 v[k] = htl;
646 col[k].j = j + 1;
647 col[k].i = i - 1;
648 k++;
649 }
650
651 if (j < my - 1) {
652 v[k] = ht;
653 col[k].j = j + 1;
654 col[k].i = i;
655 k++;
656 }
657
658 /*
659 Set matrix values using local numbering, which was defined
660 earlier, in the main routine.
661 */
662 PetscCall(MatSetValuesStencil(Hessian, 1, &row, k, col, v, INSERT_VALUES));
663 }
664 }
665
666 PetscCall(DMDAVecRestoreArray(user->dm, localX, (void **)&x));
667 PetscCall(DMRestoreLocalVector(user->dm, &localX));
668
669 PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
670 PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
671
672 PetscCall(PetscLogFlops(199.0 * xm * ym));
673 PetscFunctionReturn(PETSC_SUCCESS);
674 }
675
676 /* ------------------------------------------------------------------- */
677 /*
678 MSA_BoundaryConditions - Calculates the boundary conditions for
679 the region.
680
681 Input Parameter:
682 . user - user-defined application context
683
684 Output Parameter:
685 . user - user-defined application context
686 */
MSA_BoundaryConditions(AppCtx * user)687 static PetscErrorCode MSA_BoundaryConditions(AppCtx *user)
688 {
689 PetscInt i, j, k, limit = 0, maxits = 5;
690 PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym;
691 PetscInt mx = user->mx, my = user->my;
692 PetscInt bsize = 0, lsize = 0, tsize = 0, rsize = 0;
693 PetscReal one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
694 PetscReal fnorm, det, hx, hy, xt = 0, yt = 0;
695 PetscReal u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
696 PetscReal b = -0.5, t = 0.5, l = -0.5, r = 0.5;
697 PetscReal *boundary;
698 PetscBool flg;
699
700 PetscFunctionBegin;
701 /* Get local mesh boundaries */
702 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
703 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));
704
705 bsize = xm + 2;
706 lsize = ym + 2;
707 rsize = ym + 2;
708 tsize = xm + 2;
709
710 PetscCall(PetscMalloc1(bsize, &user->bottom));
711 PetscCall(PetscMalloc1(tsize, &user->top));
712 PetscCall(PetscMalloc1(lsize, &user->left));
713 PetscCall(PetscMalloc1(rsize, &user->right));
714
715 hx = (r - l) / (mx + 1);
716 hy = (t - b) / (my + 1);
717
718 for (j = 0; j < 4; j++) {
719 if (j == 0) {
720 yt = b;
721 xt = l + hx * xs;
722 limit = bsize;
723 boundary = user->bottom;
724 } else if (j == 1) {
725 yt = t;
726 xt = l + hx * xs;
727 limit = tsize;
728 boundary = user->top;
729 } else if (j == 2) {
730 yt = b + hy * ys;
731 xt = l;
732 limit = lsize;
733 boundary = user->left;
734 } else { /* if (j==3) */
735 yt = b + hy * ys;
736 xt = r;
737 limit = rsize;
738 boundary = user->right;
739 }
740
741 for (i = 0; i < limit; i++) {
742 u1 = xt;
743 u2 = -yt;
744 for (k = 0; k < maxits; k++) {
745 nf1 = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
746 nf2 = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
747 fnorm = PetscSqrtReal(nf1 * nf1 + nf2 * nf2);
748 if (fnorm <= tol) break;
749 njac11 = one + u2 * u2 - u1 * u1;
750 njac12 = two * u1 * u2;
751 njac21 = -two * u1 * u2;
752 njac22 = -one - u1 * u1 + u2 * u2;
753 det = njac11 * njac22 - njac21 * njac12;
754 u1 = u1 - (njac22 * nf1 - njac12 * nf2) / det;
755 u2 = u2 - (njac11 * nf2 - njac21 * nf1) / det;
756 }
757
758 boundary[i] = u1 * u1 - u2 * u2;
759 if (j == 0 || j == 1) {
760 xt = xt + hx;
761 } else { /* if (j==2 || j==3) */
762 yt = yt + hy;
763 }
764 }
765 }
766
767 /* Scale the boundary if desired */
768 if (1 == 1) {
769 PetscReal scl = 1.0;
770
771 PetscCall(PetscOptionsGetReal(NULL, NULL, "-bottom", &scl, &flg));
772 if (flg) {
773 for (i = 0; i < bsize; i++) user->bottom[i] *= scl;
774 }
775
776 PetscCall(PetscOptionsGetReal(NULL, NULL, "-top", &scl, &flg));
777 if (flg) {
778 for (i = 0; i < tsize; i++) user->top[i] *= scl;
779 }
780
781 PetscCall(PetscOptionsGetReal(NULL, NULL, "-right", &scl, &flg));
782 if (flg) {
783 for (i = 0; i < rsize; i++) user->right[i] *= scl;
784 }
785
786 PetscCall(PetscOptionsGetReal(NULL, NULL, "-left", &scl, &flg));
787 if (flg) {
788 for (i = 0; i < lsize; i++) user->left[i] *= scl;
789 }
790 }
791 PetscFunctionReturn(PETSC_SUCCESS);
792 }
793
794 /* ------------------------------------------------------------------- */
795 /*
796 MSA_InitialPoint - Calculates the initial guess in one of three ways.
797
798 Input Parameters:
799 . user - user-defined application context
800 . X - vector for initial guess
801
802 Output Parameters:
803 . X - newly computed initial guess
804 */
MSA_InitialPoint(AppCtx * user,Vec X)805 static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X)
806 {
807 PetscInt start2 = -1, i, j;
808 PetscReal start1 = 0;
809 PetscBool flg1, flg2;
810
811 PetscFunctionBegin;
812 PetscCall(PetscOptionsGetReal(NULL, NULL, "-start", &start1, &flg1));
813 PetscCall(PetscOptionsGetInt(NULL, NULL, "-random", &start2, &flg2));
814
815 if (flg1) { /* The zero vector is reasonable */
816 PetscCall(VecSet(X, start1));
817 } else if (flg2 && start2 > 0) { /* Try a random start between -0.5 and 0.5 */
818 PetscRandom rctx;
819 PetscReal np5 = -0.5;
820
821 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
822 PetscCall(VecSetRandom(X, rctx));
823 PetscCall(PetscRandomDestroy(&rctx));
824 PetscCall(VecShift(X, np5));
825 } else { /* Take an average of the boundary conditions */
826 PetscInt xs, xm, ys, ym;
827 PetscInt mx = user->mx, my = user->my;
828 PetscReal **x;
829
830 /* Get local mesh boundaries */
831 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
832
833 /* Get pointers to vector data */
834 PetscCall(DMDAVecGetArray(user->dm, X, (void **)&x));
835
836 /* Perform local computations */
837 for (j = ys; j < ys + ym; j++) {
838 for (i = xs; i < xs + xm; i++) x[j][i] = (((j + 1) * user->bottom[i - xs + 1] + (my - j + 1) * user->top[i - xs + 1]) / (my + 2) + ((i + 1) * user->left[j - ys + 1] + (mx - i + 1) * user->right[j - ys + 1]) / (mx + 2)) / 2.0;
839 }
840 PetscCall(DMDAVecRestoreArray(user->dm, X, (void **)&x));
841 PetscCall(PetscLogFlops(9.0 * xm * ym));
842 }
843 PetscFunctionReturn(PETSC_SUCCESS);
844 }
845
846 /*-----------------------------------------------------------------------*/
My_Monitor(Tao tao,PetscCtx ctx)847 PetscErrorCode My_Monitor(Tao tao, PetscCtx ctx)
848 {
849 Vec X;
850
851 PetscFunctionBegin;
852 PetscCall(TaoGetSolution(tao, &X));
853 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
854 PetscFunctionReturn(PETSC_SUCCESS);
855 }
856
857 /*TEST
858
859 build:
860 requires: !complex
861
862 test:
863 args: -tao_monitor_short -tao_type lmvm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
864 requires: !single
865
866 test:
867 suffix: 2
868 nsize: 2
869 args: -tao_monitor_short -tao_type nls -tao_nls_ksp_max_it 15 -tao_gatol 1.e-4
870 filter: grep -v "nls ksp"
871 requires: !single
872
873 test:
874 suffix: 2_snes
875 nsize: 2
876 args: -tao_monitor_short -tao_type snes -ksp_converged_maxits -ksp_max_it 15 -snes_atol 1.e-4
877 filter: grep -v "nls ksp"
878 requires: !single
879
880 test:
881 suffix: 3
882 nsize: 3
883 args: -tao_monitor_short -tao_type cg -tao_cg_type fr -da_grid_x 10 -da_grid_y 10 -tao_gatol 1.e-3
884 requires: !single
885
886 test:
887 suffix: 3_snes
888 nsize: 3
889 args: -tao_monitor_short -tao_type snes -snes_type ncg -snes_ncg_type fr -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4
890 requires: !single
891
892 test:
893 suffix: 4_snes_ngmres
894 args: -tao_type snes -snes_type ngmres -npc_snes_type ncg -da_grid_x 10 -da_grid_y 10 -snes_atol 1.e-4 -npc_snes_ncg_type fr -snes_converged_reason -snes_monitor ::ascii_info_detail -snes_ngmres_monitor -snes_ngmres_select_type {{linesearch difference}separate output}
895 requires: !single
896
897 test:
898 suffix: 5
899 nsize: 2
900 args: -tao_monitor_short -tao_type bmrm -da_grid_x 10 -da_grid_y 8 -tao_gatol 1.e-3
901 requires: !single
902
903 TEST*/
904