***CHOLESKY(CROUTS) METHOD*************************************************
*  THE PROGRAM FOR SYSTEM OF LINEAR EQUATION BY LU DECOMPOSITION METHOD   *
*                         WITH COMPLETE PIVOTING                          *
***************************************************************************
*
      program LU
      implicit doubleprecision(a-h,o-z)
      integer p,r
      parameter (nd = 40)
      doubleprecision A(nd,nd),L(nd,nd),U(nd,nd),C(nd),S(nd),X(nd)
     $,Y(nd),ORD(nd),ORDC(nd)
      parameter (epsil = 1d-9)
c
c     segment for reading the equation
      print*, 'Enter number of equations, n:'
      read(*,*) n
      m = n + 1
      print*, 'Enter coefficient matrix row-wise:'
      read(*,*) ((A(i,j), j = 1,n),i = 1, n)
      write(6,*)
c
c     establish initial ordering in order vectors 
      do i = 1, n
         ORD(i) = i
         ORDC(i) = i
      enddo
c
c     segment for complete pivoting
      nn = n - 1
      do p = 1, n
         CALL PIVOT(A,n,ORD,ORDC,p)
c
c        LU decomposition
         do i = 1, n
            L(i,1) = A(i,1)
         enddo
         do j = 1, n
            U(1,j) = A(1,j)/L(1,1)
         enddo
*
         do j = 2, n
            do i = j, n
               sum = 0
               do k = 1, j-1
                  sum = sum + L(i,k) *U(k,j)
               enddo
               L(i,j) = A(i,j) - sum
            enddo
*
            U(j,j) = 1
            do i = j+1, n
               sum = 0
               do k = 1, j-1
                  sum = sum + L(j,k) *U(k,i)
               enddo
               if (abs(L(j,j)) > epsil) then
                  U(j,i) = (A(j,i) - sum) /L(j,j)
               endif
            enddo
         enddo
      enddo
c
c     checking for the singularity of coefficient matrix    
      do i = 1, n
         if (abs(L(i,i)) < epsil) then
            write(6,*)'-----------------------------------------------'
            write(6,*)'Coefficient matrix singular. No Solution exist.'
            write(6,*)'-----------------------------------------------'
            stop
         endif
      enddo
*
      icontinue = 1
10    if(icontinue == 1) then
         print*, 'Enter coefficient Vector:'
         read(*,*) (S(i), i = 1,n)
c
c        rearrange the elements of the S vector; C is used to hold them
         do i = 1, n
            j = ORD(i)
            C(i) = S(j)
         enddo
c
c        compute B vector
         C(1) = C(1) /L(1,1)
         do i = 2, n
            sum = 0
            do k = 1, i-1
               sum = sum + L(i,k) *C(k)
            enddo
            C(i) = (C(i) - sum) /L(i,i)
         enddo
c
c        augmenting with constant vector
         do i = 1, n
            U(i,m) = C(i)
         enddo
c
c        back substitution
         nn = n - 1
         X(n) = U(n,m)
         do i = nn, 1, -1
            sum = 0
            k = i + 1
            do j = k, n
               sum = sum + U(i,j) *X(j)
            enddo
            X(i) = U(i,m) - sum
         enddo
         do i = 1, n
            j = ORDC(i)
            Y(j) = X(i)
         enddo           
c
c        outputting the solution vector
         write(6,*)'-----------------------------------------------'
         write(6,*) '                       The solution Vector is:'
         write(6,*)
         write(6,101) (Y(i), i = 1,n)
         write(6,*)'-----------------------------------------------'
         write(6,*)
         write(6,*) ' More right-hand side vector ?  < 0 for NO / 1 for
     $YES >:'
         read(*,*) icontinue
         go to 10
      endif
      write(6,*)'---------------------------------------------------'
*
101   format('',29x,f12.4)
*
      stop
      end
*
***COMPLETE PIVOTING (Loop in the calling program)********************
*
      subroutine PIVOT(A,n,ORD,ORDC,i)
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),ORD(nd),ORDC(nd)
      integer col,row,p,r
c
c     complete pivoting
      row = i
      col = i
      do p = i, n
      do r = i, n
         if (dabs(A(row,col)) < dabs(A(p,r))) then
            row = p
            col = r
         endif
      enddo
      enddo
*
      if (col.ne.i) then
         do ii = 1, n
            temp = A(ii,i)
            A(ii,i) = A(ii,col)
            A(ii,col) = temp
         enddo
         t = ORDC(i)
         ORDC(i) = ORDC(col)
         ORDC(col) = t
      endif
*
      m = n + 1 
      if (row.ne.i) then
         do jj = 1, m
            temp = A(i,jj)
            A(i,jj) = A(row,jj)
            A(row,jj) = temp
         enddo
         t = ORD(i)
         ORD(i) = ORD(row)
         ORD(row) = t
      endif
*     
      return
      end