regrid_solvers.F90

Go to the documentation of this file.

module regrid_solvers
!==============================================================================
!
! This file is part of MOM.
!
! 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.
!
!
!==============================================================================

use mom_error_handler, only : mom_error, fatal

implicit none ; private

public :: solve_linear_system, solve_tridiagonal_system

! -----------------------------------------------------------------------------
! This module contains the following routines
! -----------------------------------------------------------------------------
contains

! -----------------------------------------------------------------------------
! Solve the linear system AX = B
! -----------------------------------------------------------------------------
subroutine solve_linear_system( A, B, X, system_size )
! -----------------------------------------------------------------------------
! 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.
! -----------------------------------------------------------------------------

  ! Arguments
  real, dimension(:,:), intent(inout)   :: A
  real, dimension(:), intent(inout)     :: B
  real, dimension(:), intent(inout)     :: X
  integer                               :: system_size

  ! 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 .LE. system_size ) )

        if ( abs( a(k,i) ) .GT. 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
        end if

    end do ! 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 !' )
    end if

    ! 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 .NE. i ) then
      do j = 1,system_size
        swap_a = a(i,j)
        a(i,j) = a(k,j)
        a(k,j) = swap_a
      end do
      swap_b = b(i)
      b(i) = b(k)
      b(k) = swap_b
    end if

    ! 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
    end do
    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)
      end do
      b(k) = b(k) - factor * b(i)
    end do

  end do ! 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)
    end do
    x(i) = x(i) / a(i,i)
  end do

end subroutine solve_linear_system


! -----------------------------------------------------------------------------
! Solve the tridiagonal system AX = B
! -----------------------------------------------------------------------------
subroutine solve_tridiagonal_system( Al, Ad, Au, B, X, system_size )
! -----------------------------------------------------------------------------
! This routine uses Thomas's algorithm to solve the tridiagonal system AX = B.
! (A is made up of lower, middle and upper diagonals)
! -----------------------------------------------------------------------------
  ! Arguments
  real, dimension(:), intent(inout) :: Al, Ad, Au   ! lo., mid. and up. diagonals
  real, dimension(:), intent(inout) :: B            ! system right-hand side
  real, dimension(:), intent(inout) :: X            ! solution vector
  integer, intent(in)               :: system_size

  ! 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)
  end do

  ! Forward sweep
  do k = 2,n
    b(k) = b(k) - al(k) * b(k-1)
  end do

  ! 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)
  end do

end subroutine solve_tridiagonal_system

end module regrid_solvers