regrid_solvers Module Reference

Functions/Subroutines
subroutine, public	solve_linear_system (A, B, X, system_size)
subroutine, public	solve_tridiagonal_system (Al, Ad, Au, B, X, system_size)

Function/Subroutine Documentation

solve_linear_system()

subroutine, public regrid_solvers::solve_linear_system	(	real, dimension(:,:), intent(inout)	A,
		real, dimension(:), intent(inout)	B,
		real, dimension(:), intent(inout)	X,
		integer	system_size
	)

Definition at line 30 of file regrid_solvers.F90.

References mom_error_handler::mom_error().

Referenced by regrid_edge_slopes::edge_slopes_implicit_h3(), regrid_edge_slopes::edge_slopes_implicit_h5(), regrid_edge_values::edge_values_explicit_h4(), regrid_edge_values::edge_values_implicit_h4(), and regrid_edge_values::edge_values_implicit_h6().

! -----------------------------------------------------------------------------
! 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

Here is the call graph for this function:

Here is the caller graph for this function:

solve_tridiagonal_system()

subroutine, public regrid_solvers::solve_tridiagonal_system	(	real, dimension(:), intent(inout)	Al,
		real, dimension(:), intent(inout)	Ad,
		real, dimension(:), intent(inout)	Au,
		real, dimension(:), intent(inout)	B,
		real, dimension(:), intent(inout)	X,
		integer, intent(in)	system_size
	)

Definition at line 132 of file regrid_solvers.F90.

Referenced by regrid_edge_slopes::edge_slopes_implicit_h3(), regrid_edge_slopes::edge_slopes_implicit_h5(), regrid_edge_values::edge_values_implicit_h4(), and regrid_edge_values::edge_values_implicit_h6().

! -----------------------------------------------------------------------------
! 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

Here is the caller graph for this function: