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 "rosenbrock1f.h"

 17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 18: !                   Variable declarations
 19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: !
 21: !  See additional variable declarations in the file rosenbrock1f.h

 23:       PetscErrorCode   ierr    ! used to check for functions returning nonzeros
 24:       Vec              x       ! solution vector
 25:       Mat              H       ! hessian matrix
 26:       Tao        tao     ! TAO_SOVER context
 27:       PetscBool       flg
 28:       PetscInt         i2,i1
 29:       PetscMPIInt     size
 30:       PetscReal      zero

 32: !  Note: Any user-defined Fortran routines (such as FormGradient)
 33: !  MUST be declared as external.

 35:       external FormFunctionGradient,FormHessian

 37:       zero = 0.0d0
 38:       i2 = 2
 39:       i1 = 1

 41: !  Initialize TAO and PETSc
 42:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 43:       if (ierr .ne. 0) then
 44:          print*,'Unable to initialize PETSc'
 45:          stop
 46:       endif

 48:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 49:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif

 51: !  Initialize problem parameters
 52:       n     = 2
 53:       alpha = 99.0d0

 55: ! Check for command line arguments to override defaults
 56:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,    &
 57:      &                        '-n',n,flg,ierr)
 58:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,   &
 59:      &                         '-alpha',alpha,flg,ierr)

 61: !  Allocate vectors for the solution and gradient
 62:       call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)

 64: !  Allocate storage space for Hessian;
 65:       call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,                   &
 66:      &     PETSC_NULL_INTEGER, H,ierr)

 68:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)

 70: !  The TAO code begins here

 72: !  Create TAO solver
 73:       call TaoCreate(PETSC_COMM_SELF,tao,ierr)
 74:       CHKERRA(ierr)
 75:       call TaoSetType(tao,TAOLMVM,ierr)
 76:       CHKERRA(ierr)

 78: !  Set routines for function, gradient, and hessian evaluation
 79:       call TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,                 &
 80:      &      FormFunctionGradient,0,ierr)
 81:       CHKERRA(ierr)
 82:       call TaoSetHessian(tao,H,H,FormHessian,                             &
 83:      &     0,ierr)
 84:       CHKERRA(ierr)

 86: !  Optional: Set initial guess
 87:       call VecSet(x, zero, ierr)
 88:       call TaoSetSolution(tao, x, ierr)
 89:       CHKERRA(ierr)

 91: !  Check for TAO command line options
 92:       call TaoSetFromOptions(tao,ierr)
 93:       CHKERRA(ierr)

 95: !  SOLVE THE APPLICATION
 96:       call TaoSolve(tao,ierr)

 98: !  TaoView() prints ierr about the TAO solver; the option
 99: !      -tao_view
100: !  can alternatively be used to activate this at runtime.
101: !      call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)

103: !  Free TAO data structures
104:       call TaoDestroy(tao,ierr)

106: !  Free PETSc data structures
107:       call VecDestroy(x,ierr)
108:       call MatDestroy(H,ierr)

110:       call PetscFinalize(ierr)
111:       end

113: ! --------------------------------------------------------------------
114: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
115: !
116: !  Input Parameters:
117: !  tao - the Tao context
118: !  X   - input vector
119: !  dummy - not used
120: !
121: !  Output Parameters:
122: !  G - vector containing the newly evaluated gradient
123: !  f - function value

125:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
126: #include "rosenbrock1f.h"

128:       Tao        tao
129:       Vec              X,G
130:       PetscReal        f
131:       PetscErrorCode   ierr
132:       PetscInt         dummy

134:       PetscReal        ff,t1,t2
135:       PetscInt         i,nn

137: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
138: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
139: ! will return an array of doubles referenced by x_array offset by x_index.
140: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
141: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
142:       PetscReal        g_v(0:1),x_v(0:1)
143:       PetscOffset      g_i,x_i

145:       0
146:       nn = n/2
147:       ff = 0

149: !     Get pointers to vector data
150:       call VecGetArrayRead(X,x_v,x_i,ierr)
151:       call VecGetArray(G,g_v,g_i,ierr)

153: !     Compute G(X)
154:       do i=0,nn-1
155:          t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
156:          t2 = 1.0 - x_v(x_i + 2*i)
157:          ff = ff + alpha*t1*t1 + t2*t2
158:          g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
159:          g_v(g_i + 2*i + 1) = 2.0*alpha*t1
160:       enddo

162: !     Restore vectors
163:       call VecRestoreArrayRead(X,x_v,x_i,ierr)
164:       call VecRestoreArray(G,g_v,g_i,ierr)

166:       f = ff
167:       call PetscLogFlops(15.0d0*nn,ierr)

169:       return
170:       end

172: !
173: ! ---------------------------------------------------------------------
174: !
175: !  FormHessian - Evaluates Hessian matrix.
176: !
177: !  Input Parameters:
178: !  tao     - the Tao context
179: !  X       - input vector
180: !  dummy   - optional user-defined context, as set by SNESSetHessian()
181: !            (not used here)
182: !
183: !  Output Parameters:
184: !  H      - Hessian matrix
185: !  PrecH  - optionally different preconditioning matrix (not used here)
186: !  flag   - flag indicating matrix structure
187: !  ierr   - error code
188: !
189: !  Note: Providing the Hessian may not be necessary.  Only some solvers
190: !  require this matrix.

192:       subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
193: #include "rosenbrock1f.h"

195: !  Input/output variables:
196:       Tao        tao
197:       Vec              X
198:       Mat              H, PrecH
199:       PetscErrorCode   ierr
200:       PetscInt         dummy

202:       PetscReal        v(0:1,0:1)
203:       PetscBool assembled

205: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
206: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
207: ! will return an array of doubles referenced by x_array offset by x_index.
208: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
209: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
210:       PetscReal        x_v(0:1)
211:       PetscOffset      x_i
212:       PetscInt         i,nn,ind(0:1),i2

214:       0
215:       nn= n/2
216:       i2 = 2

218: !  Zero existing matrix entries
219:       call MatAssembled(H,assembled,ierr)
220:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)

222: !  Get a pointer to vector data

224:       call VecGetArrayRead(X,x_v,x_i,ierr)

226: !  Compute Hessian entries

228:       do i=0,nn-1
229:          v(1,1) = 2.0*alpha
230:          v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) -                          &
231:      &                3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
232:          v(1,0) = -4.0*alpha*x_v(x_i+2*i)
233:          v(0,1) = v(1,0)
234:          ind(0) = 2*i
235:          ind(1) = 2*i + 1
236:          call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
237:       enddo

239: !  Restore vector

241:       call VecRestoreArrayRead(X,x_v,x_i,ierr)

243: !  Assemble matrix

245:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
246:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)

248:       call PetscLogFlops(9.0d0*nn,ierr)

250:       return
251:       end

253: !
254: !/*TEST
255: !
256: !   build:
257: !      requires: !complex
258: !
259: !   test:
260: !      args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
261: !      requires: !single
262: !
263: !TEST*/