L-BFGS-B  3.0
Large-scale Bound-constrained Optimization
formk.f
Go to the documentation of this file.
1 c> \file formk.f
2 
3 c> \brief Forms the LEL^T factorization of the indefinite matrix K.
4 c>
5 c> This subroutine forms the LEL^T factorization of the indefinite matrix
6 c>
7 c> K = [-D -Y'ZZ'Y/theta L_a'-R_z' ]
8 c> [L_a -R_z theta*S'AA'S ]
9 c> where E = [-I 0]
10 c> [ 0 I]
11 c>
12 c> The matrix K can be shown to be equal to the matrix M^[-1]N
13 c> occurring in section 5.1 of [1], as well as to the matrix
14 c> Mbar^[-1] Nbar in section 5.3.
15 c>
16 c> @param n On entry n is the dimension of the problem.<br/>
17 c> On exit n is unchanged.
18 c>
19 c> @param nsub On entry nsub is the number of subspace variables in free set.<br/>
20 c> On exit nsub is not changed.
21 c>
22 c> @param ind On entry ind specifies the indices of subspace variables.<br/>
23 c> On exit ind is unchanged.
24 c>
25 c> @param nenter On entry nenter is the number of variables entering the
26 c> free set.<br/>
27 c> On exit nenter is unchanged.
28 c>
29 c> @param ileave On entry indx2(ileave),...,indx2(n) are the variables leaving
30 c> the free set.<br/>
31 c> On exit ileave is unchanged.
32 c>
33 c> @param indx2 On entry indx2(1),...,indx2(nenter) are the variables entering
34 c> the free set, while indx2(ileave),...,indx2(n) are the
35 c> variables leaving the free set.<br/>
36 c> On exit indx2 is unchanged.
37 c>
38 c> @param iupdat On entry iupdat is the total number of BFGS updates made so far.<br/>
39 c> On exit iupdat is unchanged.
40 c>
41 c> @param updatd On entry 'updatd' is true if the L-BFGS matrix is updated.<br/>
42 c> On exit 'updatd' is unchanged.
43 c>
44 c> @param wn On entry wn is unspecified.<br/>
45 c> On exit the upper triangle of wn stores the LEL^T factorization
46 c> of the 2*col x 2*col indefinite matrix
47 c> [-D -Y'ZZ'Y/theta L_a'-R_z' ]
48 c> [L_a -R_z theta*S'AA'S ]
49 c>
50 c> @param wn1 On entry wn1 stores the lower triangular part of
51 c> [Y' ZZ'Y L_a'+R_z']
52 c> [L_a+R_z S'AA'S ]
53 c> in the previous iteration.<br/>
54 c> On exit wn1 stores the corresponding updated matrices.<br/>
55 c> The purpose of wn1 is just to store these inner products
56 c> so they can be easily updated and inserted into wn.
57 c>
58 c> @param m On entry m is the maximum number of variable metric corrections
59 c> used to define the limited memory matrix.<br/>
60 c> On exit m is unchanged.
61 c>
62 c> @param ws On entry this stores S, a set of s-vectors, that defines the
63 c> limited memory BFGS matrix.<br/>
64 c> On exit this array is unchanged.
65 c>
66 c> @param wy On entry this stores Y, a set of y-vectors, that defines the
67 c> limited memory BFGS matrix.<br/>
68 c> On exit this array is unchanged.
69 c>
70 c> @param sy On entry this stores S'Y, that defines the
71 c> limited memory BFGS matrix.<br/>
72 c> On exit this array is unchanged.
73 c>
74 c> @param theta On entry theta is the scaling factor specifying B_0 = theta I.<br/>
75 c> On exit theta is unchanged.
76 c>
77 c> @param col On entry col is the actual number of variable metric
78 c> corrections stored so far.<br/>
79 c> On exit col is unchanged.
80 c>
81 c> @param head On entry head is the location of the first s-vector (or y-vector)
82 c> in S (or Y).<br/>
83 c> On exit col is unchanged.
84 c>
85 c> @param info On entry info is unspecified.<br/>
86 c> On exit info<ul><li>= 0 for normal return;</li>
87 c> <li>= -1 when the 1st Cholesky factorization failed;</li>
88 c> <li>= -2 when the 2st Cholesky factorization failed.</li></ul>
89  subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat,
90  + updatd, wn, wn1, m, ws, wy, sy, theta, col,
91  + head, info)
92 
93  integer n, nsub, m, col, head, nenter, ileave, iupdat,
94  + info, ind(n), indx2(n)
95  double precision theta, wn(2*m, 2*m), wn1(2*m, 2*m),
96  + ws(n, m), wy(n, m), sy(m, m)
97  logical updatd
98 c
99 c References:
100 c [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited
101 c memory algorithm for bound constrained optimization'',
102 c SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208.
103 c
104 c [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a
105 c limited memory FORTRAN code for solving bound constrained
106 c optimization problems'', Tech. Report, NAM-11, EECS Department,
107 c Northwestern University, 1994.
108 c
109 c (Postscript files of these papers are available via anonymous
110 c ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.)
111 c
112 c * * *
113 c
114 c NEOS, November 1994. (Latest revision June 1996.)
115 c Optimization Technology Center.
116 c Argonne National Laboratory and Northwestern University.
117 c Written by
118 c Ciyou Zhu
119 c in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.
120 c
121 c
122 c ************
123 
124  integer m2,ipntr,jpntr,iy,is,jy,js,is1,js1,k1,i,k,
125  + col2,pbegin,pend,dbegin,dend,upcl
126  double precision ddot,temp1,temp2,temp3,temp4
127  double precision one,zero
128  parameter(one=1.0d0,zero=0.0d0)
129 
130 c Form the lower triangular part of
131 c WN1 = [Y' ZZ'Y L_a'+R_z']
132 c [L_a+R_z S'AA'S ]
133 c where L_a is the strictly lower triangular part of S'AA'Y
134 c R_z is the upper triangular part of S'ZZ'Y.
135 
136  if (updatd) then
137  if (iupdat .gt. m) then
138 c shift old part of WN1.
139  do 10 jy = 1, m - 1
140  js = m + jy
141  call dcopy(m-jy,wn1(jy+1,jy+1),1,wn1(jy,jy),1)
142  call dcopy(m-jy,wn1(js+1,js+1),1,wn1(js,js),1)
143  call dcopy(m-1,wn1(m+2,jy+1),1,wn1(m+1,jy),1)
144  10 continue
145  endif
146 
147 c put new rows in blocks (1,1), (2,1) and (2,2).
148  pbegin = 1
149  pend = nsub
150  dbegin = nsub + 1
151  dend = n
152  iy = col
153  is = m + col
154  ipntr = head + col - 1
155  if (ipntr .gt. m) ipntr = ipntr - m
156  jpntr = head
157  do 20 jy = 1, col
158  js = m + jy
159  temp1 = zero
160  temp2 = zero
161  temp3 = zero
162 c compute element jy of row 'col' of Y'ZZ'Y
163  do 15 k = pbegin, pend
164  k1 = ind(k)
165  temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
166  15 continue
167 c compute elements jy of row 'col' of L_a and S'AA'S
168  do 16 k = dbegin, dend
169  k1 = ind(k)
170  temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
171  temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
172  16 continue
173  wn1(iy,jy) = temp1
174  wn1(is,js) = temp2
175  wn1(is,jy) = temp3
176  jpntr = mod(jpntr,m) + 1
177  20 continue
178 
179 c put new column in block (2,1).
180  jy = col
181  jpntr = head + col - 1
182  if (jpntr .gt. m) jpntr = jpntr - m
183  ipntr = head
184  do 30 i = 1, col
185  is = m + i
186  temp3 = zero
187 c compute element i of column 'col' of R_z
188  do 25 k = pbegin, pend
189  k1 = ind(k)
190  temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
191  25 continue
192  ipntr = mod(ipntr,m) + 1
193  wn1(is,jy) = temp3
194  30 continue
195  upcl = col - 1
196  else
197  upcl = col
198  endif
199 
200 c modify the old parts in blocks (1,1) and (2,2) due to changes
201 c in the set of free variables.
202  ipntr = head
203  do 45 iy = 1, upcl
204  is = m + iy
205  jpntr = head
206  do 40 jy = 1, iy
207  js = m + jy
208  temp1 = zero
209  temp2 = zero
210  temp3 = zero
211  temp4 = zero
212  do 35 k = 1, nenter
213  k1 = indx2(k)
214  temp1 = temp1 + wy(k1,ipntr)*wy(k1,jpntr)
215  temp2 = temp2 + ws(k1,ipntr)*ws(k1,jpntr)
216  35 continue
217  do 36 k = ileave, n
218  k1 = indx2(k)
219  temp3 = temp3 + wy(k1,ipntr)*wy(k1,jpntr)
220  temp4 = temp4 + ws(k1,ipntr)*ws(k1,jpntr)
221  36 continue
222  wn1(iy,jy) = wn1(iy,jy) + temp1 - temp3
223  wn1(is,js) = wn1(is,js) - temp2 + temp4
224  jpntr = mod(jpntr,m) + 1
225  40 continue
226  ipntr = mod(ipntr,m) + 1
227  45 continue
228 
229 c modify the old parts in block (2,1).
230  ipntr = head
231  do 60 is = m + 1, m + upcl
232  jpntr = head
233  do 55 jy = 1, upcl
234  temp1 = zero
235  temp3 = zero
236  do 50 k = 1, nenter
237  k1 = indx2(k)
238  temp1 = temp1 + ws(k1,ipntr)*wy(k1,jpntr)
239  50 continue
240  do 51 k = ileave, n
241  k1 = indx2(k)
242  temp3 = temp3 + ws(k1,ipntr)*wy(k1,jpntr)
243  51 continue
244  if (is .le. jy + m) then
245  wn1(is,jy) = wn1(is,jy) + temp1 - temp3
246  else
247  wn1(is,jy) = wn1(is,jy) - temp1 + temp3
248  endif
249  jpntr = mod(jpntr,m) + 1
250  55 continue
251  ipntr = mod(ipntr,m) + 1
252  60 continue
253 
254 c Form the upper triangle of WN = [D+Y' ZZ'Y/theta -L_a'+R_z' ]
255 c [-L_a +R_z S'AA'S*theta]
256 
257  m2 = 2*m
258  do 70 iy = 1, col
259  is = col + iy
260  is1 = m + iy
261  do 65 jy = 1, iy
262  js = col + jy
263  js1 = m + jy
264  wn(jy,iy) = wn1(iy,jy)/theta
265  wn(js,is) = wn1(is1,js1)*theta
266  65 continue
267  do 66 jy = 1, iy - 1
268  wn(jy,is) = -wn1(is1,jy)
269  66 continue
270  do 67 jy = iy, col
271  wn(jy,is) = wn1(is1,jy)
272  67 continue
273  wn(iy,iy) = wn(iy,iy) + sy(iy,iy)
274  70 continue
275 
276 c Form the upper triangle of WN= [ LL' L^-1(-L_a'+R_z')]
277 c [(-L_a +R_z)L'^-1 S'AA'S*theta ]
278 
279 c first Cholesky factor (1,1) block of wn to get LL'
280 c with L' stored in the upper triangle of wn.
281  !call dpofa(wn,m2,col,info)
282  call dpotrf('U',col,wn,m2,info)
283 
284  if (info .ne. 0) then
285  info = -1
286  return
287  endif
288 c then form L^-1(-L_a'+R_z') in the (1,2) block.
289  col2 = 2*col
290  do 71 js = col+1 ,col2
291  !call dtrsl(wn,m2,col,wn(1,js),11,info)
292 
293  ! TODO: can combine this loop into a single dtrsm call?
294  call dtrsm('l','u','t','n',col,1,one,wn,m2,wn(1,js),col)
295  71 continue
296 
297 c Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
298 c upper triangle of (2,2) block of wn.
299 
300 
301  do 72 is = col+1, col2
302  do 74 js = is, col2
303  wn(is,js) = wn(is,js) + ddot(col,wn(1,is),1,wn(1,js),1)
304  74 continue
305  72 continue
306 
307 c Cholesky factorization of (2,2) block of wn.
308 
309  !call dpofa(wn(col+1,col+1),m2,col,info)
310  call dpotrf('U',col,wn(col+1,col+1),m2,info)
311 
312  if (info .ne. 0) then
313  info = -2
314  return
315  endif
316 
317  return
318 
319  end
formk
subroutine formk(n, nsub, ind, nenter, ileave, indx2, iupdat, updatd, wn, wn1, m, ws, wy, sy, theta, col, head, info)
Forms the LEL^T factorization of the indefinite matrix K.
Definition: formk.f:92