Actual source code: rosenbrock1f.F

  1: !  Program usage: mpirun -np 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: !/*T
 11: !  Concepts: TAO - Solving an unconstrained minimization problem
 12: !  Routines: TaoInitialize(); TaoFinalize();
 13: !  Routines: TaoCreate();
 14: !  Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
 15: !  Routines: TaoSetHessianRoutine();
 16: !  Routines: TaoSetInitialVector();
 17: !  Routines: TaoSetFromOptions();
 18: !  Routines: TaoSolve();
 19: !  Routines: TaoGetTerminationReason(); TaoDestroy(); 
 20: !  Processors: 1
 21: !T*/ 
 22: !

 24: ! ---------------------------------------------------------------------- 
 25: !
 26:       implicit none

 28: #include "rosenbrock1f.h"

 30: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 31: !                   Variable declarations
 32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 33: !
 34: !  See additional variable declarations in the file rosenbrock1f.h

 36:       PetscErrorCode   ierr    ! used to check for functions returning nonzeros
 37:       Vec              x       ! solution vector
 38:       Mat              H       ! hessian matrix
 39:       TaoSolver        tao     ! TAO_SOVER context
 40:       PetscBool       flg     
 41:       PetscInt         i2,i1
 42:       integer          size,rank    ! number of processes running
 43:       PetscReal      zero
 44:       TaoSolverTerminationReason reason
 45:       


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

 51:       external FormFunctionGradient,FormHessian

 53:       zero = 0.0d0
 54:       i2 = 2
 55:       i1 = 1

 57: !  Initialize TAO and PETSc
 58:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 59:       call TaoInitialize(PETSC_NULL_CHARACTER,ierr)

 61:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 62:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 63:       if (size .ne. 1) then
 64:          if (rank .eq. 0) then
 65:             write(6,*) 'This is a uniprocessor example only!'
 66:          endif
 67:          SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
 68:       endif

 70: !  Initialize problem parameters
 71:       n     = 2
 72:       alpha = 99.0d0



 76: ! Check for command line arguments to override defaults        
 77:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,           &
 78:      &                        ierr)
 79:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha',            &
 80:      &                           alpha,flg,ierr)

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

 85: !  Allocate storage space for Hessian; 
 86:       call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,                   &
 87:      &     PETSC_NULL_INTEGER, H,ierr)

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


 92: !  The TAO code begins here 

 94: !  Create TAO solver
 95:       call TaoCreate(PETSC_COMM_SELF,tao,ierr)
 96:       CHKERRQ(ierr)
 97:       call TaoSetType(tao,"tao_lmvm",ierr)
 98:       CHKERRQ(ierr)

