1! Program usage: mpiexec -n 1 rosenbrock1f [-help] [all TAO options] 2! 3! Description: This example demonstrates use of the TAO package to solve an 4! unconstrained minimization problem on a single processor. We minimize the 5! extended Rosenbrock function: 6! sum_{i=0}^{n/2-1} (alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2) 7! 8! The C version of this code is rosenbrock1.c 9! 10 11! 12 13! ---------------------------------------------------------------------- 14! 15#include "rosenbrock1f.h" 16 17! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 18! Variable declarations 19! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 20! 21! See additional variable declarations in the file rosenbrock1f.h 22 23 PetscErrorCode ierr ! used to check for functions returning nonzeros 24 Vec x ! solution vector 25 Mat H ! hessian matrix 26 Tao tao ! TAO_SOVER context 27 PetscBool flg 28 PetscInt i2,i1 29 PetscMPIInt size 30 PetscReal zero 31 32! Note: Any user-defined Fortran routines (such as FormGradient) 33! MUST be declared as external. 34 35 external FormFunctionGradient,FormHessian 36 37 zero = 0.0d0 38 i2 = 2 39 i1 = 1 40 41! Initialize TAO and PETSc 42 PetscCallA(PetscInitialize(ierr)) 43 44 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)) 45 if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif 46 47! Initialize problem parameters 48 n = 2 49 alpha = 99.0d0 50 51! Check for command line arguments to override defaults 52 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr)) 53 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-alpha',alpha,flg,ierr)) 54 55! Allocate vectors for the solution and gradient 56 PetscCallA(VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)) 57 58! Allocate storage space for Hessian; 59 PetscCallA(MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,PETSC_NULL_INTEGER, H,ierr)) 60 61 PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)) 62 63! The TAO code begins here 64 65! Create TAO solver 66 PetscCallA(TaoCreate(PETSC_COMM_SELF,tao,ierr)) 67 PetscCallA(TaoSetType(tao,TAOLMVM,ierr)) 68 69! Set routines for function, gradient, and hessian evaluation 70 PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr)) 71 PetscCallA(TaoSetHessian(tao,H,H,FormHessian,0,ierr)) 72 73! Optional: Set initial guess 74 PetscCallA(VecSet(x, zero, ierr)) 75 PetscCallA(TaoSetSolution(tao, x, ierr)) 76 77! Check for TAO command line options 78 PetscCallA(TaoSetFromOptions(tao,ierr)) 79 80! SOLVE THE APPLICATION 81 PetscCallA(TaoSolve(tao,ierr)) 82 83! TaoView() prints ierr about the TAO solver; the option 84! -tao_view 85! can alternatively be used to activate this at runtime. 86! PetscCallA(TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)) 87 88! Free TAO data structures 89 PetscCallA(TaoDestroy(tao,ierr)) 90 91! Free PETSc data structures 92 PetscCallA(VecDestroy(x,ierr)) 93 PetscCallA(MatDestroy(H,ierr)) 94 95 PetscCallA(PetscFinalize(ierr)) 96 end 97 98! -------------------------------------------------------------------- 99! FormFunctionGradient - Evaluates the function f(X) and gradient G(X) 100! 101! Input Parameters: 102! tao - the Tao context 103! X - input vector 104! dummy - not used 105! 106! Output Parameters: 107! G - vector containing the newly evaluated gradient 108! f - function value 109 110 subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr) 111#include "rosenbrock1f.h" 112 113 Tao tao 114 Vec X,G 115 PetscReal f 116 PetscErrorCode ierr 117 PetscInt dummy 118 119 PetscReal ff,t1,t2 120 PetscInt i,nn 121 122! PETSc's VecGetArray acts differently in Fortran than it does in C. 123! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 124! will return an array of doubles referenced by x_array offset by x_index. 125! i.e., to reference the kth element of X, use x_array(k + x_index). 126! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 127 PetscReal g_v(0:1),x_v(0:1) 128 PetscOffset g_i,x_i 129 130 ierr = 0 131 nn = n/2 132 ff = 0 133 134! Get pointers to vector data 135 PetscCall(VecGetArrayRead(X,x_v,x_i,ierr)) 136 PetscCall(VecGetArray(G,g_v,g_i,ierr)) 137 138! Compute G(X) 139 do i=0,nn-1 140 t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i) 141 t2 = 1.0 - x_v(x_i + 2*i) 142 ff = ff + alpha*t1*t1 + t2*t2 143 g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2 144 g_v(g_i + 2*i + 1) = 2.0*alpha*t1 145 enddo 146 147! Restore vectors 148 PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr)) 149 PetscCall(VecRestoreArray(G,g_v,g_i,ierr)) 150 151 f = ff 152 PetscCall(PetscLogFlops(15.0d0*nn,ierr)) 153 154 return 155 end 156 157! 158! --------------------------------------------------------------------- 159! 160! FormHessian - Evaluates Hessian matrix. 161! 162! Input Parameters: 163! tao - the Tao context 164! X - input vector 165! dummy - optional user-defined context, as set by SNESSetHessian() 166! (not used here) 167! 168! Output Parameters: 169! H - Hessian matrix 170! PrecH - optionally different preconditioning matrix (not used here) 171! flag - flag indicating matrix structure 172! ierr - error code 173! 174! Note: Providing the Hessian may not be necessary. Only some solvers 175! require this matrix. 176 177 subroutine FormHessian(tao,X,H,PrecH,dummy,ierr) 178#include "rosenbrock1f.h" 179 180! Input/output variables: 181 Tao tao 182 Vec X 183 Mat H, PrecH 184 PetscErrorCode ierr 185 PetscInt dummy 186 187 PetscReal v(0:1,0:1) 188 PetscBool assembled 189 190! PETSc's VecGetArray acts differently in Fortran than it does in C. 191! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 192! will return an array of doubles referenced by x_array offset by x_index. 193! i.e., to reference the kth element of X, use x_array(k + x_index). 194! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 195 PetscReal x_v(0:1) 196 PetscOffset x_i 197 PetscInt i,nn,ind(0:1),i2 198 199 ierr = 0 200 nn= n/2 201 i2 = 2 202 203! Zero existing matrix entries 204 PetscCall(MatAssembled(H,assembled,ierr)) 205 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr)) 206 207! Get a pointer to vector data 208 209 PetscCall(VecGetArrayRead(X,x_v,x_i,ierr)) 210 211! Compute Hessian entries 212 213 do i=0,nn-1 214 v(1,1) = 2.0*alpha 215 v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2 216 v(1,0) = -4.0*alpha*x_v(x_i+2*i) 217 v(0,1) = v(1,0) 218 ind(0) = 2*i 219 ind(1) = 2*i + 1 220 PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)) 221 enddo 222 223! Restore vector 224 225 PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr)) 226 227! Assemble matrix 228 229 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)) 230 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)) 231 232 PetscCall(PetscLogFlops(9.0d0*nn,ierr)) 233 234 return 235 end 236 237! 238!/*TEST 239! 240! build: 241! requires: !complex 242! 243! test: 244! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5 245! requires: !single 246! 247!TEST*/ 248