regrid_solvers Module Reference

Detailed Description

Solvers of linear systems.

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.

Functions/Subroutines
subroutine, public	solve_linear_system (A, B, X, system_size)
	Solve the linear system AX = B by Gaussian elimination. More...
subroutine, public	solve_tridiagonal_system (Al, Ad, Au, B, X, system_size)
	Solve the tridiagonal system AX = B. More...

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, intent(in)	system_size
	)

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.

Parameters

[in,out]	a	The matrix being inverted
[in,out]	b	system right-hand side
[in,out]	x	solution vector
[in]	system_size	The size of the system

Definition at line 20 of file regrid_solvers.F90.

  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
 

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().

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
	)

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)

Parameters

[in,out]	ad	Maxtix center diagonal
[in,out]	al	Matrix lower diagonal
[in,out]	au	Matrix upper diagonal
[in,out]	b	system right-hand side
[in,out]	x	solution vector
[in]	system_size	The size of the system

Definition at line 114 of file regrid_solvers.F90.

  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
 

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().

Here is the caller graph for this function: