      program lutest
      parameter (np=48)
      dimension a(np,np),b(np,np),indx(np)
c
      open(1,file='matrix.prn')
      do i=1,np
       read(1,*)(a(i,j),j=1,np)
      enddo
      close(1)
      open(1,file='in_mat.dat')
      do i=1,np
       write(1,'(24(f7.2,1x))')(a(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f7.2,1x))')(a(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      call ludcmp(a,np,np,indx,d)
c
      do i=1,np
       write(2,*)indx(i)
      enddo
c
      open(1,file='out_mat.dat')
      do i=1,np
       write(1,'(24(f7.2,1x))')(a(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f7.2,1x))')(a(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      do i=1,np
       do j=1,np
        b(i,j)=0.
       enddo
      enddo
      do i=1,np
       do j=1,i
        b(np-i+1,np-j+1)=a(np-i+1,np-j+1)
       enddo
      enddo
      open(1,file='lu_upper.dat')
      do i=1,np
       write(1,'(24(f7.2,1x))')(b(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f7.2,1x))')(b(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      do i=1,np
       do j=1,np
        b(i,j)=0.
       enddo
       b(i,i)=1.
      enddo
      do i=2,np
       do j=1,i-1
        b(i,j)=a(i,j)
       enddo
      enddo
      open(1,file='lu_lower.dat')
      do i=1,np
       write(1,'(24(f7.2,1x))')(b(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f7.2,1x))')(b(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      end

      SUBROUTINE LUDCMP(A,N,NP,INDX,D)
      PARAMETER (NMAX=100,TINY=1.0E-20)
      DIMENSION A(NP,NP),INDX(N),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
        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
        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


