Defining real and complex variants of an abstract type for vector calculus

Help
3

Unlimited polymorphism certainly is an option indeed although I ain’t very familiar with it. How would it play along with the rest of the library ? Consider for instance the subroutine computing the Arnoldi factorization of a linear operator A. Simplified version of this code is shown below.

subroutine arnoldi_factorization(A, X, H)
      !! Compute the Arnoldi factorization A @ X = X @ H.
      class(abstract_linop), intent(in) :: A
      !! Linear operator to be factorized.
      class(abstract_vector), intent(inout) :: X(:)
      !! Orthogonal Krylov subspace.
      real(kind=wp), intent(inout) :: H(:, :)
      !! Upper Hessenberg matrix.

      ! Internal variables.
      integer :: i, k
      real(kind=wp) :: alpha
      
      do k = 1, size(X)-1
         ! Matrix-vector product X[k+1] = A @ X[k].
         call A%matvec(X(k), X(k+1))
         ! Update Hessenberg matrix and orthogonalize current vector.
         do i = 1, k
              alpha = X(i)%dot(X(k+1)) ; H(i, k) = alpha
              call X(k+1)%axpby(1.0_wp, X(i), -alpha)
         enddo
         ! Normalize current vector
         alpha = X(k+1)%norm() ; H(k+1, k) = alpha
         call X(k+1)%scal(1.0_wp / alpha)
      end do

      return
   end subroutine arnoldi_factorization

Here, A is an abstract linear operator with the matvec type-bound procedure satisfying the following interface

abstract interface
    subroutine matvec(self, x, y)
        !! Compute the matrix-vector product y = A @ x.
        import abstract_vector, abstract_linop
        class(abstract_linop) , intent(in)  :: self
        class(abstract_vector), intent(in)  :: x
        class(abstract_vector), intent(out) :: y
    end subroutine matvec
end interface

Would taking the unlimited polymorphism route implies that:

  • the Hessenberg matrix H needs to be defined as class(*), intent(inout) :: H(:, :),
  • the number alpha be also defined as class(*) :: alpha,
  • Pepper the code with some select class/end select every time I need to access one of these variables ?
  • If the answers to the first two questions are yes, does it also imply that, in order to use arnoldi_factorization in their code, users would have to declare H as an unlimited polymorphic variable as well ?

The idea I originally had in mind was something along the following lines

type, abstract, public :: abstract_vector
!! Base type from which all other vector types will be extended.
contains
end type abstract_vector

!----------------------------------------------------------------
!-----     ABSTRACT REAL VECTOR + TYPE-BOUND PROCEDURES     -----
!----------------------------------------------------------------

type, abstract, extends(abstract_vector), public :: abstract_real_vector
contains
private
procedure(abstract_real_scal) , pass(self), deferred, public :: scal
procedure(abstract_real_axpby), pass(self), deferred, public :: axpby
procedure(abstract_real_dot)  , pass(self), deferred, public :: dot
procedure                     , pass(self), public           :: norm
end type abstract_vector

abstract interface
    subroutine abstract_real_scal(self, alpha)
    import abstract_vector, wp
    class(abstract_vector), intent(inout) :: self
    real(kind=wp), intent(in) :: alpha
    end abstract_real_scal

    subroutine abstract_real_axpby(self, alpha, vec, beta)
    import abstract_vector, wp
    class(abstract_vector), intent(inout) :: self
    real(kind=wp), intent(in) :: alpha
    class(abstract_vector), intent(in) :: vec
    real(kind=wp), intent(in) :: beta
    end subroutine abstract_real_axpby

    real(kind=wp) function abstract_real_dot(self, vec) result(alpha)
    import abstract_vector, wp
    class(abstract_vector), intent(in) :: self, vec
    end function abstract_real_dot
end interface

!-------------------------------------------------------------------
!-----     ABSTRACT COMPLEX VECTOR + TYPE-BOUND PROCEDURES     -----
!-------------------------------------------------------------------

