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