!=======================================================================
! ch08_gaussian_elimination.f90
!
! Solve a system of linear equations A*x = b by Gaussian elimination
! with scaled row pivoting and back substitution.
!
! Origin:
!   Based on Gauss_2.f90 by Alexander Godunov, October-November 2009.
!   The original pivot choice and elimination method are retained.
!   Comments, singular-matrix reporting, and verification have been
!   revised for the Chapter 8 companion website.
!
! Example:
!   A = [ 0  1  2; 2  1  4; 2  4  6 ], b = [ 4, 3, 7 ].
!   Exact solution: x = [ -2.5, 0, 2 ].
!   The zero in A(1,1) requires a row interchange.
!
! The program prints the original system, the upper triangular system,
! the solution, and the residual computed with the original A and b.
! A small residual checks the equations, but does not by itself prove
! that a poorly conditioned system has an accurate solution.
!=======================================================================
program ch08_gaussian_elimination
implicit none

integer, parameter :: n = 3
double precision :: a(n,n), b(n), a_work(n,n), b_work(n), x(n)
double precision :: residual(n), max_residual
integer :: i, j, info

! Each implied-DO fills one row of A.
data (a(1,j),j=1,n) / 0.0d0, 1.0d0, 2.0d0 /
data (a(2,j),j=1,n) / 2.0d0, 1.0d0, 4.0d0 /
data (a(3,j),j=1,n) / 2.0d0, 4.0d0, 6.0d0 /
data (b(i),i=1,n)   / 4.0d0, 3.0d0, 7.0d0 /

write (*,'(a)') 'Gaussian elimination with scaled row pivoting'
write (*,'(/,a)') 'Original matrix A and right-hand side b:'
do i = 1,n
   write (*,'(3f12.5,3x,f12.5)') (a(i,j),j=1,n), b(i)
end do

! The elimination changes its matrix and right-hand side. Keep the
! originals for an independent check of the computed solution.
a_work = a
b_work = b
call gauss_scaled(a_work,b_work,x,n,info)

if (info /= 0) then
   write (*,'(/,a,i0)') 'No unique solution: zero pivot at step ',info
   stop 1
end if

! Show the upper triangular equations, as in the original program.
write (*,'(/,a)') 'Matrix A and right-hand side b after elimination:'
do i = 1,n
   write (*,'(3f12.5,3x,f12.5)') (a_work(i,j),j=1,n), b_work(i)
end do

write (*,'(/,a)') 'Solution x:'
do i = 1,n
   write (*,'(a,i0,a,f12.6)') 'x(',i,') = ',x(i)
end do

! Check each original equation without calling a matrix-library routine.
! The elimination has changed a_work and b_work, so use a and b here.
max_residual = 0.0d0
do i = 1,n
   residual(i) = -b(i)
   do j = 1,n
      residual(i) = residual(i) + a(i,j)*x(j)
   end do
   max_residual = max(max_residual,abs(residual(i)))
end do
write (*,'(/,a,es12.4)') 'Maximum residual |A*x-b| = ',max_residual

contains

subroutine gauss_scaled(a,b,x,n,info)
!=======================================================================
! Solve A*x = b by forward elimination and back substitution.
!
! Input:
!   n       - number of equations (n >= 1)
!   a(n,n)  - coefficient matrix
!   b(n)    - right-hand side
!
! Output:
!   x(n)    - solution when info = 0
!   info    - 0 on success; pivot step where a zero pivot was found
!
! On return, a contains the upper triangular matrix and b contains
! the transformed right-hand side; the original arrays are destroyed.
!
! Method:
!   At step k, compute s(i) = max(|a(i,j)|), j = k,...,n, for each
!   candidate row. Choose the row with largest |a(i,k)|/s(i), exchange
!   it with row k, and eliminate the coefficients below the pivot.
!   This recomputed active-row scale follows the original Gauss_2.
!   Finally solve the upper triangular equations from bottom to top.
!
! A zero row is skipped before dividing by s(i). An exactly zero
! pivot reports a singular system. A nonzero but small pivot can still
! indicate an ill-conditioned system; use a separate sensitivity check.
!=======================================================================
implicit none
integer, intent(in) :: n
double precision, intent(inout) :: a(n,n), b(n)
double precision, intent(out) :: x(n)
integer, intent(out) :: info

double precision :: scale, score, best, factor, store, sum
integer :: i, j, k, pivot_row

x = 0.0d0
info = 0

! Forward elimination: choose a scaled pivot in each active column.
do k = 1,n-1
   best = 0.0d0
   pivot_row = k
   do i = k,n
      scale = 0.0d0
      do j = k,n
         scale = max(scale,abs(a(i,j)))
      end do
      if (scale == 0.0d0) cycle
      score = abs(a(i,k))/scale
      if (score > best) then
         best = score
         pivot_row = i
      end if
   end do

   if (best == 0.0d0) then
      info = k
      return
   end if

   ! Exchange the active coefficient rows and the matching right sides.
   if (pivot_row /= k) then
      do j = k,n
         store = a(k,j)
         a(k,j) = a(pivot_row,j)
         a(pivot_row,j) = store
      end do
      store = b(k)
      b(k) = b(pivot_row)
      b(pivot_row) = store
   end if

   ! Subtract a multiple of the pivot equation from each lower row.
   do i = k+1,n
      factor = a(i,k)/a(k,k)
      a(i,k) = 0.0d0
      b(i) = b(i) - factor*b(k)
      do j = k+1,n
         a(i,j) = a(i,j) - factor*a(k,j)
      end do
   end do
end do

! The final pivot is not selected inside the loop above.
if (a(n,n) == 0.0d0) then
   info = n
   return
end if

! Back substitution through the upper triangular system.
do i = n,1,-1
   sum = 0.0d0
   do j = i+1,n
      sum = sum + a(i,j)*x(j)
   end do
   x(i) = (b(i)-sum)/a(i,i)
end do
end subroutine gauss_scaled

end program ch08_gaussian_elimination