splitting_gauss_seidel_redblack_cpu_s.f90 Source File

submodule(splitting_m) splitting_gauss_seidel_redblack_cpu_s
    !! Splitting solver with Red-Black Gauss-Seidel splitting method on CPU
    use screen_io_m, only : get_stderr
    implicit none

contains

    module subroutine create_gauss_seidel_redblack_cpu(self, mesh)
        class(splitting_gauss_seidel_redblack_cpu_t), intent(inout) :: self
        type(mesh_cart_t), intent(in), target :: mesh

        self%ndim = mesh%get_n_points()

        self%n_points_red   = mesh%get_n_points_red()
        self%n_points_black = mesh%get_n_points_black()
        self%redblack_indices => mesh%redblack_indices

    end subroutine

    module subroutine apply_gauss_seidel_redblack_cpu(self, a, dinv, x, b)
        class(splitting_gauss_seidel_redblack_cpu_t), intent(inout) :: self
        type(csrmat_t), intent(in) :: a
        real(FP), intent(in), dimension(self%ndim) :: dinv
        real(FP), intent(inout), dimension(self%ndim) :: x
        real(FP), intent(in), dimension(self%ndim) :: b

        if (.not.((a%ndim) == (self%ndim))) then
            write(get_stderr(), *) &
                'error splitting_gauss_seidel_redblack_cpu%apply): &
                &called for wrong-size matrix'
            write(get_stderr(), *) &
                'matrix ndim is', a%ndim, ', gsrb%ndim is', self%ndim
            error stop
        endif

        call compute_subloop(1, self%n_points_red)
        call compute_subloop(self%n_points_red + 1, &
                             self%n_points_red + self%n_points_black)
        call compute_subloop(self%n_points_red + self%n_points_black + 1, &
                             self%ndim)

        contains

            subroutine compute_subloop(kinit, kfinal)
                !! Contained subroutine for the copy/pasted inner loop lines
                integer, intent(in) :: kinit
                integer, intent(in) :: kfinal

                integer :: kk, ll

                ! Sum is supposed to avoid the diagonal term (which would be
                ! x(l) / d_inverse(l)). Adding x(l) to the entire right-hand
                ! side effectively removes the diagonal term. Also this avoids
                ! the divisions/multiplications required if the diagonal term
                ! would be dealt with directly by the dot_product.

                !$omp do private(kk, ll)
                do kk = kinit, kfinal
                    ll = self%redblack_indices(kk)
                    x(ll) = x(ll) + (b(ll) - &
                            dot_product(a%val(a%i(ll):a%i(ll+1)-1), &
                                        x(a%j(a%i(ll):a%i(ll+1)-1)))) * dinv(ll)
                enddo
                !$omp end do

            end subroutine compute_subloop

    end subroutine

    module subroutine destructor_gauss_seidel_redblack_cpu(self)
        type(splitting_gauss_seidel_redblack_cpu_t), intent(inout) :: self

        self%redblack_indices => null()
        self%ndim = 0
        self%n_points_red = 0
        self%n_points_black = 0

    end subroutine

end submodule