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
|
The library code originally uses LINPACK (from an embedded
copy). Since LINPACK has largely been superseded by LAPACK, this patch
replaces calls to the former with equivalent calls to the
latter. Specifically, dpofa is replaced by dpotrf, and dtrsl is
replaced by dtrtrs.
--- a/lbfgsb.f
+++ b/lbfgsb.f
@@ -1185,7 +1185,7 @@
p(i2) = v(i2) + sum
20 continue
c Solve the triangular system
- 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 +1197,7 @@
c [ 0 J' ] [ p2 ] [ p2 ].
c solve J^Tp2=p2.
- 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 +2135,7 @@
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)
+ call dpotrf('U', col, wn, m2, info)
if (info .ne. 0) then
info = -1
return
@@ -2143,7 +2143,7 @@
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)
+ 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 +2158,7 @@
c Cholesky factorization of (2,2) block of wn.
- 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 +2227,7 @@
c Cholesky factorize T to J*J' with
c J' stored in the upper triangle of wt.
- call dpofa(wt,m,col,info)
+ call dpotrf('U', col, wt, m, info)
if (info .ne. 0) then
info = -3
endif
@@ -3208,12 +3208,12 @@
m2 = 2*m
col2 = 2*col
- 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)
+ 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.
|
