L-BFGS-B  3.0
Large-scale Bound-constrained Optimization
subsm.f
Go to the documentation of this file.
1 c> \file subsm.f
2 
3 c> \brief Performs the subspace minimization.
4 c>
5 c> Given xcp, l, u, r, an index set that specifies
6 c> the active set at xcp, and an l-BFGS matrix B
7 c> (in terms of WY, WS, SY, WT, head, col, and theta),
8 c> this subroutine computes an approximate solution
9 c> of the subspace problem
10 c>
11 c> (P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)
12 c>
13 c> subject to l<=x<=u
14 c> x_i=xcp_i for all i in A(xcp)
15 c>
16 c> along the subspace unconstrained Newton direction
17 c>
18 c> d = -(Z'BZ)^(-1) r.
19 c>
20 c> The formula for the Newton direction, given the L-BFGS matrix
21 c> and the Sherman-Morrison formula, is
22 c>
23 c> d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.
24 c>
25 c> where
26 c> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
27 c> [L_a -R_z theta*S'AA'S ]
28 c>
29 c> Note that this procedure for computing d differs
30 c> from that described in [1]. One can show that the matrix K is
31 c> equal to the matrix M^[-1]N in that paper.
32 c>
33 c> @param n On entry n is the dimension of the problem.<br/>
34 c> On exit n is unchanged.
35 c>
36 c> @param m On entry m is the maximum number of variable metric corrections
37 c> used to define the limited memory matrix.<br/>
38 c> On exit m is unchanged.
39 c>
40 c> @param nsub On entry nsub is the number of free variables.<br/>
41 c> On exit nsub is unchanged.
42 c>
43 c> @param ind On entry ind specifies the coordinate indices of free variables.<br/>
44 c> On exit ind is unchanged.
45 c>
46 c> @param l On entry l is the lower bound of x.<br/>
47 c> On exit l is unchanged.
48 c>
49 c> @param u On entry u is the upper bound of x.<br/>
50 c> On exit u is unchanged.
51 c>
52 c> @param nbd On entry nbd represents the type of bounds imposed on the
53 c> variables, and must be specified as follows:
54 c> nbd(i)=<ul><li>0 if x(i) is unbounded,</li>
55 c> <li>1 if x(i) has only a lower bound,</li>
56 c> <li>2 if x(i) has both lower and upper bounds, and</li>
57 c> <li>3 if x(i) has only an upper bound.</li></ul>
58 c> On exit nbd is unchanged.
59 c>
60 c> @param x On entry x specifies the Cauchy point xcp.<br/>
61 c> On exit x(i) is the minimizer of Q over the subspace of
62 c> free variables.
63 c>
64 c> @param d On entry d is the reduced gradient of Q at xcp.<br/>
65 c> On exit d is the Newton direction of Q.
66 c>
67 c> @param xp used to safeguard the projected Newton direction<br/>
68 c>
69 c> @param xx On entry it holds the current iterate.<br/>
70 c> On output it is unchanged.
71 c>
72 c> @param gg On entry it holds the gradient at the current iterate.<br/>
73 c> On output it is unchanged.
74 c>
75 c> @param ws On entry this stores S, a set of s-vectors, that defines the
76 c> limited memory BFGS matrix.<br/>
77 c> On exit this array is unchanged.
78 c>
79 c> @param wy On entry this stores Y, a set of y-vectors, that defines the
80 c> limited memory BFGS matrix.<br/>
81 c> On exit this array is unchanged.
82 c>
83 c> @param theta On entry theta is the scaling factor specifying B_0 = theta I.<br/>
84 c> On exit theta is unchanged.
85 c>
86 c> @param col On entry col is the actual number of variable metric
87 c> corrections stored so far.<br/>
88 c> On exit col is unchanged.
89 c>
90 c> @param head On entry head is the location of the first s-vector (or y-vector)
91 c> in S (or Y).<br/>
92 c> On exit col is unchanged.
93 c>
94 c> @param iword On entry iword is unspecified.<br/>
95 c> On exit iword specifies the status of the subspace solution.
96 c> iword = <ul><li>0 if the solution is in the box,</li>
97 c> <li>1 if some bound is encountered.</li></ul>
98 c>
99 c> @param wv working array
100 c>
101 c> @param wn On entry the upper triangle of wn stores the LEL^T factorization
102 c> of the indefinite matrix<br/>
103 c> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
104 c> [L_a -R_z theta*S'AA'S ]<br/>
105 c> where E = [-I 0]
106 c> [ 0 I]<br/>
107 c> On exit wn is unchanged.
108 c>
109 c> @param iprint must be set by the user;
110 c> It controls the frequency and type of output generated:<ul>
111 c> <li>iprint<0 no output is generated;</li>
112 c> <li>iprint=0 print only one line at the last iteration;</li>
113 c> <li>0<iprint<99 print also f and |proj g| every iprint iterations;</li>
114 c> <li>iprint=99 print details of every iteration except n-vectors;</li>
115 c> <li>iprint=100 print also the changes of active set and final x;</li>
116 c> <li>iprint>100 print details of every iteration including x and g;</li></ul>
117 c> When iprint > 0, the file iterate.dat will be created to
118 c> summarize the iteration.
119 c>
120 c>
121 c> Historical note: this routine used to take an `info` output parameter
122 c> for an "ill-conditioned K" status. Since K's conditioning is checked
123 c> in `formt`/`formk` (which fail loudly with their own `info`) and the
124 c> two `dtrsm` calls inside `subsm` cannot fail on a non-singular factor,
125 c> the parameter was always 0 on exit and has been removed.
126  subroutine subsm ( n, m, nsub, ind, l, u, nbd, x, d, xp, ws, wy,
127  + theta, xx, gg,
128  + col, head, iword, wv, wn, iprint )
129  implicit none
130  integer n, m, nsub, col, head, iword, iprint,
131  + ind(nsub), nbd(n)
132  double precision theta,
133  + l(n), u(n), x(n), d(n), xp(n), xx(n), gg(n),
134  + ws(n, m), wy(n, m),
135  + wv(2*m), wn(2*m, 2*m)
136 
137 c **********************************************************************
138 c
139 c This routine contains the major changes in the updated version.
140 c The changes are described in the accompanying paper
141 c
142 c Jose Luis Morales, Jorge Nocedal
143 c "Remark On Algorithm 788: L-BFGS-B: Fortran Subroutines for Large-Scale
144 c Bound Constrained Optimization". Decemmber 27, 2010.
145 c
146 c J.L. Morales Departamento de Matematicas,
147 c Instituto Tecnologico Autonomo de Mexico
148 c Mexico D.F.
149 c
150 c J, Nocedal Department of Electrical Engineering and
151 c Computer Science.
152 c Northwestern University. Evanston, IL. USA
153 c
154 c January 17, 2011
155 c
156 c **********************************************************************
157 c
158 c References:
159 c
160 c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
161 c memory algorithm for bound constrained optimization'',
162 c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
163 c
164 c * * *
165 c
166 c NEOS, November 1994. (Latest revision June 1996.)
167 c Optimization Technology Center.
168 c Argonne National Laboratory and Northwestern University.
169 c Written by
170 c Ciyou Zhu
171 c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
172 c
173 c
174 c ************
175 
176  integer pointr,m2,col2,ibd,jy,js,i,j,k
177  double precision alpha, xk, dk, temp1, temp2
178  double precision one,zero
179  parameter(one=1.0d0,zero=0.0d0)
180 c
181  double precision dd_p
182 
183  if (nsub .le. 0) return
184  if (iprint .ge. 99) write (6,1001)
185 
186 c Compute wv = W'Zd.
187 
188  pointr = head
189  do 20 i = 1, col
190  temp1 = zero
191  temp2 = zero
192  do 10 j = 1, nsub
193  k = ind(j)
194  temp1 = temp1 + wy(k,pointr)*d(j)
195  temp2 = temp2 + ws(k,pointr)*d(j)
196  10 continue
197  wv(i) = temp1
198  wv(col + i) = theta*temp2
199  pointr = mod(pointr,m) + 1
200  20 continue
201 
202 c Compute wv:=K^(-1)wv.
203 
204  m2 = 2*m
205  col2 = 2*col
206 
207  call dtrsm('l','u','t','n',col2,1,one,wn,m2,wv,col2)
208 
209  do 25 i = 1, col
210  wv(i) = -wv(i)
211  25 continue
212 
213  call dtrsm('l','u','n','n',col2,1,one,wn,m2,wv,col2)
214 
215 c Compute d = (1/theta)d + (1/theta**2)Z'W wv.
216 
217  pointr = head
218  do 40 jy = 1, col
219  js = col + jy
220  do 30 i = 1, nsub
221  k = ind(i)
222  d(i) = d(i) + wy(k,pointr)*wv(jy)/theta
223  + + ws(k,pointr)*wv(js)
224  30 continue
225  pointr = mod(pointr,m) + 1
226  40 continue
227 
228  call dscal( nsub, one/theta, d, 1 )
229 c
230 c-----------------------------------------------------------------
231 c Let us try the projection, d is the Newton direction
232 
233  iword = 0
234 
235  call dcopy ( n, x, 1, xp, 1 )
236 c
237  do 50 i=1, nsub
238  k = ind(i)
239  dk = d(i)
240  xk = x(k)
241  if ( nbd(k) .ne. 0 ) then
242 c
243  if ( nbd(k).eq.1 ) then ! lower bounds only
244  x(k) = max( l(k), xk + dk )
245  if ( x(k).eq.l(k) ) iword = 1
246  else
247 c
248  if ( nbd(k).eq.2 ) then ! upper and lower bounds
249  xk = max( l(k), xk + dk )
250  x(k) = min( u(k), xk )
251  if ( x(k).eq.l(k) .or. x(k).eq.u(k) ) iword = 1
252  else
253 c
254  if ( nbd(k).eq.3 ) then ! upper bounds only
255  x(k) = min( u(k), xk + dk )
256  if ( x(k).eq.u(k) ) iword = 1
257  end if
258  end if
259  end if
260 c
261  else ! free variables
262  x(k) = xk + dk
263  end if
264  50 continue
265 c
266  if ( iword.eq.0 ) then
267  go to 911
268  end if
269 c
270 c check sign of the directional derivative
271 c
272  dd_p = zero
273  do 55 i=1, n
274  dd_p = dd_p + (x(i) - xx(i))*gg(i)
275  55 continue
276  if ( dd_p .gt.zero ) then
277  call dcopy( n, xp, 1, x, 1 )
278  write(6,*) ' Positive dir derivative in projection '
279  write(6,*) ' Using the backtracking step '
280  else
281  go to 911
282  endif
283 c
284 c-----------------------------------------------------------------
285 c
286  alpha = one
287  temp1 = alpha
288  ibd = 0
289  do 60 i = 1, nsub
290  k = ind(i)
291  dk = d(i)
292  if (nbd(k) .ne. 0) then
293  if (dk .lt. zero .and. nbd(k) .le. 2) then
294  temp2 = l(k) - x(k)
295  if (temp2 .ge. zero) then
296  temp1 = zero
297  else if (dk*alpha .lt. temp2) then
298  temp1 = temp2/dk
299  endif
300  else if (dk .gt. zero .and. nbd(k) .ge. 2) then
301  temp2 = u(k) - x(k)
302  if (temp2 .le. zero) then
303  temp1 = zero
304  else if (dk*alpha .gt. temp2) then
305  temp1 = temp2/dk
306  endif
307  endif
308  if (temp1 .lt. alpha) then
309  alpha = temp1
310  ibd = i
311  endif
312  endif
313  60 continue
314 
315  if (alpha .lt. one) then
316  dk = d(ibd)
317  k = ind(ibd)
318  if (dk .gt. zero) then
319  x(k) = u(k)
320  d(ibd) = zero
321  else if (dk .lt. zero) then
322  x(k) = l(k)
323  d(ibd) = zero
324  endif
325  endif
326  do 70 i = 1, nsub
327  k = ind(i)
328  x(k) = x(k) + alpha*d(i)
329  70 continue
330 cccccc
331  911 continue
332 
333  if (iprint .ge. 99) write (6,1004)
334 
335  1001 format (/,'----------------SUBSM entered-----------------',/)
336  1004 format (/,'----------------exit SUBSM --------------------',/)
337 
338  return
339 
340  end
subsm
subroutine subsm(n, m, nsub, ind, l, u, nbd, x, d, xp, ws, wy, theta, xx, gg, col, head, iword, wv, wn, iprint)
Performs the subspace minimization.
Definition: subsm.f:129