Patch against upstream lbfgsb (needs packaging) that replaces its bundled
LINPACK with LAPACK.
--- a/lbfgsb.f
+++ b/lbfgsb.f
@@ -1185,7 +1185,8 @@
p(i2) = v(i2) + sum
20 continue
c Solve the triangular system
- call dtrsl(wt,m,col,p(col+1),11,info)
+c call dtrsl(wt,m,col,p(col+1),11,info)
+ call dtrtrs('U', 'T', 'N', col, 1, wt, m, p(col+1), col, info)
if (info .ne. 0) return
c solve D^(1/2)p1=v1.
@@ -1197,7 +1198,8 @@
c [ 0 J' ] [ p2 ] [ p2 ].
c solve J^Tp2=p2.
- call dtrsl(wt,m,col,p(col+1),01,info)
+c call dtrsl(wt,m,col,p(col+1),01,info)
+ call dtrtrs('U', 'N', 'N', col, 1, wt, m, p(col+1), col, info)
if (info .ne. 0) return
c compute p1=-D^(-1/2)(p1-D^(-1/2)L'p2)
@@ -2135,7 +2137,8 @@
c first Cholesky factor (1,1) block of wn to get LL'
c with L' stored in the upper triangle of wn.
- call dpofa(wn,m2,col,info)
+c call dpofa(wn,m2,col,info)
+ call dpotrf('U', col, wn, m2, info)
if (info .ne. 0) then
info = -1
return
@@ -2143,7 +2146,8 @@
c then form L^-1(-L_a'+R_z') in the (1,2) block.
col2 = 2*col
do 71 js = col+1 ,col2
- call dtrsl(wn,m2,col,wn(1,js),11,info)
+c call dtrsl(wn,m2,col,wn(1,js),11,info)
+ call dtrtrs('U', 'T', 'N', col, 1, wn, m2, wn(1,js), col, info)
71 continue
c Form S'AA'S*theta + (L^-1(-L_a'+R_z'))'L^-1(-L_a'+R_z') in the
@@ -2158,7 +2162,8 @@
c Cholesky factorization of (2,2) block of wn.
- call dpofa(wn(col+1,col+1),m2,col,info)
+c call dpofa(wn(col+1,col+1),m2,col,info)
+ call dpotrf('U', col, wn(col+1,col+1), m2, info)
if (info .ne. 0) then
info = -2
return
@@ -2227,7 +2232,8 @@
c Cholesky factorize T to J*J' with
c J' stored in the upper triangle of wt.
- call dpofa(wt,m,col,info)
+c call dpofa(wt,m,col,info)
+ call dpotrf('U', col, wt, m, info)
if (info .ne. 0) then
info = -3
endif
@@ -3208,12 +3214,14 @@
m2 = 2*m
col2 = 2*col
- call dtrsl(wn,m2,col2,wv,11,info)
+c call dtrsl(wn,m2,col2,wv,11,info)
+ call dtrtrs('U', 'T', 'N', col2, 1, wn, m2, wv, col2, info)
if (info .ne. 0) return
do 25 i = 1, col
wv(i) = -wv(i)
25 continue
- call dtrsl(wn,m2,col2,wv,01,info)
+c call dtrsl(wn,m2,col2,wv,01,info)
+ call dtrtrs('U', 'N', 'N', col2, 1, wn, m2, wv, col2, info)
if (info .ne. 0) return
c Compute d = (1/theta)d + (1/theta**2)Z'W wv.
