# linpack.f Source File

L-BFGS-B submodule

## Source Code

!! L-BFGS-B submodule
c
c  L-BFGS-B is released under the “New BSD License” (aka “Modified BSD License”
c  or “3-clause license”)
c  Please read attached file License.txt
c
subroutine dpofa(a,lda,n,info)
integer lda,n,info
double precision a(lda,*)
c
c     dpofa factors a double precision symmetric positive definite
c     matrix.
c
c     dpofa is usually called by dpoco, but it can be called
c     directly with a saving in time if  rcond  is not needed.
c     (time for dpoco) = (1 + 18/n)*(time for dpofa) .
c
c     on entry
c
c        a       double precision(lda, n)
c                the symmetric matrix to be factored.  only the
c                diagonal and upper triangle are used.
c
c        lda     integer
c                the leading dimension of the array  a .
c
c        n       integer
c                the order of the matrix  a .
c
c     on return
c
c        a       an upper triangular matrix  r  so that  a = trans(r)*r
c                where  trans(r)  is the transpose.
c                the strict lower triangle is unaltered.
c                if  info .ne. 0 , the factorization is not complete.
c
c        info    integer
c                = 0  for normal return.
c                = k  signals an error condition.  the leading minor
c                     of order  k  is not positive definite.
c
c     linpack.  this version dated 08/14/78 .
c     cleve moler, university of new mexico, argonne national lab.
c
c     subroutines and functions
c
c     blas ddot
c     fortran sqrt
c
c     internal variables
c
double precision ddot,t
double precision s
integer j,jm1,k
c     begin block with ...exits to 40
c
c
do 30 j = 1, n
info = j
s = 0.0d0
jm1 = j - 1
if (jm1 .lt. 1) go to 20
do 10 k = 1, jm1
t = a(k,j) - ddot(k-1,a(1,k),1,a(1,j),1)
t = t/a(k,k)
a(k,j) = t
s = s + t*t
10       continue
20       continue
s = a(j,j) - s
c     ......exit
if (s .le. 0.0d0) go to 40
a(j,j) = sqrt(s)
30    continue
info = 0
40 continue
return
end

c====================== The end of dpofa ===============================

subroutine dtrsl(t,ldt,n,b,job,info)
integer ldt,n,job,info
double precision t(ldt,*),b(*)
c
c
c     dtrsl solves systems of the form
c
c                   t * x = b
c     or
c                   trans(t) * x = b
c
c     where t is a triangular matrix of order n. here trans(t)
c     denotes the transpose of the matrix t.
c
c     on entry
c
c         t         double precision(ldt,n)
c                   t contains the matrix of the system. the zero
c                   elements of the matrix are not referenced, and
c                   the corresponding elements of the array can be
c                   used to store other information.
c
c         ldt       integer
c                   ldt is the leading dimension of the array t.
c
c         n         integer
c                   n is the order of the system.
c
c         b         double precision(n).
c                   b contains the right hand side of the system.
c
c         job       integer
c                   job specifies what kind of system is to be solved.
c                   if job is
c
c                        00   solve t*x=b, t lower triangular,
c                        01   solve t*x=b, t upper triangular,
c                        10   solve trans(t)*x=b, t lower triangular,
c                        11   solve trans(t)*x=b, t upper triangular.
c
c     on return
c
c         b         b contains the solution, if info .eq. 0.
c                   otherwise b is unaltered.
c
c         info      integer
c                   info contains zero if the system is nonsingular.
c                   otherwise info contains the index of
c                   the first zero diagonal element of t.
c
c     linpack. this version dated 08/14/78 .
c     g. w. stewart, university of maryland, argonne national lab.
c
c     subroutines and functions
c
c     blas daxpy,ddot
c     fortran mod
c
c     internal variables
c
double precision ddot,temp
integer case,j,jj
c
c     begin block permitting ...exits to 150
c
c        check for zero diagonal elements.
c
do 10 info = 1, n
c     ......exit
if (t(info,info) .eq. 0.0d0) go to 150
10    continue
info = 0
c
c        determine the task and go to it.
c
case = 1
if (mod(job,10) .ne. 0) case = 2
if (mod(job,100)/10 .ne. 0) case = case + 2
go to (20,50,80,110), case
c
c        solve t*x=b for t lower triangular
c
20    continue
b(1) = b(1)/t(1,1)
if (n .lt. 2) go to 40
do 30 j = 2, n
temp = -b(j-1)
call daxpy(n-j+1,temp,t(j,j-1),1,b(j),1)
b(j) = b(j)/t(j,j)
30       continue
40       continue
go to 140
c
c        solve t*x=b for t upper triangular.
c
50    continue
b(n) = b(n)/t(n,n)
if (n .lt. 2) go to 70
do 60 jj = 2, n
j = n - jj + 1
temp = -b(j+1)
call daxpy(j,temp,t(1,j+1),1,b(1),1)
b(j) = b(j)/t(j,j)
60       continue
70       continue
go to 140
c
c        solve trans(t)*x=b for t lower triangular.
c
80    continue
b(n) = b(n)/t(n,n)
if (n .lt. 2) go to 100
do 90 jj = 2, n
j = n - jj + 1
b(j) = b(j) - ddot(jj-1,t(j+1,j),1,b(j+1),1)
b(j) = b(j)/t(j,j)
90       continue
100       continue
go to 140
c
c        solve trans(t)*x=b for t upper triangular.
c
110    continue
b(1) = b(1)/t(1,1)
if (n .lt. 2) go to 130
do 120 j = 2, n
b(j) = b(j) - ddot(j-1,t(1,j),1,b(1),1)
b(j) = b(j)/t(j,j)
120       continue
130       continue
140    continue
150 continue
return
end

c====================== The end of dtrsl ===============================