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:       Vec             x       ! solution vector
 27:       Mat             H       ! hessian matrix
 28:       Tao             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:       if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif

 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:       Tao              tao
120:       Vec              X,G
121:       PetscReal        f
122:       PetscErrorCode   ierr
123:       PetscInt         dummy

125:       PetscReal        ff,t1,t2
126:       PetscInt         i,nn

128: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
129: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
130: ! will return an array of doubles referenced by x_array offset by x_index.
131: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
132: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
133:       PetscReal        g_v(0:1),x_v(0:1)
134:       PetscOffset      g_i,x_i
135:       PetscReal        alpha
136:       PetscInt         n
137:       common /params/ alpha, n

139:       0
140:       nn = n/2
141:       ff = 0

143: !     Get pointers to vector data
144:       PetscCall(VecGetArrayRead(X,x_v,x_i,ierr))
145:       PetscCall(VecGetArray(G,g_v,g_i,ierr))

147: !     Compute G(X)
148:       do i=0,nn-1
149:          t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
150:          t2 = 1.0 - x_v(x_i + 2*i)
151:          ff = ff + alpha*t1*t1 + t2*t2
152:          g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
153:          g_v(g_i + 2*i + 1) = 2.0*alpha*t1
154:       enddo

156: !     Restore vectors
157:       PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr))
158:       PetscCall(VecRestoreArray(G,g_v,g_i,ierr))

160:       f = ff
161:       PetscCall(PetscLogFlops(15.0d0*nn,ierr))

163:       return
164:       end

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

186:       subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
187:       use petsctao
188:       implicit none

190: !  Input/output variables:
191:       Tao              tao
192:       Vec              X
193:       Mat              H, PrecH
194:       PetscErrorCode   ierr
195:       PetscInt         dummy

197:       PetscReal        v(0:1,0:1)
198:       PetscBool        assembled

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

212:       0
213:       nn= n/2
214:       i2 = 2

216: !  Zero existing matrix entries
217:       PetscCall(MatAssembled(H,assembled,ierr))
218:       if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr))

220: !  Get a pointer to vector data

222:       PetscCall(VecGetArrayRead(X,x_v,x_i,ierr))

224: !  Compute Hessian entries

226:       do i=0,nn-1
227:          v(1,1) = 2.0*alpha
228:          v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
229:          v(1,0) = -4.0*alpha*x_v(x_i+2*i)
230:          v(0,1) = v(1,0)
231:          ind(0) = 2*i
232:          ind(1) = 2*i + 1
233:          PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr))
234:       enddo

236: !  Restore vector

238:       PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr))

240: !  Assemble matrix

242:       PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
243:       PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))

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

247:       return
248:       end

250: !
251: !/*TEST
252: !
253: !   build:
254: !      requires: !complex
255: !
256: !   test:
257: !      args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
258: !      requires: !single
259: !
260: !TEST*/