Actual source code: rosenbrock1f.F

  1: ! "$Id$";
  2: !
  3: !  Program usage: mpirun -np 1 rosenbrock1f [-help] [all TAO options] 
  4: !
  5: !  Description:  This example demonstrates use of the TAO package to solve an
  6: !  unconstrained minimization problem on a single processor.  We minimize the 
  7: !  extended Rosenbrock function: 
  8: !       sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2 )
  9: !
 10: !  The C version of this code is rosenbrock1.c
 11: !
 12: !/*T
 13: !  Concepts: TAO - Solving an unconstrained minimization problem
 14: !  Routines: TaoInitialize(); TaoFinalize(); TaoSetOptions();
 15: !  Routines: TaoApplicationCreate();
 16: !  Routines: TaoCreate(); TaoAppSetObjectiveAndGradientRoutine(); 
 17: !  Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
 18: !  Routines: TaoAppSetInitialSolutionVec(); 
 19: !  Routines: TaoSolveApplication(); TaoDestroy(); TaoAppDestroy();
 20: !  Routines: TaoView(); TaoGetTerminationReason();                              
 21: !  Processors: 1
 22: !T*/ 
 23: !
 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:       integer          info    ! used to check for functions returning nonzeros
 37:       Vec              x       ! solution vector
 38:       Mat              H       ! hessian matrix
 39:       TAO_SOLVER       tao     ! TAO_SOVER context
 40:       TAO_APPLICATION  taoapp   ! TAO application context
 41:       PetscTruth       flg     
 42:       PetscInt         i2,i1
 43:       integer          size,rank    ! number of processes running
 44:       PetscScalar      zero
 45:       TaoTerminateReason reason

 47:       


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

 53:       external FormFunctionGradient,FormHessian

 55:       zero = 0.0d0
 56:       i2 = 2
 57:       i1 = 1

 59: !  Initialize TAO and PETSc
 60:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
 61:       call TaoInitialize(PETSC_NULL_CHARACTER,info)

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

 72: !  Initialize problem parameters
 73:       n     = 2
 74:       alpha = 99.0d0



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

 84: !  Allocate vectors for the solution and gradient
 85:       call VecCreateSeq(PETSC_COMM_SELF,n,x,info)

 87: !  Allocate storage space for Hessian; 
 88:       call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,                   &
 89:      &     PETSC_NULL_INTEGER, H,info)

 91:       call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)


 94: !  The TAO code begins here 

 96: !  Create TAO solver
 97:       call TaoCreate(PETSC_COMM_SELF,TAO_NULL_CHARACTER,tao,info)
 98:       call TaoApplicationCreate(PETSC_COMM_SELF,taoapp,info)
 99: !  Set routines for function, gradient, and hessian evaluation 

101: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
102:       call TaoAppSetObjectiveAndGradientRo(taoapp,                      &
103:      &     FormFunctionGradient,TAO_NULL_OBJECT,info)
104:       call TaoAppSetHessianMat(taoapp,H,H,info)
105:       call TaoAppSetHessianRoutine(taoapp,FormHessian,TAO_NULL_OBJECT,   &
106:      &     info)

108: !  Optional: Set initial guess
109:       call VecSet(x, zero, info)
110:       call TaoAppSetInitialSolutionVec(taoapp, x, info)


113: !  Check for TAO command line options
114:       call TaoSetTolerances(tao,1.0d-3,1.0d-3,1.0d-3,1.0d-3,info)
115:       call TaoSetOptions(taoapp,tao,info)

117: !  SOLVE THE APPLICATION
118:       call TaoSolveApplication(taoapp,tao,info)

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

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

132: !  Free TAO data structures
133:       call TaoDestroy(tao,info)
134:       call TaoAppDestroy(taoapp,info)

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

140: !  Finalize TAO 
141:       call TaoFinalize(info)
142:       call PetscFinalize(info)

144:       end


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

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

165:       TAO_APPLICATION  taoapp
166:       Vec              X,G
167:       PetscScalar      f
168:       integer          info
169:       PetscInt         dummy


172:       PetscScalar      ff,t1,t2
173:       PetscInt         i,nn

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

183:       info = 0
184:       nn = n/2
185:       ff = 0

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


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

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

205:       f = ff
206:       call PetscLogFlops(nn*15,info)

208:       return
209:       end

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

231:       subroutine FormHessian(taoapp,X,H,PrecH,flag,dummy,info)
232:       implicit none

234: #include "rosenbrock1f.h"

236: !  Input/output variables:
237:       TAO_APPLICATION  taoapp
238:       Vec              X
239:       Mat              H, PrecH
240:       MatStructure     flag
241:       integer          info
242:       PetscInt         dummy
243:       
244:       PetscScalar      v(0:1,0:1)
245:       PetscTruth assembled

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


257:       info = 0
258:       nn= n/2
259:       i2 = 2

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

265: !  Get a pointer to vector data

267:       call VecGetArray(X,x_v,x_i,info)

269: !  Compute Hessian entries

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

282: !  Restore vector

284:       call VecRestoreArray(X,x_v,x_i,info)

286: !  Assemble matrix

288:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
289:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)

291:       call PetscLogFlops(nn*9,info)

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

309:       flag = SAME_NONZERO_PATTERN

311:       return
312:       end