type, abstract, extends(abstract_vector), public :: abstract_cmplx_vector
contains
private
procedure(abstract_cmplx_scal) , pass(self), deferred, public :: scal
procedure(abstract_cmplx_axpby), pass(self), deferred, public :: axpby
procedure(abstract_cmplx_dot)  , pass(self), deferred, public :: dot
procedure                      , pass(self), public           :: norm
end type abstract_vector

abstract interface
    subroutine abstract_cmplx_scal(self, alpha)
    import abstract_vector, wp
    class(abstract_vector), intent(inout) :: self
    cmplx(kind=wp), intent(in) :: alpha
    end abstract_real_scal

    subroutine abstract_cmplx_axpby(self, alpha, vec, beta)
    import abstract_vector, wp
    class(abstract_vector), intent(inout) :: self
    cmplx(kind=wp), intent(in) :: alpha
    class(abstract_vector), intent(in) :: vec
    cmplx(kind=wp), intent(in) :: beta
    end subroutine abstract_cmplx_axpby

    cmplx(kind=wp) function abstract_cmplx_dot(self, vec) result(alpha)
    import abstract_vector, wp
    class(abstract_vector), intent(in) :: self, vec
    end function abstract_cmplx_dot
end interface

I would also need some subroutines to do cmplx_vector = cmplx_vector + real_vector although the opposite is not the case. Then have something like what follows for Arnoldi

interface arnoldi_factorization
     module procedure real_arnoldi_factorization
     module procedure cmplx_arnoldi_factorization
end interface

subroutine real_arnoldi_factorization(A, X, H)
      !! Compute the Arnoldi factorization A @ X = X @ H.
      class(abstract_linop), intent(in) :: A
      !! Linear operator to be factorized.
      class(abstract_real_vector), intent(inout) :: X(:)
      !! Orthogonal Krylov subspace.
      real(kind=wp), intent(inout) :: H(:, :)
      !! Upper Hessenberg matrix.

      ! Internal variables.
      integer :: i, k
      real(kind=wp) :: alpha
      
      do k = 1, size(X)-1
         ! Matrix-vector product X[k+1] = A @ X[k].
         call A%matvec(X(k), X(k+1))
         ! Update Hessenberg matrix and orthogonalize current vector.
         do i = 1, k
              alpha = X(i)%dot(X(k+1)) ; H(i, k) = alpha
              call X(k+1)%axpby(1.0_wp, X(i), -alpha)
         enddo
         ! Normalize current vector.
         alpha = X(k+1)%norm() ; H(k+1, k) = alpha
         call X(k+1)%scal(1.0_wp / alpha)
      end do

      return
   end subroutine arnoldi_factorization

subroutine cmplx_arnoldi_factorization(A, X, H)
      !! Compute the Arnoldi factorization A @ X = X @ H.
      class(abstract_linop), intent(in) :: A
      !! Linear operator to be factorized.
      class(abstract_cmplx_vector), intent(inout) :: X(:)
      !! Orthogonal Krylov subspace.
      cmplx(kind=wp), intent(inout) :: H(:, :)
      !! Upper Hessenberg matrix.

      ! Internal variables.
      integer :: i, k
      cmplx(kind=wp) :: alpha
      
      do k = 1, size(X)-1
         ! Matrix-vector product X[k+1] = A @ X[k].
         call A%matvec(X(k), X(k+1))
         ! Update Hessenberg matrix and orthogonalized current vector.
         do i = 1, k
              alpha = X(i)%dot(X(k+1)) ; H(i, k) = alpha
              call X(k+1)%axpby(1.0_wp, X(i), -alpha)
         enddo
         ! Normalize current vector.
         alpha = X(k+1)%norm() ; H(k+1, k) = alpha
         call X(k+1)%scal(1.0_wp / alpha)
      end do

      return
   end subroutine arnoldi_factorization

What would be the benefit of unlimited polymorphism compared to these real/complex definitions of the different procedures ? And to what extent would unlimited polymorphism impact the overall performances ?

The reason I’m asking these details is in relation with this discussion. I basically want to make this library as user-friendly as possible by exposing a scipy-like API in the end for user-defined notions of vectors and linear operators while being as efficient as possible. And cherry on the cake, I’d like the source code to be easily readable for someone with a basic knowledge of Krylov methods and Fortran.

1 Like