      program LUbackSUB
      integer n1,m1,L1
      integer NP,MP,SR
      INTEGER i,j,k
      INTEGER fp
      PARAMETER(NP=50,MP=50)
      INTEGER indx(NP)
      real ma1(NP,MP),TA1(NP,MP),XMM(NP,MP), b(NP),x(NP)
      real AL(NP,MP),AU(NP,MP),d

      SR=1
      open (unit=8,file="linsys1.dat")
      fp=8

      read (fp,*) 
      read (fp,*) n1,m1   
      call readmatrix(MA1,n1,m1,fp)

      read (fp,*) 
      read (fp,*) L1 
      call readvektor(b,L1,fp)
      do 13 i=1,n1
         do 11 j=1,m1
            TA1(i,j)=MA1(i,j)
 11      continue

 13   continue

      do  k=1,L1
         x(k)=b(k)
      End Do

      write(*,*) ' Matrix A : '
      call displaymatrix(TA1,N1,M1)
      call ludcmp(MA1,N1,NP,indx,d)
      print *, "Composed LU "
      call displaymatrix(MA1,N1,M1)
c     solve equation
      call lubksb(MA1,n1,NP,indx,x)
C     test results with original matrix
      write(*,*) 'Right-hand side vector:' 
      call displayvektor(x,n1)
      
      end program

      subroutine getLU(am,al,au,n)
      integer n,k,l
      real am(50,50),al(50,50),au(50,50)
        do  k=1,n
         do  l=1,n
            if (l.gt.k) then
               au(k,l)=am(k,l)
               al(k,l)=0.0
            else if (l.lt.k) then
               au(k,l)=0.0
               al(k,l)=am(k,l)
            else
               au(k,l)=am(k,l)
               al(k,l)=1.0
            endif
         end do
      end do
      end subroutine

      subroutine genindex(indx,jndx,n)
      integer k,dum,n
      integer indx(50),jndx(50)
      do  k=1,n
         jndx(k)=k
      End do
      do 17 k=1,n
        dum=jndx(indx(k))
        jndx(indx(k))=jndx(k)
        jndx(k)=dum
17    continue
      end subroutine

      subroutine permuteMat(a,aout,n,indx)
      integer n,k,j,l
      integer indx(50),jndx(50)
      real a(50,50), aout(50,50)
      call genindex(indx,jndx,n)
      do 19 k=1,n
         do 18 j=1,n
            if (jndx(j).eq.k) then
               do l=1,n
                  aout(k,l)=a(j,l)
               End Do
            endif
 18      continue
 19   continue
      end subroutine
