Compiler Explorer

Source code

! kiss_mkl.F90 -- Procedural interface to oneMKL RCI ISS
!
! This module provides a procedural interface around the
! reverse-communication-based solvers provided in MKL. To use
! these, the caller has to implement a matvec subroutine, which
! is passed as a procedure pointer.
!
! Only double precision is available.
!

module kiss_mkl
implicit none
private

public :: dp
public :: rci_iss_cache
public :: matvec_fun
public :: pcg

integer, parameter :: dp = kind(1.0d0)

type, abstract :: rci_iss_cache
    integer :: n
    integer :: ipar(128)
    real(dp) :: dpar(128)
    real(dp), allocatable :: tmp(:)
    integer :: nrestart = 0
    procedure(matvec_fun), pointer :: psolve => null()
contains
    procedure(matvec_fun), deferred :: matvec
end type

abstract interface
    subroutine matvec_fun(cache,x,y)
        import rci_iss_cache, dp
        class(rci_iss_cache) :: cache
        real(dp), intent(in) :: x(cache%n)
        real(dp), intent(out) :: y(cache%n)
    end subroutine
end interface

real(dp), parameter :: DEFAULT_RTOL = 1.0d-6
real(dp), parameter :: DEFAULT_ATOL = 0

INTERFACE
       SUBROUTINE DCG(n,x,b,rci_request,ipar,dpar,tmp)
       INTEGER            n, rci_request
       DOUBLE PRECISION   x( * ), b( * ), dpar( * )
       DOUBLE PRECISION   tmp( n, * )
       INTEGER            ipar( * )
       END SUBROUTINE
      END INTERFACE

INTERFACE
       SUBROUTINE DCG_INIT(n,x,b,rci_request,ipar,dpar,tmp)
       INTEGER            n, rci_request
       DOUBLE PRECISION   x( * ), b( * ), dpar( * )
       DOUBLE PRECISION   tmp( n, * )
       INTEGER            ipar( * )
       END SUBROUTINE
      END INTERFACE

INTERFACE
       SUBROUTINE DCG_CHECK(n,x,b,rci_request,ipar,dpar,tmp)
       INTEGER            n, rci_request
       DOUBLE PRECISION   x( * ), b( * ), dpar( * )
       DOUBLE PRECISION   tmp( n, * )
       INTEGER            ipar( * )
       END SUBROUTINE
      END INTERFACE

INTERFACE
       SUBROUTINE DCG_GET(n,x,b,rci_request,ipar,dpar,tmp,itercount)
       INTEGER            n, rci_request, itercount
       DOUBLE PRECISION   x( * ), b( * ), dpar( * )
       DOUBLE PRECISION   tmp( n, * )
       INTEGER            ipar( * )
       END SUBROUTINE
      END INTERFACE

contains

subroutine pcg(cache,x,b,info,rtol,atol,maxiter,silent)
        class(rci_iss_cache), intent(inout) :: cache
        real(dp), intent(inout) :: x(cache%n)
        real(dp), intent(in) :: b(cache%n)
        integer, intent(out) :: info
        real(dp), intent(in), optional :: rtol, atol
        integer, intent(in), optional :: maxiter
        logical, intent(in), optional :: silent

integer :: rci_request, itercount

associate(n=>cache%n,ipar=>cache%ipar,dpar=>cache%dpar,tmp=>cache%tmp)

call dcg_init(n,x,b,rci_request,ipar,dpar,tmp)
        if (rci_request /= 0) error stop "Failure in dcg_init"

! Silence warnings and errors
        if (present(silent)) then
            if (silent) ipar(6:7) = 0
        end if

if (present(rtol)) dpar(1) = rtol
        if (present(atol)) dpar(2) = atol

if (present(maxiter)) ipar(5) = maxiter

! Let MKL do the convergence test
        ipar(8) = 1
        ipar(9) = 1
        ipar(10) = 0

! Preconditioning
        if (associated(cache%psolve)) ipar(11) = 1

! Let MKL check zero norm of generated vector
        ipar(12) = 1

call dcg_check(n,x,b,rci_request,ipar,dpar,tmp)
        if (rci_request /= 0) error stop "Failure in dcg_check"

! FIXME: if ipar(6) == 0 or ipar(7) == 0, we must handles
        !        errors ourselves

revcom: do
            call dcg(n,x,b,rci_request,ipar,dpar,tmp)
            if (rci_request == 0) then
                info = 0
                return
            end if
            select case(rci_request)
            case(1)
                ! Matvec
                call cache%matvec(tmp(1),tmp(n+1))
            case(3)
                ! Preconditioner
                call cache%psolve(tmp(2*n+1),tmp(3*n+1))
            case default
                error stop "in iterate loop, unknown request"
            end select
        end do revcom

call dcg_get(n,x,b,rci_request,ipar,dpar,tmp,itercount)
        info = itercount

end associate
    end subroutine pcg

end module kiss_mkl