100: !  Set routines for function, gradient, and hessian evaluation 
101:       call TaoSetObjectiveAndGradientRoutine(tao,                       &
102:      &      FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
103:       CHKERRQ(ierr)
104:       call TaoSetHessianRoutine(tao,H,H,FormHessian,                    &
105:      &     PETSC_NULL_OBJECT,ierr)
106:       CHKERRQ(ierr)


109: !  Optional: Set initial guess
110:       call VecSet(x, zero, ierr)
111:       call TaoSetInitialVector(tao, x, ierr)
112:       CHKERRQ(ierr)


115: !  Check for TAO command line options
116:       call TaoSetFromOptions(tao,ierr)
117:       CHKERRQ(ierr)

119: !  SOLVE THE APPLICATION
120:       call TaoSolve(tao,ierr)

122:       call TaoGetTerminationReason(tao, reason, ierr)
123:       if (reason .le. 0) then
124:          print *,'Try a different TAO method, adjust some parameters,'
125:          print *,'or check the function evaluation routines.'
126:       endif

128: !  TaoView() prints ierr about the TAO solver; the option
129: !      -tao_view 
130: !  can alternatively be used to activate this at runtime.
131: !      call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
132:       

134: !  Free TAO data structures
135:       call TaoDestroy(tao,ierr)

137: !  Free PETSc data structures 
138:       call VecDestroy(x,ierr)
139:       call MatDestroy(H,ierr)

141: !  Finalize TAO 
142:       call TaoFinalize(ierr)
143:       call PetscFinalize(ierr)

145:       end


148: ! --------------------------------------------------------------------
149: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
150: !
151: !  Input Parameters:
152: !  tao - the TaoSolver context
153: !  X   - input vector
154: !  dummy - not used
155: !
156: !  Output Parameters:
157: !  G - vector containing the newly evaluated gradient
158: !  f - function value
159:       
160:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
161:       implicit none

163: ! n,alpha defined in rosenbrock1f.h
164: #include "rosenbrock1f.h"

166:       TaoSolver        tao
167:       Vec              X,G
168:       PetscReal        f
169:       PetscErrorCode   ierr
170:       PetscInt         dummy


173:       PetscReal        ff,t1,t2
174:       PetscInt         i,nn

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

184:       0
185:       nn = n/2
186:       ff = 0

188: !     Get pointers to vector data
189:       call VecGetArray(X,x_v,x_i,ierr)
190:       call VecGetArray(G,g_v,g_i,ierr)


193: !     Compute G(X)
194:       do i=0,nn-1
195:          t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
196:          t2 = 1.0 - x_v(x_i + 2*i)
197:          ff = ff + alpha*t1*t1 + t2*t2
198:          g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
199:          g_v(g_i + 2*i + 1) = 2.0*alpha*t1
200:       enddo

202: !     Restore vectors
203:       call VecRestoreArray(X,x_v,x_i,ierr)
204:       call VecRestoreArray(G,g_v,g_i,ierr)

206:       f = ff
207:       call PetscLogFlops(nn*15.0d0,ierr)
208:       
209:       return
210:       end

212: !  
213: ! ---------------------------------------------------------------------
214: !
215: !  FormHessian - Evaluates Hessian matrix.
216: !
217: !  Input Parameters:
218: !  tao     - the TaoSolver context
219: !  X       - input vector
220: !  dummy   - optional user-defined context, as set by SNESSetHessian()
221: !            (not used here)
222: !
223: !  Output Parameters:
224: !  H      - Hessian matrix
225: !  PrecH  - optionally different preconditioning matrix (not used here)
226: !  flag   - flag indicating matrix structure
227: !  ierr   - error code
228: !
229: !  Note: Providing the Hessian may not be necessary.  Only some solvers
230: !  require this matrix.

232:       subroutine FormHessian(tao,X,H,PrecH,flag,dummy,ierr)
233:       implicit none

235: #include "rosenbrock1f.h"

237: !  Input/output variables:
238:       TaoSolver        tao
239:       Vec              X
240:       Mat              H, PrecH
241:       MatStructure     flag
242:       PetscErrorCode   ierr
243:       PetscInt         dummy
244:       
245:       PetscReal        v(0:1,0:1)
246:       PetscBool assembled

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


258:       0
259:       nn= n/2
260:       i2 = 2

262: !  Zero existing matrix entries
263:       call MatAssembled(H,assembled,ierr)
264:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)

266: !  Get a pointer to vector data

268:       call VecGetArray(X,x_v,x_i,ierr)

270: !  Compute Hessian entries

272:       do i=0,nn-1
273:          v(1,1) = 2.0*alpha
274:          v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) -                          &
275:      &                3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
276:          v(1,0) = -4.0*alpha*x_v(x_i+2*i)
277:          v(0,1) = v(1,0)
278:          ind(0) = 2*i
279:          ind(1) = 2*i + 1
280:          call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
281:       enddo

283: !  Restore vector

285:       call VecRestoreArray(X,x_v,x_i,ierr)

287: !  Assemble matrix

289:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
290:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)

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

294: !  Set flag to indicate that the Hessian matrix retains an identical
295: !  nonzero structure throughout all nonlinear iterations (although the
296: !  values of the entries change). Thus, we can save some work in setting
297: !  up the preconditioner (e.g., no need to redo symbolic factorization for
298: !  ICC preconditioners).
299: !   - If the nonzero structure of the matrix is different during
300: !     successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
301: !     must be used instead.  If you are unsure whether the matrix
302: !     structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
303: !   - Caution:  If you specify SAME_NONZERO_PATTERN, the software
304: !     believes your assertion and does not check the structure
305: !     of the matrix.  If you erroneously claim that the structure
306: !     is the same when it actually is not, the new preconditioner
307: !     will not function correctly.  Thus, use this optimization
308: !     feature with caution!

310:       flag = SAME_NONZERO_PATTERN

312:       return
313:       end