preconditioners.f90

! preconditioners.f90 --- preconditioners for symmetric 5-point stencil system
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver preconditioners ilu0 icc jacobi ssor
! Author: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! Maintainer: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! URL: https://www.dkrz.de/redmine/projects/show/scales-ppm
!
! Redistribution and use in source and binary forms, with or without
! modification, are  permitted provided that the following conditions are
! met:
!
! Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! Neither the name of the DKRZ GmbH nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Commentary:
!
! Code:
!

MODULE preconditioners
    USE solver_public
    USE solver_internal
    USE mo_kind, ONLY: wp, dp
    USE ppm_base, ONLY: abort_ppm
    USE ppm_extents, ONLY: extent_start, extent_end, extent_size, extent_type => extent

    IMPLICIT NONE

    PRIVATE

    !-------------------------------------------------------------------
    ! Internal State Variables
    !-------------------------------------------------------------------
    ! Auxiliary variable for jacobi preconditioner
    REAL(wp), ALLOCATABLE :: jacobi_diag(:,:)
    ! Auxiliary variables for ilu0 preconditioner
    REAL(wp), ALLOCATABLE :: ilu0_diag(:,:)
    ! Auxiliary variables for icc(p) preconditioner
    REAL(wp), ALLOCATABLE :: ICC_C(:,:), ICC_W(:,:,:), ICC_S(:,:,:)

    !-------------------------------------------------------------------
    ! Module Function & Procedure Interface
    !-------------------------------------------------------------------
    PUBLIC :: NONE_PRECOND, JACOBI_PRECOND, ILU0_PRECOND, SSOR_PRECOND, ICC_PRECOND
    PUBLIC :: identity, prep_jacobi, jacobi, prep_ilu0, ilu0, ssor, prep_icc, icc, precond_prepared &
              , jacobi_precond_stencil, ilu0_precond_stencil, ssor_precond_stencil, icc_precond_stencil &
              , jacobi_precond_shifted_stencil, ilu0_precond_shifted_stencil, ssor_precond_shifted_stencil &
              , icc_precond_shifted_stencil

    CONTAINS

    ! Identity function, use this as preconditioner to get the non-preconditioned CG-method
    SUBROUTINE identity(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
    END SUBROUTINE identity

    ! Preparation for Jacobi preconditioner
    SUBROUTINE prep_jacobi()
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ALLOCATE(jacobi_diag(ib:il,jb:jl))

        DO j=jb,jl
            DO i=ib,il
                jacobi_diag(i,j) = 1._wp/stencil%central(i,j)
            ENDDO
        ENDDO

    END SUBROUTINE prep_jacobi

    ! Jacobi preconditioner
    SUBROUTINE jacobi(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: ib, jb, il, jl, i, j

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        DO j=jb,jl
            DO i=ib,il
                r(i,j) = r(i,j)*jacobi_diag(i,j)
            ENDDO
        ENDDO

    END SUBROUTINE jacobi

     ! Preparation of the ILU(0) preconditioner
    SUBROUTINE prep_ilu0()
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ALLOCATE(ilu0_diag(ib:il,jb:jl))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! Calculate diagonal elements of ILU(0) factorization
        ilu0_diag = 0.
        ilu0_diag(ib,jb) = ff(ib,jb)
        DO i=ib+1,il
            ilu0_diag(i,jb) = ff(i,jb) - uf(i-1,jb)**2/ilu0_diag(i-1,jb)
        ENDDO
        DO j=jb+1,jl
            ilu0_diag(ib,j) = ff(ib,j) - vf(ib,j-1)**2/ilu0_diag(ib,j-1)
        ENDDO
        DO j=jb+1,jl
           DO i=ib+1,il
                ilu0_diag(i,j) = ff(i,j) - uf(i-1,j)**2/ilu0_diag(i-1,j) &
                                         - vf(i,j-1)**2/ilu0_diag(i,j-1)
           ENDDO
        ENDDO

    END SUBROUTINE prep_ilu0

    ! ILU(0) preconditioner
    SUBROUTINE ilu0(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: ib, jb, il, jl, i, j
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! The overall system to solve is As = r <-> (LU)s = r

        ! 1.) Solve system Ly = r
        DO i=ib+1,il
            r(i,jb) = r(i,jb) + uf(i-1,jb)/ilu0_diag(i-1,jb)*r(i-1,jb)
        ENDDO
        DO j=jb+1,jl
            r(ib,j) = r(ib,j) + vf(ib,j-1)/ilu0_diag(ib,j-1)*r(ib,j-1)
        ENDDO
        DO j = jb+1,jl
            DO i = ib+1,il
                r(i,j) = r(i,j) + uf(i-1,j)/ilu0_diag(i-1,j)*r(i-1,j) &
                                + vf(i,j-1)/ilu0_diag(i,j-1)*r(i,j-1)
            ENDDO
        ENDDO

        ! 2.) Solve system Us = y
        r(il,jl) = 1._wp/ilu0_diag(il,jl)*r(il,jl)
        DO i=il-1,ib,-1
            r(i,jl) = 1._wp/ilu0_diag(i,jl)*( r(i,jl) + uf(i,jl)*r(i+1,jl) )
        ENDDO
        DO j=jl-1,jb,-1
            r(il,j) = 1._wp/ilu0_diag(il,j)*( r(il,j) + vf(il,j)*r(il,j+1) )
        ENDDO
        DO j = jl-1,jb,-1
            DO i = il-1,ib,-1
                r(i,j) = 1._wp/ilu0_diag(i,j)*( r(i,j) + uf(i,j)*r(i+1,j)   &
                                                       + vf(i,j)*r(i,j+1) )
            ENDDO
        ENDDO

    END SUBROUTINE ilu0

    ! Symmetric SOR preconditioner
    SUBROUTINE ssor(r)
        USE solver_config, ONLY: get_ssor_param

        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: i, j, ib, jb, il, jl
        REAL(wp):: ssorpar
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Get SSOR parameter from configuration
        ssorpar = get_ssor_param()

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        !r(ib,jb) = r(ib,jb)
        DO i=ib+1,il
           r(i,jb) = r(i,jb) + ssorpar*uf(i-1,jb)/ff(i-1,jb)*r(i-1,jb)
        ENDDO
        DO j=jb+1,jl
           r(ib,j) = r(ib,j) + ssorpar*vf(ib,j-1)/ff(ib,j-1)*r(ib,j-1)
        ENDDO
        DO j=jb+1,jl
          DO i=ib+1,il
             r(i,j) = r(i,j) + ssorpar*( uf(i-1,j)/ff(i-1,j)*r(i-1,j)  &
                                       + vf(i,j-1)/ff(i,j-1)*r(i,j-1) )
          ENDDO
        ENDDO

        r(il,jl) = r(il,jl)/ff(il,jl)
        DO i=il-1,ib,-1
           r(i,jl) = ( r(i,jl) + ssorpar*uf(i,jl)*r(i+1,jl) )/ff(i,jl)
        ENDDO
        DO j=jl-1,jb,-1
           r(il,j) = ( r(il,j) + ssorpar*vf(il,j)*r(il,j+1) )/ff(il,j)
        ENDDO
        DO j = jl-1,jb,-1
           DO i = il-1,ib,-1
              r(i,j) = ( r(i,j) + ssorpar*( uf(i,j)*r(i+1,j)           &
                                          + vf(i,j)*r(i,j+1) ) )/ff(i,j)
           ENDDO
        ENDDO

    END SUBROUTINE ssor

    ! Calculates Matrix for Incomplete Cholesky Decomposition (ICC) with fillin p
    ! This routine makes use of a transformation from the barotropic stencil to a
    ! real matrix A with indices m_i, m_j. Variables s_i, s_j refer to the stencil
    ! coordinates, i, j are used according the actual context.
    ! The result L of L*L' ~= A is saved in ICC_C, ICC_W, ICC_S
    SUBROUTINE prep_icc(p)
        INTEGER, INTENT(IN) :: p  ! FILL-IN
        INTEGER :: ib, il, jb, jl, m_i, m_j, m_n, m_m, k
        REAL(wp) :: tmp
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! MAX(1,p-1) handles structurel exception for p=1
        ! ICC_C = Center, _W = West, _S = South in stencil
        ALLOCATE(ICC_C(ib:il,jb:jl), ICC_W(MAX(1,p-1),ib:il,jb:jl), ICC_S(1+p,ib:il,jb:jl))

        ! Initialize
        ICC_C = 1.
        ICC_W = 0.
        ICC_S = 0.

        m_n = extent_size(stencil%extent)    ! dim of barotropic Matrix, m_ = Matrix
        m_m = extent_size(stencil%extent(1)) ! m_m is distance between diagonal and outer diagonal

        DO m_i = 1,m_n

            ! start treating diagonal
            ! L(i,i) = sqrt( A(i,i) - sum(L(i,1:i-1).^2) ) in OCTAVE
            tmp = get_value_of_A(m_i,m_i)

            ! outer diagonal
            DO k = MAX(1,m_i-m_m),m_i-m_m+p
                tmp = tmp - get_value_of_L(m_i,k)**2
            ENDDO

            ! secondary diagonals
            DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                tmp = tmp - get_value_of_L(m_i,k)**2
            ENDDO

            tmp = SQRT(tmp)
            CALL set_value_of_L(m_i, m_i, tmp)
            ! end treating diagonal

            ! start treating secondary diagonals
            DO m_j = m_i+1,m_i+MAX(1,p-1)
                IF ( m_j > m_n ) EXIT
                ! L(j,i) = ( A(j,i) - sum(L(j,1:i-1).*L(i,1:i-1)) )/L(i,i) in OCTAVE
                tmp = get_value_of_A(m_j, m_i)

                ! outer diagonal
                DO k = MAX(1,m_i-m_m),m_i-m_m+p
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                ! secondary diagonals
                DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                tmp = tmp / get_value_of_L(m_i,m_i)
                CALL set_value_of_L(m_j, m_i, tmp)
            ENDDO
            ! end treating secondary diagonals

            ! start treating outer diagonals
            DO m_j = m_i+m_m-p,m_i+m_m
                IF ( m_j > m_n ) EXIT
                ! L(j,i) = ( A(j,i) - sum(L(j,1:i-1).*L(i,1:i-1)) )/L(i,i) in OCTAVE
                tmp = get_value_of_A(m_j, m_i)

                ! outer diagonal
                DO k = MAX(1,m_i-m_m),m_i-m_m+p
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                ! secondary diagonals
                DO k = MAX(1,m_i-MAX(1,p-1)),m_i-1
                    tmp = tmp - get_value_of_L(m_j,k) * get_value_of_L(m_i,k)
                ENDDO

                tmp = tmp / get_value_of_L(m_i,m_i)
                CALL set_value_of_L(m_j, m_i, tmp)
            ENDDO
            ! end treating outer diagonals

        ENDDO

        ! Calculate inverse of ICC_C to avoid division afterwards
        ICC_C = 1._wp/ICC_C

        CONTAINS

        ! Calculates stencil coordinates given a row in matrix A
        SUBROUTINE coord_by_row(row, i, j)
            INTEGER, INTENT(IN) :: row
            INTEGER, INTENT(OUT) :: i, j

            i = MOD(row - 1, m_m) + ib
            j = (row - 1) / m_m + jb

        END SUBROUTINE coord_by_row

        ! Get value in matrix L
        FUNCTION get_value_of_L(i,j) RESULT(res)
            INTEGER, INTENT(IN) :: i, j
            INTEGER :: d, s_i, s_j
            REAL(wp) :: res

            res = 0.

            ! diagonal
            IF ( i == j) THEN
                CALL coord_by_row(i, s_i, s_j)
                res = ICC_C(s_i, s_j)
                RETURN
            ENDIF

            ! secondary diagonal
            d = i - j
            IF ( d >= 1 .AND. d <= MAX(1,p-1) ) THEN
                CALL coord_by_row(i, s_i, s_j)
                res = ICC_W(d, s_i, s_j)
                RETURN
            ENDIF

            ! outer diagonal
            IF ( i - m_m > 0 ) THEN
                d = j - (i - m_m) + 1
                IF ( d >= 1 .AND. d <= p+1 ) THEN
                    CALL coord_by_row(i, s_i, s_j)
                    res = ICC_S(d, s_i, s_j)
                    RETURN
                ENDIF
            ENDIF

        END FUNCTION get_value_of_L

        ! Set value in matrix L
        SUBROUTINE set_value_of_L(i,j,v)
            INTEGER, INTENT(IN) :: i, j
            REAL(wp), INTENT(IN) :: v
            INTEGER :: d, s_i, s_j

            ! diagonal
            IF ( i == j) THEN
                CALL coord_by_row(i, s_i, s_j)
                ICC_C(s_i, s_j) = v
                RETURN
            ENDIF

            ! secondary diagonal
            d = i - j
            IF ( d >= 1 .AND. d <= MAX(1,p-1) ) THEN
                CALL coord_by_row(i, s_i, s_j)
                ICC_W(d, s_i, s_j) = v
                RETURN
            ENDIF

            ! outer diagonal
            IF ( i - m_m > 0 ) THEN
                d = j - (i - m_m) + 1
                IF ( d >= 1 .AND. d <= p+1 ) THEN
                    CALL coord_by_row(i, s_i, s_j)
                    ICC_S(d, s_i, s_j) = v
                    RETURN
                ENDIF
            ENDIF

        END SUBROUTINE set_value_of_L

        ! Get value of virtual barotropic matrix A
        FUNCTION get_value_of_A(i,j) RESULT(res)
            INTEGER, INTENT(IN) :: i, j
            INTEGER :: s_i, s_j
            REAL(wp) :: res

            res = 0.

            IF ( i == j) THEN ! diagonal
                CALL coord_by_row(i, s_i, s_j)
                res = ff(s_i, s_j)
            ELSEIF ( j-i == 1 ) THEN ! right secondary diagonal
                CALL coord_by_row(i, s_i, s_j)
                IF (s_i >= il) RETURN ! hit boundary
                res = -uf(s_i,s_j)
            ELSEIF ( j-i == -1 ) THEN ! left secondary diagonal
                CALL coord_by_row(i, s_i, s_j)
                IF (s_i <= ib) RETURN ! hit boundary
                res = -uf(s_i-1,s_j)
            ELSEIF ( j-i == m_m ) THEN ! right outer diagonal
                CALL coord_by_row(i, s_i, s_j)
                IF (s_j >= jl) RETURN ! hit boundary
                res = -vf(s_i,s_j)
            ELSEIF ( j-i == -m_m ) THEN ! left outer diagonal
                CALL coord_by_row(i, s_i, s_j)
                IF (s_j <= jb) RETURN ! hit boundary
                res = -vf(s_i,s_j-1)
            ENDIF

        END FUNCTION get_value_of_A

    END SUBROUTINE prep_icc

    ! ICC(p) preconditioner
    SUBROUTINE icc(r)
        REAL(wp), INTENT(INOUT) :: r(:,:)
        INTEGER :: ib, il, jb, jl, i, j, k, p, s_l, w_l
        REAL(wp), DIMENSION(:,:), POINTER :: uf, vf, ff

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Define some aliases
        uf => stencil%zonal
        vf => stencil%meridional
        ff => stencil%central

        ! Set halos to zero
        CALL clear_halos(r, stencil%extent)

        ! get the fill-in p of ICC
        p = SIZE(ICC_S,1) - 1
        s_l = p + 1 ! length of south diagonal
        w_l = MAX(1,p-1) ! length of west diagonal

        ! The overall system to solve is As = r <-> (LL')s = r

        ! 1.) Solve system Ly = r

        ! first column
        r(ib,jb) = r(ib,jb)*ICC_C(ib,jb)
        DO i=ib+1,il
           r(i,jb) = ( r(i,jb) - ICC_W(1,i,jb)*r(i-1,jb) )*ICC_C(i,jb)
        ENDDO

        ! all other columns
        DO j = jb+1,jl
            DO i = ib,il
                DO k = 1,s_l
                    r(i,j) = r(i,j) - ICC_S(k,i,j)*r(i+k-1,j-1)
                ENDDO
                DO k = 1,w_l
                    r(i,j) = r(i,j) - ICC_W(k,i,j)*r(i-k,j)
                ENDDO
                r(i,j) = r(i,j)*ICC_C(i,j)
            ENDDO
        ENDDO

        ! 2.) Solve system L's = y
        ! Backward substitution is done on a column level instead of row-wise

        ! all but first columns
        DO j = jl,jb+1,-1
            DO i = il,ib,-1
                r(i,j) = r(i,j)*ICC_C(i,j)
                DO k = w_l,1,-1
                    r(i-k,j) = r(i-k,j) - r(i,j)*ICC_W(k,i,j) ! saxpy
                ENDDO
                DO k = s_l,1,-1
                    r(i+k-1,j-1) = r(i+k-1,j-1) - r(i,j)*ICC_S(k,i,j) ! saxpy
                ENDDO
            ENDDO
        ENDDO

        ! first column
        DO i = il,ib+1,-1
            r(i,jb) = r(i,jb)*ICC_C(i,jb)
            r(i-1,jb) = r(i-1,jb) - r(i,jb)*ICC_W(1,i,jb) ! saxpy
        ENDDO
        r(ib,jb) = r(ib,jb)*ICC_C(ib,jb)

    END SUBROUTINE icc

    ! Check if a certain preconditioner is already prepared
    FUNCTION precond_prepared(preconditioner) RESULT(ans)
        INTEGER, INTENT(IN) :: preconditioner
        LOGICAL :: ans

        SELECT CASE (preconditioner)
        CASE (NONE_PRECOND)
            ans = .TRUE.
        CASE (JACOBI_PRECOND)
            ans = ALLOCATED(jacobi_diag)
        CASE (ILU0_PRECOND)
            ans = ALLOCATED(ilu0_diag)
        CASE (SSOR_PRECOND)
            ans = config%ssor_param /= 0._wp
        CASE (ICC_PRECOND)
            ans = ALLOCATED(ICC_C) .AND. ALLOCATED(ICC_W) .AND. ALLOCATED(ICC_S)
        CASE DEFAULT
            CALL abort_ppm("Selected preconditioner [" &
                 // int2str(preconditioner) // "] does not exist!", &
                 __FILE__, &
                 __LINE__)
        END SELECT

    END FUNCTION precond_prepared

    ! Jacobi preconditioned matrix stencil operation
    SUBROUTINE jacobi_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        IF (.NOT. precond_prepared(JACOBI_PRECOND)) THEN
            CALL prep_jacobi()
        ENDIF
        CALL apply_stencil(field, res_field)

        ! Apply Jacobi preconditioner
        CALL jacobi(res_field)

    END SUBROUTINE jacobi_precond_stencil

    ! Jacobi preconditioned matrix stencil operation shifted by a value
    SUBROUTINE jacobi_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        IF (.NOT. precond_prepared(JACOBI_PRECOND)) THEN
            CALL prep_jacobi()
        ENDIF
        CALL apply_stencil(field, res_field)

        ! Apply Jacobi preconditioner
        CALL jacobi(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)

    END SUBROUTINE jacobi_precond_shifted_stencil

    ! ILU(0) preconditioned matrix stencil operation
    SUBROUTINE ilu0_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply ILU0 preconditioner
        IF (.NOT. precond_prepared(ILU0_PRECOND)) THEN
            CALL prep_ilu0()
        ENDIF
        CALL ilu0(res_field)

    END SUBROUTINE ilu0_precond_stencil

    ! ILU(0) preconditioned matrix stencil operation shifted by a value
    SUBROUTINE ilu0_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply ILU0 preconditioner
        IF (.NOT. precond_prepared(ILU0_PRECOND)) THEN
            CALL prep_ilu0()
        ENDIF
        CALL ilu0(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)

    END SUBROUTINE ilu0_precond_shifted_stencil

    ! SSOR preconditioned matrix stencil operation
    SUBROUTINE ssor_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply SSOR preconditioner
        IF (.NOT. precond_prepared(SSOR_PRECOND)) THEN
            CALL abort_ppm("No SSOR parameter provided!", &
                 __FILE__, &
                 __LINE__)
        ENDIF
        CALL ssor(res_field)

    END SUBROUTINE ssor_precond_stencil

    ! SSOR preconditioned matrix stencil operation shifted by a value
    SUBROUTINE ssor_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply SSOR preconditioner
        IF (.NOT. precond_prepared(SSOR_PRECOND)) THEN
            CALL abort_ppm("No SSOR parameter provided!", &
                 __FILE__, &
                 __LINE__)
        ENDIF
        CALL ssor(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)

    END SUBROUTINE ssor_precond_shifted_stencil

    ! ICC(p) preconditioned matrix stencil operation
    SUBROUTINE icc_precond_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply ICC preconditioner
        IF (.NOT. precond_prepared(ICC_PRECOND)) THEN
            CALL prep_icc(config%icc_param)
        ENDIF
        CALL icc(res_field)

    END SUBROUTINE icc_precond_stencil

    ! ICC(p) preconditioned matrix stencil operation shifted by a value
    SUBROUTINE icc_precond_shifted_stencil(field, res_field)
        USE solver_config, ONLY: apply_stencil

        REAL(wp), INTENT(IN) :: field(:,:)
        REAL(wp), INTENT(OUT) :: res_field(:,:)
        INTEGER :: i, j , ib, jb, il, jl

        ! Define some variables for the ranges of the fields
        ib = extent_start(stencil%extent(1))
        jb = extent_start(stencil%extent(2))
        il = extent_end(stencil%extent(1))
        jl = extent_end(stencil%extent(2))

        ! Apply stencil first
        CALL apply_stencil(field, res_field)

        ! Apply ICC preconditioner
        IF (.NOT. precond_prepared(ICC_PRECOND)) THEN
            CALL prep_icc(config%icc_param)
        ENDIF
        CALL icc(res_field)

        ! Perform the shift
        FORALL(i = ib:il, j = jb:jl) res_field(i,j) = stencil%shift*field(i,j) - res_field(i,j)

    END SUBROUTINE icc_precond_shifted_stencil

END MODULE preconditioners