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 "petsc/finclude/petsctao.h"
16: use petsctao
17: implicit none
19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: ! Variable declarations
21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
22: !
23: ! See additional variable declarations in the file rosenbrock1f.h
25: PetscErrorCode ierr ! used to check for functions returning nonzeros
26: type(tVec) x ! solution vector
27: type(tMat) H ! hessian matrix
28: type(tTao) tao ! TAO_SOVER context
29: PetscBool flg
30: PetscInt i2,i1
31: PetscMPIInt size
32: PetscReal zero
33: PetscReal alpha
34: PetscInt n
35: common /params/ alpha, n
37: ! Note: Any user-defined Fortran routines (such as FormGradient)
38: ! MUST be declared as external.
40: external FormFunctionGradient,FormHessian
42: zero = 0.0d0
43: i2 = 2
44: i1 = 1
46: ! Initialize TAO and PETSc
47: PetscCallA(PetscInitialize(ierr))
49: PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
50: PetscCheckA(size .eq. 1,PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only')
52: ! Initialize problem parameters
53: n = 2
54: alpha = 99.0d0
56: ! Check for command line arguments to override defaults
57: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
58: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-alpha',alpha,flg,ierr))
60: ! Allocate vectors for the solution and gradient
61: PetscCallA(VecCreateSeq(PETSC_COMM_SELF,n,x,ierr))
63: ! Allocate storage space for Hessian;
64: PetscCallA(MatCreateSeqBAIJ(PETSC_COMM_SELF,i2,n,n,i1,PETSC_NULL_INTEGER_ARRAY, H,ierr))
66: PetscCallA(MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr))
68: ! The TAO code begins here
70: ! Create TAO solver
71: PetscCallA(TaoCreate(PETSC_COMM_SELF,tao,ierr))
72: PetscCallA(TaoSetType(tao,TAOLMVM,ierr))
74: ! Set routines for function, gradient, and hessian evaluation
75: PetscCallA(TaoSetObjectiveAndGradient(tao,PETSC_NULL_VEC,FormFunctionGradient,0,ierr))
76: PetscCallA(TaoSetHessian(tao,H,H,FormHessian,0,ierr))
78: ! Optional: Set initial guess
79: PetscCallA(VecSet(x, zero, ierr))
80: PetscCallA(TaoSetSolution(tao, x, ierr))
82: ! Check for TAO command line options
83: PetscCallA(TaoSetFromOptions(tao,ierr))
85: ! SOLVE THE APPLICATION
86: PetscCallA(TaoSolve(tao,ierr))
88: ! TaoView() prints ierr about the TAO solver; the option
89: ! -tao_view
90: ! can alternatively be used to activate this at runtime.
91: ! PetscCallA(TaoView(tao,PETSC_VIEWER_STDOUT_SELF,ierr))
93: ! Free TAO data structures
94: PetscCallA(TaoDestroy(tao,ierr))
96: ! Free PETSc data structures
97: PetscCallA(VecDestroy(x,ierr))
98: PetscCallA(MatDestroy(H,ierr))
100: PetscCallA(PetscFinalize(ierr))
101: end
103: ! --------------------------------------------------------------------
104: ! FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
105: !
106: ! Input Parameters:
107: ! tao - the Tao context
108: ! X - input vector
109: ! dummy - not used
110: !
111: ! Output Parameters:
112: ! G - vector containing the newly evaluated gradient
113: ! f - function value
115: subroutine FormFunctionGradient(tao, X, f, G, dummy, ierr)
116: use petsctao
117: implicit none
119: type(tTao) tao
120: type(tVec) X,G
121: PetscReal f
122: PetscErrorCode ierr
123: PetscInt dummy
125: PetscReal ff,t1,t2
126: PetscInt i,nn
127: PetscReal, pointer :: g_v(:),x_v(:)
128: PetscReal alpha
129: PetscInt n
130: common /params/ alpha, n
132: ierr = 0
133: nn = n/2
134: ff = 0
136: ! Get pointers to vector data
137: PetscCall(VecGetArrayReadF90(X,x_v,ierr))
138: PetscCall(VecGetArrayF90(G,g_v,ierr))
140: ! Compute G(X)
141: do i=0,nn-1
142: t1 = x_v(1+2*i+1) - x_v(1+2*i)*x_v(1+2*i)
143: t2 = 1.0 - x_v(1 + 2*i)
144: ff = ff + alpha*t1*t1 + t2*t2
145: g_v(1 + 2*i) = -4*alpha*t1*x_v(1 + 2*i) - 2.0*t2
146: g_v(1 + 2*i + 1) = 2.0*alpha*t1
147: enddo
149: ! Restore vectors
150: PetscCall(VecRestoreArrayReadF90(X,x_v,ierr))
151: PetscCall(VecRestoreArrayF90(G,g_v,ierr))
153: f = ff
154: PetscCall(PetscLogFlops(15.0d0*nn,ierr))
156: end
158: !
159: ! ---------------------------------------------------------------------
160: !
161: ! FormHessian - Evaluates Hessian matrix.
162: !
163: ! Input Parameters:
164: ! tao - the Tao context
165: ! X - input vector
166: ! dummy - optional user-defined context, as set by SNESSetHessian()
167: ! (not used here)
168: !
169: ! Output Parameters:
170: ! H - Hessian matrix
171: ! PrecH - optionally different preconditioning matrix (not used here)
172: ! flag - flag indicating matrix structure
173: ! ierr - error code
174: !
175: ! Note: Providing the Hessian may not be necessary. Only some solvers
176: ! require this matrix.
178: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
179: use petsctao
180: implicit none
182: ! Input/output variables:
183: type(tTao) tao
184: type(tVec) X
185: type(tMat) H, PrecH
186: PetscErrorCode ierr
187: PetscInt dummy
189: PetscReal v(0:1,0:1)
190: PetscBool assembled
192: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
193: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
194: ! will return an array of doubles referenced by x_array offset by x_index.
195: ! i.e., to reference the kth element of X, use x_array(k + x_index).
196: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
197: PetscReal, pointer :: x_v(:)
198: PetscInt i,nn,ind(0:1),i2
199: PetscReal alpha
200: PetscInt n
201: common /params/ alpha, n
203: ierr = 0
204: nn= n/2
205: i2 = 2
207: ! Zero existing matrix entries
208: PetscCall(MatAssembled(H,assembled,ierr))
209: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr))
211: ! Get a pointer to vector data
213: PetscCall(VecGetArrayReadF90(X,x_v,ierr))
215: ! Compute Hessian entries
217: do i=0,nn-1
218: v(1,1) = 2.0*alpha
219: v(0,0) = -4.0*alpha*(x_v(1+2*i+1) - 3*x_v(1+2*i)*x_v(1+2*i))+2
220: v(1,0) = -4.0*alpha*x_v(1+2*i)
221: v(0,1) = v(1,0)
222: ind(0) = 2*i
223: ind(1) = 2*i + 1
224: PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr))
225: enddo
227: ! Restore vector
229: PetscCall(VecRestoreArrayReadF90(X,x_v,ierr))
231: ! Assemble matrix
233: PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
234: PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))
236: PetscCall(PetscLogFlops(9.0d0*nn,ierr))
238: end
240: !
241: !/*TEST
242: !
243: ! build:
244: ! requires: !complex
245: !
246: ! test:
247: ! args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5
248: ! requires: !single
249: !
250: !TEST*/