Gaussian Elimination (no pivoting)

      program gauss
      implicit real*4 (a-h,o-z)
      parameter(n=10)
      parameter(tol=1.0e-20)
      dimension a(n,n+1), u(n), b(n)
      dimension c(n,n)
 
c     the following algorithm uses gaussian elimination to solve a
c     system of n linear equations of the form Au = b.
c     it is assumed that no pivoting is required
c     input:
c     integer n - - order of the system
c     array a(i,j) - - n by n coefficient matrix
c     array b(i) - - right side of system
c     output
c     solution vector u or a message that a may be singular
c
c  step1 define linear system for hilbert matrix example
      xi = 0.0
      do 10 i = 1,n
      xi = xi + 1.0
      sum = 0.0
      xj = 0.0
      do 11 j = 1,n
      xj = xj + 1.0
      c(i,j) = 1.0/(i + j - 1) 
      a(i,j) = 1.0/(xi+xj-1.0)
      sum = sum + a(i,j)
11    continue
      b(i) = sum
      a(i,n+1) = sum
c  the following print statement is for investigative 
c  purposes only -- take it out after you see the effects
c  of promoting integers to floats in constructing the
c  hilbert matrix  -- also remove the matrix c when finished   
c  investigating
      print*,'i=',i,'  a(i,i)= ',a(i,i),'         all floats'
      print*,'i=',i,'  c(i,i)= ',c(i,i),'         mixed mode'

10    continue
c  step2   begin forward substitution
      do 1 i = 1,n-1
      if(abs(a(i,i)).lt.tol) call mesg()
c  step 3   calculate multiplier
      do 2 k = i+1,n
      r = a(k,i)/a(i,i)
c  step 4  subtract r times eq i fron eq k
      do 3 j = i,n+1
      a(k,j) = a(k,j) - r*a(i,j)
 3    continue
 2    continue
 1    continue
c  step 5  check for singular matrix
      if(abs(a(n,n)).lt.tol) call mesg()
c  step 6  initialize backward substitution
      u(n) = a(n,n+1)/a(n,n)
c  step 7  back substitute and store solution in u
      do 4 i = n-1,1,-1
      sum = 0
      do 5 j = i+1,n
      sum = sum + a(i,j)*u(j)
 5    continue

      u(i) = (a(i,n+1) - sum)/a(i,i)
 4    continue
c  step 8  output solution
      print*,'Solution Vector'
      do 6 i = 1,n
      print*,'i=',i,'  u(i)= ',u(i)
 6    continue
c  step 9  output residuals and errors
      print*,'Residual and Error Vectors'
      xi = 0.0
      do 7 i = 1,n
      xi = xi + 1.0
      sum = 0.0
      xk = 0.0
      do 77 k = 1,n
      xk = xk + 1.0
      hik = 1.0/(xi+xk-1.0)
      sum = sum + hik*u(k)
77    continue
      resid = b(i) - sum
      error = u(i) - 1.0
      print*,'i=',i,'   resid(i)= ',resid
      print*,'i=',i,'   error(i)= ',error
7     continue
      end
 
      subroutine mesg()
      print*,'matrix may be singular or a pivot is required'
      stop
      end