regrid_solvers.F90

!> Solvers of linear systems.
module regrid_solvers
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_error_handler, only : mom_error, fatal
 
implicit none ; private
 
public :: solve_linear_system, solve_tridiagonal_system
 
contains
 
!> Solve the linear system AX = B by Gaussian elimination
!!
!! This routine uses Gauss's algorithm to transform the system's original
!! matrix into an upper triangular matrix. Back substitution yields the answer.
!! The matrix A must be square and its size must be that of the vectors B and X.
subroutine solve_linear_system( A, B, X, system_size )
  real, dimension(:,:), intent(inout) :: a  !< The matrix being inverted
  real, dimension(:),   intent(inout) :: b  !< system right-hand side
  real, dimension(:),   intent(inout) :: x  !< solution vector
  integer, intent(in)                 :: system_size !< The size of the system
  ! Local variables
  integer               :: i, j, k
  real, parameter       :: eps = 0.0        ! Minimum pivot magnitude allowed
  real                  :: factor
  real                  :: pivot
  real                  :: swap_a, swap_b
  logical               :: found_pivot      ! boolean indicating whether
                                            ! a pivot has been found
  ! Loop on rows
  do i = 1,system_size-1
 
    found_pivot = .false.
 
    ! Start to look for a pivot in row i. If the pivot
    ! in row i -- which is the current row -- is not valid,
    ! we keep looking for a valid pivot by searching the
    ! entries of column i in rows below row i. Once a valid
    ! pivot is found (say in row k), rows i and k are swaped.
    k = i
    do while ( ( .NOT. found_pivot ) .AND. ( k <= system_size ) )
 
        if ( abs( a(k,i) ) > eps ) then  ! a valid pivot is found
          found_pivot = .true.
        else                                ! Go to the next row to see
                                            ! if there is a valid pivot there
          k = k + 1
        endif
 
    enddo ! end loop to find pivot
 
    ! If no pivot could be found, the system is singular and we need
    ! to end the execution
    if ( .NOT. found_pivot ) then
      write(0,*) ' A=',a
      call mom_error( fatal, 'The linear system is singular !' )
    endif
 
    ! If the pivot is in a row that is different than row i, that is if
    ! k is different than i, we need to swap those two rows
    if ( k /= i ) then
      do j = 1,system_size
        swap_a = a(i,j)
        a(i,j) = a(k,j)
        a(k,j) = swap_a
      enddo
      swap_b = b(i)
      b(i) = b(k)
      b(k) = swap_b
    endif
 
    ! Transform pivot to 1 by dividing the entire row
    ! (right-hand side included) by the pivot
    pivot = a(i,i)
    do j = i,system_size
      a(i,j) = a(i,j) / pivot
    enddo
    b(i) = b(i) / pivot
 
    ! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1
 
    ! Put zeros in column for all rows below that containing
    ! pivot (which is row i)
    do k = (i+1),system_size    ! k is the row index
      factor = a(k,i)
      do j = (i+1),system_size      ! j is the column index
        a(k,j) = a(k,j) - factor * a(i,j)
      enddo
      b(k) = b(k) - factor * b(i)
    enddo
 
  enddo ! end loop on i
 
 
  ! Solve system by back substituting
  x(system_size) = b(system_size) / a(system_size,system_size)
  do i = system_size-1,1,-1 ! loop on rows, starting from second to last row
    x(i) = b(i)
    do j = (i+1),system_size
      x(i) = x(i) - a(i,j) * x(j)
    enddo
    x(i) = x(i) / a(i,i)
  enddo
 
end subroutine solve_linear_system
 
!> Solve the tridiagonal system AX = B
!!
!! This routine uses Thomas's algorithm to solve the tridiagonal system AX = B.
!! (A is made up of lower, middle and upper diagonals)
subroutine solve_tridiagonal_system( Al, Ad, Au, B, X, system_size )
  real, dimension(:), intent(inout) :: ad  !< Maxtix center diagonal
  real, dimension(:), intent(inout) :: al  !< Matrix lower diagonal
  real, dimension(:), intent(inout) :: au  !< Matrix upper diagonal
  real, dimension(:), intent(inout) :: b   !< system right-hand side
  real, dimension(:), intent(inout) :: x   !< solution vector
  integer, intent(in)               :: system_size !< The size of the system
  ! Local variables
  integer                               :: k        ! Loop index
  integer                               :: n        ! system size
 
  n = system_size
 
  ! Factorization
  do k = 1,n-1
    al(k+1) = al(k+1) / ad(k)
    ad(k+1) = ad(k+1) - al(k+1) * au(k)
  enddo
 
  ! Forward sweep
  do k = 2,n
    b(k) = b(k) - al(k) * b(k-1)
  enddo
 
  ! Backward sweep
  x(n) = b(n) / ad(n)
  do k = n-1,1,-1
    x(k) = ( b(k) - au(k)*x(k+1) ) / ad(k)
  enddo
 
end subroutine solve_tridiagonal_system
 
!> \namespace regrid_solvers
!!
!! Date of creation: 2008.06.12
!! L. White
!!
!! This module contains solvers of linear systems.
!! These routines could (should ?) be replaced later by more efficient ones.
 
end module regrid_solvers