c=====================================================================
      SUBROUTINE dlsx2(b,m1,z,icon)
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
c input
      integer m1,icon
      real*8 b(120),z(16)
c.....................................................................
c  Reference: Press, Flannery, Teukoisky, and Vetterling,
c       Numerical Recipes in C, Cambridge Univ. Press, (1991).
c.....................................................................
c const
      integer NP
      parameter( NP=16 )
c var
      integer indx(NP)
      real*8 a(NP,NP)
c begin
      l=0
      do 10 i=1, m1
         do 20 j=1, i
            l = l+1
            a(i,j) = b(l)
            a(j,i) = b(l)
 20      continue
 10   continue
      call ludcmp(a,m1,NP,indx,d,icon)
      call lubksb(a,m1,NP,indx,z)

      return
      end
c=====================================================================
      SUBROUTINE ludcmp(a,n,np,indx,d,icon)
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      INTEGER n,np,indx(n),NMAX
      integer icon
      REAL*8 d,a(np,np),TINY
      PARAMETER (NMAX=500,TINY=1.0e-20)
      INTEGER i,imax,j,k
      REAL*8 aamax,dum,sum,vv(NMAX)
      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
ce      if (aamax.eq.0.) pause 'singular matrix in ludcmp'
        if (aamax.eq.0.) then
           icon=20000
           return
        endif
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.)a(j,j)=TINY
        if(j.ne.n)then
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      END
c=====================================================================
      SUBROUTINE lubksb(a,n,np,indx,b)
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      INTEGER n,np,indx(n)
      REAL*8 a(np,np),b(n)
      INTEGER i,ii,j,ll
      REAL*8 sum
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      END
