! simeq.f90  solve simultaneous equations AX=Y 

subroutine simeq(n, sz, A, Y, X)
  implicit none
  integer, intent(in) :: n  ! number of equations
  integer, intent(in) :: sz ! dimension of arrays
  double precision, dimension(sz,sz), intent(in) :: A
  double precision, dimension(sz), intent(in) :: Y
  double precision, dimension(sz), intent(out) :: X

!      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH REAL     
!                COEFFICIENTS   (A) * |X| = |Y|                     
!                                                                   
!      INPUT  : THE NUMBER OF EQUATIONS  n                          
!               THE REAL MATRIX  A   should be A(i)(j) but A(i,j) 
!               THE REAL VECTOR  Y                                  
!      OUTPUT : THE REAL VECTOR  X                                  
!                                                                   
!      METHOD : GAUSS-JORDAN ELIMINATION USING MAXIMUM ELEMENT      
!               FOR PIVOT.                                          
!                                                                   
!      USAGE  :     call simeq(n,A,X,Y)                                 
!                                                                   
!                                                                   
!    WRITTEN BY : JON SQUIRE , 28 MAY 1983                          
!    ORIGONAL DEC 1959 for IBM 650, TRANSLATED TO OTHER LANGUAGES   
!    e.g. FORTRAN converted to Ada converted to C to f90                   

    double precision, dimension(n,n+1) :: B  ! (n,n+1)  WORKING MATRIX 
    integer, dimension(n) :: ROW             ! ROW INTERCHANGE INDICIES 
    integer :: HOLD , I_PIVOT                ! PIVOT INDICIES 
    double precision :: PIVOT                ! PIVOT ELEMENT VALUE 
    double precision :: ABS_PIVOT
    integer :: i, j, k, m

    m = n+1

    ! BUILD WORKING DATA STRUCTURE 
    do i=1,n
      do j=1,n
        B(i,j) = A(i,j)
      end do ! j
      B(i,n+1) = Y(i)
    end do ! i
    ! SET UP ROW  INTERCHANGE VECTORS 
    do k=1,n
      ROW(k) = k
    end do ! k

    ! BEGIN MAIN REDUCTION LOOP 
    do k=1,n

      ! FIND LARGEST ELEMENT FOR PIVOT 
      PIVOT = B(ROW(k),k)
      ABS_PIVOT = abs(PIVOT)
      I_PIVOT = k
      do i=k,n
        if( abs(B(ROW(i),k)) > ABS_PIVOT) then
          I_PIVOT = i
          PIVOT = B(ROW(i),k)
          ABS_PIVOT = abs ( PIVOT )
        end if
      end do ! i

      ! HAVE PIVOT, INTERCHANGE ROW POINTERS 
      HOLD = ROW(k)
      ROW(k) = ROW(I_PIVOT)
      ROW(I_PIVOT) = HOLD

      ! CHECK FOR NEAR SINGULAR 
      if( ABS_PIVOT < 1.0D-52 ) then
        do j=k+1,n+1
          B(ROW(k),j) = 0.0D0
        end do ! j
        print *, 'redundant row (singular) ', ROW(k)
      else

        ! REDUCE ABOUT PIVOT 
        do j=k+1,n+1
          B(ROW(k),j) = B(ROW(k),j) / B(ROW(k),k)
        end do ! j

        ! INNER REDUCTION LOOP 
        do i=1,n
          if( i .ne. k) then
            do j=k+1,n+1
              B(ROW(i),j) = B(ROW(i),j) - B(ROW(i),k) * B(ROW(k),j)
            end do ! j
          end if
        end do ! i
      end if 
      ! FINISHED INNER REDUCTION 
    end do ! k

    ! END OF MAIN REDUCTION LOOP 
    ! BUILD  X  FOR RETURN, UNSCRAMBLING ROWS 
    do i=1,n
      X(i) = B(ROW(i),n+1)
    end do ! i
end subroutine simeq