Actual source code: rosenbrock1f.F90

  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: !

 11: !

 13: ! ----------------------------------------------------------------------
 14: !
 15: #include "petsc/finclude/petsctao.h"
 16: module rosenbrock1fmodule
 17:   use petsctao

 19:   implicit none
 20:   PetscReal alpha
 21:   PetscInt n

 23: contains
 24: ! --------------------------------------------------------------------
 25: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
 26: !
 27: !  Input Parameters:
 28: !  tao - the Tao context
 29: !  X   - input vector
 30: !  dummy - not used
 31: !
 32: !  Output Parameters:
 33: !  G - vector containing the newly evaluated gradient
 34: !  f - function value

 36:   subroutine FormFunctionGradient(ta, X, f, G, dummy, ierr)
 37:     type(tTao) ta
 38:     type(tVec) X, G
 39:     PetscReal, intent(out) :: f
 40:     PetscErrorCode, intent(out) :: ierr
 41:     PetscInt dummy

 43:     PetscReal t1, t2
 44:     PetscInt i, nn
 45:     PetscReal, pointer :: g_v(:), x_v(:)

 47: !   Get pointers to vector data
 48:     PetscCall(VecGetArrayRead(X, x_v, ierr))
 49:     PetscCall(VecGetArray(G, g_v, ierr))

 51:     nn = n/2
 52:     f = 0.0_PETSC_REAL_KIND
 53: !   Compute G(X)
 54:     do i = 0, nn - 1
 55:       t1 = x_v(1 + 2*i + 1) - x_v(1 + 2*i)*x_v(1 + 2*i)
 56:       t2 = 1.0 - x_v(1 + 2*i)
 57:       f = f + alpha*t1**2 + t2**2
 58:       g_v(1 + 2*i) = -4*alpha*t1*x_v(1 + 2*i) - 2.0*t2
 59:       g_v(1 + 2*i + 1) = 2.0*alpha*t1
 60:     end do

 62: !   Restore vectors
 63:     PetscCall(VecRestoreArrayRead(X, x_v, ierr))
 64:     PetscCall(VecRestoreArray(G, g_v, ierr))

 66:     PetscCall(PetscLogFlops(15.0d0*nn, ierr))

 68:   end

 70: !
 71: ! ---------------------------------------------------------------------
 72: !
 73: !  FormHessian - Evaluates Hessian matrix.
 74: !
 75: !  Input Parameters:
 76: !  tao     - the Tao context
 77: !  X       - input vector
 78: !  dummy   - optional user-defined context, as set by SNESSetHessian()
 79: !            (not used here)
 80: !
 81: !  Output Parameters:
 82: !  H      - Hessian matrix
 83: !  PrecH  - optionally different matrix used to compute the preconditioner (not used here)
 84: !  ierr   - error code
 85: !
 86: !  Note: Providing the Hessian may not be necessary.  Only some solvers
 87: !  require this matrix.

 89:   subroutine FormHessian(ta, X, H, PrecH, dummy, ierr)
 90: !  Input/output variables:
 91:     type(tTao) ta
 92:     type(tVec) X
 93:     type(tMat) H, PrecH
 94:     PetscErrorCode, intent(out) :: ierr
 95:     PetscInt dummy

 97:     PetscReal v(0:1, 0:1)
 98:     PetscBool assembled

100: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
101:     PetscReal, pointer :: x_v(:)
102:     PetscInt i, nn, ind(0:1)

104: !   Zero existing matrix entries
105:     PetscCall(MatAssembled(H, assembled, ierr))
106:     if (assembled) PetscCall(MatZeroEntries(H, ierr))

108: !   Get a pointer to vector data
109:     PetscCall(VecGetArrayRead(X, x_v, ierr))

