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: !
11: !
13: ! ----------------------------------------------------------------------
14: !
15: #include "rosenbrock1f.h"
17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18: ! Variable declarations
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: !
21: ! See additional variable declarations in the file rosenbrock1f.h
23: PetscErrorCode ierr ! used to check for functions returning nonzeros
24: Vec x ! solution vector
25: Mat H ! hessian matrix
26: Tao tao ! TAO_SOVER context
27: PetscBool flg
28: PetscInt i2,i1
29: PetscMPIInt size
30: PetscReal zero
32: ! Note: Any user-defined Fortran routines (such as FormGradient)
33: ! MUST be declared as external.
35: external FormFunctionGradient,FormHessian
37: zero = 0.0d0
38: i2 = 2
39: i1 = 1
41: ! Initialize TAO and PETSc
42: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
43: if (ierr .ne. 0) then
44: print*,'Unable to initialize PETSc'
45: stop
46: endif
48: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
49: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
51: ! Initialize problem parameters
52: n = 2
53: alpha = 99.0d0
55: ! Check for command line arguments to override defaults
56: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
57: & '-n',n,flg,ierr)
58: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
59: & '-alpha',alpha,flg,ierr)
61: ! Allocate vectors for the solution and gradient
62: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
64: ! Allocate storage space for Hessian;
65: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
66: & PETSC_NULL_INTEGER, H,ierr)
68: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
70: ! The TAO code begins here
72: ! Create TAO solver
73: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
74: CHKERRA(ierr)
75: call TaoSetType(tao,TAOLMVM,ierr)
76: CHKERRA(ierr)
78: ! Set routines for function, gradient, and hessian evaluation
79: call TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC, &
80: & FormFunctionGradient,0,ierr)
81: CHKERRA(ierr)
82: call TaoSetHessian(tao,H,H,FormHessian, &
83: & 0,ierr)
84: CHKERRA(ierr)
86: ! Optional: Set initial guess
87: call VecSet(x, zero, ierr)
88: call TaoSetSolution(tao, x, ierr)
89: CHKERRA(ierr)
91: ! Check for TAO command line options
92: call TaoSetFromOptions(tao,ierr)
93: CHKERRA(ierr)
95: ! SOLVE THE APPLICATION
96: call TaoSolve(tao,ierr)
98: ! TaoView() prints ierr about the TAO solver; the option
99: ! -tao_view
100: ! can alternatively be used to activate this at runtime.
101: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
103: ! Free TAO data structures
104: call TaoDestroy(tao,ierr)
106: ! Free PETSc data structures
107: call VecDestroy(x,ierr)
108: call MatDestroy(H,ierr)
110: call PetscFinalize(ierr)
111: end
113: ! --------------------------------------------------------------------
114: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
115: !
116: ! Input Parameters:
117: ! tao - the Tao context
118: ! X - input vector
119: ! dummy - not used
120: !
121: ! Output Parameters:
122: ! G - vector containing the newly evaluated gradient
123: ! f - function value
125: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
126: #include "rosenbrock1f.h"
128: Tao tao
129: Vec X,G
130: PetscReal f
131: PetscErrorCode ierr
132: PetscInt dummy
134: PetscReal ff,t1,t2
135: PetscInt i,nn
137: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
138: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
139: ! will return an array of doubles referenced by x_array offset by x_index.
140: ! i.e., to reference the kth element of X, use x_array(k + x_index).
141: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
142: PetscReal g_v(0:1),x_v(0:1)
143: PetscOffset g_i,x_i
145: 0
146: nn = n/2
147: ff = 0
149: ! Get pointers to vector data
150: call VecGetArrayRead(X,x_v,x_i,ierr)
151: call VecGetArray(G,g_v,g_i,ierr)
153: ! Compute G(X)
154: do i=0,nn-1
155: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
156: t2 = 1.0 - x_v(x_i + 2*i)
157: ff = ff + alpha*t1*t1 + t2*t2
158: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
159: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
160: enddo
162: ! Restore vectors
163: call VecRestoreArrayRead(X,x_v,x_i,ierr)
164: call VecRestoreArray(G,g_v,g_i,ierr)
166: f = ff
167: call PetscLogFlops(15.0d0*nn,ierr)
169: return
170: end
172: !
173: ! ---------------------------------------------------------------------
174: !
175: ! FormHessian - Evaluates Hessian matrix.
176: !
177: ! Input Parameters:
178: ! tao - the Tao context
179: ! X - input vector
180: ! dummy - optional user-defined context, as set by SNESSetHessian()
181: ! (not used here)
182: !
183: ! Output Parameters:
184: ! H - Hessian matrix
185: ! PrecH - optionally different preconditioning matrix (not used here)
186: ! flag - flag indicating matrix structure
187: ! ierr - error code
188: !
189: ! Note: Providing the Hessian may not be necessary. Only some solvers
190: ! require this matrix.
192: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
193: #include "rosenbrock1f.h"
195: ! Input/output variables:
196: Tao tao
197: Vec X
198: Mat H, PrecH
199: PetscErrorCode ierr
200: PetscInt dummy
202: PetscReal v(0:1,0:1)
203: PetscBool assembled
205: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
206: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
207: ! will return an array of doubles referenced by x_array offset by x_index.
208: ! i.e., to reference the kth element of X, use x_array(k + x_index).
209: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
210: PetscReal x_v(0:1)
211: PetscOffset x_i
212: PetscInt i,nn,ind(0:1),i2
214: 0
215: nn= n/2
216: i2 = 2
218: ! Zero existing matrix entries
219: call MatAssembled(H,assembled,ierr)
220: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
222: ! Get a pointer to vector data
224: call VecGetArrayRead(X,x_v,x_i,ierr)
226: ! Compute Hessian entries
228: do i=0,nn-1
229: v(1,1) = 2.0*alpha
230: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
231: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
232: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
233: v(0,1) = v(1,0)
234: ind(0) = 2*i
235: ind(1) = 2*i + 1
236: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
237: enddo
239: ! Restore vector
241: call VecRestoreArrayRead(X,x_v,x_i,ierr)
243: ! Assemble matrix
245: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
246: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
248: call PetscLogFlops(9.0d0*nn,ierr)
250: return
251: end
253: !
254: !/*TEST
255: !
256: ! build:
257: ! requires: !complex
258: !
259: ! test:
260: ! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
261: ! requires: !single
262: !
263: !TEST*/