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