Actual source code: rosenbrock1f.F
1: ! "$Id$";
2: !
3: ! Program usage: mpirun -np 1 rosenbrock1f [-help] [all TAO options]
4: !
5: ! Description: This example demonstrates use of the TAO package to solve an
6: ! unconstrained minimization problem on a single processor. We minimize the
7: ! extended Rosenbrock function:
8: ! sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2 )
9: !
10: ! The C version of this code is rosenbrock1.c
11: !
12: !/*T
13: ! Concepts: TAO - Solving an unconstrained minimization problem
14: ! Routines: TaoInitialize(); TaoFinalize(); TaoSetOptions();
15: ! Routines: TaoApplicationCreate();
16: ! Routines: TaoCreate(); TaoAppSetObjectiveAndGradientRoutine();
17: ! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
18: ! Routines: TaoAppSetInitialSolutionVec();
19: ! Routines: TaoSolveApplication(); TaoDestroy(); TaoAppDestroy();
20: ! Routines: TaoView(); TaoGetTerminationReason();
21: ! Processors: 1
22: !T*/
23: !
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: integer info ! used to check for functions returning nonzeros
37: Vec x ! solution vector
38: Mat H ! hessian matrix
39: TAO_SOLVER tao ! TAO_SOVER context
40: TAO_APPLICATION taoapp ! TAO application context
41: PetscTruth flg
42: PetscInt i2,i1
43: integer size,rank ! number of processes running
44: PetscScalar zero
45: TaoTerminateReason reason
47:
50: ! Note: Any user-defined Fortran routines (such as FormGradient)
51: ! MUST be declared as external.
53: external FormFunctionGradient,FormHessian
55: zero = 0.0d0
56: i2 = 2
57: i1 = 1
59: ! Initialize TAO and PETSc
60: call PetscInitialize(PETSC_NULL_CHARACTER,info)
61: call TaoInitialize(PETSC_NULL_CHARACTER,info)
63: call MPI_Comm_size(PETSC_COMM_WORLD,size,info)
64: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,info)
65: if (size .ne. 1) then
66: if (rank .eq. 0) then
67: write(6,*) 'This is a uniprocessor example only!'
68: endif
69: SETERRQ(1,' ',info)
70: endif
72: ! Initialize problem parameters
73: n = 2
74: alpha = 99.0d0
78: ! Check for command line arguments to override defaults
79: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg, &
80: & info)
81: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha', &
82: & alpha,flg,info)
84: ! Allocate vectors for the solution and gradient
85: call VecCreateSeq(PETSC_COMM_SELF,n,x,info)
87: ! Allocate storage space for Hessian;
88: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
89: & PETSC_NULL_INTEGER, H,info)
91: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,info)
94: ! The TAO code begins here
96: ! Create TAO solver
97: call TaoCreate(PETSC_COMM_SELF,TAO_NULL_CHARACTER,tao,info)
98: call TaoApplicationCreate(PETSC_COMM_SELF,taoapp,info)
99: ! Set routines for function, gradient, and hessian evaluation
101: ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
102: call TaoAppSetObjectiveAndGradientRo(taoapp, &
103: & FormFunctionGradient,TAO_NULL_OBJECT,info)
104: call TaoAppSetHessianMat(taoapp,H,H,info)
105: call TaoAppSetHessianRoutine(taoapp,FormHessian,TAO_NULL_OBJECT, &
106: & info)
108: ! Optional: Set initial guess
109: call VecSet(x, zero, info)
110: call TaoAppSetInitialSolutionVec(taoapp, x, info)
113: ! Check for TAO command line options
114: call TaoSetTolerances(tao,1.0d-3,1.0d-3,1.0d-3,1.0d-3,info)
115: call TaoSetOptions(taoapp,tao,info)
117: ! SOLVE THE APPLICATION
118: call TaoSolveApplication(taoapp,tao,info)
120: call TaoGetTerminationReason(tao, reason, info)
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 info about the TAO solver; the option
127: ! -tao_view
128: ! can alternatively be used to activate this at runtime.
129: ! call TaoView(tao,info)
130:
132: ! Free TAO data structures
133: call TaoDestroy(tao,info)
134: call TaoAppDestroy(taoapp,info)
136: ! Free PETSc data structures
137: call VecDestroy(x,info)
138: call MatDestroy(H,info)
140: ! Finalize TAO
141: call TaoFinalize(info)
142: call PetscFinalize(info)
144: end
147: ! --------------------------------------------------------------------
148: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
149: !
150: ! Input Parameters:
151: ! tao - the TAO_SOLVER context
152: ! X - input vector
153: ! dummy - not used
154: !
155: ! Output Parameters:
156: ! G - vector containing the newly evaluated gradient
157: ! f - function value
158:
159: subroutine FormFunctionGradient(taoapp, X, f, G, dummy, info)
160: implicit none
162: ! n,alpha defined in rosenbrock1f.h
163: #include "rosenbrock1f.h"
165: TAO_APPLICATION taoapp
166: Vec X,G
167: PetscScalar f
168: integer info
169: PetscInt dummy
172: PetscScalar ff,t1,t2
173: PetscInt i,nn
175: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
176: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
177: ! will return an array of doubles referenced by x_array offset by x_index.
178: ! i.e., to reference the kth element of X, use x_array(k + x_index).
179: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
180: PetscScalar g_v(0:1),x_v(0:1)
181: PetscOffset g_i,x_i
183: info = 0
184: nn = n/2
185: ff = 0
187: ! Get pointers to vector data
188: call VecGetArray(X,x_v,x_i,info)
189: call VecGetArray(G,g_v,g_i,info)
192: ! Compute G(X)
193: do i=0,nn-1
194: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
195: t2 = 1.0 - x_v(x_i + 2*i)
196: ff = ff + alpha*t1*t1 + t2*t2
197: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
198: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
199: enddo
201: ! Restore vectors
202: call VecRestoreArray(X,x_v,x_i,info)
203: call VecRestoreArray(G,g_v,g_i,info)
205: f = ff
206: call PetscLogFlops(nn*15,info)
208: return
209: end
211: !
212: ! ---------------------------------------------------------------------
213: !
214: ! FormHessian - Evaluates Hessian matrix.
215: !
216: ! Input Parameters:
217: ! tao - the TAO_SOLVER context
218: ! X - input vector
219: ! dummy - optional user-defined context, as set by SNESSetHessian()
220: ! (not used here)
221: !
222: ! Output Parameters:
223: ! H - Hessian matrix
224: ! PrecH - optionally different preconditioning matrix (not used here)
225: ! flag - flag indicating matrix structure
226: ! info - error code
227: !
228: ! Note: Providing the Hessian may not be necessary. Only some solvers
229: ! require this matrix.
231: subroutine FormHessian(taoapp,X,H,PrecH,flag,dummy,info)
232: implicit none
234: #include "rosenbrock1f.h"
236: ! Input/output variables:
237: TAO_APPLICATION taoapp
238: Vec X
239: Mat H, PrecH
240: MatStructure flag
241: integer info
242: PetscInt dummy
243:
244: PetscScalar v(0:1,0:1)
245: PetscTruth assembled
247: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
248: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
249: ! will return an array of doubles referenced by x_array offset by x_index.
250: ! i.e., to reference the kth element of X, use x_array(k + x_index).
251: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
252: PetscScalar x_v(0:1)
253: PetscOffset x_i
254: PetscInt i,nn,ind(0:1),i2
257: info = 0
258: nn= n/2
259: i2 = 2
261: ! Zero existing matrix entries
262: call MatAssembled(H,assembled,info)
263: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,info)
265: ! Get a pointer to vector data
267: call VecGetArray(X,x_v,x_i,info)
269: ! Compute Hessian entries
271: do i=0,nn-1
272: v(1,1) = 2.0*alpha
273: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
274: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
275: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
276: v(0,1) = v(1,0)
277: ind(0) = 2*i
278: ind(1) = 2*i + 1
279: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,info)
280: enddo
282: ! Restore vector
284: call VecRestoreArray(X,x_v,x_i,info)
286: ! Assemble matrix
288: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
289: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)
291: call PetscLogFlops(nn*9,info)
293: ! Set flag to indicate that the Hessian matrix retains an identical
294: ! nonzero structure throughout all nonlinear iterations (although the
295: ! values of the entries change). Thus, we can save some work in setting
296: ! up the preconditioner (e.g., no need to redo symbolic factorization for
297: ! ICC preconditioners).
298: ! - If the nonzero structure of the matrix is different during
299: ! successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
300: ! must be used instead. If you are unsure whether the matrix
301: ! structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
302: ! - Caution: If you specify SAME_NONZERO_PATTERN, the software
303: ! believes your assertion and does not check the structure
304: ! of the matrix. If you erroneously claim that the structure
305: ! is the same when it actually is not, the new preconditioner
306: ! will not function correctly. Thus, use this optimization
307: ! feature with caution!
309: flag = SAME_NONZERO_PATTERN
311: return
312: end