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.
