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 type(tVec) x ! solution vector 27 type(tMat) H ! hessian matrix 28 type(tTao) 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 PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only') 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 type(tTao) tao 120 type(tVec) X,G 121 PetscReal f 122 PetscErrorCode ierr 123 PetscInt dummy 124 125 PetscReal ff,t1,t2 126 PetscInt i,nn 127 PetscReal, pointer :: g_v(:),x_v(:) 128 PetscReal alpha 129 PetscInt n 130 common /params/ alpha, n 131 132 ierr = 0 133 nn = n/2 134 ff = 0 135 136! Get pointers to vector data 137 PetscCall(VecGetArrayReadF90(X,x_v,ierr)) 138 PetscCall(VecGetArrayF90(G,g_v,ierr)) 139 140! Compute G(X) 141 do i=0,nn-1 142 t1 = x_v(1+2*i+1) - x_v(1+2*i)*x_v(1+2*i) 143 t2 = 1.0 - x_v(1 + 2*i) 144 ff = ff + alpha*t1*t1 + t2*t2 145 g_v(1 + 2*i) = -4*alpha*t1*x_v(1 + 2*i) - 2.0*t2 146 g_v(1 + 2*i + 1) = 2.0*alpha*t1 147 enddo 148 149! Restore vectors 150 PetscCall(VecRestoreArrayReadF90(X,x_v,ierr)) 151 PetscCall(VecRestoreArrayF90(G,g_v,ierr)) 152 153 f = ff 154 PetscCall(PetscLogFlops(15.0d0*nn,ierr)) 155 156 return 157 end 158 159! 160! --------------------------------------------------------------------- 161! 162! FormHessian - Evaluates Hessian matrix. 163! 164! Input Parameters: 165! tao - the Tao context 166! X - input vector 167! dummy - optional user-defined context, as set by SNESSetHessian() 168! (not used here) 169! 170! Output Parameters: 171! H - Hessian matrix 172! PrecH - optionally different preconditioning matrix (not used here) 173! flag - flag indicating matrix structure 174! ierr - error code 175! 176! Note: Providing the Hessian may not be necessary. Only some solvers 177! require this matrix. 178 179 subroutine FormHessian(tao,X,H,PrecH,dummy,ierr) 180 use petsctao 181 implicit none 182 183! Input/output variables: 184 type(tTao) tao 185 type(tVec) X 186 type(tMat) H, PrecH 187 PetscErrorCode ierr 188 PetscInt dummy 189 190 PetscReal v(0:1,0:1) 191 PetscBool assembled 192 193! PETSc's VecGetArray acts differently in Fortran than it does in C. 194! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 195! will return an array of doubles referenced by x_array offset by x_index. 196! i.e., to reference the kth element of X, use x_array(k + x_index). 197! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 198 PetscReal, pointer :: x_v(:) 199 PetscInt i,nn,ind(0:1),i2 200 PetscReal alpha 201 PetscInt n 202 common /params/ alpha, n 203 204 ierr = 0 205 nn= n/2 206 i2 = 2 207 208! Zero existing matrix entries 209 PetscCall(MatAssembled(H,assembled,ierr)) 210 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr)) 211 212! Get a pointer to vector data 213 214 PetscCall(VecGetArrayReadF90(X,x_v,ierr)) 215 216! Compute Hessian entries 217 218 do i=0,nn-1 219 v(1,1) = 2.0*alpha 220 v(0,0) = -4.0*alpha*(x_v(1+2*i+1) - 3*x_v(1+2*i)*x_v(1+2*i))+2 221 v(1,0) = -4.0*alpha*x_v(1+2*i) 222 v(0,1) = v(1,0) 223 ind(0) = 2*i 224 ind(1) = 2*i + 1 225 PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)) 226 enddo 227 228! Restore vector 229 230 PetscCall(VecRestoreArrayReadF90(X,x_v,ierr)) 231 232! Assemble matrix 233 234 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)) 235 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)) 236 237 PetscCall(PetscLogFlops(9.0d0*nn,ierr)) 238 239 return 240 end 241 242! 243!/*TEST 244! 245! build: 246! requires: !complex 247! 248! test: 249! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5 250! requires: !single 251! 252!TEST*/ 253