Actual source code: rosenbrock1f.F
petsc-3.6.1 2015-08-06
1: ! Program usage: mpirun -np 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: TaoGetConvergedReason(); 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
43: TaoConvergedReason reason
47: ! Note: Any user-defined Fortran routines (such as FormGradient)
48: ! MUST be declared as external.
50: external FormFunctionGradient,FormHessian
52: zero = 0.0d0
53: i2 = 2
54: i1 = 1
56: ! Initialize TAO and PETSc
57: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
59: call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
60: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
61: if (size .ne. 1) then
62: if (rank .eq. 0) then
63: write(6,*) 'This is a uniprocessor example only!'
64: endif
65: SETERRQ(PETSC_COMM_SELF,1,' ',ierr)
66: endif
68: ! Initialize problem parameters
69: n = 2
70: alpha = 99.0d0
74: ! Check for command line arguments to override defaults
75: call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg, &
76: & ierr)
77: call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha', &
78: & alpha,flg,ierr)
80: ! Allocate vectors for the solution and gradient
81: call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
83: ! Allocate storage space for Hessian;
84: call MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1, &
85: & PETSC_NULL_INTEGER, H,ierr)
87: call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr)
90: ! The TAO code begins here
92: ! Create TAO solver
93: call TaoCreate(PETSC_COMM_SELF,tao,ierr)
94: CHKERRQ(ierr)
95: call TaoSetType(tao,TAOLMVM,ierr)
96: CHKERRQ(ierr)
98: ! Set routines for function, gradient, and hessian evaluation
99: call TaoSetObjectiveAndGradientRoutine(tao, &
100: & FormFunctionGradient,PETSC_NULL_OBJECT,ierr)
101: CHKERRQ(ierr)
102: call TaoSetHessianRoutine(tao,H,H,FormHessian, &
103: & PETSC_NULL_OBJECT,ierr)
104: CHKERRQ(ierr)
107: ! Optional: Set initial guess
108: call VecSet(x, zero, ierr)
109: call TaoSetInitialVector(tao, x, ierr)
110: CHKERRQ(ierr)
113: ! Check for TAO command line options
114: call TaoSetFromOptions(tao,ierr)
115: CHKERRQ(ierr)
117: ! SOLVE THE APPLICATION
118: call TaoSolve(tao,ierr)
120: call TaoGetConvergedReason(tao, reason, ierr)
121: if (reason .le. 0) then
122: print *,'Try a different TAO method, adjust some parameters,'
123: print *,'or check the function evaluation routines.'
124: endif
126: ! TaoView() prints ierr about the TAO solver; the option
127: ! -tao_view
128: ! can alternatively be used to activate this at runtime.
129: ! call TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr)
132: ! Free TAO data structures
133: call TaoDestroy(tao,ierr)
135: ! Free PETSc data structures
136: call VecDestroy(x,ierr)
137: call MatDestroy(H,ierr)
139: call PetscFinalize(ierr)
140: end
143: ! --------------------------------------------------------------------
144: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
145: !
146: ! Input Parameters:
147: ! tao - the Tao context
148: ! X - input vector
149: ! dummy - not used
150: !
151: ! Output Parameters:
152: ! G - vector containing the newly evaluated gradient
153: ! f - function value
155: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
156: implicit none
158: ! n,alpha defined in rosenbrock1f.h
159: #include "rosenbrock1f.h"
161: Tao tao
162: Vec X,G
163: PetscReal f
164: PetscErrorCode ierr
165: PetscInt dummy
168: PetscReal ff,t1,t2
169: PetscInt i,nn
171: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
172: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
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: PetscReal g_v(0:1),x_v(0:1)
177: PetscOffset g_i,x_i
179: 0
180: nn = n/2
181: ff = 0
183: ! Get pointers to vector data
184: call VecGetArray(X,x_v,x_i,ierr)
185: call VecGetArray(G,g_v,g_i,ierr)
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,ierr)
199: call VecRestoreArray(G,g_v,g_i,ierr)
201: f = ff
202: call PetscLogFlops(15.0d0*nn,ierr)
204: return
205: end
207: !
208: ! ---------------------------------------------------------------------
209: !
210: ! FormHessian - Evaluates Hessian matrix.
211: !
212: ! Input Parameters:
213: ! tao - the Tao 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: ! ierr - error code
223: !
224: ! Note: Providing the Hessian may not be necessary. Only some solvers
225: ! require this matrix.
227: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
228: implicit none
230: #include "rosenbrock1f.h"
232: ! Input/output variables:
233: Tao tao
234: Vec X
235: Mat H, PrecH
236: PetscErrorCode ierr
237: PetscInt dummy
239: PetscReal v(0:1,0:1)
240: PetscBool assembled
242: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
243: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
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: PetscReal x_v(0:1)
248: PetscOffset x_i
249: PetscInt i,nn,ind(0:1),i2
252: 0
253: nn= n/2
254: i2 = 2
256: ! Zero existing matrix entries
257: call MatAssembled(H,assembled,ierr)
258: if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(H,ierr)
260: ! Get a pointer to vector data
262: call VecGetArray(X,x_v,x_i,ierr)
264: ! Compute Hessian entries
266: do i=0,nn-1
267: v(1,1) = 2.0*alpha
268: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - &
269: & 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
270: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
271: v(0,1) = v(1,0)
272: ind(0) = 2*i
273: ind(1) = 2*i + 1
274: call MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr)
275: enddo
277: ! Restore vector
279: call VecRestoreArray(X,x_v,x_i,ierr)
281: ! Assemble matrix
283: call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr)
284: call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr)
286: call PetscLogFlops(9.0d0*nn,ierr)
288: return
289: end