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
|
From: Gard Spreemann <gspr@nonempty.org>
Date: Thu, 30 Jul 2020 09:57:45 +0200
Subject: replace-linpack-with-lapack
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.
---
lbfgsb.f | 16 ++++++++--------
1 file changed, 8 insertions(+), 8 deletions(-)
diff --git a/lbfgsb.f b/lbfgsb.f
index 9c9e7d9..54fd1f9 100644
--- a/lbfgsb.f
+++ b/lbfgsb.f
@@ -1185,7 +1185,7 @@ c solve Jp2=v2+LD^(-1)v1.
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 PART II: solve [ -D^(1/2) D^(-1/2)*L' ] [ p1 ] = [ p1 ]
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 [(-L_a +R_z)L'^-1 S'AA'S*theta ]
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 with L' stored in the upper triangle of wn.
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 upper triangle of (2,2) block of wn.
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 store T in the upper triangle of the array wt.
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 @@ c Compute wv:=K^(-1)wv.
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.
|
