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: !
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*/
22: !
24: ! ----------------------------------------------------------------------
25: !
26: #include "rosenbrock1f.h"
28: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: ! Variable declarations
30: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: !
32: ! See additional variable declarations in the file rosenbrock1f.h
34: PetscErrorCode ierr ! used to check for functions returning nonzeros
35: Vec x ! solution vector
36: Mat H ! hessian matrix
37: Tao tao ! TAO_SOVER context
38: PetscBool flg
39: PetscInt i2,i1
40: PetscMPIInt size
41: PetscReal zero
43: ! Note: Any user-defined Fortran routines (such as FormGradient)
44: ! MUST be declared as external.
46: external FormFunctionGradient,FormHessian
48: zero = 0.0d0
49: i2 = 2
50: i1 = 1
52: ! Initialize TAO and PETSc
53: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
54: if (ierr .ne. 0) then
55: print*,'Unable to initialize PETSc'
56: stop
57: endif
59: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
60: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
62: ! Initialize problem parameters
63: n = 2
64: alpha = 99.0d0
66: ! Check for command line arguments to override defaults
67: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
68: & '-n',n,flg,ierr)
69: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
70: & '-alpha',alpha,flg,ierr)
72: ! Allocate vectors for the solution and gradient
73: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
75: ! Allocate storage space for Hessian;
76: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
77: & PETSC_NULL_INTEGER, H,ierr)
79: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
81: ! The TAO code begins here
83: ! Create TAO solver
84: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
85: CHKERRA(ierr)
86: call TaoSetType(tao,TAOLMVM,ierr)
87: CHKERRA(ierr)
89: ! Set routines for function, gradient, and hessian evaluation
90: call TaoSetObjectiveAndGradientRoutine(tao, &
91: & FormFunctionGradient,0,ierr)
92: CHKERRA(ierr)
93: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
94: & 0,ierr)
95: CHKERRA(ierr)
97: ! Optional: Set initial guess
98: call VecSet(x, zero, ierr)
99: call TaoSetInitialVector(tao, x, ierr)
100: CHKERRA(ierr)
102: ! Check for TAO command line options
103: call TaoSetFromOptions(tao,ierr)
104: CHKERRA(ierr)
106: ! SOLVE THE APPLICATION
107: call TaoSolve(tao,ierr)
109: ! TaoView() prints ierr about the TAO solver; the option
110: ! -tao_view
111: ! can alternatively be used to activate this at runtime.
112: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
114: ! Free TAO data structures
115: call TaoDestroy(tao,ierr)
117: ! Free PETSc data structures
118: call VecDestroy(x,ierr)
119: call MatDestroy(H,ierr)
121: call PetscFinalize(ierr)
122: end
124: ! --------------------------------------------------------------------
125: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
126: !
127: ! Input Parameters:
128: ! tao - the Tao context
129: ! X - input vector
130: ! dummy - not used
131: !
132: ! Output Parameters:
133: ! G - vector containing the newly evaluated gradient
134: ! f - function value
136: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
137: #include "rosenbrock1f.h"
139: Tao tao
140: Vec X,G
141: PetscReal f
142: PetscErrorCode ierr
143: PetscInt dummy
145: PetscReal ff,t1,t2
146: PetscInt i,nn
148: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
149: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
150: ! will return an array of doubles referenced by x_array offset by x_index.
151: ! i.e., to reference the kth element of X, use x_array(k + x_index).
152: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
153: PetscReal g_v(0:1),x_v(0:1)
154: PetscOffset g_i,x_i
156: 0
157: nn = n/2
158: ff = 0
160: ! Get pointers to vector data
161: call VecGetArrayRead(X,x_v,x_i,ierr)
162: call VecGetArray(G,g_v,g_i,ierr)
164: ! Compute G(X)
165: do i=0,nn-1
166: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
167: t2 = 1.0 - x_v(x_i + 2*i)
168: ff = ff + alpha*t1*t1 + t2*t2
169: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
170: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
171: enddo
173: ! Restore vectors
174: call VecRestoreArrayRead(X,x_v,x_i,ierr)
175: call VecRestoreArray(G,g_v,g_i,ierr)
177: f = ff
178: call PetscLogFlops(15.0d0*nn,ierr)
180: return
181: end
183: !
184: ! ---------------------------------------------------------------------
185: !
186: ! FormHessian - Evaluates Hessian matrix.
187: !
188: ! Input Parameters:
189: ! tao - the Tao context
190: ! X - input vector
191: ! dummy - optional user-defined context, as set by SNESSetHessian()
192: ! (not used here)
193: !
194: ! Output Parameters:
195: ! H - Hessian matrix
196: ! PrecH - optionally different preconditioning matrix (not used here)
197: ! flag - flag indicating matrix structure
198: ! ierr - error code
199: !
200: ! Note: Providing the Hessian may not be necessary. Only some solvers
201: ! require this matrix.
203: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
204: #include "rosenbrock1f.h"
206: ! Input/output variables:
207: Tao tao
208: Vec X
209: Mat H, PrecH
210: PetscErrorCode ierr
211: PetscInt dummy
213: PetscReal v(0:1,0:1)
214: PetscBool assembled
216: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
217: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
218: ! will return an array of doubles referenced by x_array offset by x_index.
219: ! i.e., to reference the kth element of X, use x_array(k + x_index).
220: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
221: PetscReal x_v(0:1)
222: PetscOffset x_i
223: PetscInt i,nn,ind(0:1),i2
225: 0
226: nn= n/2
227: i2 = 2
229: ! Zero existing matrix entries
230: call MatAssembled(H,assembled,ierr)
231: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
233: ! Get a pointer to vector data
235: call VecGetArrayRead(X,x_v,x_i,ierr)
237: ! Compute Hessian entries
239: do i=0,nn-1
240: v(1,1) = 2.0*alpha
241: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
242: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
243: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
244: v(0,1) = v(1,0)
245: ind(0) = 2*i
246: ind(1) = 2*i + 1
247: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
248: enddo
250: ! Restore vector
252: call VecRestoreArrayRead(X,x_v,x_i,ierr)
254: ! Assemble matrix
256: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
257: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
259: call PetscLogFlops(9.0d0*nn,ierr)
261: return
262: end
264: !
265: !/*TEST
266: !
267: ! build:
268: ! requires: !complex
269: !
270: ! test:
271: ! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
272: ! requires: !single
273: !
274: !TEST*/