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.