Actual source code: rosenbrock1f.F90
petsc-3.9.4 2018-09-11
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: !
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: TaoDestroy();
19: ! Processors: 1
20: !T*/
23: !
25: ! ----------------------------------------------------------------------
26: !
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
42: PetscReal zero
44: ! Note: Any user-defined Fortran routines (such as FormGradient)
45: ! MUST be declared as external.
47: external FormFunctionGradient,FormHessian
49: zero = 0.0d0
50: i2 = 2
51: i1 = 1
53: ! Initialize TAO and PETSc
54: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
55: if (ierr .ne. 0) then
56: print*,'Unable to initialize PETSc'
57: stop
58: endif
60: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
61: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,1,'This is a uniprocessor example only'); endif
63: ! Initialize problem parameters
64: n = 2
65: alpha = 99.0d0
68: ! Check for command line arguments to override defaults
69: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
70: & '-n',n,flg,ierr)
71: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
72: & '-alpha',alpha,flg,ierr)
74: ! Allocate vectors for the solution and gradient
75: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
77: ! Allocate storage space for Hessian;
78: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
79: & PETSC_NULL_INTEGER, H,ierr)
81: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
84: ! The TAO code begins here
86: ! Create TAO solver
87: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
88: CHKERRA(ierr)
89: call TaoSetType(tao,TAOLMVM,ierr)
90: CHKERRA(ierr)
92: ! Set routines for function, gradient, and hessian evaluation
93: call TaoSetObjectiveAndGradientRoutine(tao, &
94: & FormFunctionGradient,0,ierr)
95: CHKERRA(ierr)
96: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
97: & 0,ierr)
98: CHKERRA(ierr)
101: ! Optional: Set initial guess
102: call VecSet(x, zero, ierr)
103: call TaoSetInitialVector(tao, x, ierr)
104: CHKERRA(ierr)
107: ! Check for TAO command line options
108: call TaoSetFromOptions(tao,ierr)
109: CHKERRA(ierr)
111: ! SOLVE THE APPLICATION
112: call TaoSolve(tao,ierr)
114: ! TaoView() prints ierr about the TAO solver; the option
115: ! -tao_view
116: ! can alternatively be used to activate this at runtime.
117: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
120: ! Free TAO data structures
121: call TaoDestroy(tao,ierr)
123: ! Free PETSc data structures
124: call VecDestroy(x,ierr)
125: call MatDestroy(H,ierr)
127: call PetscFinalize(ierr)
128: end
131: ! --------------------------------------------------------------------
132: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
133: !
134: ! Input Parameters:
135: ! tao - the Tao context
136: ! X - input vector
137: ! dummy - not used
138: !
139: ! Output Parameters:
140: ! G - vector containing the newly evaluated gradient
141: ! f - function value
143: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
144: #include "rosenbrock1f.h"
146: Tao tao
147: Vec X,G
148: PetscReal f
149: PetscErrorCode ierr
150: PetscInt dummy
153: PetscReal ff,t1,t2
154: PetscInt i,nn
156: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
157: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
158: ! will return an array of doubles referenced by x_array offset by x_index.
159: ! i.e., to reference the kth element of X, use x_array(k + x_index).
160: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
161: PetscReal g_v(0:1),x_v(0:1)
162: PetscOffset g_i,x_i
164: 0
165: nn = n/2
166: ff = 0
168: ! Get pointers to vector data
169: call VecGetArrayRead(X,x_v,x_i,ierr)
170: call VecGetArray(G,g_v,g_i,ierr)
173: ! Compute G(X)
174: do i=0,nn-1
175: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
176: t2 = 1.0 - x_v(x_i + 2*i)
177: ff = ff + alpha*t1*t1 + t2*t2
178: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
179: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
180: enddo
182: ! Restore vectors
183: call VecRestoreArrayRead(X,x_v,x_i,ierr)
184: call VecRestoreArray(G,g_v,g_i,ierr)
186: f = ff
187: call PetscLogFlops(15.0d0*nn,ierr)
189: return
190: end
192: !
193: ! ---------------------------------------------------------------------
194: !
195: ! FormHessian - Evaluates Hessian matrix.
196: !
197: ! Input Parameters:
198: ! tao - the Tao context
199: ! X - input vector
200: ! dummy - optional user-defined context, as set by SNESSetHessian()
201: ! (not used here)
202: !
203: ! Output Parameters:
204: ! H - Hessian matrix
205: ! PrecH - optionally different preconditioning matrix (not used here)
206: ! flag - flag indicating matrix structure
207: ! ierr - error code
208: !
209: ! Note: Providing the Hessian may not be necessary. Only some solvers
210: ! require this matrix.
212: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
213: #include "rosenbrock1f.h"
215: ! Input/output variables:
216: Tao tao
217: Vec X
218: Mat H, PrecH
219: PetscErrorCode ierr
220: PetscInt dummy
222: PetscReal v(0:1,0:1)
223: PetscBool assembled
225: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
226: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
227: ! will return an array of doubles referenced by x_array offset by x_index.
228: ! i.e., to reference the kth element of X, use x_array(k + x_index).
229: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
230: PetscReal x_v(0:1)
231: PetscOffset x_i
232: PetscInt i,nn,ind(0:1),i2
235: 0
236: nn= n/2
237: i2 = 2
239: ! Zero existing matrix entries
240: call MatAssembled(H,assembled,ierr)
241: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
243: ! Get a pointer to vector data
245: call VecGetArrayRead(X,x_v,x_i,ierr)
247: ! Compute Hessian entries
249: do i=0,nn-1
250: v(1,1) = 2.0*alpha
251: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
252: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
253: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
254: v(0,1) = v(1,0)
255: ind(0) = 2*i
256: ind(1) = 2*i + 1
257: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
258: enddo
260: ! Restore vector
262: call VecRestoreArrayRead(X,x_v,x_i,ierr)
264: ! Assemble matrix
266: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
267: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
269: call PetscLogFlops(9.0d0*nn,ierr)
271: return
272: end
278: !
279: !/*TEST
280: !
281: ! build:
282: ! requires: !complex
283: !
284: ! test:
285: ! args: -tao_smonitor -tao_type ntr
286: ! requires: !single
287: !
288: !TEST*/