Actual source code: rosenbrock1f.F

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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: TaoCreate();
 13: !  Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
 14: !  Routines: TaoSetHessianRoutine();
 15: !  Routines: TaoSetInitialVector();
 16: !  Routines: TaoSetFromOptions();
 17: !  Routines: TaoSolve();
 18: !  Routines: TaoGetConvergedReason(); TaoDestroy();
 19: !  Processors: 1
 20: !T*/
 21: !

 23: ! ----------------------------------------------------------------------
 24: !
 25:       implicit none

 27: #include "rosenbrock1f.h"

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

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



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

 50:       external FormFunctionGradient,FormHessian

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

 56: !  Initialize TAO and PETSc
 57:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

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

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



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

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

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

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


 90: !  The TAO code begins here

 92: !  Create TAO solver
 93:       call TaoCreate(PETSC_COMM_SELF,tao,ierr)
 94:       CHKERRQ(ierr)
 95:       call TaoSetType(tao,TAOLMVM,ierr)
 96:       CHKERRQ(ierr)

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


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


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

117: !  SOLVE THE APPLICATION
118:       call TaoSolve(tao,ierr)

120:       call TaoGetConvergedReason(tao, reason, ierr)
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 ierr about the TAO solver; the option
127: !      -tao_view
128: !  can alternatively be used to activate this at runtime.
129: !      call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)


132: !  Free TAO data structures
133:       call TaoDestroy(tao,ierr)

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

139:       call PetscFinalize(ierr)
140:       end


143: ! --------------------------------------------------------------------
144: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
145: !
146: !  Input Parameters:
147: !  tao - the Tao context
148: !  X   - input vector
149: !  dummy - not used
150: !
151: !  Output Parameters:
152: !  G - vector containing the newly evaluated gradient
153: !  f - function value

155:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
156:       implicit none

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

161:       Tao        tao
162:       Vec              X,G
163:       PetscReal        f
164:       PetscErrorCode   ierr
165:       PetscInt         dummy


168:       PetscReal        ff,t1,t2
169:       PetscInt         i,nn

171: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
172: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
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:       PetscReal        g_v(0:1),x_v(0:1)
177:       PetscOffset      g_i,x_i

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

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


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,ierr)
199:       call VecRestoreArray(G,g_v,g_i,ierr)

201:       f = ff
202:       call PetscLogFlops(15.0d0*nn,ierr)

204:       return
205:       end

207: !
208: ! ---------------------------------------------------------------------
209: !
210: !  FormHessian - Evaluates Hessian matrix.
211: !
212: !  Input Parameters:
213: !  tao     - the Tao 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: !  ierr   - error code
223: !
224: !  Note: Providing the Hessian may not be necessary.  Only some solvers
225: !  require this matrix.

227:       subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
228:       implicit none

230: #include "rosenbrock1f.h"

232: !  Input/output variables:
233:       Tao        tao
234:       Vec              X
235:       Mat              H, PrecH
236:       PetscErrorCode   ierr
237:       PetscInt         dummy

239:       PetscReal        v(0:1,0:1)
240:       PetscBool assembled

242: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
243: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
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:       PetscReal        x_v(0:1)
248:       PetscOffset      x_i
249:       PetscInt         i,nn,ind(0:1),i2


252:       0
253:       nn= n/2
254:       i2 = 2

256: !  Zero existing matrix entries
257:       call MatAssembled(H,assembled,ierr)
258:       if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)

260: !  Get a pointer to vector data

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

264: !  Compute Hessian entries

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

277: !  Restore vector

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

281: !  Assemble matrix

283:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
284:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)

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

288:       return
289:       end