Actual source code: rosenbrock1f.F
1: ! "$Id: rosenbrock1f.F 1.52 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
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(); TaoSetApplication();
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: integer size,rank ! number of processes running
43: PetscScalar zero
44: TaoTerminateReason reason
48: ! Note: Any user-defined Fortran routines (such as FormGradient)
49: ! MUST be declared as external.
51: external FormFunctionGradient,FormHessian
54: zero = 0.0d0
56: ! Initialize TAO and PETSc
57: call PetscInitialize(PETSC_NULL_CHARACTER,info)
58: call TaoInitialize(PETSC_NULL_CHARACTER,info)
60: call MPI_Comm_size(PETSC_COMM_WORLD,size,info)
61: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,info)
62: if (size .ne. 1) then
63: if (rank .eq. 0) then
64: write(6,*) 'This is a uniprocessor example only!'
65: endif
66: SETERRQ(1,' ',info)
67: endif
69: ! Initialize problem parameters
70: n = 2
71: alpha = 99.0d0
75: ! Check for command line arguments to override defaults
76: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg, &
77: & info)
78: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha', &
79: & alpha,flg,info)
81: ! Allocate vectors for the solution and gradient
82: call VecCreateSeq(PETSC_COMM_SELF,n,x,info)
84: ! Allocate storage space for Hessian;
85: call MatCreateSeqBDiag(PETSC_COMM_SELF,n,n,0,2,0,0,H,info)
86: call MatSetOption(H,MAT_SYMMETRIC,info)
89: ! The TAO code begins here
91: ! Create TAO solver
92: call TaoCreate(PETSC_COMM_SELF,TAO_NULL_CHARACTER,tao,info)
93: call TaoApplicationCreate(PETSC_COMM_SELF,taoapp,info)
95: ! Set routines for function, gradient, and hessian evaluation
97: ! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
98: call TaoAppSetObjectiveAndGradientRo(taoapp, &
99: & FormFunctionGradient,TAO_NULL_OBJECT,info)
100: call TaoAppSetHessianMat(taoapp,H,H,info)
101: call TaoAppSetHessianRoutine(taoapp,FormHessian,TAO_NULL_OBJECT, &
102: & info)
104: ! Optional: Set initial guess
105: call VecSet(x, zero, info)
106: call TaoAppSetInitialSolutionVec(taoapp, x, info)
108: ! Now that the PETSc application is all set, attach to tao context
109: call TaoSetApplication(tao,taoapp,info)
111: ! Check for TAO command line options
112: call TaoSetOptions(taoapp,tao,info)
114: ! SOLVE THE APPLICATION
115: call TaoSolveApplication(taoapp,tao,info)
117: call TaoGetTerminationReason(tao, reason, info)
118: if (reason .le. 0) then
119: print *,'Try a different TAO method, adjust some parameters,'
120: print *,'or check the function evaluation routines.'
121: endif
123: ! TaoView() prints info about the TAO solver; the option
124: ! -tao_view
125: ! can alternatively be used to activate this at runtime.
126: ! call TaoView(tao,info)
127:
129: ! Free TAO data structures
130: call TaoDestroy(tao,info)
131: call TaoAppDestroy(taoapp,info)
133: ! Free PETSc data structures
134: call VecDestroy(x,info)
135: call MatDestroy(H,info)
137: ! Finalize TAO
138: call TaoFinalize(info)
139: call PetscFinalize(info)
141: end
144: ! --------------------------------------------------------------------
145: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
146: !
147: ! Input Parameters:
148: ! tao - the TAO_SOLVER context
149: ! X - input vector
150: ! dummy - not used
151: !
152: ! Output Parameters:
153: ! G - vector containing the newly evaluated gradient
154: ! f - function value
155:
156: subroutine FormFunctionGradient(taoapp, X, f, G, dummy, info)
157: implicit none
159: ! n,alpha defined in rosenbrock1f.h
160: #include "rosenbrock1f.h"
162: TAO_APPLICATION taoapp
163: Vec X,G
164: PetscScalar f
165: integer dummy, info
168: PetscScalar ff,t1,t2
169: integer i,nn
171: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
172: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
173: ! will return an array of doubles referenced by x_array offset by x_index.
174: ! i.e., to reference the kth element of X, use x_array(k + x_index).
175: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
176: PetscScalar g_v(0:1),x_v(0:1)
177: PetscOffset g_i,x_i
179: info = 0
180: nn = n/2
181: ff = 0
183: ! Get pointers to vector data
184: call VecGetArray(X,x_v,x_i,info)
185: call VecGetArray(G,g_v,g_i,info)
188: ! Compute G(X)
189: do i=0,nn-1
190: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
191: t2 = 1.0 - x_v(x_i + 2*i)
192: ff = ff + alpha*t1*t1 + t2*t2
193: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
194: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
195: enddo
197: ! Restore vectors
198: call VecRestoreArray(X,x_v,x_i,info)
199: call VecRestoreArray(G,g_v,g_i,info)
201: f = ff
202: call PetscLogFlops(nn*15,info)
204: return
205: end
207: !
208: ! ---------------------------------------------------------------------
209: !
210: ! FormHessian - Evaluates Hessian matrix.
211: !
212: ! Input Parameters:
213: ! tao - the TAO_SOLVER context
214: ! X - input vector
215: ! dummy - optional user-defined context, as set by SNESSetHessian()
216: ! (not used here)
217: !
218: ! Output Parameters:
219: ! H - Hessian matrix
220: ! PrecH - optionally different preconditioning matrix (not used here)
221: ! flag - flag indicating matrix structure
222: ! info - error code
223: !
224: ! Note: Providing the Hessian may not be necessary. Only some solvers
225: ! require this matrix.
227: subroutine FormHessian(taoapp,X,H,PrecH,flag,dummy,info)
228: implicit none
230: #include "rosenbrock1f.h"
232: ! Input/output variables:
233: TAO_APPLICATION taoapp
234: Vec X
235: Mat H, PrecH
236: MatStructure flag
237: integer dummy,info
238:
239: PetscScalar v(0:1,0:1)
240: PetscTruth assembled
242: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
243: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
244: ! will return an array of doubles referenced by x_array offset by x_index.
245: ! i.e., to reference the kth element of X, use x_array(k + x_index).
246: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
247: PetscScalar x_v(0:1)
248: PetscOffset x_i
249: integer i,nn,ind(0:1)
252: info = 0
253: nn= n/2
255: ! Zero existing matrix entries
256: call MatAssembled(H,assembled,info)
257: if (assembled .eq. PETSC_TRUE) call MatZeroEntries(H,info)
259: ! Get a pointer to vector data
261: call VecGetArray(X,x_v,x_i,info)
263: ! Compute Hessian entries
265: do i=0,nn-1
266: v(1,1) = 2.0*alpha
267: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
268: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
269: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
270: v(0,1) = v(1,0)
271: ind(0) = 2*i
272: ind(1) = 2*i + 1
273: call MatSetValues(H,2,ind,2,ind,v,INSERT_VALUES,info)
274: enddo
276: ! Restore vector
278: call VecRestoreArray(X,x_v,x_i,info)
280: ! Assemble matrix
282: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
283: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)
285: call PetscLogFlops(nn*9,info)
287: ! Set flag to indicate that the Hessian matrix retains an identical
288: ! nonzero structure throughout all nonlinear iterations (although the
289: ! values of the entries change). Thus, we can save some work in setting
290: ! up the preconditioner (e.g., no need to redo symbolic factorization for
291: ! ICC preconditioners).
292: ! - If the nonzero structure of the matrix is different during
293: ! successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
294: ! must be used instead. If you are unsure whether the matrix
295: ! structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
296: ! - Caution: If you specify SAME_NONZERO_PATTERN, the software
297: ! believes your assertion and does not check the structure
298: ! of the matrix. If you erroneously claim that the structure
299: ! is the same when it actually is not, the new preconditioner
300: ! will not function correctly. Thus, use this optimization
301: ! feature with caution!
303: flag = SAME_NONZERO_PATTERN
305: return
306: end