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 "petsc/finclude/petsctao.h"
 16:       use petsctao
 17:       implicit none

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

 25:       PetscErrorCode  ierr    ! used to check for functions returning nonzeros
 26:       type(tVec)      x       ! solution vector
 27:       type(tMat)      H       ! hessian matrix
 28:       type(tTao)      tao     ! TAO_SOVER context
 29:       PetscBool       flg
 30:       PetscInt        i2,i1
 31:       PetscMPIInt     size
 32:       PetscReal       zero
 33:       PetscReal       alpha
 34:       PetscInt        n
 35:       common /params/ alpha, n

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

 40:       external FormFunctionGradient,FormHessian

 42:       zero = 0.0d0
 43:       i2 = 2
 44:       i1 = 1

 46: !  Initialize TAO and PETSc
 47:       PetscCallA(PetscInitialize(ierr))

 49:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 50:       PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')

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

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

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

 63: !  Allocate storage space for Hessian;
 64:       PetscCallA(MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,PETSC_NULL_INTEGER, H,ierr))

 66:       PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))

 68: !  The TAO code begins here

 70: !  Create TAO solver
 71:       PetscCallA(TaoCreate(PETSC_COMM_SELF,tao,ierr))
 72:       PetscCallA(TaoSetType(tao,TAOLMVM,ierr))

 74: !  Set routines for function, gradient, and hessian evaluation
 75:       PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
 76:       PetscCallA(TaoSetHessian(tao,H,H,FormHessian,0,ierr))

 78: !  Optional: Set initial guess
 79:       PetscCallA(VecSet(x, zero, ierr))
 80:       PetscCallA(TaoSetSolution(tao, x, ierr))

 82: !  Check for TAO command line options
 83:       PetscCallA(TaoSetFromOptions(tao,ierr))

 85: !  SOLVE THE APPLICATION
 86:       PetscCallA(TaoSolve(tao,ierr))

 88: !  TaoView() prints ierr about the TAO solver; the option
 89: !      -tao_view
 90: !  can alternatively be used to activate this at runtime.
 91: !      PetscCallA(TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr))

 93: !  Free TAO data structures
 94:       PetscCallA(TaoDestroy(tao,ierr))

 96: !  Free PETSc data structures
 97:       PetscCallA(VecDestroy(x,ierr))
 98:       PetscCallA(MatDestroy(H,ierr))

100:       PetscCallA(PetscFinalize(ierr))
101:       end

103: ! --------------------------------------------------------------------
104: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
105: !
106: !  Input Parameters:
107: !  tao - the Tao context
108: !  X   - input vector
109: !  dummy - not used
110: !
111: !  Output Parameters:
112: !  G - vector containing the newly evaluated gradient
113: !  f - function value

115:       subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
116:       use petsctao
117:       implicit none

119:       type(tTao)       tao
120:       type(tVec)       X,G
121:       PetscReal        f
122:       PetscErrorCode   ierr
123:       PetscInt         dummy

125:       PetscReal        ff,t1,t2
126:       PetscInt         i,nn
127:       PetscReal, pointer :: g_v(:),x_v(:)
128:       PetscReal        alpha
129:       PetscInt         n
130:       common /params/ alpha, n

132:       ierr = 0
133:       nn = n/2
134:       ff = 0

136: !     Get pointers to vector data
137:       PetscCall(VecGetArrayReadF90(X,x_v,ierr))
138:       PetscCall(VecGetArrayF90(G,g_v,ierr))

140: !     Compute G(X)
141:       do i=0,nn-1
142:          t1 = x_v(1+2*i+1) - x_v(1+2*i)*x_v(1+2*i)
143:          t2 = 1.0 - x_v(1 + 2*i)
144:          ff = ff + alpha*t1*t1 + t2*t2
145:          g_v(1 + 2*i) = -4*alpha*t1*x_v(1 + 2*i) - 2.0*t2
146:          g_v(1 + 2*i + 1) = 2.0*alpha*t1
147:       enddo

149: !     Restore vectors
150:       PetscCall(VecRestoreArrayReadF90(X,x_v,ierr))
151:       PetscCall(VecRestoreArrayF90(G,g_v,ierr))

153:       f = ff
154:       PetscCall(PetscLogFlops(15.0d0*nn,ierr))

156:       return
157:       end

159: !
160: ! ---------------------------------------------------------------------
161: !
162: !  FormHessian - Evaluates Hessian matrix.
163: !
164: !  Input Parameters:
165: !  tao     - the Tao context
166: !  X       - input vector
167: !  dummy   - optional user-defined context, as set by SNESSetHessian()
168: !            (not used here)
169: !
170: !  Output Parameters:
171: !  H      - Hessian matrix
172: !  PrecH  - optionally different preconditioning matrix (not used here)
173: !  flag   - flag indicating matrix structure
174: !  ierr   - error code
175: !
176: !  Note: Providing the Hessian may not be necessary.  Only some solvers
177: !  require this matrix.

179:       subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
180:       use petsctao
181:       implicit none

183: !  Input/output variables:
184:       type(tTao)       tao
185:       type(tVec)       X
186:       type(tMat)       H, PrecH
187:       PetscErrorCode   ierr
188:       PetscInt         dummy

190:       PetscReal        v(0:1,0:1)
191:       PetscBool        assembled

193: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
194: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
195: ! will return an array of doubles referenced by x_array offset by x_index.
196: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
197: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
198:       PetscReal, pointer :: x_v(:)
199:       PetscInt         i,nn,ind(0:1),i2
200:       PetscReal        alpha
201:       PetscInt         n
202:       common /params/ alpha, n

204:       ierr = 0
205:       nn= n/2
206:       i2 = 2

208: !  Zero existing matrix entries
209:       PetscCall(MatAssembled(H,assembled,ierr))
210:       if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr))

212: !  Get a pointer to vector data

214:       PetscCall(VecGetArrayReadF90(X,x_v,ierr))

216: !  Compute Hessian entries

218:       do i=0,nn-1
219:          v(1,1) = 2.0*alpha
220:          v(0,0) = -4.0*alpha*(x_v(1+2*i+1) - 3*x_v(1+2*i)*x_v(1+2*i))+2
221:          v(1,0) = -4.0*alpha*x_v(1+2*i)
222:          v(0,1) = v(1,0)
223:          ind(0) = 2*i
224:          ind(1) = 2*i + 1
225:          PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr))
226:       enddo

228: !  Restore vector

230:       PetscCall(VecRestoreArrayReadF90(X,x_v,ierr))

232: !  Assemble matrix

234:       PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
235:       PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))

237:       PetscCall(PetscLogFlops(9.0d0*nn,ierr))

239:       return
240:       end

242: !
243: !/*TEST
244: !
245: !   build:
246: !      requires: !complex
247: !
248: !   test:
249: !      args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
250: !      requires: !single
251: !
252: !TEST*/