Actual source code: rosenbrock1f.F

  1: ! "$Id: rosenbrock1f.F 1.52 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
  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(); TaoSetApplication();
 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:       integer          size,rank    ! number of processes running
 43:       PetscScalar      zero
 44:       TaoTerminateReason reason



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

 51:       external FormFunctionGradient,FormHessian


 54:       zero = 0.0d0

 56: !  Initialize TAO and PETSc
 57:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
 58:       call TaoInitialize(PETSC_NULL_CHARACTER,info)

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

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



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

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

 84: !  Allocate storage space for Hessian;
 85:       call MatCreateSeqBDiag(PETSC_COMM_SELF,n,n,0,2,0,0,H,info)
 86:       call MatSetOption(H,MAT_SYMMETRIC,info)


 89: !  The TAO code begins here

 91: !  Create TAO solver
 92:       call TaoCreate(PETSC_COMM_SELF,TAO_NULL_CHARACTER,tao,info)
 93:       call TaoApplicationCreate(PETSC_COMM_SELF,taoapp,info)

 95: !  Set routines for function, gradient, and hessian evaluation

 97: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
 98:       call TaoAppSetObjectiveAndGradientRo(taoapp,                      &
 99:      &     FormFunctionGradient,TAO_NULL_OBJECT,info)
100:       call TaoAppSetHessianMat(taoapp,H,H,info)
101:       call TaoAppSetHessianRoutine(taoapp,FormHessian,TAO_NULL_OBJECT,   &
102:      &     info)

104: !  Optional: Set initial guess
105:       call VecSet(x, zero, info)
106:       call TaoAppSetInitialSolutionVec(taoapp, x, info)

108: !  Now that the PETSc application is all set, attach to tao context
109:       call TaoSetApplication(tao,taoapp,info)

111: !  Check for TAO command line options
112:       call TaoSetOptions(taoapp,tao,info)

114: !  SOLVE THE APPLICATION
115:       call TaoSolveApplication(taoapp,tao,info)

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

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

129: !  Free TAO data structures
130:       call TaoDestroy(tao,info)
131:       call TaoAppDestroy(taoapp,info)

133: !  Free PETSc data structures
134:       call VecDestroy(x,info)
135:       call MatDestroy(H,info)

137: !  Finalize TAO
138:       call TaoFinalize(info)
139:       call PetscFinalize(info)

141:       end


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

159: ! n,alpha defined in rosenbrock1f.h
160: #include "rosenbrock1f.h"

162:       TAO_APPLICATION  taoapp
163:       Vec              X,G
164:       PetscScalar      f
165:       integer          dummy, info


168:       PetscScalar      ff,t1,t2
169:       integer          i,nn

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

179:       info = 0
180:       nn = n/2
181:       ff = 0

183: !     Get pointers to vector data
184:       call VecGetArray(X,x_v,x_i,info)
185:       call VecGetArray(G,g_v,g_i,info)


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

197: !     Restore vectors
198:       call VecRestoreArray(X,x_v,x_i,info)
199:       call VecRestoreArray(G,g_v,g_i,info)

201:       f = ff
202:       call PetscLogFlops(nn*15,info)

204:       return
205:       end

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

227:       subroutine FormHessian(taoapp,X,H,PrecH,flag,dummy,info)
228:       implicit none

230: #include "rosenbrock1f.h"

232: !  Input/output variables:
233:       TAO_APPLICATION  taoapp
234:       Vec              X
235:       Mat              H, PrecH
236:       MatStructure     flag
237:       integer          dummy,info
238: 
239:       PetscScalar      v(0:1,0:1)
240:       PetscTruth assembled

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


252:       info = 0
253:       nn= n/2

255: !  Zero existing matrix entries
256:       call MatAssembled(H,assembled,info)
257:       if (assembled .eq. PETSC_TRUE) call MatZeroEntries(H,info)

259: !  Get a pointer to vector data

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

263: !  Compute Hessian entries

265:       do i=0,nn-1
266:          v(1,1) = 2.0*alpha
267:          v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) -                          &
268:      &                3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
269:          v(1,0) = -4.0*alpha*x_v(x_i+2*i)
270:          v(0,1) = v(1,0)
271:          ind(0) = 2*i
272:          ind(1) = 2*i + 1
273:          call MatSetValues(H,2,ind,2,ind,v,INSERT_VALUES,info)
274:       enddo

276: !  Restore vector

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

280: !  Assemble matrix

282:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
283:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)

285:       call PetscLogFlops(nn*9,info)

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

303:       flag = SAME_NONZERO_PATTERN

305:       return
306:       end