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, 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: return
157: end
159: !
160: ! ---------------------------------------------------------------------
161: !
162: ! FormHessian - Evaluates Hessian matrix.
163: !
164: ! Input Parameters:
165: ! tao - the Tao context
166: ! X - input vector
167: ! dummy - optional user-defined context, as set by SNESSetHessian()
168: ! (not used here)
169: !
170: ! Output Parameters:
171: ! H - Hessian matrix
172: ! PrecH - optionally different preconditioning matrix (not used here)
173: ! flag - flag indicating matrix structure
174: ! ierr - error code
175: !
176: ! Note: Providing the Hessian may not be necessary. Only some solvers
177: ! require this matrix.
179: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
180: use petsctao
181: implicit none
183: ! Input/output variables:
184: type(tTao) tao
185: type(tVec) X
186: type(tMat) H, PrecH
187: PetscErrorCode ierr
188: PetscInt dummy
190: PetscReal v(0:1,0:1)
191: PetscBool assembled
193: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
194: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
195: ! will return an array of doubles referenced by x_array offset by x_index.
196: ! i.e., to reference the kth element of X, use x_array(k + x_index).
197: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
198: PetscReal, pointer :: x_v(:)
199: PetscInt i,nn,ind(0:1),i2
200: PetscReal alpha
201: PetscInt n
202: common /params/ alpha, n
204: ierr = 0
205: nn= n/2
206: i2 = 2
208: ! Zero existing matrix entries
209: PetscCall(MatAssembled(H,assembled,ierr))
210: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr))
212: ! Get a pointer to vector data
214: PetscCall(VecGetArrayReadF90(X,x_v,ierr))
216: ! Compute Hessian entries
218: do i=0,nn-1
219: v(1,1) = 2.0*alpha
220: v(0,0) = -4.0*alpha*(x_v(1+2*i+1) - 3*x_v(1+2*i)*x_v(1+2*i))+2
221: v(1,0) = -4.0*alpha*x_v(1+2*i)
222: v(0,1) = v(1,0)
223: ind(0) = 2*i
224: ind(1) = 2*i + 1
225: PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr))
226: enddo
228: ! Restore vector
230: PetscCall(VecRestoreArrayReadF90(X,x_v,ierr))
232: ! Assemble matrix
234: PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
235: PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))
237: PetscCall(PetscLogFlops(9.0d0*nn,ierr))
239: return
240: end
242: !
243: !/*TEST
244: !
245: ! build:
246: ! requires: !complex
247: !
248: ! test:
249: ! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
250: ! requires: !single
251: !
252: !TEST*/