xref: /petsc/src/tao/unconstrained/tutorials/minsurf1.c (revision 9371c9d470a9602b6d10a8bf50c9b2280a79e45a)
1c4762a1bSJed Brown /* Program usage: mpiexec -n 1 minsurf1 [-help] [all TAO options] */
2c4762a1bSJed Brown 
3c4762a1bSJed Brown /*  Include "petsctao.h" so we can use TAO solvers.  */
4c4762a1bSJed Brown #include <petsctao.h>
5c4762a1bSJed Brown 
6*9371c9d4SSatish Balay static char help[] = "This example demonstrates use of the TAO package to \n\
7c4762a1bSJed Brown solve an unconstrained minimization problem.  This example is based on a \n\
8c4762a1bSJed Brown problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and \n\
9c4762a1bSJed Brown boundary values along the edges of the domain, the objective is to find the\n\
10c4762a1bSJed Brown surface with the minimal area that satisfies the boundary conditions.\n\
11c4762a1bSJed Brown The command line options are:\n\
12c4762a1bSJed Brown   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
13c4762a1bSJed Brown   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
14c4762a1bSJed Brown   -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise \n\n";
15c4762a1bSJed Brown 
16c4762a1bSJed Brown /*
17c4762a1bSJed Brown    User-defined application context - contains data needed by the
18c4762a1bSJed Brown    application-provided call-back routines, FormFunctionGradient()
19c4762a1bSJed Brown    and FormHessian().
20c4762a1bSJed Brown */
21c4762a1bSJed Brown typedef struct {
22c4762a1bSJed Brown   PetscInt   mx, my;                      /* discretization in x, y directions */
23c4762a1bSJed Brown   PetscReal *bottom, *top, *left, *right; /* boundary values */
24c4762a1bSJed Brown   Mat        H;
25c4762a1bSJed Brown } AppCtx;
26c4762a1bSJed Brown 
27c4762a1bSJed Brown /* -------- User-defined Routines --------- */
28c4762a1bSJed Brown 
29c4762a1bSJed Brown static PetscErrorCode MSA_BoundaryConditions(AppCtx *);
30c4762a1bSJed Brown static PetscErrorCode MSA_InitialPoint(AppCtx *, Vec);
31c4762a1bSJed Brown static PetscErrorCode QuadraticH(AppCtx *, Vec, Mat);
32c4762a1bSJed Brown PetscErrorCode        FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
33c4762a1bSJed Brown PetscErrorCode        FormHessian(Tao, Vec, Mat, Mat, void *);
34c4762a1bSJed Brown 
35*9371c9d4SSatish Balay int main(int argc, char **argv) {
36c4762a1bSJed Brown   PetscInt    N;    /* Size of vector */
37c4762a1bSJed Brown   PetscMPIInt size; /* Number of processors */
38c4762a1bSJed Brown   Vec         x;    /* solution, gradient vectors */
39c4762a1bSJed Brown   KSP         ksp;  /*  PETSc Krylov subspace method */
40c4762a1bSJed Brown   PetscBool   flg;  /* A return value when checking for user options */
41c4762a1bSJed Brown   Tao         tao;  /* Tao solver context */
42c4762a1bSJed Brown   AppCtx      user; /* user-defined work context */
43c4762a1bSJed Brown 
44c4762a1bSJed Brown   /* Initialize TAO,PETSc */
45327415f7SBarry Smith   PetscFunctionBeginUser;
469566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
47c4762a1bSJed Brown 
489566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
493c859ba3SBarry Smith   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Incorrect number of processors");
50c4762a1bSJed Brown 
51c4762a1bSJed Brown   /* Specify default dimension of the problem */
52*9371c9d4SSatish Balay   user.mx = 4;
53*9371c9d4SSatish Balay   user.my = 4;
54c4762a1bSJed Brown 
55c4762a1bSJed Brown   /* Check for any command line arguments that override defaults */
569566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));
579566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
58c4762a1bSJed Brown 
599566063dSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n---- Minimum Surface Area Problem -----\n"));
6063a3b9bcSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", user.mx, user.my));
61c4762a1bSJed Brown 
62c4762a1bSJed Brown   /* Calculate any derived values from parameters */
63c4762a1bSJed Brown   N = user.mx * user.my;
64c4762a1bSJed Brown 
65c4762a1bSJed Brown   /* Create TAO solver and set desired solution method  */
669566063dSJacob Faibussowitsch   PetscCall(TaoCreate(PETSC_COMM_SELF, &tao));
679566063dSJacob Faibussowitsch   PetscCall(TaoSetType(tao, TAOLMVM));
68c4762a1bSJed Brown 
69c4762a1bSJed Brown   /* Initialize minsurf application data structure for use in the function evaluations  */
709566063dSJacob Faibussowitsch   PetscCall(MSA_BoundaryConditions(&user)); /* Application specific routine */
71c4762a1bSJed Brown 
72c4762a1bSJed Brown   /*
73c4762a1bSJed Brown      Create a vector to hold the variables.  Compute an initial solution.
74c4762a1bSJed Brown      Set this vector, which will be used by TAO.
75c4762a1bSJed Brown   */
769566063dSJacob Faibussowitsch   PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &x));
779566063dSJacob Faibussowitsch   PetscCall(MSA_InitialPoint(&user, x)); /* Application specific routine */
789566063dSJacob Faibussowitsch   PetscCall(TaoSetSolution(tao, x));     /* A TAO routine                */
79c4762a1bSJed Brown 
80c4762a1bSJed Brown   /* Provide TAO routines for function, gradient, and Hessian evaluation */
819566063dSJacob Faibussowitsch   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
82c4762a1bSJed Brown 
83c4762a1bSJed Brown   /* Create a matrix data structure to store the Hessian.  This structure will be used by TAO */
849566063dSJacob Faibussowitsch   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, N, N, 7, NULL, &(user.H)));
859566063dSJacob Faibussowitsch   PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
869566063dSJacob Faibussowitsch   PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
87c4762a1bSJed Brown 
88c4762a1bSJed Brown   /* Check for any TAO command line options */
899566063dSJacob Faibussowitsch   PetscCall(TaoSetFromOptions(tao));
90c4762a1bSJed Brown 
91c4762a1bSJed Brown   /* Limit the number of iterations in the KSP linear solver */
929566063dSJacob Faibussowitsch   PetscCall(TaoGetKSP(tao, &ksp));
931baa6e33SBarry Smith   if (ksp) PetscCall(KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, user.mx * user.my));
94c4762a1bSJed Brown 
95c4762a1bSJed Brown   /* SOLVE THE APPLICATION */
969566063dSJacob Faibussowitsch   PetscCall(TaoSolve(tao));
97c4762a1bSJed Brown 
989566063dSJacob Faibussowitsch   PetscCall(TaoDestroy(&tao));
999566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&x));
1009566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&user.H));
1019566063dSJacob Faibussowitsch   PetscCall(PetscFree(user.bottom));
1029566063dSJacob Faibussowitsch   PetscCall(PetscFree(user.top));
1039566063dSJacob Faibussowitsch   PetscCall(PetscFree(user.left));
1049566063dSJacob Faibussowitsch   PetscCall(PetscFree(user.right));
105c4762a1bSJed Brown 
1069566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
107b122ec5aSJacob Faibussowitsch   return 0;
108c4762a1bSJed Brown }
109c4762a1bSJed Brown 
110c4762a1bSJed Brown /* -------------------------------------------------------------------- */
111c4762a1bSJed Brown 
112c4762a1bSJed Brown /*  FormFunctionGradient - Evaluates function and corresponding gradient.
113c4762a1bSJed Brown 
114c4762a1bSJed Brown     Input Parameters:
115c4762a1bSJed Brown .   tao     - the Tao context
116c4762a1bSJed Brown .   X       - input vector
117c4762a1bSJed Brown .   userCtx - optional user-defined context, as set by TaoSetFunctionGradient()
118c4762a1bSJed Brown 
119c4762a1bSJed Brown     Output Parameters:
120c4762a1bSJed Brown .   fcn     - the newly evaluated function
121c4762a1bSJed Brown .   G       - vector containing the newly evaluated gradient
122c4762a1bSJed Brown */
123*9371c9d4SSatish Balay PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn, Vec G, void *userCtx) {
124c4762a1bSJed Brown   AppCtx            *user = (AppCtx *)userCtx;
125c4762a1bSJed Brown   PetscInt           i, j, row;
126c4762a1bSJed Brown   PetscInt           mx = user->mx, my = user->my;
127c4762a1bSJed Brown   PetscReal          rhx = mx + 1, rhy = my + 1;
128c4762a1bSJed Brown   PetscReal          hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy, area = 0.5 * hx * hy, ft = 0;
129c4762a1bSJed Brown   PetscReal          f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
130c4762a1bSJed Brown   PetscReal          df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
131c4762a1bSJed Brown   PetscReal          zero = 0.0;
132c4762a1bSJed Brown   PetscScalar       *g;
133c4762a1bSJed Brown   const PetscScalar *x;
134c4762a1bSJed Brown 
135c4762a1bSJed Brown   PetscFunctionBeginUser;
1369566063dSJacob Faibussowitsch   PetscCall(VecSet(G, zero));
137c4762a1bSJed Brown 
1389566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(X, &x));
1399566063dSJacob Faibussowitsch   PetscCall(VecGetArray(G, &g));
140c4762a1bSJed Brown 
141c4762a1bSJed Brown   /* Compute function over the locally owned part of the mesh */
142c4762a1bSJed Brown   for (j = 0; j < my; j++) {
143c4762a1bSJed Brown     for (i = 0; i < mx; i++) {
144c4762a1bSJed Brown       row = (j)*mx + (i);
145c4762a1bSJed Brown       xc  = x[row];
146c4762a1bSJed Brown       xlt = xrb = xl = xr = xb = xt = xc;
147c4762a1bSJed Brown       if (i == 0) { /* left side */
148c4762a1bSJed Brown         xl  = user->left[j + 1];
149c4762a1bSJed Brown         xlt = user->left[j + 2];
150c4762a1bSJed Brown       } else {
151c4762a1bSJed Brown         xl = x[row - 1];
152c4762a1bSJed Brown       }
153c4762a1bSJed Brown 
154c4762a1bSJed Brown       if (j == 0) { /* bottom side */
155c4762a1bSJed Brown         xb  = user->bottom[i + 1];
156c4762a1bSJed Brown         xrb = user->bottom[i + 2];
157c4762a1bSJed Brown       } else {
158c4762a1bSJed Brown         xb = x[row - mx];
159c4762a1bSJed Brown       }
160c4762a1bSJed Brown 
161c4762a1bSJed Brown       if (i + 1 == mx) { /* right side */
162c4762a1bSJed Brown         xr  = user->right[j + 1];
163c4762a1bSJed Brown         xrb = user->right[j];
164c4762a1bSJed Brown       } else {
165c4762a1bSJed Brown         xr = x[row + 1];
166c4762a1bSJed Brown       }
167c4762a1bSJed Brown 
168c4762a1bSJed Brown       if (j + 1 == 0 + my) { /* top side */
169c4762a1bSJed Brown         xt  = user->top[i + 1];
170c4762a1bSJed Brown         xlt = user->top[i];
171c4762a1bSJed Brown       } else {
172c4762a1bSJed Brown         xt = x[row + mx];
173c4762a1bSJed Brown       }
174c4762a1bSJed Brown 
175*9371c9d4SSatish Balay       if (i > 0 && j + 1 < my) { xlt = x[row - 1 + mx]; }
176*9371c9d4SSatish Balay       if (j > 0 && i + 1 < mx) { xrb = x[row + 1 - mx]; }
177c4762a1bSJed Brown 
178c4762a1bSJed Brown       d1 = (xc - xl);
179c4762a1bSJed Brown       d2 = (xc - xr);
180c4762a1bSJed Brown       d3 = (xc - xt);
181c4762a1bSJed Brown       d4 = (xc - xb);
182c4762a1bSJed Brown       d5 = (xr - xrb);
183c4762a1bSJed Brown       d6 = (xrb - xb);
184c4762a1bSJed Brown       d7 = (xlt - xl);
185c4762a1bSJed Brown       d8 = (xt - xlt);
186c4762a1bSJed Brown 
187c4762a1bSJed Brown       df1dxc = d1 * hydhx;
188c4762a1bSJed Brown       df2dxc = (d1 * hydhx + d4 * hxdhy);
189c4762a1bSJed Brown       df3dxc = d3 * hxdhy;
190c4762a1bSJed Brown       df4dxc = (d2 * hydhx + d3 * hxdhy);
191c4762a1bSJed Brown       df5dxc = d2 * hydhx;
192c4762a1bSJed Brown       df6dxc = d4 * hxdhy;
193c4762a1bSJed Brown 
194c4762a1bSJed Brown       d1 *= rhx;
195c4762a1bSJed Brown       d2 *= rhx;
196c4762a1bSJed Brown       d3 *= rhy;
197c4762a1bSJed Brown       d4 *= rhy;
198c4762a1bSJed Brown       d5 *= rhy;
199c4762a1bSJed Brown       d6 *= rhx;
200c4762a1bSJed Brown       d7 *= rhy;
201c4762a1bSJed Brown       d8 *= rhx;
202c4762a1bSJed Brown 
203c4762a1bSJed Brown       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
204c4762a1bSJed Brown       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
205c4762a1bSJed Brown       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
206c4762a1bSJed Brown       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
207c4762a1bSJed Brown       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
208c4762a1bSJed Brown       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
209c4762a1bSJed Brown 
210c4762a1bSJed Brown       ft = ft + (f2 + f4);
211c4762a1bSJed Brown 
212c4762a1bSJed Brown       df1dxc /= f1;
213c4762a1bSJed Brown       df2dxc /= f2;
214c4762a1bSJed Brown       df3dxc /= f3;
215c4762a1bSJed Brown       df4dxc /= f4;
216c4762a1bSJed Brown       df5dxc /= f5;
217c4762a1bSJed Brown       df6dxc /= f6;
218c4762a1bSJed Brown 
219c4762a1bSJed Brown       g[row] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
220c4762a1bSJed Brown     }
221c4762a1bSJed Brown   }
222c4762a1bSJed Brown 
223c4762a1bSJed Brown   for (j = 0; j < my; j++) { /* left side */
224c4762a1bSJed Brown     d3 = (user->left[j + 1] - user->left[j + 2]) * rhy;
225c4762a1bSJed Brown     d2 = (user->left[j + 1] - x[j * mx]) * rhx;
226c4762a1bSJed Brown     ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
227c4762a1bSJed Brown   }
228c4762a1bSJed Brown 
229c4762a1bSJed Brown   for (i = 0; i < mx; i++) { /* bottom */
230c4762a1bSJed Brown     d2 = (user->bottom[i + 1] - user->bottom[i + 2]) * rhx;
231c4762a1bSJed Brown     d3 = (user->bottom[i + 1] - x[i]) * rhy;
232c4762a1bSJed Brown     ft = ft + PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
233c4762a1bSJed Brown   }
234c4762a1bSJed Brown 
235c4762a1bSJed Brown   for (j = 0; j < my; j++) { /* right side */
236c4762a1bSJed Brown     d1 = (x[(j + 1) * mx - 1] - user->right[j + 1]) * rhx;
237c4762a1bSJed Brown     d4 = (user->right[j] - user->right[j + 1]) * rhy;
238c4762a1bSJed Brown     ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
239c4762a1bSJed Brown   }
240c4762a1bSJed Brown 
241c4762a1bSJed Brown   for (i = 0; i < mx; i++) { /* top side */
242c4762a1bSJed Brown     d1 = (x[(my - 1) * mx + i] - user->top[i + 1]) * rhy;
243c4762a1bSJed Brown     d4 = (user->top[i + 1] - user->top[i]) * rhx;
244c4762a1bSJed Brown     ft = ft + PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
245c4762a1bSJed Brown   }
246c4762a1bSJed Brown 
247c4762a1bSJed Brown   /* Bottom left corner */
248c4762a1bSJed Brown   d1 = (user->left[0] - user->left[1]) * rhy;
249c4762a1bSJed Brown   d2 = (user->bottom[0] - user->bottom[1]) * rhx;
250c4762a1bSJed Brown   ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
251c4762a1bSJed Brown 
252c4762a1bSJed Brown   /* Top right corner */
253c4762a1bSJed Brown   d1 = (user->right[my + 1] - user->right[my]) * rhy;
254c4762a1bSJed Brown   d2 = (user->top[mx + 1] - user->top[mx]) * rhx;
255c4762a1bSJed Brown   ft += PetscSqrtScalar(1.0 + d1 * d1 + d2 * d2);
256c4762a1bSJed Brown 
257c4762a1bSJed Brown   (*fcn) = ft * area;
258c4762a1bSJed Brown 
259c4762a1bSJed Brown   /* Restore vectors */
2609566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(X, &x));
2619566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(G, &g));
2629566063dSJacob Faibussowitsch   PetscCall(PetscLogFlops(67.0 * mx * my));
263c4762a1bSJed Brown   PetscFunctionReturn(0);
264c4762a1bSJed Brown }
265c4762a1bSJed Brown 
266c4762a1bSJed Brown /* ------------------------------------------------------------------- */
267c4762a1bSJed Brown /*
268c4762a1bSJed Brown    FormHessian - Evaluates the Hessian matrix.
269c4762a1bSJed Brown 
270c4762a1bSJed Brown    Input Parameters:
271c4762a1bSJed Brown .  tao  - the Tao context
272c4762a1bSJed Brown .  x    - input vector
273c4762a1bSJed Brown .  ptr  - optional user-defined context, as set by TaoSetHessian()
274c4762a1bSJed Brown 
275c4762a1bSJed Brown    Output Parameters:
276c4762a1bSJed Brown .  H    - Hessian matrix
277c4762a1bSJed Brown .  Hpre - optionally different preconditioning matrix
278c4762a1bSJed Brown .  flg  - flag indicating matrix structure
279c4762a1bSJed Brown 
280c4762a1bSJed Brown */
281*9371c9d4SSatish Balay PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr) {
282c4762a1bSJed Brown   AppCtx *user = (AppCtx *)ptr;
283c4762a1bSJed Brown 
284c4762a1bSJed Brown   PetscFunctionBeginUser;
285c4762a1bSJed Brown   /* Evaluate the Hessian entries*/
2869566063dSJacob Faibussowitsch   PetscCall(QuadraticH(user, X, H));
287c4762a1bSJed Brown   PetscFunctionReturn(0);
288c4762a1bSJed Brown }
289c4762a1bSJed Brown 
290c4762a1bSJed Brown /* ------------------------------------------------------------------- */
291c4762a1bSJed Brown /*
292c4762a1bSJed Brown    QuadraticH - Evaluates the Hessian matrix.
293c4762a1bSJed Brown 
294c4762a1bSJed Brown    Input Parameters:
295c4762a1bSJed Brown .  user - user-defined context, as set by TaoSetHessian()
296c4762a1bSJed Brown .  X    - input vector
297c4762a1bSJed Brown 
298c4762a1bSJed Brown    Output Parameter:
299c4762a1bSJed Brown .  H    - Hessian matrix
300c4762a1bSJed Brown */
301*9371c9d4SSatish Balay PetscErrorCode QuadraticH(AppCtx *user, Vec X, Mat Hessian) {
302c4762a1bSJed Brown   PetscInt           i, j, k, row;
303c4762a1bSJed Brown   PetscInt           mx = user->mx, my = user->my;
304c4762a1bSJed Brown   PetscInt           col[7];
305c4762a1bSJed Brown   PetscReal          hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), hydhx = hy / hx, hxdhy = hx / hy;
306c4762a1bSJed Brown   PetscReal          rhx = mx + 1, rhy = my + 1;
307c4762a1bSJed Brown   PetscReal          f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
308c4762a1bSJed Brown   PetscReal          hl, hr, ht, hb, hc, htl, hbr;
309c4762a1bSJed Brown   const PetscScalar *x;
310c4762a1bSJed Brown   PetscReal          v[7];
311c4762a1bSJed Brown 
312362febeeSStefano Zampini   PetscFunctionBeginUser;
313c4762a1bSJed Brown   /* Get pointers to vector data */
3149566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(X, &x));
315c4762a1bSJed Brown 
316c4762a1bSJed Brown   /* Initialize matrix entries to zero */
3179566063dSJacob Faibussowitsch   PetscCall(MatZeroEntries(Hessian));
318c4762a1bSJed Brown 
319c4762a1bSJed Brown   /* Set various matrix options */
3209566063dSJacob Faibussowitsch   PetscCall(MatSetOption(Hessian, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
321c4762a1bSJed Brown 
322c4762a1bSJed Brown   /* Compute Hessian over the locally owned part of the mesh */
323c4762a1bSJed Brown   for (i = 0; i < mx; i++) {
324c4762a1bSJed Brown     for (j = 0; j < my; j++) {
325c4762a1bSJed Brown       row = (j)*mx + (i);
326c4762a1bSJed Brown 
327c4762a1bSJed Brown       xc  = x[row];
328c4762a1bSJed Brown       xlt = xrb = xl = xr = xb = xt = xc;
329c4762a1bSJed Brown 
330c4762a1bSJed Brown       /* Left side */
331c4762a1bSJed Brown       if (i == 0) {
332c4762a1bSJed Brown         xl  = user->left[j + 1];
333c4762a1bSJed Brown         xlt = user->left[j + 2];
334c4762a1bSJed Brown       } else {
335c4762a1bSJed Brown         xl = x[row - 1];
336c4762a1bSJed Brown       }
337c4762a1bSJed Brown 
338c4762a1bSJed Brown       if (j == 0) {
339c4762a1bSJed Brown         xb  = user->bottom[i + 1];
340c4762a1bSJed Brown         xrb = user->bottom[i + 2];
341c4762a1bSJed Brown       } else {
342c4762a1bSJed Brown         xb = x[row - mx];
343c4762a1bSJed Brown       }
344c4762a1bSJed Brown 
345c4762a1bSJed Brown       if (i + 1 == mx) {
346c4762a1bSJed Brown         xr  = user->right[j + 1];
347c4762a1bSJed Brown         xrb = user->right[j];
348c4762a1bSJed Brown       } else {
349c4762a1bSJed Brown         xr = x[row + 1];
350c4762a1bSJed Brown       }
351c4762a1bSJed Brown 
352c4762a1bSJed Brown       if (j + 1 == my) {
353c4762a1bSJed Brown         xt  = user->top[i + 1];
354c4762a1bSJed Brown         xlt = user->top[i];
355c4762a1bSJed Brown       } else {
356c4762a1bSJed Brown         xt = x[row + mx];
357c4762a1bSJed Brown       }
358c4762a1bSJed Brown 
359*9371c9d4SSatish Balay       if (i > 0 && j + 1 < my) { xlt = x[row - 1 + mx]; }
360*9371c9d4SSatish Balay       if (j > 0 && i + 1 < mx) { xrb = x[row + 1 - mx]; }
361c4762a1bSJed Brown 
362c4762a1bSJed Brown       d1 = (xc - xl) * rhx;
363c4762a1bSJed Brown       d2 = (xc - xr) * rhx;
364c4762a1bSJed Brown       d3 = (xc - xt) * rhy;
365c4762a1bSJed Brown       d4 = (xc - xb) * rhy;
366c4762a1bSJed Brown       d5 = (xrb - xr) * rhy;
367c4762a1bSJed Brown       d6 = (xrb - xb) * rhx;
368c4762a1bSJed Brown       d7 = (xlt - xl) * rhy;
369c4762a1bSJed Brown       d8 = (xlt - xt) * rhx;
370c4762a1bSJed Brown 
371c4762a1bSJed Brown       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
372c4762a1bSJed Brown       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
373c4762a1bSJed Brown       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
374c4762a1bSJed Brown       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
375c4762a1bSJed Brown       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
376c4762a1bSJed Brown       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);
377c4762a1bSJed Brown 
378c4762a1bSJed Brown       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
379c4762a1bSJed Brown       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
380c4762a1bSJed Brown       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
381c4762a1bSJed Brown       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);
382c4762a1bSJed Brown 
383c4762a1bSJed Brown       hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
384c4762a1bSJed Brown       htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);
385c4762a1bSJed Brown 
386*9371c9d4SSatish Balay       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);
387c4762a1bSJed Brown 
388*9371c9d4SSatish Balay       hl *= 0.5;
389*9371c9d4SSatish Balay       hr *= 0.5;
390*9371c9d4SSatish Balay       ht *= 0.5;
391*9371c9d4SSatish Balay       hb *= 0.5;
392*9371c9d4SSatish Balay       hbr *= 0.5;
393*9371c9d4SSatish Balay       htl *= 0.5;
394*9371c9d4SSatish Balay       hc *= 0.5;
395c4762a1bSJed Brown 
396c4762a1bSJed Brown       k = 0;
397c4762a1bSJed Brown       if (j > 0) {
398*9371c9d4SSatish Balay         v[k]   = hb;
399*9371c9d4SSatish Balay         col[k] = row - mx;
400*9371c9d4SSatish Balay         k++;
401c4762a1bSJed Brown       }
402c4762a1bSJed Brown 
403c4762a1bSJed Brown       if (j > 0 && i < mx - 1) {
404*9371c9d4SSatish Balay         v[k]   = hbr;
405*9371c9d4SSatish Balay         col[k] = row - mx + 1;
406*9371c9d4SSatish Balay         k++;
407c4762a1bSJed Brown       }
408c4762a1bSJed Brown 
409c4762a1bSJed Brown       if (i > 0) {
410*9371c9d4SSatish Balay         v[k]   = hl;
411*9371c9d4SSatish Balay         col[k] = row - 1;
412*9371c9d4SSatish Balay         k++;
413c4762a1bSJed Brown       }
414c4762a1bSJed Brown 
415*9371c9d4SSatish Balay       v[k]   = hc;
416*9371c9d4SSatish Balay       col[k] = row;
417*9371c9d4SSatish Balay       k++;
418c4762a1bSJed Brown 
419c4762a1bSJed Brown       if (i < mx - 1) {
420*9371c9d4SSatish Balay         v[k]   = hr;
421*9371c9d4SSatish Balay         col[k] = row + 1;
422*9371c9d4SSatish Balay         k++;
423c4762a1bSJed Brown       }
424c4762a1bSJed Brown 
425c4762a1bSJed Brown       if (i > 0 && j < my - 1) {
426*9371c9d4SSatish Balay         v[k]   = htl;
427*9371c9d4SSatish Balay         col[k] = row + mx - 1;
428*9371c9d4SSatish Balay         k++;
429c4762a1bSJed Brown       }
430c4762a1bSJed Brown 
431c4762a1bSJed Brown       if (j < my - 1) {
432*9371c9d4SSatish Balay         v[k]   = ht;
433*9371c9d4SSatish Balay         col[k] = row + mx;
434*9371c9d4SSatish Balay         k++;
435c4762a1bSJed Brown       }
436c4762a1bSJed Brown 
437c4762a1bSJed Brown       /*
438c4762a1bSJed Brown          Set matrix values using local numbering, which was defined
439c4762a1bSJed Brown          earlier, in the main routine.
440c4762a1bSJed Brown       */
4419566063dSJacob Faibussowitsch       PetscCall(MatSetValues(Hessian, 1, &row, k, col, v, INSERT_VALUES));
442c4762a1bSJed Brown     }
443c4762a1bSJed Brown   }
444c4762a1bSJed Brown 
445c4762a1bSJed Brown   /* Restore vectors */
4469566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(X, &x));
447c4762a1bSJed Brown 
448c4762a1bSJed Brown   /* Assemble the matrix */
4499566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(Hessian, MAT_FINAL_ASSEMBLY));
4509566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(Hessian, MAT_FINAL_ASSEMBLY));
451c4762a1bSJed Brown 
4529566063dSJacob Faibussowitsch   PetscCall(PetscLogFlops(199.0 * mx * my));
453c4762a1bSJed Brown   PetscFunctionReturn(0);
454c4762a1bSJed Brown }
455c4762a1bSJed Brown 
456c4762a1bSJed Brown /* ------------------------------------------------------------------- */
457c4762a1bSJed Brown /*
458c4762a1bSJed Brown    MSA_BoundaryConditions -  Calculates the boundary conditions for
459c4762a1bSJed Brown    the region.
460c4762a1bSJed Brown 
461c4762a1bSJed Brown    Input Parameter:
462c4762a1bSJed Brown .  user - user-defined application context
463c4762a1bSJed Brown 
464c4762a1bSJed Brown    Output Parameter:
465c4762a1bSJed Brown .  user - user-defined application context
466c4762a1bSJed Brown */
467*9371c9d4SSatish Balay static PetscErrorCode MSA_BoundaryConditions(AppCtx *user) {
468c4762a1bSJed Brown   PetscInt   i, j, k, limit = 0;
469c4762a1bSJed Brown   PetscInt   maxits = 5;
470c4762a1bSJed Brown   PetscInt   mx = user->mx, my = user->my;
471c4762a1bSJed Brown   PetscInt   bsize = 0, lsize = 0, tsize = 0, rsize = 0;
472c4762a1bSJed Brown   PetscReal  one = 1.0, two = 2.0, three = 3.0, tol = 1e-10;
473c4762a1bSJed Brown   PetscReal  fnorm, det, hx, hy, xt = 0, yt = 0;
474c4762a1bSJed Brown   PetscReal  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
475c4762a1bSJed Brown   PetscReal  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
476c4762a1bSJed Brown   PetscReal *boundary;
477c4762a1bSJed Brown 
478c4762a1bSJed Brown   PetscFunctionBeginUser;
479*9371c9d4SSatish Balay   bsize = mx + 2;
480*9371c9d4SSatish Balay   lsize = my + 2;
481*9371c9d4SSatish Balay   rsize = my + 2;
482*9371c9d4SSatish Balay   tsize = mx + 2;
483c4762a1bSJed Brown 
4849566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(bsize, &user->bottom));
4859566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(tsize, &user->top));
4869566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(lsize, &user->left));
4879566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(rsize, &user->right));
488c4762a1bSJed Brown 
489*9371c9d4SSatish Balay   hx = (r - l) / (mx + 1);
490*9371c9d4SSatish Balay   hy = (t - b) / (my + 1);
491c4762a1bSJed Brown 
492c4762a1bSJed Brown   for (j = 0; j < 4; j++) {
493c4762a1bSJed Brown     if (j == 0) {
494c4762a1bSJed Brown       yt       = b;
495c4762a1bSJed Brown       xt       = l;
496c4762a1bSJed Brown       limit    = bsize;
497c4762a1bSJed Brown       boundary = user->bottom;
498c4762a1bSJed Brown     } else if (j == 1) {
499c4762a1bSJed Brown       yt       = t;
500c4762a1bSJed Brown       xt       = l;
501c4762a1bSJed Brown       limit    = tsize;
502c4762a1bSJed Brown       boundary = user->top;
503c4762a1bSJed Brown     } else if (j == 2) {
504c4762a1bSJed Brown       yt       = b;
505c4762a1bSJed Brown       xt       = l;
506c4762a1bSJed Brown       limit    = lsize;
507c4762a1bSJed Brown       boundary = user->left;
508c4762a1bSJed Brown     } else { /*  if (j==3) */
509c4762a1bSJed Brown       yt       = b;
510c4762a1bSJed Brown       xt       = r;
511c4762a1bSJed Brown       limit    = rsize;
512c4762a1bSJed Brown       boundary = user->right;
513c4762a1bSJed Brown     }
514c4762a1bSJed Brown 
515c4762a1bSJed Brown     for (i = 0; i < limit; i++) {
516c4762a1bSJed Brown       u1 = xt;
517c4762a1bSJed Brown       u2 = -yt;
518c4762a1bSJed Brown       for (k = 0; k < maxits; k++) {
519c4762a1bSJed Brown         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
520c4762a1bSJed Brown         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
521c4762a1bSJed Brown         fnorm = PetscSqrtScalar(nf1 * nf1 + nf2 * nf2);
522c4762a1bSJed Brown         if (fnorm <= tol) break;
523c4762a1bSJed Brown         njac11 = one + u2 * u2 - u1 * u1;
524c4762a1bSJed Brown         njac12 = two * u1 * u2;
525c4762a1bSJed Brown         njac21 = -two * u1 * u2;
526c4762a1bSJed Brown         njac22 = -one - u1 * u1 + u2 * u2;
527c4762a1bSJed Brown         det    = njac11 * njac22 - njac21 * njac12;
528c4762a1bSJed Brown         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
529c4762a1bSJed Brown         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
530c4762a1bSJed Brown       }
531c4762a1bSJed Brown 
532c4762a1bSJed Brown       boundary[i] = u1 * u1 - u2 * u2;
533c4762a1bSJed Brown       if (j == 0 || j == 1) {
534c4762a1bSJed Brown         xt = xt + hx;
535c4762a1bSJed Brown       } else { /*  if (j==2 || j==3) */
536c4762a1bSJed Brown         yt = yt + hy;
537c4762a1bSJed Brown       }
538c4762a1bSJed Brown     }
539c4762a1bSJed Brown   }
540c4762a1bSJed Brown   PetscFunctionReturn(0);
541c4762a1bSJed Brown }
542c4762a1bSJed Brown 
543c4762a1bSJed Brown /* ------------------------------------------------------------------- */
544c4762a1bSJed Brown /*
545c4762a1bSJed Brown    MSA_InitialPoint - Calculates the initial guess in one of three ways.
546c4762a1bSJed Brown 
547c4762a1bSJed Brown    Input Parameters:
548c4762a1bSJed Brown .  user - user-defined application context
549c4762a1bSJed Brown .  X - vector for initial guess
550c4762a1bSJed Brown 
551c4762a1bSJed Brown    Output Parameters:
552c4762a1bSJed Brown .  X - newly computed initial guess
553c4762a1bSJed Brown */
554*9371c9d4SSatish Balay static PetscErrorCode MSA_InitialPoint(AppCtx *user, Vec X) {
555c4762a1bSJed Brown   PetscInt  start = -1, i, j;
556c4762a1bSJed Brown   PetscReal zero  = 0.0;
557c4762a1bSJed Brown   PetscBool flg;
558c4762a1bSJed Brown 
559c4762a1bSJed Brown   PetscFunctionBeginUser;
5609566063dSJacob Faibussowitsch   PetscCall(VecSet(X, zero));
5619566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetInt(NULL, NULL, "-start", &start, &flg));
562c4762a1bSJed Brown 
563c4762a1bSJed Brown   if (flg && start == 0) { /* The zero vector is reasonable */
5649566063dSJacob Faibussowitsch     PetscCall(VecSet(X, zero));
565c4762a1bSJed Brown   } else { /* Take an average of the boundary conditions */
566c4762a1bSJed Brown     PetscInt   row;
567c4762a1bSJed Brown     PetscInt   mx = user->mx, my = user->my;
568c4762a1bSJed Brown     PetscReal *x;
569c4762a1bSJed Brown 
570c4762a1bSJed Brown     /* Get pointers to vector data */
5719566063dSJacob Faibussowitsch     PetscCall(VecGetArray(X, &x));
572c4762a1bSJed Brown     /* Perform local computations */
573c4762a1bSJed Brown     for (j = 0; j < my; j++) {
574c4762a1bSJed Brown       for (i = 0; i < mx; i++) {
575c4762a1bSJed Brown         row    = (j)*mx + (i);
576c4762a1bSJed Brown         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;
577c4762a1bSJed Brown       }
578c4762a1bSJed Brown     }
579c4762a1bSJed Brown     /* Restore vectors */
5809566063dSJacob Faibussowitsch     PetscCall(VecRestoreArray(X, &x));
581c4762a1bSJed Brown   }
582c4762a1bSJed Brown   PetscFunctionReturn(0);
583c4762a1bSJed Brown }
584c4762a1bSJed Brown 
585c4762a1bSJed Brown /*TEST
586c4762a1bSJed Brown 
587c4762a1bSJed Brown    build:
588c4762a1bSJed Brown       requires: !complex
589c4762a1bSJed Brown 
590c4762a1bSJed Brown    test:
591c4762a1bSJed Brown       args: -tao_smonitor -tao_type nls -mx 10 -my 8 -tao_gatol 1.e-4
592c4762a1bSJed Brown       requires: !single
593c4762a1bSJed Brown 
594c4762a1bSJed Brown    test:
595c4762a1bSJed Brown       suffix: 2
596c4762a1bSJed Brown       args: -tao_smonitor -tao_type bmrm -mx 10 -my 8 -tao_gatol 1.e-3
597c4762a1bSJed Brown       requires: !single
598c4762a1bSJed Brown 
599c4762a1bSJed Brown    test:
600c4762a1bSJed Brown       suffix: 3
601c4762a1bSJed Brown       args: -tao_smonitor -tao_type lmvm -mx 10 -my 8 -tao_gatol 1.e-3
602c4762a1bSJed Brown       requires: !single
603c4762a1bSJed Brown 
604c4762a1bSJed Brown    test:
605c4762a1bSJed Brown       suffix: 4
606c4762a1bSJed Brown       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4
607c4762a1bSJed Brown       requires: !single
608c4762a1bSJed Brown 
609c4762a1bSJed Brown    test:
6100f0abf79SStefano Zampini       suffix: 4_ew
6110f0abf79SStefano Zampini       args: -tao_smonitor -tao_type bntr -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
6120f0abf79SStefano Zampini       requires: !single
6130f0abf79SStefano Zampini 
6140f0abf79SStefano Zampini    test:
615c4762a1bSJed Brown       suffix: 5
616c4762a1bSJed Brown       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4
617c4762a1bSJed Brown       requires: !single
618c4762a1bSJed Brown 
619c4762a1bSJed Brown    test:
6200f0abf79SStefano Zampini       suffix: 5_ew
6210f0abf79SStefano Zampini       args: -tao_smonitor -tao_type bntl -tao_ksp_ew -mx 10 -my 8 -tao_gatol 1.e-4
6220f0abf79SStefano Zampini       requires: !single
6230f0abf79SStefano Zampini 
6240f0abf79SStefano Zampini    test:
625c4762a1bSJed Brown       suffix: 6
626c4762a1bSJed Brown       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_gatol 1.e-4
627c4762a1bSJed Brown       requires: !single
628c4762a1bSJed Brown 
629c4762a1bSJed Brown    test:
6300f0abf79SStefano Zampini       suffix: 6_ew
6310f0abf79SStefano Zampini       args: -tao_smonitor -tao_type bnls -tao_ksp_ew -tao_bnk_ksp_ew_version 3 -mx 10 -my 8 -tao_gatol 1.e-4
6320f0abf79SStefano Zampini       requires: !single
6330f0abf79SStefano Zampini 
6340f0abf79SStefano Zampini    test:
635c4762a1bSJed Brown       suffix: 7
636c4762a1bSJed Brown       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
637c4762a1bSJed Brown       requires: !single
638c4762a1bSJed Brown 
639c4762a1bSJed Brown    test:
640c4762a1bSJed Brown       suffix: 8
641c4762a1bSJed Brown       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
642c4762a1bSJed Brown       requires: !single
643c4762a1bSJed Brown 
644c4762a1bSJed Brown    test:
645c4762a1bSJed Brown       suffix: 9
646c4762a1bSJed Brown       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4
647c4762a1bSJed Brown       requires: !single
648c4762a1bSJed Brown 
64934ad9904SAlp Dener    test:
65034ad9904SAlp Dener       suffix: 10
65134ad9904SAlp Dener       args: -tao_smonitor -tao_type bnls -mx 10 -my 8 -tao_bnk_max_cg_its 3 -tao_gatol 1.e-4 -tao_mf_hessian
65234ad9904SAlp Dener 
65334ad9904SAlp Dener    test:
65434ad9904SAlp Dener       suffix: 11
65534ad9904SAlp Dener       args: -tao_smonitor -tao_type bntr -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
65634ad9904SAlp Dener       requires: !single
65734ad9904SAlp Dener 
65834ad9904SAlp Dener    test:
65934ad9904SAlp Dener       suffix: 12
66034ad9904SAlp Dener       args: -tao_smonitor -tao_type bntl -mx 10 -my 8 -tao_gatol 1.e-4 -tao_mf_hessian
66134ad9904SAlp Dener       requires: !single
66234ad9904SAlp Dener 
663c4762a1bSJed Brown TEST*/
664