      program lutest
      parameter (np=48)
      real*8 a,b,v
      dimension a(np,np),b(np,np),indx(np),v(np*(np+1)/2)
 
      open(1,file='matrix.prn')
      do i=1,np
       read(1,*)(a(i,j),j=1,np)
      enddo
      close(1)
 
c     a(1,1)=2.
c     a(2,1)=1.
c     a(3,1)=1.
c     a(1,2)=1.
c     a(2,2)=3.
c     a(3,2)=2.
c     a(1,3)=1.
c     a(2,3)=2.
c     a(3,3)=4.

      open(1,file='in_mat.dat')
      do i=1,np
       write(1,'(24(f8.3,1x))')(a(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f8.3,1x))')(a(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      jj=1
      do i=1,np
       do j=i,np
        v(jj)=a(i,j)
        jj=jj+1
       enddo
      enddo
c
      ia=1
      do ic=1,np
       in=np-ic+1
       call cholesky(v(ia),in)
       ia=ia+in
      enddo
c
      jj=1
      do i=1,np
       do j=i,np
        a(i,j)=v(jj)
        a(j,i)=a(i,j)
        jj=jj+1
       enddo
      enddo
c      
      open(1,file='out_mat.dat')
      do i=1,np
       write(1,'(24(f8.3,1x))')(a(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f8.3,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(f8.3,1x))')(b(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f8.3,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
      enddo
      do i=1,np
       do j=1,i
        b(i,j)=a(i,j)
       enddo
      enddo
      open(1,file='lu_lower.dat')
      do i=1,np
       write(1,'(24(f8.3,1x))')(b(i,j),j=1,np/2)
      enddo
      write(1,*)
      do i=1,np
       write(1,'(24(f8.3,1x))')(b(i,j),j=np/2+1,np)
      enddo
      close(1)
c
      end

      subroutine cholesky(a,n)
c
      real*8 a,b
      dimension a(*),b(48,48)
c
      a(1)=sqrt(a(1))
      do i=2,n
       a(i)=a(i)/a(1)
      enddo
c
      ii=2
      jj=2
      do i=n+1,n*(n+1)/2
       a(i)=a(i)-a(ii)*a(jj)
       jj=jj+1
       if(jj.gt.n)then
        ii=ii+1
        jj=ii
       endif
      enddo

      jj=1
      do i=1,n
       do j=i,n
        b(i,j)=a(jj)
        b(j,i)=b(i,j)
        jj=jj+1
       enddo
       write(12,'(48(f7.2,1x))')(b(i,j),j=1,n)
      enddo
      write(12,*)
c
      return
      end

