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 end 157 158! 159! --------------------------------------------------------------------- 160! 161! FormHessian - Evaluates Hessian matrix. 162! 163! Input Parameters: 164! tao - the Tao context 165! X - input vector 166! dummy - optional user-defined context, as set by SNESSetHessian() 167! (not used here) 168! 169! Output Parameters: 170! H - Hessian matrix 171! PrecH - optionally different preconditioning matrix (not used here) 172! flag - flag indicating matrix structure 173! ierr - error code 174! 175! Note: Providing the Hessian may not be necessary. Only some solvers 176! require this matrix. 177 178 subroutine FormHessian(tao,X,H,PrecH,dummy,ierr) 179 use petsctao 180 implicit none 181 182! Input/output variables: 183 type(tTao) tao 184 type(tVec) X 185 type(tMat) H, PrecH 186 PetscErrorCode ierr 187 PetscInt dummy 188 189 PetscReal v(0:1,0:1) 190 PetscBool assembled 191 192! PETSc's VecGetArray acts differently in Fortran than it does in C. 193! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 194! will return an array of doubles referenced by x_array offset by x_index. 195! i.e., to reference the kth element of X, use x_array(k + x_index). 196! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 197 PetscReal, pointer :: x_v(:) 198 PetscInt i,nn,ind(0:1),i2 199 PetscReal alpha 200 PetscInt n 201 common /params/ alpha, n 202 203 ierr = 0 204 nn= n/2 205 i2 = 2 206 207! Zero existing matrix entries 208 PetscCall(MatAssembled(H,assembled,ierr)) 209 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr)) 210 211! Get a pointer to vector data 212 213 PetscCall(VecGetArrayReadF90(X,x_v,ierr)) 214 215! Compute Hessian entries 216 217 do i=0,nn-1 218 v(1,1) = 2.0*alpha 219 v(0,0) = -4.0*alpha*(x_v(1+2*i+1) - 3*x_v(1+2*i)*x_v(1+2*i))+2 220 v(1,0) = -4.0*alpha*x_v(1+2*i) 221 v(0,1) = v(1,0) 222 ind(0) = 2*i 223 ind(1) = 2*i + 1 224 PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)) 225 enddo 226 227! Restore vector 228 229 PetscCall(VecRestoreArrayReadF90(X,x_v,ierr)) 230 231! Assemble matrix 232 233 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)) 234 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)) 235 236 PetscCall(PetscLogFlops(9.0d0*nn,ierr)) 237 238 end 239 240! 241!/*TEST 242! 243! build: 244! requires: !complex 245! 246! test: 247! args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5 248! requires: !single 249! 250!TEST*/ 251