-
Notifications
You must be signed in to change notification settings - Fork 1
/
DriverRosenbrockp.f90~
327 lines (265 loc) · 12.6 KB
/
DriverRosenbrockp.f90~
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
!
! L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
! or “3-clause license”)
! Please read attached file License.txt
!
! DRIVER 3 in Fortran 90
! --------------------------------------------------------------
! TIME-CONTROLLED DRIVER FOR L-BFGS-B
! --------------------------------------------------------------
!
! L-BFGS-B is a code for solving large nonlinear optimization
! problems with simple bounds on the variables.
!
! The code can also be used for unconstrained problems and is
! as efficient for these problems as the earlier limited memory
! code L-BFGS.
!
! This driver shows how to terminate a run after some prescribed
! CPU time has elapsed, and how to print the desired information
! before exiting.
!
! References:
!
! [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
! memory algorithm for bound constrained optimization'',
! SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
!
! [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: FORTRAN
! Subroutines for Large Scale Bound Constrained Optimization''
! Tech. Report, NAM-11, EECS Department, Northwestern University,
! 1994.
!
!
! (Postscript files of these papers are available via anonymous
! ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
!
! * * *
!
! February 2011 (latest revision)
! Optimization Center at Northwestern University
! Instituto Tecnologico Autonomo de Mexico
!
! Jorge Nocedal and Jose Luis Morales
!
! **************
program driver
! This time-controlled driver shows that it is possible to terminate
! a run by elapsed CPU time, and yet be able to print all desired
! information. This driver also illustrates the use of two
! stopping criteria that may be used in conjunction with a limit
! on execution time.
implicit none
! We specify a limit on the CPU time (tlimit = 10 seconds)
!
! We suppress the default output (iprint = -1). The user could
! also elect to use the default output by choosing iprint >= 0.)
! We suppress the code-supplied stopping tests because we will
! provide our own termination conditions
! We specify the dimension n of the sample problem and the number
! m of limited memory corrections stored.
integer, parameter :: iprint = -10, numargs = 0, jmax = 20d0
integer, parameter :: dp = kind(1.0d0)
real(dp), parameter :: factr = 0.0d0, pgtol = 0.0d0, &
tlimit = 1000.0d0, taux = 1.0d-3
!
character(len=120) :: task, csave
logical :: lsave(4)
integer :: isave(47), xindex
real(dp) :: f, r, taud
real(dp) :: dsave(30)
integer, allocatable :: nbd(:), iwa(:)
real(dp), allocatable :: x(:), l(:), u(:), g(:), wa(:)
!
real(dp) :: t1, t2, time1, time2, p, z, r1
integer :: i, j, m, n, nfg, startx, startg, nbisect = 0
! Get outside parameters
CHARACTER(len=32) :: arg
CALL getarg(1 , arg)
read (arg,*) m
call getarg(2, arg)
read(arg,*) n
call getarg(3, arg)
read(arg,*) p
call getarg(4, arg)
read(arg,*) taud
allocate ( nbd(n), x(n), l(n), u(n), g(n) )
allocate ( iwa(3*n) )
allocate ( wa(2*m*n + 5*n + 11*m*m + 8*m + 2*jmax*n) )
wa = 0.0d0
! This time-controlled driver shows that it is possible to terminate
! a run by elapsed CPU time, and yet be able to print all desired
! information. This driver also illustrates the use of two
! stopping criteria that may be used in conjunction with a limit
! on execution time. The sample problem used here is the same as in
! driver1 and driver2 (the extended Rosenbrock function with bounds
! on the variables).
! We now specify nbd which defines the bounds on the variables:
! l specifies the lower bounds,
! u specifies the upper bounds.
! First set bounds on the odd-numbered variables.
do 10 i=1, n,2
nbd(i)=2
l(i)=1.0d1
u(i)=1.0d2
10 continue
! Next set bounds on the even-numbered variables.
do 12 i=2, n,2
nbd(i)=2
l(i)=-1.0d2
u(i)=1.0d2
12 continue
! We now define the starting point.
do i=1, n
x(i) = (u(i) + l(i)) / 2.0d0
enddo
x(1) = x(1) - 1d0
do 14 i=2, n
x(i)=x(i) - (1-2d0**(1-i))
14 continue
x = 3d0
! do 1454 i=1, n
! x(i)=3.0d0
!1454 continue
! x(1) = 0.5d0
!x(2) = 0.6d0
!x(3) = 0.7d0
! We now write the heading of the output.
! write (6,16)
!16 format(/,5x, 'Solving sample problem.',&
! /,5x, ' (f = 0.0 at the optimal solution.)',/)
! We start the iteration by initializing task.
task = 'START'
! ------- the beginning of the loop ----------
! We begin counting the CPU time.
call timer(time1)
nfg = 0
do while( task(1:2).eq.'FG'.or.task.eq.'NEW_X'.or. &
task.eq.'START')
! This is the call to the L-BFGS-B code.
call setulb(n,m,x,l,u,nbd,f,g,factr,pgtol,wa,iwa, &
task,iprint, csave,lsave,isave,dsave,taux, nfg, jmax, taud, nbisect)
!write (*,*) 'final results rosenbrock (new) run:', m, n, p, isave(30), isave(34), f, dsave(30), task
!write(*,*) 'number of bisections: ', nbisect
!write (*,*) 'partial results iteration, f:', isave(30), f
if (task(1:2) .eq. 'FG') then
! the minimization routine has returned to request the
! function f and gradient g values at the current x.
! Before evaluating f and g we check the CPU time spent.
call timer(time2)
! if (time2-time1 .gt. tlimit) then
if (.false.) then
task='STOP: CPU EXCEEDING THE TIME LIMIT.'
! Note: Assigning task(1:4)='STOP' will terminate the run;
! setting task(7:9)='CPU' will restore the information at
! the latest iterate generated by the code so that it can
! be correctly printed by the driver.
! In this driver we have chosen to disable the
! printing options of the code (we set iprint=-1);
! instead we are using customized output: we print the
! latest value of x, the corresponding function value f and
! the norm of the projected gradient |proj g|.
! We print out the information contained in task.
write (6,*) task
! We print the latest iterate contained in wa(j+1:j+n), where
j = 3*n+2*m*n+11*m**2
write (6,*) 'Latest iterate X ='
write (6,'((1x,1p, 6(1x,d11.4)))') (wa(i),i = j+1,j+n)
! We print the function value f and the norm of the projected
! gradient |proj g| at the last iterate; they are stored in
! dsave(2) and dsave(13) respectively.
write (6,'(a,1p,d12.5,4x,a,1p,d12.5)') &
'At latest iterate f =',dsave(2),'|proj g| =',dsave(13)
else
! The time limit has not been reached and we compute
! the function value f for the sample problem.
do 1823 i=1,n
g(i)=0d0
1823 continue
f=((x(1)-1d0)**2)
g(1)=2d0*(x(1)-1d0)
do 20 i=1, (n-1)
z=x(i+1)-x(i)**2
f=f+abs(z)**p
! r1=p*(z**2)**(p/2)/z
r1=p * abs(z)**(p - 1)
if (z < 0) then
r1 = -r1
endif
g(i+1)=r1
g(i)=g(i)-2d0*x(i)*r1
20 continue
nfg = nfg + 1
! Write f on the wa matrix for later stopping condition
startx = isave(46) + mod(nfg - 1, jmax) * n
startg = isave(45) + mod(nfg - 1, jmax) * n
xindex = 0
do i = startx, (startx + n - 1)
xindex = xindex + 1
wa(i) = x(xindex)
end do
xindex = 0
do i = startg, (startg + n - 1)
xindex = xindex + 1
wa(i) = g(xindex)
end do
!write(*,*) 'x', x
!write(*, *) 'wa', wa(isave(46):(isave(46)+(jmax*n)-1))
! write (6,*) 'Current X for debugging ='
! write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
! write (6,*) 'Current g for debugging ='
! write (6,'((1x,1p, 6(1x,d11.4)))') (g(i),i = 1,n)
! write (6,*) 'Current f value ='
! write (*,'(A, F8.3)') 'f=', f
endif
! go back to the minimization routine.
else
if (task(1:5) .eq. 'NEW_X') then
! THE ONLY STOPPING CONDITION THAT WE ARE GOING TO USE
if (isave(30) .eq. 10000) &
task= 'STOP: TOTAL NUMBER OF ITERATIONS REACHED 10000'
! the minimization routine has returned with a new iterate.
! The time limit has not been reached, and we test whether
! the following two stopping tests are satisfied:
! 1) Terminate if the total number of f and g evaluations
! exceeds 900.
!if (isave(34) .ge. 90000) &
! task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'
! 2) Terminate if |proj g|/(1+|f|) < 1.0d-10.
!if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) &
! task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
! We wish to print the following information at each iteration:
! 1) the current iteration number, isave(30),
! 2) the total number of f and g evaluations, isave(34),
! 3) the value of the objective function f,
! 4) the norm of the projected gradient, dsve(13)
!
! See the comments at the end of driver1 for a description
! of the variables isave and dsave.
!write (6,'(2(a,i5,4x),a,1p,d12.5,4x,a,1p,d12.5)') 'Iterate' &
! ,isave(30),'nfg =',isave(34),'f =',f,'|proj g| =',dsave(13)
! If the run is to be terminated, we print also the information
! contained in task as well as the final value of x.
!if (task(1:4) .eq. 'STOP') then
! write (6,*) task
! write (6,*) 'Final X='
! write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
!endif
endif
end if
end do
! If task is neither FG nor NEW_X we terminate execution.
call timer(time2)
write (*,*) 'final results rosenbrock (new) run:', m, n, p, isave(30), isave(34), f, dsave(30), task
!write (6,*) task
!write (6,*) 'Final X='
!write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
!write (6,*) 'Final G='
!write (6,'((1x,1p, 6(1x,d11.4)))') (g(i),i = 1,n)
!write (*,*) 'mnp = ', m, n, p
!write (*,*) 'Value of f = (NEW CODE)', f
!write (*,*) 'Iterate ', isave(30)
!write (6,*) 'Final X='
!write (6,'((1x,1p, 6(1x,d11.4)))') (x(i),i = 1,n)
end program driver
!======================= The end of driver3 ============================