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) ta ! 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 == 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, ta, ierr))
72: PetscCallA(TaoSetType(ta, TAOLMVM, ierr))
74: ! Set routines for function, gradient, and hessian evaluation
75: PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
76: PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))
78: ! Optional: Set initial guess
79: PetscCallA(VecSet(x, zero, ierr))
80: PetscCallA(TaoSetSolution(ta, x, ierr))
82: ! Check for TAO command line options
83: PetscCallA(TaoSetFromOptions(ta, ierr))
85: ! SOLVE THE APPLICATION
86: PetscCallA(TaoSolve(ta, 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(ta,PETSC_VIEWER_STDOUT_SELF,ierr))
93: ! Free TAO data structures
94: PetscCallA(TaoDestroy(ta, 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(ta, X, f, G, dummy, ierr)
116: use petsctao
117: implicit none
119: type(tTao) ta
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(VecGetArrayRead(X, x_v, ierr))
138: PetscCall(VecGetArray(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: end do
149: ! Restore vectors
150: PetscCall(VecRestoreArrayRead(X, x_v, ierr))
151: PetscCall(VecRestoreArray(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 matrix used to compute the preconditioner (not used here)
172: ! ierr - error code
173: !
174: ! Note: Providing the Hessian may not be necessary. Only some solvers
175: ! require this matrix.
177: subroutine FormHessian(ta, X, H, PrecH, dummy, ierr)
178: use petsctao
179: implicit none
181: ! Input/output variables:
182: type(tTao) ta
183: type(tVec) X
184: type(tMat) H, PrecH
185: PetscErrorCode ierr
186: PetscInt dummy
188: PetscReal v(0:1, 0:1)
189: PetscBool assembled
191: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
192: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
193: ! will return an array of doubles referenced by x_array offset by x_index.
194: ! i.e., to reference the kth element of X, use x_array(k + x_index).
195: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
196: PetscReal, pointer :: x_v(:)
197: PetscInt i, nn, ind(0:1), i2
198: PetscReal alpha
199: PetscInt n
200: common/params/alpha, n
202: ierr = 0
203: nn = n/2
204: i2 = 2
206: ! Zero existing matrix entries
207: PetscCall(MatAssembled(H, assembled, ierr))
208: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H, ierr))
210: ! Get a pointer to vector data
212: PetscCall(VecGetArrayRead(X, x_v, ierr))
214: ! Compute Hessian entries
216: do i = 0, nn - 1
217: v(1, 1) = 2.0*alpha
218: v(0, 0) = -4.0*alpha*(x_v(1 + 2*i + 1) - 3*x_v(1 + 2*i)*x_v(1 + 2*i)) + 2
219: v(1, 0) = -4.0*alpha*x_v(1 + 2*i)
220: v(0, 1) = v(1, 0)
221: ind(0) = 2*i
222: ind(1) = 2*i + 1
223: PetscCall(MatSetValues(H, i2, ind, i2, ind, reshape(v, [i2*i2]), INSERT_VALUES, ierr))
224: end do
226: ! Restore vector
228: PetscCall(VecRestoreArrayRead(X, x_v, ierr))
230: ! Assemble matrix
232: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
233: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))
235: PetscCall(PetscLogFlops(9.0d0*nn, ierr))
237: end
239: !
240: !/*TEST
241: !
242: ! build:
243: ! requires: !complex
244: !
245: ! test:
246: ! args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5
247: ! requires: !single
248: !
249: !TEST*/