program GaussianElimination
! Solve a linear system of equations with Gaussian Elimination
!   and Back Substitution

! SUBROUTINES: mtxrd,mtxwrt,mtxmlt

! Always declare ALL variables
  implicit none
  INTEGER :: indx,jndx,kndx,lndx ! loop counters
  INTEGER :: nsize               ! # of equations (size of the A matrix)
  REAL :: amtx(10,10)            ! A matrix in Ax=B
  REAL :: bvct(10)               ! B vector in Ax=B
  REAL :: xvct(10)               ! x vector in Ax=B
  REAL :: scl(10)                ! Scale factors
  INTEGER :: irow(10)            ! Row numbers (for swapping rows)
  INTEGER :: ipivot              ! Pivot row
  INTEGER :: itemp               ! Temporary value while swapping
  REAL :: tmpval                 ! Temporary value while pivoting
  REAL :: pival                  ! Pivot value
  REAL :: fctr                   ! Multiplication factor for row reduction
  CHARACTER(len = 64) :: filein  ! Input file name
  INTEGER :: ieof                ! Flag to ensure the file opens correctly

! Initially, no row swaps
  do indx = 1, 10
    irow(indx) = indx
  end do

! Introduce the program
  write (*, *) &
    'Program to solve a linear system of equations using'
  write (*, *) &
    '  Gaussian Elimination and Back Substitution'
  write (*, *) ! A Blank Line

