Actual source code: rosenbrock1f.F90
petsc-3.8.4 2018-03-24
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: #include "rosenbrock1f.h"
27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28: ! Variable declarations
29: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
30: !
31: ! See additional variable declarations in the file rosenbrock1f.h
33: PetscErrorCode ierr ! used to check for functions returning nonzeros
34: Vec x ! solution vector
35: Mat H ! hessian matrix
36: Tao tao ! TAO_SOVER context
37: PetscBool flg
38: PetscInt i2,i1
39: integer size,rank ! number of processes running
40: PetscReal zero
42: ! Note: Any user-defined Fortran routines (such as FormGradient)
43: ! MUST be declared as external.
45: external FormFunctionGradient,FormHessian
47: zero = 0.0d0
48: i2 = 2
49: i1 = 1
51: ! Initialize TAO and PETSc
52: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
53: if (ierr .ne. 0) then
54: print*,'Unable to initialize PETSc'
55: stop
56: endif
58: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
59: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
60: if (size .ne. 1) then
61: if (rank .eq. 0) then
62: write(6,*) 'This is a uniprocessor example only!'
63: endif
64: SETERRA(PETSC_COMM_SELF,1,' ')
65: endif
67: ! Initialize problem parameters
68: n = 2
69: alpha = 99.0d0
73: ! Check for command line arguments to override defaults
74: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
75: & '-n',n,flg,ierr)
76: call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
77: & '-alpha',alpha,flg,ierr)
79: ! Allocate vectors for the solution and gradient
80: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
82: ! Allocate storage space for Hessian;
83: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
84: & PETSC_NULL_INTEGER, H,ierr)
86: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
89: ! The TAO code begins here
91: ! Create TAO solver
92: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
93: CHKERRA(ierr)
94: call TaoSetType(tao,TAOLMVM,ierr)
95: CHKERRA(ierr)
97: ! Set routines for function, gradient, and hessian evaluation
98: call TaoSetObjectiveAndGradientRoutine(tao, &
99: & FormFunctionGradient,0,ierr)
100: CHKERRA(ierr)
101: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
102: & 0,ierr)
103: CHKERRA(ierr)
106: ! Optional: Set initial guess
107: call VecSet(x, zero, ierr)
108: call TaoSetInitialVector(tao, x, ierr)
109: CHKERRA(ierr)
112: ! Check for TAO command line options
113: call TaoSetFromOptions(tao,ierr)
114: CHKERRA(ierr)
116: ! SOLVE THE APPLICATION
117: call TaoSolve(tao,ierr)
119: ! TaoView() prints ierr about the TAO solver; the option
120: ! -tao_view
121: ! can alternatively be used to activate this at runtime.
122: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
125: ! Free TAO data structures
126: call TaoDestroy(tao,ierr)
128: ! Free PETSc data structures
129: call VecDestroy(x,ierr)
130: call MatDestroy(H,ierr)
132: call PetscFinalize(ierr)
133: end
136: ! --------------------------------------------------------------------
137: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
138: !
139: ! Input Parameters:
140: ! tao - the Tao context
141: ! X - input vector
142: ! dummy - not used
143: !
144: ! Output Parameters:
145: ! G - vector containing the newly evaluated gradient
146: ! f - function value
148: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
149: #include "rosenbrock1f.h"
151: Tao tao
152: Vec X,G
153: PetscReal f
154: PetscErrorCode ierr
155: PetscInt dummy
158: PetscReal ff,t1,t2
159: PetscInt i,nn
161: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
162: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
163: ! will return an array of doubles referenced by x_array offset by x_index.
164: ! i.e., to reference the kth element of X, use x_array(k + x_index).
165: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
166: PetscReal g_v(0:1),x_v(0:1)
167: PetscOffset g_i,x_i
169: 0
170: nn = n/2
171: ff = 0
173: ! Get pointers to vector data
174: call VecGetArray(X,x_v,x_i,ierr)
175: call VecGetArray(G,g_v,g_i,ierr)
178: ! Compute G(X)
179: do i=0,nn-1
180: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
181: t2 = 1.0 - x_v(x_i + 2*i)
182: ff = ff + alpha*t1*t1 + t2*t2
183: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
184: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
185: enddo
187: ! Restore vectors
188: call VecRestoreArray(X,x_v,x_i,ierr)
189: call VecRestoreArray(G,g_v,g_i,ierr)
191: f = ff
192: call PetscLogFlops(15.0d0*nn,ierr)
194: return
195: end
197: !
198: ! ---------------------------------------------------------------------
199: !
200: ! FormHessian - Evaluates Hessian matrix.
201: !
202: ! Input Parameters:
203: ! tao - the Tao context
204: ! X - input vector
205: ! dummy - optional user-defined context, as set by SNESSetHessian()
206: ! (not used here)
207: !
208: ! Output Parameters:
209: ! H - Hessian matrix
210: ! PrecH - optionally different preconditioning matrix (not used here)
211: ! flag - flag indicating matrix structure
212: ! ierr - error code
213: !
214: ! Note: Providing the Hessian may not be necessary. Only some solvers
215: ! require this matrix.
217: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
218: #include "rosenbrock1f.h"
220: ! Input/output variables:
221: Tao tao
222: Vec X
223: Mat H, PrecH
224: PetscErrorCode ierr
225: PetscInt dummy
227: PetscReal v(0:1,0:1)
228: PetscBool assembled
230: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
231: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
232: ! will return an array of doubles referenced by x_array offset by x_index.
233: ! i.e., to reference the kth element of X, use x_array(k + x_index).
234: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
235: PetscReal x_v(0:1)
236: PetscOffset x_i
237: PetscInt i,nn,ind(0:1),i2
240: 0
241: nn= n/2
242: i2 = 2
244: ! Zero existing matrix entries
245: call MatAssembled(H,assembled,ierr)
246: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
248: ! Get a pointer to vector data
250: call VecGetArray(X,x_v,x_i,ierr)
252: ! Compute Hessian entries
254: do i=0,nn-1
255: v(1,1) = 2.0*alpha
256: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
257: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
258: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
259: v(0,1) = v(1,0)
260: ind(0) = 2*i
261: ind(1) = 2*i + 1
262: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
263: enddo
265: ! Restore vector
267: call VecRestoreArray(X,x_v,x_i,ierr)
269: ! Assemble matrix
271: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
272: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
274: call PetscLogFlops(9.0d0*nn,ierr)
276: return
277: end