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) ta ! 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 == 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_ARRAY, 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, ta, ierr)) 72 PetscCallA(TaoSetType(ta, TAOLMVM, ierr)) 73 74! Set routines for function, gradient, and hessian evaluation 75 PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr)) 76 PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr)) 77 78! Optional: Set initial guess 79 PetscCallA(VecSet(x, zero, ierr)) 80 PetscCallA(TaoSetSolution(ta, x, ierr)) 81 82! Check for TAO command line options 83 PetscCallA(TaoSetFromOptions(ta, ierr)) 84 85! SOLVE THE APPLICATION 86 PetscCallA(TaoSolve(ta, 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(ta,PETSC_VIEWER_STDOUT_SELF,ierr)) 92 93! Free TAO data structures 94 PetscCallA(TaoDestroy(ta, 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(ta, X, f, G, dummy, ierr) 116 use petsctao 117 implicit none 118 119 type(tTao) ta 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(VecGetArrayRead(X, x_v, ierr)) 138 PetscCall(VecGetArray(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 end do 148 149! Restore vectors 150 PetscCall(VecRestoreArrayRead(X, x_v, ierr)) 151 PetscCall(VecRestoreArray(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 matrix used to compute the preconditioner (not used here) 172! ierr - error code 173! 174! Note: Providing the Hessian may not be necessary. Only some solvers 175! require this matrix. 176 177 subroutine FormHessian(ta, X, H, PrecH, dummy, ierr) 178 use petsctao 179 implicit none 180 181! Input/output variables: 182 type(tTao) ta 183 type(tVec) X 184 type(tMat) H, PrecH 185 PetscErrorCode ierr 186 PetscInt dummy 187 188 PetscReal v(0:1, 0:1) 189 PetscBool assembled 190 191! PETSc's VecGetArray acts differently in Fortran than it does in C. 192! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 193! will return an array of doubles referenced by x_array offset by x_index. 194! i.e., to reference the kth element of X, use x_array(k + x_index). 195! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 196 PetscReal, pointer :: x_v(:) 197 PetscInt i, nn, ind(0:1), i2 198 PetscReal alpha 199 PetscInt n 200 common/params/alpha, n 201 202 ierr = 0 203 nn = n/2 204 i2 = 2 205 206! Zero existing matrix entries 207 PetscCall(MatAssembled(H, assembled, ierr)) 208 if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H, ierr)) 209 210! Get a pointer to vector data 211 212 PetscCall(VecGetArrayRead(X, x_v, ierr)) 213 214! Compute Hessian entries 215 216 do i = 0, nn - 1 217 v(1, 1) = 2.0*alpha 218 v(0, 0) = -4.0*alpha*(x_v(1 + 2*i + 1) - 3*x_v(1 + 2*i)*x_v(1 + 2*i)) + 2 219 v(1, 0) = -4.0*alpha*x_v(1 + 2*i) 220 v(0, 1) = v(1, 0) 221 ind(0) = 2*i 222 ind(1) = 2*i + 1 223 PetscCall(MatSetValues(H, i2, ind, i2, ind, reshape(v, [i2*i2]), INSERT_VALUES, ierr)) 224 end do 225 226! Restore vector 227 228 PetscCall(VecRestoreArrayRead(X, x_v, ierr)) 229 230! Assemble matrix 231 232 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr)) 233 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr)) 234 235 PetscCall(PetscLogFlops(9.0d0*nn, ierr)) 236 237 end 238 239! 240!/*TEST 241! 242! build: 243! requires: !complex 244! 245! test: 246! args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5 247! requires: !single 248! 249!TEST*/ 250