Actual source code: rosenbrock1f.F
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: TaoInitialize(); TaoFinalize();
13: ! Routines: TaoCreate();
14: ! Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
15: ! Routines: TaoSetHessianRoutine();
16: ! Routines: TaoSetInitialVector();
17: ! Routines: TaoSetFromOptions();
18: ! Routines: TaoSolve();
19: ! Routines: TaoGetTerminationReason(); TaoDestroy();
20: ! Processors: 1
21: !T*/
22: !
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: PetscErrorCode ierr ! used to check for functions returning nonzeros
37: Vec x ! solution vector
38: Mat H ! hessian matrix
39: TaoSolver tao ! TAO_SOVER context
40: PetscBool flg
41: PetscInt i2,i1
42: integer size,rank ! number of processes running
43: PetscReal zero
44: TaoSolverTerminationReason reason
45:
48: ! Note: Any user-defined Fortran routines (such as FormGradient)
49: ! MUST be declared as external.
51: external FormFunctionGradient,FormHessian
53: zero = 0.0d0
54: i2 = 2
55: i1 = 1
57: ! Initialize TAO and PETSc
58: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
59: call TaoInitialize(PETSC_NULL_CHARACTER,ierr)
61: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
62: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
63: if (size .ne. 1) then
64: if (rank .eq. 0) then
65: write(6,*) 'This is a uniprocessor example only!'
66: endif
67: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
68: endif
70: ! Initialize problem parameters
71: n = 2
72: alpha = 99.0d0
76: ! Check for command line arguments to override defaults
77: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg, &
78: & ierr)
79: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha', &
80: & alpha,flg,ierr)
82: ! Allocate vectors for the solution and gradient
83: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
85: ! Allocate storage space for Hessian;
86: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
87: & PETSC_NULL_INTEGER, H,ierr)
89: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
92: ! The TAO code begins here
94: ! Create TAO solver
95: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
96: CHKERRQ(ierr)
97: call TaoSetType(tao,"tao_lmvm",ierr)
98: CHKERRQ(ierr)
100: ! Set routines for function, gradient, and hessian evaluation
101: call TaoSetObjectiveAndGradientRoutine(tao, &
102: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
103: CHKERRQ(ierr)
104: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
105: & PETSC_NULL_OBJECT,ierr)
106: CHKERRQ(ierr)
109: ! Optional: Set initial guess
110: call VecSet(x, zero, ierr)
111: call TaoSetInitialVector(tao, x, ierr)
112: CHKERRQ(ierr)
115: ! Check for TAO command line options
116: call TaoSetFromOptions(tao,ierr)
117: CHKERRQ(ierr)
119: ! SOLVE THE APPLICATION
120: call TaoSolve(tao,ierr)
122: call TaoGetTerminationReason(tao, reason, ierr)
123: if (reason .le. 0) then
124: print *,'Try a different TAO method, adjust some parameters,'
125: print *,'or check the function evaluation routines.'
126: endif
128: ! TaoView() prints ierr about the TAO solver; the option
129: ! -tao_view
130: ! can alternatively be used to activate this at runtime.
131: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
132:
134: ! Free TAO data structures
135: call TaoDestroy(tao,ierr)
137: ! Free PETSc data structures
138: call VecDestroy(x,ierr)
139: call MatDestroy(H,ierr)
141: ! Finalize TAO
142: call TaoFinalize(ierr)
143: call PetscFinalize(ierr)
145: end
148: ! --------------------------------------------------------------------
149: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
150: !
151: ! Input Parameters:
152: ! tao - the TaoSolver context
153: ! X - input vector
154: ! dummy - not used
155: !
156: ! Output Parameters:
157: ! G - vector containing the newly evaluated gradient
158: ! f - function value
159:
160: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
161: implicit none
163: ! n,alpha defined in rosenbrock1f.h
164: #include "rosenbrock1f.h"
166: TaoSolver tao
167: Vec X,G
168: PetscReal f
169: PetscErrorCode ierr
170: PetscInt dummy
173: PetscReal ff,t1,t2
174: PetscInt i,nn
176: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
177: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
178: ! will return an array of doubles referenced by x_array offset by x_index.
179: ! i.e., to reference the kth element of X, use x_array(k + x_index).
180: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
181: PetscReal g_v(0:1),x_v(0:1)
182: PetscOffset g_i,x_i
184: 0
185: nn = n/2
186: ff = 0
188: ! Get pointers to vector data
189: call VecGetArray(X,x_v,x_i,ierr)
190: call VecGetArray(G,g_v,g_i,ierr)
193: ! Compute G(X)
194: do i=0,nn-1
195: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
196: t2 = 1.0 - x_v(x_i + 2*i)
197: ff = ff + alpha*t1*t1 + t2*t2
198: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
199: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
200: enddo
202: ! Restore vectors
203: call VecRestoreArray(X,x_v,x_i,ierr)
204: call VecRestoreArray(G,g_v,g_i,ierr)
206: f = ff
207: call PetscLogFlops(nn*15.0d0,ierr)
208:
209: return
210: end
212: !
213: ! ---------------------------------------------------------------------
214: !
215: ! FormHessian - Evaluates Hessian matrix.
216: !
217: ! Input Parameters:
218: ! tao - the TaoSolver context
219: ! X - input vector
220: ! dummy - optional user-defined context, as set by SNESSetHessian()
221: ! (not used here)
222: !
223: ! Output Parameters:
224: ! H - Hessian matrix
225: ! PrecH - optionally different preconditioning matrix (not used here)
226: ! flag - flag indicating matrix structure
227: ! ierr - error code
228: !
229: ! Note: Providing the Hessian may not be necessary. Only some solvers
230: ! require this matrix.
232: subroutine FormHessian(tao,X,H,PrecH,flag,dummy,ierr)
233: implicit none
235: #include "rosenbrock1f.h"
237: ! Input/output variables:
238: TaoSolver tao
239: Vec X
240: Mat H, PrecH
241: MatStructure flag
242: PetscErrorCode ierr
243: PetscInt dummy
244:
245: PetscReal v(0:1,0:1)
246: PetscBool assembled
248: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
249: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
250: ! will return an array of doubles referenced by x_array offset by x_index.
251: ! i.e., to reference the kth element of X, use x_array(k + x_index).
252: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
253: PetscReal x_v(0:1)
254: PetscOffset x_i
255: PetscInt i,nn,ind(0:1),i2
258: 0
259: nn= n/2
260: i2 = 2
262: ! Zero existing matrix entries
263: call MatAssembled(H,assembled,ierr)
264: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
266: ! Get a pointer to vector data
268: call VecGetArray(X,x_v,x_i,ierr)
270: ! Compute Hessian entries
272: do i=0,nn-1
273: v(1,1) = 2.0*alpha
274: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
275: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
276: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
277: v(0,1) = v(1,0)
278: ind(0) = 2*i
279: ind(1) = 2*i + 1
280: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
281: enddo
283: ! Restore vector
285: call VecRestoreArray(X,x_v,x_i,ierr)
287: ! Assemble matrix
289: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
290: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
292: call PetscLogFlops(nn*9.0d0,ierr)
294: ! Set flag to indicate that the Hessian matrix retains an identical
295: ! nonzero structure throughout all nonlinear iterations (although the
296: ! values of the entries change). Thus, we can save some work in setting
297: ! up the preconditioner (e.g., no need to redo symbolic factorization for
298: ! ICC preconditioners).
299: ! - If the nonzero structure of the matrix is different during
300: ! successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
301: ! must be used instead. If you are unsure whether the matrix
302: ! structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
303: ! - Caution: If you specify SAME_NONZERO_PATTERN, the software
304: ! believes your assertion and does not check the structure
305: ! of the matrix. If you erroneously claim that the structure
306: ! is the same when it actually is not, the new preconditioner
307: ! will not function correctly. Thus, use this optimization
308: ! feature with caution!
310: flag = SAME_NONZERO_PATTERN
312: return
313: end