! MPI Tutorial
! Dr. Andrew C. Pineda, HPCERC/AHPCC
! Dr. Brian Smith, HPCERC/AHPCC
! The University of New Mexico
! November 17, 1997
! Last Revised: September 18, 1998

module GaussianSolver
  implicit none

  ! The default value for the smallest pivot that will be accepted
  ! using the GaussianSolver subroutines. Pivots smaller than this
  ! threshold will cause premature termination of the linear equation
  ! solver and return false as the return value of the function.
  real, parameter :: DEFAULT_SMALLEST_PIVOT = 1.0e-6

contains

  ! Use Gaussian elimination to calculate the solution to the linear
  ! system, A x = b. No partial pivoting is done. If the threshold
  ! argument is present, it is used as the smallest allowable pivot
  ! encountered in the computation; otherwise, DEFAULT_SMALLEST_PIVOT,
  ! defined in this module, is used as the default threshold. The status
  ! of the computation is a logical returned by the function indicating
  ! the existence of a unique solution (.true.), or the nonexistence of
  ! a unique solution or threshold passed (.false.).
  ! Note that this is an inappropriate method for some linear systems.
  ! In particular, the linear system, M x = b, where M = 10e-12 I, will
  ! cause this routine to fail due to the presence of small pivots.
  ! However, this system is perfectly conditioned, with solution x = b.
  function gaussianElimination( A, b, x, threshold )
    implicit none
    logical gaussianElimination
    real, dimension( :, : ), intent( in ) :: A ! Assume the shape of A.
    real, dimension( : ), intent( in ) :: b ! Assume the shape of b.
    real, dimension( : ), intent( out ) :: x ! Assume the shape of x.
    ! The optional attribute specifies that the indicated argument 40
    ! is not required to be present in a call to the function. The
    ! presence of optional arguments, such as threshold, may be checked
    ! using the intrinsic logical function, present (see below).
    real, optional, intent( in ) :: threshold
    integer i, j ! Local index variables.
    integer N ! Order of the linear system.
    real m ! Multiplier.
    real :: smallestPivot = DEFAULT_SMALLEST_PIVOT
    ! Pointers to the appropriate rows of the matrix during the elmination.
    real, dimension( : ), pointer :: pivotRow
    real, dimension( : ), pointer :: currentRow
    ! Copies of the input arguments. These copies are modified during
    ! the computation. The target attribute is used to indicate that
    ! the specified variable may be the target of a pointer. Rows of
    ! ACopy are targets of pivotRow and currentRow, defined above.
    real, dimension( size( A, 1 ), size( A, 2) ), target :: ACopy
    real, dimension( size( b ) ) :: bCopy
    !
    ! Status of the computation. The return value of the function.
    !
    logical successful

    !
    ! Change the smallestPivot if the threshold argument was included.
    !
    if ( present( threshold ) ) smallestPivot = abs( threshold )
    !
    ! Setup the order of the system by using the intrinsic function size.
    ! size returns the number of elements in the specified dimension of
    ! an array or the total number of elements if the dimension is not
    ! specified. Also assume that a unique solution exists initially.
    !
    N = size( b )
    ACopy = A
    bCopy = b
    successful = .true.
    !
    ! Begin the Gaussian elimination algorithm. Note the use of array
    ! sections in the following loops. These eliminate the need for
    ! many do loops that are common in Fortran 77 code. Pointers are
    ! also used below and enhance the readability of the elimination
    ! process. Begin with the first row.
    !
    i = 1
    ! Reduce the system to upper triangular.
    do while ( ( successful ) .and. ( i < N ) )
       !
       ! The following statement is called pointer assignment and uses
       ! the pointer assignment operator '=>'. This causes pivotRow
       ! to be an alias for the ith row of ACopy. Note that this does
       ! not cause any movement of data.
       ! Assign the pivot row.
       !
       pivotRow => ACopy( i, : )
       !
       ! Verify that the current pivot is not smaller than smallestPivot.
       !
       successful = abs( pivotRow( i ) ) >= smallestPivot
       if ( successful ) then
          !
          ! Eliminate the entries in the pivot column below the pivot row.
          !
          do j = i+1, N
             ! Assign the current row.
             currentRow => ACopy( j, : )
             ! Calculate the multiplier.
             m = currentRow( i ) / pivotRow( i )
             ! Perform the elimination step on currentRow and right
             ! hand side, bCopy.
             currentRow = currentRow - m * pivotRow
             bCopy( j ) = bCopy( j ) - m * bCopy( i )
          enddo
       endif
       ! Move to the next row.
       i = i + 1
    end do
    ! Check the last pivot.
    pivotRow => ACopy( N, : )
    if ( successful ) successful = abs( pivotRow( N ) ) >= smallestPivot
    if ( successful ) then
       do i = N, 2, -1 ! Backward substitution.
          ! Determine the ith unknown, x( i ).
          x( i ) = bCopy( i ) / ACopy( i, i )
          ! Substitute the now known value of x( i ), reducing the order of
          ! the system by 1.
          bCopy = bCopy - x( i ) * ACopy( :, i )
       enddo
    endif
    ! Determine the value of x( 1 ) as a special case.
    if ( successful ) x( 1 ) = bCopy( 1 ) / ACopy( 1, 1 )
    ! Prepare the return value of the function.
    gaussianElimination = successful
  end function gaussianElimination

  ! Output A in Matlab format, using name in the Matlab assignment statement.
  subroutine printMatrix( A, name )
    implicit none
    real, dimension( :, : ) :: A ! Assume the shape of A.
    character name ! Name for use in assignment, ie, name =
    ! ......
    integer n, m, i, j
    n = size( A, 1 )
    m = size( A, 2 )
    write( *, fmt="(a1,a5)", advance = "no" ) name, ' = [ '
    ! Output the matrix, except for the last row, which needs no `;'.
    do i = 1, n-1
       ! Output current row.
       do j = 1, m-1
          write( *, fmt="(f10.6,a2)", advance = "no" ) A( i, j ), ', '
       enddo
       ! Output last element in row and end current row.
       write( *, fmt="(f10.6,a1)" ) A( i, m ), ';'
    enddo
    ! Output the last row.
    do j = 1, m-1
       write( *, fmt="(f10.6,a2)", advance = "no" ) A( i, j ), ', '
    enddo
    ! Output last element in row and end.
    write( *, fmt="(f10.6,a1)" ) A( i, m ), ']'
  end subroutine printMatrix

  ! Output b in Matlab format, using name in the Matlab assignment statement.
  subroutine printVector( b, name )
    implicit none
    real, dimension( : ) :: b ! Assume the shape of b.
    character name ! Name for use in assignment, ie, name = ......
    integer n, i
    n = size( b )
    write( *, fmt="(a1,a5)", advance = "no" ) name, ' = [ '
    do i = 1, n-1
       write( *, fmt = "(f10.6,a2)", advance = "no" ) b( i ), ', '
    enddo
    write( *, fmt = "(f10.6,a2)" ) b( n ), ']'
  end subroutine printVector
end module GaussianSolver



! A program to solve linear systems using the GaussianSolver module.
program SolveLinearSystem
  ! Include the module for the various linear solvers.
  use GaussianSolver
  implicit none
  integer, parameter :: N = 5 ! Order of the linear system.
  real, parameter :: TOO_SMALL = 1.0e-7 ! Threshold for pivots.
  ! Declare the necessary arrays and vectors to solve the linear system
  ! A x = b.
  real, dimension( N, N ) :: A ! Coefficient matrix.
  real, dimension( N ) :: x, b ! Vector of unknowns, and right hand side.
  real, dimension( N, N ) :: LU ! Matrix for LU factorization of A.
  logical successful ! Status of computations.

  ! The intrinsic subroutine, random_number, fills a real array or scalar,
  ! with uniformly distributed random variates in the interval [0,1).
  call random_number( A ) ! Initialize the coefficient matrix.
  call random_number( b ) ! Initialize the right-hand side.

  ! Output the matrix in Matlab format for ease of checking the solution.
  call printMatrix( A, 'A' )
  call printVector( b, 'b')

  ! Use Gaussian elmination to calcuate the solution of the linear system.
  ! The call below uses the default threshold specified in the
  ! GaussianSolver module by omitting the optional argument.
  successful = gaussianElimination( A, b, x )
  print *, '===================================='
  print *, 'Gaussian Elimination:'
  print *, '------------------------------------'
  if ( successful ) then
     call printVector( x, 'x' )
     print *, 'Infinity Norm of Difference = ', &
          maxval( abs ( matmul( A, x ) - b ) )
  else
     print *, 'No unique solution or threshold passed.'
  endif
end program SolveLinearSystem

