Program lintest
       implicit none
!  Test program for linear solver
!  JWP, updated November 1998
       integer :: n,i,j,ifg
       integer, parameter :: max=50
       real, dimension (1:max,1:max) :: a
       real, dimension (1:max) :: x, c
      
       write (6,*) 'Solver for small sets of linear equations...'
       write (6,*) 'Enter number of equations, less than ',max
       read  (5,*) n
       write (6,*)  &
        'Now enter coefficients for each row and corresponding constant'
       write(6,*) 'i.e. ',n+1, ' numbers per row...'
       do i=1,n
        write(6,*) 'Elements for row ', i, ':'
        read(5,*) a(i,1:n), c(i)
       end do
!
       write(6,*) 'Matrix and constants are:'
       do  i=1,n
!      check print
       if (n>=4) then
         write(6,'(7f10.4)') a(i,1:n), c(i) ! To make larger printouts tidier
       else
         write(6,*) a(i,1:n), c(i)
       end if

       end do
!
! Solve the equations...
       call gauss  (a, max, x, n, c, ifg)

!
! Print the solution ...
       write(6,*) 'flag =',ifg
       write(6,*) 'solution..'
       write(6,*)  x(1:n)
       stop
!
end program lintest


subroutine gauss (a, max, x, n, c, iflag)
!  solve n linear equations  a.x=c
!  max is upper dimension of arrays
!  iflag is nonzero if solution fails
! *** Restriction!  the diagonal elements are always used as pivots ! ***
!
!  a(,) is the supplied coefficient matrix
!  max is the dimensioned of this array (i.e the max possible no. of equations)
!  x() is where the solution will go
!  n is the actual, not maximum, no. of eqns and unknowns
!  c() is the supplied constant vector
!  iflag must be a variable. It is set to zero if the solution procedure worked
! 
       implicit none
       integer, intent (in) :: max,n
       integer, intent(out) :: iflag
       real, intent (inout) :: a(1:max,1:max),c(1:max)
       real, intent (out) :: x(1:max) 
       real :: factor, pivot, sum
       integer :: i,j       
!
! Set iflag=0. If solution is successful it will return with this value.
! If failure occurs, iflag will be nonzero.
       iflag=0
!
! Start the elimination loop.....
      do  i=1,n-1  ! ... go down the matrix by rows...
!
! Check that the pivot on row i is not zero. Fail if it is
       pivot = a(i,i)
       if (abs(pivot) <= 1e-10) then 
        iflag=1 
        return 
       end if
! 
!
!  ... now reduce all the other rows..
       do  j=i+1,n
! 
           factor = a(j,i) / pivot
!      This is the number to multiply each row by before subtracting it
           a(j,1:n) = a(j,1:n) - factor*a(i,1:n)
!  ..  and also the constant.
           c(j)=c(j) - factor*c(i)
!
       end do  ! .. go on to next row with current elimination
!
      end do    ! .. and then eliminate next equation
!
!     Finally check the last pivot..
      if (abs(a(n,n)) <= 1e-10) then ;
        iflag=1 ;
        return ;
      end if
       
! ..  and divide last constant for back substitution..
      x(n)=c(n)/a(n,n)
! ------------------------ Elimination Complete --------------------------
!
!  --------------------     Back substitute ------------------------------
      do i=n-1,1,-1  ! each equation row 
         sum = c(i)
!
         do  j=i+1,n   ! each term
           sum = sum - a(i,j)*x(j)
         end do
!
         x(i) = sum / a(i,i)
!
      end do
!
       
      return
end