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