! Read from a file (we don't have to type all those numbers)
  write (*, *) 'Enter the input file name: '
  read (*, *) filein
  open (10, file = filein, status = 'old', iostat = ieof)

! If I have trouble opening the file, blow me out of the program!
  if (ieof .ne. 0) then
    write (*, *) 'File error!'
    stop
  end if

! Input the A matrix - write it out for the user
  write (*, *) 'A Matrix:'
  call mtxrd(amtx,nsize,nsize)
  call mtxwrt(amtx,irow,nsize,nsize)
! Input the B vector - write it out for the user
  write (*, *) ! Blank line
  write (*, *) 'B Vector:'
  call vctrd(bvct,nsize)
  call vctwrt(bvct,nsize)
  
! *** FIND SCALE FACTORS *** !

  do indx = 1,nsize
    scl(indx) = abs(amtx(indx,1))
    do jndx = 2, nsize
      if (abs(amtx(indx,jndx)) > scl(indx)) then
        scl(indx) = abs(amtx(indx,jndx))
      end if
    end do
  end do
  
  write (*, *) ! Blank line
  write (*, *) 'Scale factors:'
  call vctwrt(scl,nsize)
  
! *** BEGIN GAUSSIAN ELIMINATION (ROW REDUCTION) *** !

  do kndx = 1, nsize-1  ! For each column k...

    ! Choose the pivot element
    ipivot = kndx   ! WAS ipivot = irow(kndx)
    pival = abs(amtx(irow(ipivot), kndx) / scl(irow(ipivot)))
    do lndx = kndx+1, nsize
      tmpval = abs(amtx(irow(lndx), kndx) / scl(irow(lndx)))
      if (tmpval > pival) then
        pival = tmpval
        ipivot = lndx   ! WAS ipivot = irow(lndx)
      end if
    end do
    
    !write (*, *) ! Blank line
    !write (*, *) 'Swapping rows ', kndx, ' and ', ipivot, ' -- pival is ', pival
    
    ! Row swap
    itemp = irow(kndx)
    irow(kndx) = irow(ipivot)
    irow(ipivot) = itemp

    !call augwrt(amtx,irow,bvct,nsize)
    
    ! Row reduce
    do indx = kndx+1, nsize  ! For each row i > k...
      fctr = amtx(irow(indx),kndx) / amtx(irow(kndx),kndx)
      amtx(irow(indx),kndx) = 0.0      ! We're zeroing out a(i,k)
      do jndx = kndx+1,nsize     ! Scale & subtract the rest of the row
        amtx(irow(indx),jndx) = amtx(irow(indx),jndx) - fctr * amtx(irow(kndx),jndx)
      end do
      bvct(irow(indx)) = bvct(irow(indx)) - fctr * bvct(irow(kndx))
    end do
    
    !write (*, *) 'Reduction:'
    !call augwrt(amtx,irow,bvct,nsize)
    
  end do

! *** END GAUSSIAN ELIMINATION (ROW REDUCTION) *** !


  write (*, *) ! Blank line
  write (*, *) 'Reduced, Augmented Matrix:'
  call augwrt(amtx,irow,bvct,nsize)

! Back Substitute
  do indx = nsize, 1, -1  ! Start at lower right and work upwards
    xvct(indx) = bvct(irow(indx))
! Subtract the previous x-values times their respective coefficients
    do jndx = indx+1, nsize
      xvct(indx) = xvct(indx) - amtx(irow(indx),jndx) * xvct(jndx)
    end do
! Divide by the current coefficient
    xvct(indx) = xvct(indx) / amtx(irow(indx),indx)
  end do

  write (*, *) ! Blank line
  write (*, *) 'Solution Vector:'
  call vctwrt(xvct,nsize)

  stop

contains

  subroutine mtxrd(amtx,mrow,ncol)

! ABSTRACT: Subroutine to read an m x n matrix from a file
!   FORMAT: first line has the dimensions
!           each following line has a row of the matrix

! Always declare ALL variables
    implicit none
    INTEGER :: irow,jcol               ! row, column loop counters
    INTEGER, INTENT(OUT) :: mrow,ncol  ! # rows, columns
    REAL, INTENT(OUT) :: amtx(:,:)     ! matrix

! Read the size of the matrix
    read (10, *) mrow,ncol

! Read the matrix, one row at a time
    do irow = 1,mrow
      read (10, *) (amtx(irow,jcol),jcol=1,ncol)
    end do

    return
  end subroutine mtxrd

  subroutine mtxwrt(amtx,irow,mrow,ncol)

! ABSTRACT: Subroutine to write an m x n matrix to the screen

! Always declare ALL variables
    implicit none
    INTEGER :: indx,jcol              ! row, column loop counters
    INTEGER, INTENT(IN) :: mrow,ncol  ! # rows, columns
    INTEGER, INTENT(IN) :: irow(:)    ! For row swapping
    REAL, INTENT(IN) :: amtx(:,:)     ! matrix

! Write the matrix, one row at a time
    do indx = 1,mrow
      write (*, '('' '',10g10.4)') (amtx(irow(indx),jcol),jcol=1,ncol)
    end do
    return
  end subroutine mtxwrt

  subroutine augwrt(amtx,irow,bvct,nsize)

! ABSTRACT: Subroutine to write an n x n matrix with a vector augmented to it

! Always declare ALL variables
    implicit none
    INTEGER :: indx,jcol            ! row, column loop counters
    INTEGER, INTENT(IN) :: nsize   ! # rows
    INTEGER, INTENT(IN) :: irow(:)    ! For row swapping
    REAL, INTENT(IN) :: amtx(:,:)  ! matrix
    REAL, INTENT(IN) :: bvct(:)    ! vector

! Write the matrix, one row at a time
    do indx = 1,nsize
      write (*, '('' '',10g10.4)') &
        (amtx(irow(indx),jcol),jcol=1,nsize),bvct(irow(indx))
    end do
    return
  end subroutine augwrt

  subroutine vctrd(vctr,nsize)

! ABSTRACT: Subroutine to read a vector from a file
!   FORMAT: first line has the dimensions
!           each following line has a row of the vector

! Always declare ALL variables
    implicit none
    INTEGER :: irow               ! row loop counters
    INTEGER, INTENT(OUT) :: nsize ! # rows
    REAL, INTENT(OUT) :: vctr(:)  ! vector

! Read the size of the matrix
    read (10, *) nsize

! Read the matrix, one row at a time
    do irow = 1,nsize
      read (10, *) vctr(irow)
    end do

    return
  end subroutine vctrd

  subroutine vctwrt(vctr,nsize)

! ABSTRACT: Subroutine to write a vector to the screen

! Always declare ALL variables
    implicit none
    INTEGER :: irow               ! row loop counters
    INTEGER, INTENT(IN) :: nsize ! # rows
    REAL, INTENT(IN) :: vctr(:)  ! vector

! Write the matrix, one row at a time
    do irow = 1,nsize
      write (*, *) vctr(irow)
    end do
    return
  end subroutine vctwrt

end program GaussianElimination