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