111: !   Compute Hessian entries
112:     nn = n/2
113:     do i = 0, nn - 1
114:       v(1, 1) = 2.0*alpha
115:       v(0, 0) = -4.0*alpha*(x_v(1 + 2*i + 1) - 3*x_v(1 + 2*i)*x_v(1 + 2*i)) + 2
116:       v(1, 0) = -4.0*alpha*x_v(1 + 2*i)
117:       v(0, 1) = v(1, 0)
118:       ind(0) = 2*i
119:       ind(1) = 2*i + 1
120:       PetscCall(MatSetValues(H, 2_PETSC_INT_KIND, ind, 2_PETSC_INT_KIND, ind, reshape(v, [2_PETSC_INT_KIND**2]), INSERT_VALUES, ierr))
121:     end do

123: !  Restore vector

125:     PetscCall(VecRestoreArrayRead(X, x_v, ierr))

127: !  Assemble matrix

129:     PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
130:     PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))

132:     PetscCall(PetscLogFlops(9.0d0*nn, ierr))

134:   end
135: end module rosenbrock1fmodule

137: program rosenbrock1f
138:   use petsctao
139:   use rosenbrock1fmodule
140:   implicit none

142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: !                   Variable declarations
144: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: !
146: !  See additional variable declarations in the file rosenbrock1f.h

148:   PetscErrorCode ierr    ! used to check for functions returning nonzeros
149:   type(tVec) x       ! solution vector
150:   type(tMat) H       ! hessian matrix
151:   type(tTao) ta     ! TAO_SOVER context
152:   PetscBool flg
153:   PetscMPIInt size

155: !  Initialize TAO and PETSc
156:   PetscCallA(PetscInitialize(ierr))

158:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
159:   PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')

161: ! Initialize problem parameters and
162: ! check for command line arguments to override defaults
163:   n = 2
164:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
165:   alpha = 99.0d0
166:   PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-alpha', alpha, flg, ierr))

168: ! Allocate vectors for the solution and gradient
169:   PetscCallA(VecCreateSeq(PETSC_COMM_SELF, n, x, ierr))

171: ! Allocate storage space for Hessian
172:   PetscCallA(MatCreateSeqBAIJ(PETSC_COMM_SELF, 2_PETSC_INT_KIND, n, n, 1_PETSC_INT_KIND, PETSC_NULL_INTEGER_ARRAY, H, ierr))

174:   PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))

176: !  The TAO code begins here

178: ! Create TAO solver
179:   PetscCallA(TaoCreate(PETSC_COMM_SELF, ta, ierr))
180:   PetscCallA(TaoSetType(ta, TAOLMVM, ierr))

182: ! Set routines for function, gradient, and hessian evaluation
183:   PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
184:   PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))

186: ! Optional: Set initial guess
187:   PetscCallA(VecSet(x, 0.0_PETSC_REAL_KIND, ierr))
188:   PetscCallA(TaoSetSolution(ta, x, ierr))

190: ! Check for TAO command line options
191:   PetscCallA(TaoSetFromOptions(ta, ierr))

193: ! SOLVE THE APPLICATION
194:   PetscCallA(TaoSolve(ta, ierr))

196: !  TaoView() prints ierr about the TAO solver; the option
197: !      -tao_view
198: !  can alternatively be used to activate this at runtime.
199: !      PetscCallA(TaoView(ta,PETSC_VIEWER_STDOUT_SELF,ierr))

201: ! Free TAO data structures
202:   PetscCallA(TaoDestroy(ta, ierr))

204: ! Free PETSc data structures
205:   PetscCallA(VecDestroy(x, ierr))
206:   PetscCallA(MatDestroy(H, ierr))

208:   PetscCallA(PetscFinalize(ierr))
209: end program rosenbrock1f
210: !
211: !/*TEST
212: !
213: !   build:
214: !      requires: !complex
215: !
216: !   test:
217: !      args: -tao_monitor -tao_type ntr -tao_gatol 1.e-5
218: !      requires: !single
219: !
220: !TEST*/