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