xref: /petsc/src/tao/unconstrained/tutorials/eptorsion1.c (revision b122ec5aa1bd4469eb4e0673542fb7de3f411254)
1c4762a1bSJed Brown /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */
2c4762a1bSJed Brown 
3c4762a1bSJed Brown /* ----------------------------------------------------------------------
4c4762a1bSJed Brown 
5c4762a1bSJed Brown   Elastic-plastic torsion problem.
6c4762a1bSJed Brown 
7c4762a1bSJed Brown   The elastic plastic torsion problem arises from the determination
8c4762a1bSJed Brown   of the stress field on an infinitely long cylindrical bar, which is
9c4762a1bSJed Brown   equivalent to the solution of the following problem:
10c4762a1bSJed Brown 
11c4762a1bSJed Brown   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
12c4762a1bSJed Brown 
13c4762a1bSJed Brown   where C is the torsion angle per unit length.
14c4762a1bSJed Brown 
15c4762a1bSJed Brown   The multiprocessor version of this code is eptorsion2.c.
16c4762a1bSJed Brown 
17c4762a1bSJed Brown ---------------------------------------------------------------------- */
18c4762a1bSJed Brown 
19c4762a1bSJed Brown /*
20c4762a1bSJed Brown   Include "petsctao.h" so that we can use TAO solvers.  Note that this
21c4762a1bSJed Brown   file automatically includes files for lower-level support, such as those
22c4762a1bSJed Brown   provided by the PETSc library:
23c4762a1bSJed Brown      petsc.h       - base PETSc routines   petscvec.h - vectors
24a5b23f4aSJose E. Roman      petscsys.h    - system routines        petscmat.h - matrices
25c4762a1bSJed Brown      petscis.h     - index sets            petscksp.h - Krylov subspace methods
26c4762a1bSJed Brown      petscviewer.h - viewers               petscpc.h  - preconditioners
27c4762a1bSJed Brown */
28c4762a1bSJed Brown 
29c4762a1bSJed Brown #include <petsctao.h>
30c4762a1bSJed Brown 
31c4762a1bSJed Brown static  char help[]=
32c4762a1bSJed Brown "Demonstrates use of the TAO package to solve \n\
33c4762a1bSJed Brown unconstrained minimization problems on a single processor.  This example \n\
34c4762a1bSJed Brown is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
35c4762a1bSJed Brown test suite.\n\
36c4762a1bSJed Brown The command line options are:\n\
37c4762a1bSJed Brown   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
38c4762a1bSJed Brown   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
39c4762a1bSJed Brown   -par <param>, where <param> = angle of twist per unit length\n\n";
40c4762a1bSJed Brown 
41c4762a1bSJed Brown /*T
42c4762a1bSJed Brown    Concepts: TAO^Solving an unconstrained minimization problem
43c4762a1bSJed Brown    Routines: TaoCreate(); TaoSetType();
44a82e8c82SStefano Zampini    Routines: TaoSetSolution();
45a82e8c82SStefano Zampini    Routines: TaoSetObjectiveAndGradient();
46a82e8c82SStefano Zampini    Routines: TaoSetHessian(); TaoSetFromOptions();
47c4762a1bSJed Brown    Routines: TaoGetKSP(); TaoSolve();
48c4762a1bSJed Brown    Routines: TaoDestroy();
49c4762a1bSJed Brown    Processors: 1
50c4762a1bSJed Brown T*/
51c4762a1bSJed Brown 
52c4762a1bSJed Brown /*
53c4762a1bSJed Brown    User-defined application context - contains data needed by the
54c4762a1bSJed Brown    application-provided call-back routines, FormFunction(),
55c4762a1bSJed Brown    FormGradient(), and FormHessian().
56c4762a1bSJed Brown */
57c4762a1bSJed Brown 
58c4762a1bSJed Brown typedef struct {
59c4762a1bSJed Brown    PetscReal  param;      /* nonlinearity parameter */
60c4762a1bSJed Brown    PetscInt   mx, my;     /* discretization in x- and y-directions */
61c4762a1bSJed Brown    PetscInt   ndim;       /* problem dimension */
62c4762a1bSJed Brown    Vec        s, y, xvec; /* work space for computing Hessian */
63c4762a1bSJed Brown    PetscReal  hx, hy;     /* mesh spacing in x- and y-directions */
64c4762a1bSJed Brown } AppCtx;
65c4762a1bSJed Brown 
66c4762a1bSJed Brown /* -------- User-defined Routines --------- */
67c4762a1bSJed Brown 
68c4762a1bSJed Brown PetscErrorCode FormInitialGuess(AppCtx*,Vec);
69c4762a1bSJed Brown PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
70c4762a1bSJed Brown PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
71c4762a1bSJed Brown PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
72c4762a1bSJed Brown PetscErrorCode HessianProductMat(Mat,Vec,Vec);
73c4762a1bSJed Brown PetscErrorCode HessianProduct(void*,Vec,Vec);
74c4762a1bSJed Brown PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
75c4762a1bSJed Brown PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);
76c4762a1bSJed Brown 
77c4762a1bSJed Brown PetscErrorCode main(int argc,char **argv)
78c4762a1bSJed Brown {
79c4762a1bSJed Brown   PetscInt           mx=10;               /* discretization in x-direction */
80c4762a1bSJed Brown   PetscInt           my=10;               /* discretization in y-direction */
81c4762a1bSJed Brown   Vec                x;                   /* solution, gradient vectors */
82c4762a1bSJed Brown   PetscBool          flg;                 /* A return value when checking for use options */
83c4762a1bSJed Brown   Tao                tao;                 /* Tao solver context */
84c4762a1bSJed Brown   Mat                H;                   /* Hessian matrix */
85c4762a1bSJed Brown   AppCtx             user;                /* application context */
86c4762a1bSJed Brown   PetscMPIInt        size;                /* number of processes */
87c4762a1bSJed Brown   PetscReal          one=1.0;
88c4762a1bSJed Brown 
89c4762a1bSJed Brown   PetscBool          test_lmvm = PETSC_FALSE;
90c4762a1bSJed Brown   KSP                ksp;
91c4762a1bSJed Brown   PC                 pc;
92c4762a1bSJed Brown   Mat                M;
93c4762a1bSJed Brown   Vec                in, out, out2;
94c4762a1bSJed Brown   PetscReal          mult_solve_dist;
95c4762a1bSJed Brown 
96c4762a1bSJed Brown   /* Initialize TAO,PETSc */
97*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscInitialize(&argc,&argv,(char *)0,help));
985f80ce2aSJacob Faibussowitsch   CHKERRMPI(MPI_Comm_size(MPI_COMM_WORLD,&size));
993c859ba3SBarry Smith   PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Incorrect number of processors");
100c4762a1bSJed Brown 
101c4762a1bSJed Brown   /* Specify default parameters for the problem, check for command-line overrides */
102c4762a1bSJed Brown   user.param = 5.0;
1035f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg));
1045f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg));
1055f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg));
1065f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_lmvm",&test_lmvm,&flg));
107c4762a1bSJed Brown 
1085f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n"));
1095f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",mx,my));
110c4762a1bSJed Brown   user.ndim = mx * my; user.mx = mx; user.my = my;
111c4762a1bSJed Brown   user.hx = one/(mx+1); user.hy = one/(my+1);
112c4762a1bSJed Brown 
113c4762a1bSJed Brown   /* Allocate vectors */
1145f80ce2aSJacob Faibussowitsch   CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y));
1155f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user.y,&user.s));
1165f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDuplicate(user.y,&x));
117c4762a1bSJed Brown 
118c4762a1bSJed Brown   /* Create TAO solver and set desired solution method */
1195f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoCreate(PETSC_COMM_SELF,&tao));
1205f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetType(tao,TAOLMVM));
121c4762a1bSJed Brown 
122c4762a1bSJed Brown   /* Set solution vector with an initial guess */
1235f80ce2aSJacob Faibussowitsch   CHKERRQ(FormInitialGuess(&user,x));
1245f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetSolution(tao,x));
125c4762a1bSJed Brown 
126c4762a1bSJed Brown   /* Set routine for function and gradient evaluation */
1275f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user));
128c4762a1bSJed Brown 
129c4762a1bSJed Brown   /* From command line options, determine if using matrix-free hessian */
1305f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-my_tao_mf",&flg));
131c4762a1bSJed Brown   if (flg) {
1325f80ce2aSJacob Faibussowitsch     CHKERRQ(MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H));
1335f80ce2aSJacob Faibussowitsch     CHKERRQ(MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat));
1345f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE));
135c4762a1bSJed Brown 
1365f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetHessian(tao,H,H,MatrixFreeHessian,(void *)&user));
137c4762a1bSJed Brown   } else {
1385f80ce2aSJacob Faibussowitsch     CHKERRQ(MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H));
1395f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE));
1405f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoSetHessian(tao,H,H,FormHessian,(void *)&user));
141c4762a1bSJed Brown   }
142c4762a1bSJed Brown 
143c4762a1bSJed Brown   /* Test the LMVM matrix */
144c4762a1bSJed Brown   if (test_lmvm) {
1455f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsSetValue(NULL, "-tao_type", "bntr"));
1465f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm"));
147c4762a1bSJed Brown   }
148c4762a1bSJed Brown 
149c4762a1bSJed Brown   /* Check for any TAO command line options */
1505f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSetFromOptions(tao));
151c4762a1bSJed Brown 
152c4762a1bSJed Brown   /* SOLVE THE APPLICATION */
1535f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoSolve(tao));
154c4762a1bSJed Brown 
155c4762a1bSJed Brown   /* Test the LMVM matrix */
156c4762a1bSJed Brown   if (test_lmvm) {
1575f80ce2aSJacob Faibussowitsch     CHKERRQ(TaoGetKSP(tao, &ksp));
1585f80ce2aSJacob Faibussowitsch     CHKERRQ(KSPGetPC(ksp, &pc));
1595f80ce2aSJacob Faibussowitsch     CHKERRQ(PCLMVMGetMatLMVM(pc, &M));
1605f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDuplicate(x, &in));
1615f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDuplicate(x, &out));
1625f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDuplicate(x, &out2));
1635f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSet(in, 5.0));
1645f80ce2aSJacob Faibussowitsch     CHKERRQ(MatMult(M, in, out));
1655f80ce2aSJacob Faibussowitsch     CHKERRQ(MatSolve(M, out, out2));
1665f80ce2aSJacob Faibussowitsch     CHKERRQ(VecAXPY(out2, -1.0, in));
1675f80ce2aSJacob Faibussowitsch     CHKERRQ(VecNorm(out2, NORM_2, &mult_solve_dist));
1685f80ce2aSJacob Faibussowitsch     CHKERRQ(PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", mult_solve_dist));
1695f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&in));
1705f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&out));
1715f80ce2aSJacob Faibussowitsch     CHKERRQ(VecDestroy(&out2));
172c4762a1bSJed Brown   }
173c4762a1bSJed Brown 
1745f80ce2aSJacob Faibussowitsch   CHKERRQ(TaoDestroy(&tao));
1755f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.s));
1765f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&user.y));
1775f80ce2aSJacob Faibussowitsch   CHKERRQ(VecDestroy(&x));
1785f80ce2aSJacob Faibussowitsch   CHKERRQ(MatDestroy(&H));
179c4762a1bSJed Brown 
180*b122ec5aSJacob Faibussowitsch   CHKERRQ(PetscFinalize());
181*b122ec5aSJacob Faibussowitsch   return 0;
182c4762a1bSJed Brown }
183c4762a1bSJed Brown 
184c4762a1bSJed Brown /* ------------------------------------------------------------------- */
185c4762a1bSJed Brown /*
186c4762a1bSJed Brown     FormInitialGuess - Computes an initial approximation to the solution.
187c4762a1bSJed Brown 
188c4762a1bSJed Brown     Input Parameters:
189c4762a1bSJed Brown .   user - user-defined application context
190c4762a1bSJed Brown .   X    - vector
191c4762a1bSJed Brown 
192c4762a1bSJed Brown     Output Parameters:
193c4762a1bSJed Brown .   X    - vector
194c4762a1bSJed Brown */
195c4762a1bSJed Brown PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
196c4762a1bSJed Brown {
197c4762a1bSJed Brown   PetscReal      hx = user->hx, hy = user->hy, temp;
198c4762a1bSJed Brown   PetscReal      val;
199c4762a1bSJed Brown   PetscInt       i, j, k, nx = user->mx, ny = user->my;
200c4762a1bSJed Brown 
201c4762a1bSJed Brown   /* Compute initial guess */
202c4762a1bSJed Brown   PetscFunctionBeginUser;
203c4762a1bSJed Brown   for (j=0; j<ny; j++) {
204c4762a1bSJed Brown     temp = PetscMin(j+1,ny-j)*hy;
205c4762a1bSJed Brown     for (i=0; i<nx; i++) {
206c4762a1bSJed Brown       k   = nx*j + i;
207c4762a1bSJed Brown       val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
2085f80ce2aSJacob Faibussowitsch       CHKERRQ(VecSetValues(X,1,&k,&val,ADD_VALUES));
209c4762a1bSJed Brown     }
210c4762a1bSJed Brown   }
2115f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(X));
2125f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(X));
213c4762a1bSJed Brown   PetscFunctionReturn(0);
214c4762a1bSJed Brown }
215c4762a1bSJed Brown 
216c4762a1bSJed Brown /* ------------------------------------------------------------------- */
217c4762a1bSJed Brown /*
218c4762a1bSJed Brown    FormFunctionGradient - Evaluates the function and corresponding gradient.
219c4762a1bSJed Brown 
220c4762a1bSJed Brown    Input Parameters:
221c4762a1bSJed Brown    tao - the Tao context
222c4762a1bSJed Brown    X   - the input vector
223c4762a1bSJed Brown    ptr - optional user-defined context, as set by TaoSetFunction()
224c4762a1bSJed Brown 
225c4762a1bSJed Brown    Output Parameters:
226c4762a1bSJed Brown    f   - the newly evaluated function
227c4762a1bSJed Brown    G   - the newly evaluated gradient
228c4762a1bSJed Brown */
229c4762a1bSJed Brown PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
230c4762a1bSJed Brown {
231c4762a1bSJed Brown   PetscFunctionBeginUser;
2325f80ce2aSJacob Faibussowitsch   CHKERRQ(FormFunction(tao,X,f,ptr));
2335f80ce2aSJacob Faibussowitsch   CHKERRQ(FormGradient(tao,X,G,ptr));
234c4762a1bSJed Brown   PetscFunctionReturn(0);
235c4762a1bSJed Brown }
236c4762a1bSJed Brown 
237c4762a1bSJed Brown /* ------------------------------------------------------------------- */
238c4762a1bSJed Brown /*
239c4762a1bSJed Brown    FormFunction - Evaluates the function, f(X).
240c4762a1bSJed Brown 
241c4762a1bSJed Brown    Input Parameters:
242c4762a1bSJed Brown .  tao - the Tao context
243c4762a1bSJed Brown .  X   - the input vector
244c4762a1bSJed Brown .  ptr - optional user-defined context, as set by TaoSetFunction()
245c4762a1bSJed Brown 
246c4762a1bSJed Brown    Output Parameters:
247c4762a1bSJed Brown .  f    - the newly evaluated function
248c4762a1bSJed Brown */
249c4762a1bSJed Brown PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
250c4762a1bSJed Brown {
251c4762a1bSJed Brown   AppCtx            *user = (AppCtx *) ptr;
252c4762a1bSJed Brown   PetscReal         hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
253c4762a1bSJed Brown   PetscReal         zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
254c4762a1bSJed Brown   PetscReal         v, cdiv3 = user->param/three;
255c4762a1bSJed Brown   const PetscScalar *x;
256c4762a1bSJed Brown   PetscInt          nx = user->mx, ny = user->my, i, j, k;
257c4762a1bSJed Brown 
258c4762a1bSJed Brown   PetscFunctionBeginUser;
259c4762a1bSJed Brown   /* Get pointer to vector data */
2605f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(X,&x));
261c4762a1bSJed Brown 
262c4762a1bSJed Brown   /* Compute function contributions over the lower triangular elements */
263c4762a1bSJed Brown   for (j=-1; j<ny; j++) {
264c4762a1bSJed Brown     for (i=-1; i<nx; i++) {
265c4762a1bSJed Brown       k = nx*j + i;
266c4762a1bSJed Brown       v = zero;
267c4762a1bSJed Brown       vr = zero;
268c4762a1bSJed Brown       vt = zero;
269c4762a1bSJed Brown       if (i >= 0 && j >= 0) v = x[k];
270c4762a1bSJed Brown       if (i < nx-1 && j > -1) vr = x[k+1];
271c4762a1bSJed Brown       if (i > -1 && j < ny-1) vt = x[k+nx];
272c4762a1bSJed Brown       dvdx = (vr-v)/hx;
273c4762a1bSJed Brown       dvdy = (vt-v)/hy;
274c4762a1bSJed Brown       fquad += dvdx*dvdx + dvdy*dvdy;
275c4762a1bSJed Brown       flin -= cdiv3*(v+vr+vt);
276c4762a1bSJed Brown     }
277c4762a1bSJed Brown   }
278c4762a1bSJed Brown 
279c4762a1bSJed Brown   /* Compute function contributions over the upper triangular elements */
280c4762a1bSJed Brown   for (j=0; j<=ny; j++) {
281c4762a1bSJed Brown     for (i=0; i<=nx; i++) {
282c4762a1bSJed Brown       k = nx*j + i;
283c4762a1bSJed Brown       vb = zero;
284c4762a1bSJed Brown       vl = zero;
285c4762a1bSJed Brown       v = zero;
286c4762a1bSJed Brown       if (i < nx && j > 0) vb = x[k-nx];
287c4762a1bSJed Brown       if (i > 0 && j < ny) vl = x[k-1];
288c4762a1bSJed Brown       if (i < nx && j < ny) v = x[k];
289c4762a1bSJed Brown       dvdx = (v-vl)/hx;
290c4762a1bSJed Brown       dvdy = (v-vb)/hy;
291c4762a1bSJed Brown       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
292c4762a1bSJed Brown       flin = flin - cdiv3*(vb+vl+v);
293c4762a1bSJed Brown     }
294c4762a1bSJed Brown   }
295c4762a1bSJed Brown 
296c4762a1bSJed Brown   /* Restore vector */
2975f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(X,&x));
298c4762a1bSJed Brown 
299c4762a1bSJed Brown   /* Scale the function */
300c4762a1bSJed Brown   area = p5*hx*hy;
301c4762a1bSJed Brown   *f = area*(p5*fquad+flin);
302c4762a1bSJed Brown 
3035f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogFlops(24.0*nx*ny));
304c4762a1bSJed Brown   PetscFunctionReturn(0);
305c4762a1bSJed Brown }
306c4762a1bSJed Brown 
307c4762a1bSJed Brown /* ------------------------------------------------------------------- */
308c4762a1bSJed Brown /*
309c4762a1bSJed Brown     FormGradient - Evaluates the gradient, G(X).
310c4762a1bSJed Brown 
311c4762a1bSJed Brown     Input Parameters:
312c4762a1bSJed Brown .   tao  - the Tao context
313c4762a1bSJed Brown .   X    - input vector
314c4762a1bSJed Brown .   ptr  - optional user-defined context
315c4762a1bSJed Brown 
316c4762a1bSJed Brown     Output Parameters:
317c4762a1bSJed Brown .   G - vector containing the newly evaluated gradient
318c4762a1bSJed Brown */
319c4762a1bSJed Brown PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
320c4762a1bSJed Brown {
321c4762a1bSJed Brown   AppCtx            *user = (AppCtx *) ptr;
322c4762a1bSJed Brown   PetscReal         zero=0.0, p5=0.5, three = 3.0, area, val;
323c4762a1bSJed Brown   PetscInt          nx = user->mx, ny = user->my, ind, i, j, k;
324c4762a1bSJed Brown   PetscReal         hx = user->hx, hy = user->hy;
325c4762a1bSJed Brown   PetscReal         vb, vl, vr, vt, dvdx, dvdy;
326c4762a1bSJed Brown   PetscReal         v, cdiv3 = user->param/three;
327c4762a1bSJed Brown   const PetscScalar *x;
328c4762a1bSJed Brown 
329c4762a1bSJed Brown   PetscFunctionBeginUser;
330c4762a1bSJed Brown   /* Initialize gradient to zero */
3315f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(G, zero));
332c4762a1bSJed Brown 
333c4762a1bSJed Brown   /* Get pointer to vector data */
3345f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(X,&x));
335c4762a1bSJed Brown 
336c4762a1bSJed Brown   /* Compute gradient contributions over the lower triangular elements */
337c4762a1bSJed Brown   for (j=-1; j<ny; j++) {
338c4762a1bSJed Brown     for (i=-1; i<nx; i++) {
339c4762a1bSJed Brown       k  = nx*j + i;
340c4762a1bSJed Brown       v  = zero;
341c4762a1bSJed Brown       vr = zero;
342c4762a1bSJed Brown       vt = zero;
343c4762a1bSJed Brown       if (i >= 0 && j >= 0)    v = x[k];
344c4762a1bSJed Brown       if (i < nx-1 && j > -1) vr = x[k+1];
345c4762a1bSJed Brown       if (i > -1 && j < ny-1) vt = x[k+nx];
346c4762a1bSJed Brown       dvdx = (vr-v)/hx;
347c4762a1bSJed Brown       dvdy = (vt-v)/hy;
348c4762a1bSJed Brown       if (i != -1 && j != -1) {
349c4762a1bSJed Brown         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
3505f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
351c4762a1bSJed Brown       }
352c4762a1bSJed Brown       if (i != nx-1 && j != -1) {
353c4762a1bSJed Brown         ind = k+1; val =  dvdx/hx - cdiv3;
3545f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
355c4762a1bSJed Brown       }
356c4762a1bSJed Brown       if (i != -1 && j != ny-1) {
357c4762a1bSJed Brown         ind = k+nx; val = dvdy/hy - cdiv3;
3585f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
359c4762a1bSJed Brown       }
360c4762a1bSJed Brown     }
361c4762a1bSJed Brown   }
362c4762a1bSJed Brown 
363c4762a1bSJed Brown   /* Compute gradient contributions over the upper triangular elements */
364c4762a1bSJed Brown   for (j=0; j<=ny; j++) {
365c4762a1bSJed Brown     for (i=0; i<=nx; i++) {
366c4762a1bSJed Brown       k = nx*j + i;
367c4762a1bSJed Brown       vb = zero;
368c4762a1bSJed Brown       vl = zero;
369c4762a1bSJed Brown       v  = zero;
370c4762a1bSJed Brown       if (i < nx && j > 0) vb = x[k-nx];
371c4762a1bSJed Brown       if (i > 0 && j < ny) vl = x[k-1];
372c4762a1bSJed Brown       if (i < nx && j < ny) v = x[k];
373c4762a1bSJed Brown       dvdx = (v-vl)/hx;
374c4762a1bSJed Brown       dvdy = (v-vb)/hy;
375c4762a1bSJed Brown       if (i != nx && j != 0) {
376c4762a1bSJed Brown         ind = k-nx; val = - dvdy/hy - cdiv3;
3775f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
378c4762a1bSJed Brown       }
379c4762a1bSJed Brown       if (i != 0 && j != ny) {
380c4762a1bSJed Brown         ind = k-1; val =  - dvdx/hx - cdiv3;
3815f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
382c4762a1bSJed Brown       }
383c4762a1bSJed Brown       if (i != nx && j != ny) {
384c4762a1bSJed Brown         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
3855f80ce2aSJacob Faibussowitsch         CHKERRQ(VecSetValues(G,1,&ind,&val,ADD_VALUES));
386c4762a1bSJed Brown       }
387c4762a1bSJed Brown     }
388c4762a1bSJed Brown   }
3895f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(X,&x));
390c4762a1bSJed Brown 
391c4762a1bSJed Brown   /* Assemble gradient vector */
3925f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(G));
3935f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(G));
394c4762a1bSJed Brown 
395c4762a1bSJed Brown   /* Scale the gradient */
396c4762a1bSJed Brown   area = p5*hx*hy;
3975f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScale(G, area));
3985f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogFlops(24.0*nx*ny));
399c4762a1bSJed Brown   PetscFunctionReturn(0);
400c4762a1bSJed Brown }
401c4762a1bSJed Brown 
402c4762a1bSJed Brown /* ------------------------------------------------------------------- */
403c4762a1bSJed Brown /*
404c4762a1bSJed Brown    FormHessian - Forms the Hessian matrix.
405c4762a1bSJed Brown 
406c4762a1bSJed Brown    Input Parameters:
407c4762a1bSJed Brown .  tao - the Tao context
408c4762a1bSJed Brown .  X    - the input vector
409c4762a1bSJed Brown .  ptr  - optional user-defined context, as set by TaoSetHessian()
410c4762a1bSJed Brown 
411c4762a1bSJed Brown    Output Parameters:
412c4762a1bSJed Brown .  H     - Hessian matrix
413c4762a1bSJed Brown .  PrecH - optionally different preconditioning Hessian
414c4762a1bSJed Brown .  flag  - flag indicating matrix structure
415c4762a1bSJed Brown 
416c4762a1bSJed Brown    Notes:
417c4762a1bSJed Brown    This routine is intended simply as an example of the interface
418c4762a1bSJed Brown    to a Hessian evaluation routine.  Since this example compute the
419c4762a1bSJed Brown    Hessian a column at a time, it is not particularly efficient and
420c4762a1bSJed Brown    is not recommended.
421c4762a1bSJed Brown */
422c4762a1bSJed Brown PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
423c4762a1bSJed Brown {
424c4762a1bSJed Brown   AppCtx         *user = (AppCtx *) ptr;
425c4762a1bSJed Brown   PetscInt       i,j, ndim = user->ndim;
426c4762a1bSJed Brown   PetscReal      *y, zero = 0.0, one = 1.0;
427c4762a1bSJed Brown   PetscBool      assembled;
428c4762a1bSJed Brown 
429c4762a1bSJed Brown   PetscFunctionBeginUser;
430c4762a1bSJed Brown   user->xvec = X;
431c4762a1bSJed Brown 
432c4762a1bSJed Brown   /* Initialize Hessian entries and work vector to zero */
4335f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssembled(H,&assembled));
4345f80ce2aSJacob Faibussowitsch   if (assembled)CHKERRQ(MatZeroEntries(H));
435c4762a1bSJed Brown 
4365f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(user->s, zero));
437c4762a1bSJed Brown 
438c4762a1bSJed Brown   /* Loop over matrix columns to compute entries of the Hessian */
439c4762a1bSJed Brown   for (j=0; j<ndim; j++) {
4405f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValues(user->s,1,&j,&one,INSERT_VALUES));
4415f80ce2aSJacob Faibussowitsch     CHKERRQ(VecAssemblyBegin(user->s));
4425f80ce2aSJacob Faibussowitsch     CHKERRQ(VecAssemblyEnd(user->s));
443c4762a1bSJed Brown 
4445f80ce2aSJacob Faibussowitsch     CHKERRQ(HessianProduct(ptr,user->s,user->y));
445c4762a1bSJed Brown 
4465f80ce2aSJacob Faibussowitsch     CHKERRQ(VecSetValues(user->s,1,&j,&zero,INSERT_VALUES));
4475f80ce2aSJacob Faibussowitsch     CHKERRQ(VecAssemblyBegin(user->s));
4485f80ce2aSJacob Faibussowitsch     CHKERRQ(VecAssemblyEnd(user->s));
449c4762a1bSJed Brown 
4505f80ce2aSJacob Faibussowitsch     CHKERRQ(VecGetArray(user->y,&y));
451c4762a1bSJed Brown     for (i=0; i<ndim; i++) {
452c4762a1bSJed Brown       if (y[i] != zero) {
4535f80ce2aSJacob Faibussowitsch         CHKERRQ(MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES));
454c4762a1bSJed Brown       }
455c4762a1bSJed Brown     }
4565f80ce2aSJacob Faibussowitsch     CHKERRQ(VecRestoreArray(user->y,&y));
457c4762a1bSJed Brown   }
4585f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY));
4595f80ce2aSJacob Faibussowitsch   CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY));
460c4762a1bSJed Brown   PetscFunctionReturn(0);
461c4762a1bSJed Brown }
462c4762a1bSJed Brown 
463c4762a1bSJed Brown /* ------------------------------------------------------------------- */
464c4762a1bSJed Brown /*
465c4762a1bSJed Brown    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
466c4762a1bSJed Brown    products.
467c4762a1bSJed Brown 
468c4762a1bSJed Brown    Input Parameters:
469c4762a1bSJed Brown .  tao - the Tao context
470c4762a1bSJed Brown .  X    - the input vector
471c4762a1bSJed Brown .  ptr  - optional user-defined context, as set by TaoSetHessian()
472c4762a1bSJed Brown 
473c4762a1bSJed Brown    Output Parameters:
474c4762a1bSJed Brown .  H     - Hessian matrix
475c4762a1bSJed Brown .  PrecH - optionally different preconditioning Hessian
476c4762a1bSJed Brown .  flag  - flag indicating matrix structure
477c4762a1bSJed Brown */
478c4762a1bSJed Brown PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
479c4762a1bSJed Brown {
480c4762a1bSJed Brown   AppCtx *user = (AppCtx *) ptr;
481c4762a1bSJed Brown 
482c4762a1bSJed Brown   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
483362febeeSStefano Zampini   PetscFunctionBeginUser;
484c4762a1bSJed Brown   user->xvec = X;
485c4762a1bSJed Brown   PetscFunctionReturn(0);
486c4762a1bSJed Brown }
487c4762a1bSJed Brown 
488c4762a1bSJed Brown /* ------------------------------------------------------------------- */
489c4762a1bSJed Brown /*
490c4762a1bSJed Brown    HessianProductMat - Computes the matrix-vector product
491c4762a1bSJed Brown    y = mat*svec.
492c4762a1bSJed Brown 
493c4762a1bSJed Brown    Input Parameters:
494c4762a1bSJed Brown .  mat  - input matrix
495c4762a1bSJed Brown .  svec - input vector
496c4762a1bSJed Brown 
497c4762a1bSJed Brown    Output Parameters:
498c4762a1bSJed Brown .  y    - solution vector
499c4762a1bSJed Brown */
500c4762a1bSJed Brown PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
501c4762a1bSJed Brown {
502c4762a1bSJed Brown   void           *ptr;
503c4762a1bSJed Brown 
504c4762a1bSJed Brown   PetscFunctionBeginUser;
5055f80ce2aSJacob Faibussowitsch   CHKERRQ(MatShellGetContext(mat,&ptr));
5065f80ce2aSJacob Faibussowitsch   CHKERRQ(HessianProduct(ptr,svec,y));
507c4762a1bSJed Brown   PetscFunctionReturn(0);
508c4762a1bSJed Brown }
509c4762a1bSJed Brown 
510c4762a1bSJed Brown /* ------------------------------------------------------------------- */
511c4762a1bSJed Brown /*
512c4762a1bSJed Brown    Hessian Product - Computes the matrix-vector product:
513c4762a1bSJed Brown    y = f''(x)*svec.
514c4762a1bSJed Brown 
5157a7aea1fSJed Brown    Input Parameters:
516c4762a1bSJed Brown .  ptr  - pointer to the user-defined context
517c4762a1bSJed Brown .  svec - input vector
518c4762a1bSJed Brown 
519c4762a1bSJed Brown    Output Parameters:
520c4762a1bSJed Brown .  y    - product vector
521c4762a1bSJed Brown */
522c4762a1bSJed Brown PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
523c4762a1bSJed Brown {
524c4762a1bSJed Brown   AppCtx            *user = (AppCtx *)ptr;
525c4762a1bSJed Brown   PetscReal         p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
526c4762a1bSJed Brown   const PetscScalar *x, *s;
527c4762a1bSJed Brown   PetscReal         v, vb, vl, vr, vt, hxhx, hyhy;
528c4762a1bSJed Brown   PetscInt          nx, ny, i, j, k, ind;
529c4762a1bSJed Brown 
530c4762a1bSJed Brown   PetscFunctionBeginUser;
531c4762a1bSJed Brown   nx   = user->mx;
532c4762a1bSJed Brown   ny   = user->my;
533c4762a1bSJed Brown   hx   = user->hx;
534c4762a1bSJed Brown   hy   = user->hy;
535c4762a1bSJed Brown   hxhx = one/(hx*hx);
536c4762a1bSJed Brown   hyhy = one/(hy*hy);
537c4762a1bSJed Brown 
538c4762a1bSJed Brown   /* Get pointers to vector data */
5395f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(user->xvec,&x));
5405f80ce2aSJacob Faibussowitsch   CHKERRQ(VecGetArrayRead(svec,&s));
541c4762a1bSJed Brown 
542c4762a1bSJed Brown   /* Initialize product vector to zero */
5435f80ce2aSJacob Faibussowitsch   CHKERRQ(VecSet(y, zero));
544c4762a1bSJed Brown 
545c4762a1bSJed Brown   /* Compute f''(x)*s over the lower triangular elements */
546c4762a1bSJed Brown   for (j=-1; j<ny; j++) {
547c4762a1bSJed Brown     for (i=-1; i<nx; i++) {
548c4762a1bSJed Brown        k = nx*j + i;
549c4762a1bSJed Brown        v = zero;
550c4762a1bSJed Brown        vr = zero;
551c4762a1bSJed Brown        vt = zero;
552c4762a1bSJed Brown        if (i != -1 && j != -1) v = s[k];
553c4762a1bSJed Brown        if (i != nx-1 && j != -1) {
554c4762a1bSJed Brown          vr = s[k+1];
555c4762a1bSJed Brown          ind = k+1; val = hxhx*(vr-v);
5565f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
557c4762a1bSJed Brown        }
558c4762a1bSJed Brown        if (i != -1 && j != ny-1) {
559c4762a1bSJed Brown          vt = s[k+nx];
560c4762a1bSJed Brown          ind = k+nx; val = hyhy*(vt-v);
5615f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
562c4762a1bSJed Brown        }
563c4762a1bSJed Brown        if (i != -1 && j != -1) {
564c4762a1bSJed Brown          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
5655f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
566c4762a1bSJed Brown        }
567c4762a1bSJed Brown      }
568c4762a1bSJed Brown    }
569c4762a1bSJed Brown 
570c4762a1bSJed Brown   /* Compute f''(x)*s over the upper triangular elements */
571c4762a1bSJed Brown   for (j=0; j<=ny; j++) {
572c4762a1bSJed Brown     for (i=0; i<=nx; i++) {
573c4762a1bSJed Brown        k = nx*j + i;
574c4762a1bSJed Brown        v = zero;
575c4762a1bSJed Brown        vl = zero;
576c4762a1bSJed Brown        vb = zero;
577c4762a1bSJed Brown        if (i != nx && j != ny) v = s[k];
578c4762a1bSJed Brown        if (i != nx && j != 0) {
579c4762a1bSJed Brown          vb = s[k-nx];
580c4762a1bSJed Brown          ind = k-nx; val = hyhy*(vb-v);
5815f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
582c4762a1bSJed Brown        }
583c4762a1bSJed Brown        if (i != 0 && j != ny) {
584c4762a1bSJed Brown          vl = s[k-1];
585c4762a1bSJed Brown          ind = k-1; val = hxhx*(vl-v);
5865f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
587c4762a1bSJed Brown        }
588c4762a1bSJed Brown        if (i != nx && j != ny) {
589c4762a1bSJed Brown          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
5905f80ce2aSJacob Faibussowitsch          CHKERRQ(VecSetValues(y,1,&ind,&val,ADD_VALUES));
591c4762a1bSJed Brown        }
592c4762a1bSJed Brown     }
593c4762a1bSJed Brown   }
594c4762a1bSJed Brown   /* Restore vector data */
5955f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(svec,&s));
5965f80ce2aSJacob Faibussowitsch   CHKERRQ(VecRestoreArrayRead(user->xvec,&x));
597c4762a1bSJed Brown 
598c4762a1bSJed Brown   /* Assemble vector */
5995f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyBegin(y));
6005f80ce2aSJacob Faibussowitsch   CHKERRQ(VecAssemblyEnd(y));
601c4762a1bSJed Brown 
602c4762a1bSJed Brown   /* Scale resulting vector by area */
603c4762a1bSJed Brown   area = p5*hx*hy;
6045f80ce2aSJacob Faibussowitsch   CHKERRQ(VecScale(y, area));
6055f80ce2aSJacob Faibussowitsch   CHKERRQ(PetscLogFlops(18.0*nx*ny));
606c4762a1bSJed Brown   PetscFunctionReturn(0);
607c4762a1bSJed Brown }
608c4762a1bSJed Brown 
609c4762a1bSJed Brown /*TEST
610c4762a1bSJed Brown 
611c4762a1bSJed Brown    build:
612c4762a1bSJed Brown       requires: !complex
613c4762a1bSJed Brown 
614c4762a1bSJed Brown    test:
615c4762a1bSJed Brown       suffix: 1
616c4762a1bSJed Brown       args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4
617c4762a1bSJed Brown 
618c4762a1bSJed Brown    test:
619c4762a1bSJed Brown       suffix: 2
620c4762a1bSJed Brown       args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4
621c4762a1bSJed Brown 
622c4762a1bSJed Brown    test:
623c4762a1bSJed Brown       suffix: 3
624c4762a1bSJed Brown       args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian
625c4762a1bSJed Brown 
626c4762a1bSJed Brown    test:
627c4762a1bSJed Brown      suffix: 4
628c4762a1bSJed Brown      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls
629c4762a1bSJed Brown 
630c4762a1bSJed Brown    test:
631c4762a1bSJed Brown      suffix: 5
632c4762a1bSJed Brown      args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm
633c4762a1bSJed Brown 
634c4762a1bSJed Brown    test:
635c4762a1bSJed Brown      suffix: 6
636c4762a1bSJed Brown      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
637c4762a1bSJed Brown 
638c4762a1bSJed Brown TEST*/
639