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