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: Vec x ! solution vector
27: Mat H ! hessian matrix
28: Tao 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: if (size .ne. 1) then; SETERRA(PETSC_COMM_SELF,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
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: Tao tao
120: Vec X,G
121: PetscReal f
122: PetscErrorCode ierr
123: PetscInt dummy
125: PetscReal ff,t1,t2
126: PetscInt i,nn
128: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
129: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
130: ! will return an array of doubles referenced by x_array offset by x_index.
131: ! i.e., to reference the kth element of X, use x_array(k + x_index).
132: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
133: PetscReal g_v(0:1),x_v(0:1)
134: PetscOffset g_i,x_i
135: PetscReal alpha
136: PetscInt n
137: common /params/ alpha, n
139: 0
140: nn = n/2
141: ff = 0
143: ! Get pointers to vector data
144: PetscCall(VecGetArrayRead(X,x_v,x_i,ierr))
145: PetscCall(VecGetArray(G,g_v,g_i,ierr))
147: ! Compute G(X)
148: do i=0,nn-1
149: t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
150: t2 = 1.0 - x_v(x_i + 2*i)
151: ff = ff + alpha*t1*t1 + t2*t2
152: g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
153: g_v(g_i + 2*i + 1) = 2.0*alpha*t1
154: enddo
156: ! Restore vectors
157: PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr))
158: PetscCall(VecRestoreArray(G,g_v,g_i,ierr))
160: f = ff
161: PetscCall(PetscLogFlops(15.0d0*nn,ierr))
163: return
164: end
166: !
167: ! ---------------------------------------------------------------------
168: !
169: ! FormHessian - Evaluates Hessian matrix.
170: !
171: ! Input Parameters:
172: ! tao - the Tao context
173: ! X - input vector
174: ! dummy - optional user-defined context, as set by SNESSetHessian()
175: ! (not used here)
176: !
177: ! Output Parameters:
178: ! H - Hessian matrix
179: ! PrecH - optionally different preconditioning matrix (not used here)
180: ! flag - flag indicating matrix structure
181: ! ierr - error code
182: !
183: ! Note: Providing the Hessian may not be necessary. Only some solvers
184: ! require this matrix.
186: subroutine FormHessian(tao,X,H,PrecH,dummy,ierr)
187: use petsctao
188: implicit none
190: ! Input/output variables:
191: Tao tao
192: Vec X
193: Mat H, PrecH
194: PetscErrorCode ierr
195: PetscInt dummy
197: PetscReal v(0:1,0:1)
198: PetscBool assembled
200: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
201: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
202: ! will return an array of doubles referenced by x_array offset by x_index.
203: ! i.e., to reference the kth element of X, use x_array(k + x_index).
204: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
205: PetscReal x_v(0:1)
206: PetscOffset x_i
207: PetscInt i,nn,ind(0:1),i2
208: PetscReal alpha
209: PetscInt n
210: common /params/ alpha, n
212: 0
213: nn= n/2
214: i2 = 2
216: ! Zero existing matrix entries
217: PetscCall(MatAssembled(H,assembled,ierr))
218: if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H,ierr))
220: ! Get a pointer to vector data
222: PetscCall(VecGetArrayRead(X,x_v,x_i,ierr))
224: ! Compute Hessian entries
226: do i=0,nn-1
227: v(1,1) = 2.0*alpha
228: v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) - 3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
229: v(1,0) = -4.0*alpha*x_v(x_i+2*i)
230: v(0,1) = v(1,0)
231: ind(0) = 2*i
232: ind(1) = 2*i + 1
233: PetscCall(MatSetValues(H,i2,ind,i2,ind,v,INSERT_VALUES,ierr))
234: enddo
236: ! Restore vector
238: PetscCall(VecRestoreArrayRead(X,x_v,x_i,ierr))
240: ! Assemble matrix
242: PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,ierr))
243: PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,ierr))
245: PetscCall(PetscLogFlops(9.0d0*nn,ierr))
247: return
248: end
250: !
251: !/*TEST
252: !
253: ! build:
254: ! requires: !complex
255: !
256: ! test:
257: ! args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-5
258: ! requires: !single
259: !
260: !TEST*/