Actual source code: rosenbrock1f.F
petsc-3.7.5 2017-01-01
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*/
21: !
23: ! ----------------------------------------------------------------------
24: !
25: implicit none
27: #include "rosenbrock1f.h"
29: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: ! Variable declarations
31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
32: !
33: ! See additional variable declarations in the file rosenbrock1f.h
35: PetscErrorCode ierr ! used to check for functions returning nonzeros
36: Vec x ! solution vector
37: Mat H ! hessian matrix
38: Tao tao ! TAO_SOVER context
39: PetscBool flg
40: PetscInt i2,i1
41: integer size,rank ! number of processes running
42: PetscReal zero
44: ! Note: Any user-defined Fortran routines (such as FormGradient)
45: ! MUST be declared as external.
47: external FormFunctionGradient,FormHessian
49: zero = 0.0d0
50: i2 = 2
51: i1 = 1
53: ! Initialize TAO and PETSc
54: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
56: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
57: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
58: if (size .ne. 1) then
59: if (rank .eq. 0) then
60: write(6,*) 'This is a uniprocessor example only!'
61: endif
62: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
63: endif
65: ! Initialize problem parameters
66: n = 2
67: alpha = 99.0d0
71: ! Check for command line arguments to override defaults
72: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
73: & '-n',n,flg,ierr)
74: call PetscOptionsGetReal(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
75: & '-alpha',alpha,flg,ierr)
77: ! Allocate vectors for the solution and gradient
78: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
80: ! Allocate storage space for Hessian;
81: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
82: & PETSC_NULL_INTEGER, H,ierr)
84: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
87: ! The TAO code begins here
89: ! Create TAO solver
90: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
91: CHKERRQ(ierr)
92: call TaoSetType(tao,TAOLMVM,ierr)
93: CHKERRQ(ierr)
95: ! Set routines for function, gradient, and hessian evaluation
96: call TaoSetObjectiveAndGradientRoutine(tao, &
97: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
98: CHKERRQ(ierr)
99: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
100: & PETSC_NULL_OBJECT,ierr)
101: CHKERRQ(ierr)
104: ! Optional: Set initial guess
105: call VecSet(x, zero, ierr)
106: call TaoSetInitialVector(tao, x, ierr)
107: CHKERRQ(ierr)
110: ! Check for TAO command line options
111: call TaoSetFromOptions(tao,ierr)
112: CHKERRQ(ierr)
114: ! SOLVE THE APPLICATION
115: call TaoSolve(tao,ierr)
117: ! TaoView() prints ierr about the TAO solver; the option
118: ! -tao_view
119: ! can alternatively be used to activate this at runtime.
120: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
123: ! Free TAO data structures
124: call TaoDestroy(tao,ierr)
126: ! Free PETSc data structures
127: call VecDestroy(x,ierr)
128: call MatDestroy(H,ierr)
130: call PetscFinalize(ierr)
131: end
134: ! --------------------------------------------------------------------
135: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
136: !
137: ! Input Parameters:
138: ! tao - the Tao context
139: ! X - input vector
140: ! dummy - not used
141: !
142: ! Output Parameters:
143: ! G - vector containing the newly evaluated gradient
144: ! f - function value
146: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
147: implicit none
149: ! n,alpha defined in rosenbrock1f.h
150: #include "rosenbrock1f.h"
152: Tao tao
153: Vec X,G
154: PetscReal f
155: PetscErrorCode ierr
156: PetscInt dummy
159: PetscReal ff,t1,t2
160: PetscInt i,nn
162: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
163: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
164: ! will return an array of doubles referenced by x_array offset by x_index.
165: ! i.e., to reference the kth element of X, use x_array(k + x_index).
166: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
167: PetscReal g_v(0:1),x_v(0:1)
168: PetscOffset g_i,x_i
170: 0
171: nn = n/2
172: ff = 0
174: ! Get pointers to vector data
175: call VecGetArray(X,x_v,x_i,ierr)
176: call VecGetArray(G,g_v,g_i,ierr)
179: ! Compute G(X)
180: do i=0,nn-1
181: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
182: t2 = 1.0 - x_v(x_i + 2*i)
183: ff = ff + alpha*t1*t1 + t2*t2
184: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
185: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
186: enddo
188: ! Restore vectors
189: call VecRestoreArray(X,x_v,x_i,ierr)
190: call VecRestoreArray(G,g_v,g_i,ierr)
192: f = ff
193: call PetscLogFlops(15.0d0*nn,ierr)
195: return
196: end
198: !
199: ! ---------------------------------------------------------------------
200: !
201: ! FormHessian - Evaluates Hessian matrix.
202: !
203: ! Input Parameters:
204: ! tao - the Tao context
205: ! X - input vector
206: ! dummy - optional user-defined context, as set by SNESSetHessian()
207: ! (not used here)
208: !
209: ! Output Parameters:
210: ! H - Hessian matrix
211: ! PrecH - optionally different preconditioning matrix (not used here)
212: ! flag - flag indicating matrix structure
213: ! ierr - error code
214: !
215: ! Note: Providing the Hessian may not be necessary. Only some solvers
216: ! require this matrix.
218: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
219: implicit none
221: #include "rosenbrock1f.h"
223: ! Input/output variables:
224: Tao tao
225: Vec X
226: Mat H, PrecH
227: PetscErrorCode ierr
228: PetscInt dummy
230: PetscReal v(0:1,0:1)
231: PetscBool assembled
233: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
234: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
235: ! will return an array of doubles referenced by x_array offset by x_index.
236: ! i.e., to reference the kth element of X, use x_array(k + x_index).
237: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
238: PetscReal x_v(0:1)
239: PetscOffset x_i
240: PetscInt i,nn,ind(0:1),i2
243: 0
244: nn= n/2
245: i2 = 2
247: ! Zero existing matrix entries
248: call MatAssembled(H,assembled,ierr)
249: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
251: ! Get a pointer to vector data
253: call VecGetArray(X,x_v,x_i,ierr)
255: ! Compute Hessian entries
257: do i=0,nn-1
258: v(1,1) = 2.0*alpha
259: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
260: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
261: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
262: v(0,1) = v(1,0)
263: ind(0) = 2*i
264: ind(1) = 2*i + 1
265: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
266: enddo
268: ! Restore vector
270: call VecRestoreArray(X,x_v,x_i,ierr)
272: ! Assemble matrix
274: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
275: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
277: call PetscLogFlops(9.0d0*nn,ierr)
279: return
280: end