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.