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: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
101: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
102: ! will return an array of doubles referenced by x_array offset by x_index.
103: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
104: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
105:     PetscReal, pointer :: x_v(:)
106:     PetscInt i, nn, ind(0:1)

108: !   Zero existing matrix entries
109:     PetscCall(MatAssembled(H, assembled, ierr))
110:     if (assembled) PetscCall(MatZeroEntries(H, ierr))

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

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

127: !  Restore vector

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

131: !  Assemble matrix

133:     PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
134:     PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))

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

138:   end
139: end module rosenbrock1fmodule

141: program rosenbrock1f
142:   use petsctao
143:   use rosenbrock1fmodule
144:   implicit none

146: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: !                   Variable declarations
148: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: !
150: !  See additional variable declarations in the file rosenbrock1f.h

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

159: !  Initialize TAO and PETSc
160:   PetscCallA(PetscInitialize(ierr))

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

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

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

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

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

180: !  The TAO code begins here

182: ! Create TAO solver
183:   PetscCallA(TaoCreate(PETSC_COMM_SELF, ta, ierr))
184:   PetscCallA(TaoSetType(ta, TAOLMVM, ierr))

186: ! Set routines for function, gradient, and hessian evaluation
187:   PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
188:   PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))

190: ! Optional: Set initial guess
191:   PetscCallA(VecSet(x, 0.0_PETSC_REAL_KIND, ierr))
192:   PetscCallA(TaoSetSolution(ta, x, ierr))

194: ! Check for TAO command line options
195:   PetscCallA(TaoSetFromOptions(ta, ierr))

197: ! SOLVE THE APPLICATION
198:   PetscCallA(TaoSolve(ta, ierr))

200: !  TaoView() prints ierr about the TAO solver; the option
201: !      -tao_view
202: !  can alternatively be used to activate this at runtime.
203: !      PetscCallA(TaoView(ta,PETSC_VIEWER_STDOUT_SELF,ierr))

205: ! Free TAO data structures
206:   PetscCallA(TaoDestroy(ta, ierr))

208: ! Free PETSc data structures
209:   PetscCallA(VecDestroy(x, ierr))
210:   PetscCallA(MatDestroy(H, ierr))

212:   PetscCallA(PetscFinalize(ierr))
213: end program rosenbrock1f
214: !
215: !/*TEST
216: !
217: !   build:
218: !      requires: !complex
219: !
220: !   test:
221: !      args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5
222: !      requires: !single
223: !
224: !TEST*/