***NAIVE*****************************************************************
*   THE PROGRAM FOR THE SOLUTION OF SYSTEM OF LINEAR EQUATION BY NAIVE  *
*                      GAUSS ELIMINATION METHOD                         *
*************************************************************************
*
      program NGAUSS
      implicit doubleprecision(a-h,o-z)
      parameter (nd = 40)
      doubleprecision A(nd,nd),X(nd)
      parameter (epsil = 1d-6)
c
c     segment for reading the equation
      print*, 'Enter number of equations, n:'
      read(*,*) n
      m = n + 1
      print*, 'Enter augmented coefficient matrix row-wise:'
      read(*,*) ((A(i,j), j = 1,m), i = 1,n)
      write(6,*)'------------------------------------------------------'
c
c     triangularization by eliminating the variables
      nn = n - 1
      do k = 1, nn
         kk = k + 1
         if (dabs(A(k,k)) <= epsil) then
            write(6,*) '*During triangularization diagonal element becom
     $es extremely small. Hence the Naive Gauss elimination is not suita
     $ble. Try Gauss elimination with partial pivoting.*'
            stop      ! terminating the program
         else  
            do i = kk, n
               qt = A(i,k)/A(k,k)
               do j = k, m
                  A(i,j) = A(i,j) - qt *A(k,j)
               enddo
            enddo
         endif
      enddo
c
c     checking for singularity of matrix
      do i = 1, n
         if (dabs(A(i,i)) == 0d0) then
            write(6,*)'-----------------------------------------------'
            write(6,*)'Coefficient matrix singular. No Solution exist.'
            write(6,*)'-----------------------------------------------'
            stop      ! terminating the program
         endif
      enddo
c
c     back substitution
      X(n) = A(n,m)/A(n,n)
      do i = nn, 1, -1
         sum = 0
         k = i + 1
         do j = k, n
            sum = sum + A(i,j) *X(j)
         enddo
         X(i) = (A(i,m) - sum)/A(i,i)
      enddo
c
c     outputting the solution vector
      write(6,*) '                      The solution Vector is:'
      write(6,*)
      write(6,101) (X(i), i = 1,n)
101   format('',30x,f12.4)
      write(6,*)'-----------------------------------------------------'
c
c     check the value of the determinant
      d = 1
      do i = 1, n
         d = d *A(i,i)
      enddo
      write(6,*) '      The determinant of coefficient Matrix is  =',d
      write(6,*)'-----------------------------------------------------'
*
      stop
      end