module LightKrylov_AbstractLinops !! This module provides the base classes `abtract_linop_rsp`, `abstract_linop_rdp`, !! `abstract_linop_csp` and `abstract_linop_cdp` which can be used to define your own !! linear operators. To do so, you simply need to provide two type-bound procedures: !! !! - `matvec(self, vec_in, vec_out)` : Computes the matrix-vector product. !! - `rmatvec(self, vec_in, vec_out) : Computes the transpose matrix-vector product. !! !! It also provides extended types to define the identity operator, symmetric linear !! operators, scalar-multiplication of a linear multiplication, as well as addition !! of two linear operators. use stdlib_optval, only: optval use stdlib_linalg_blas, only: gemv use LightKrylov_Logger use LightKrylov_Constants use LightKrylov_Timer_Utils, only: lightkrylov_timer use LightKrylov_AbstractVectors implicit none(type, external) private character(len=*), parameter :: this_module = 'LK_Linops' character(len=*), parameter :: this_module_long = 'Lightkrylov_AbstractLinops' type, abstract, public :: abstract_linop !! Base type to define an abstract linear operator. All other operator types defined !! in `LightKrylov` derive from this fundamental one. !! !! @warning !! Users should not extend this abstract class to define their own types. !! @endwarning integer, private :: matvec_counter = 0 integer, private :: rmatvec_counter = 0 type(lightkrylov_timer) :: matvec_timer = lightkrylov_timer('matvec timer') type(lightkrylov_timer) :: rmatvec_timer = lightkrylov_timer('rmatvec timer') contains procedure, pass(self), public :: get_counter !! Return matvec/rmatvec counter value procedure, pass(self), public :: reset_counter !! Reset matvec/rmatvec counter procedure, pass(self), public :: print_timer_info !! Print current timing data procedure, pass(self), public :: reset_timer => reset_linop_timer !! Reset current timing data procedure, pass(self), public :: finalize_timer => finalize_linop_timer !! Finalize timers and print complete history_info end type abstract_linop !------------------------------------------------------------------------------ !----- Definition of an abstract real(sp) operator with kind=sp ----- !------------------------------------------------------------------------------ type, abstract, extends(abstract_linop), public :: abstract_linop_rsp !! Base type to extend in order to define a real(sp)-valued linear operator. contains private ! User defined procedures procedure(abstract_matvec_rsp), pass(self), deferred, public :: matvec !! Procedure to compute the matrix-vector product \( \mathbf{y} = \mathbf{Ax} \). procedure(abstract_matvec_rsp), pass(self), deferred, public :: rmatvec !! Procedure to compute the reversed matrix-vector product \( \mathbf{y} = \mathbf{A}^H \mathbf{x} \). ! Wrappers including counter increment procedure, pass(self), public :: apply_matvec => apply_matvec_rsp !! Wrapper for matvec including the counter increment procedure, pass(self), public :: apply_rmatvec => apply_rmatvec_rsp !! Wrapper for rmatvec including the counter increment end type abstract_linop_rsp abstract interface subroutine abstract_matvec_rsp(self, vec_in, vec_out) !! Interface for the matrix-vector product. use lightkrylov_AbstractVectors import abstract_linop_rsp implicit none(type, external) class(abstract_linop_rsp) , intent(inout) :: self !! Linear operator \(\mathbf{A}\). class(abstract_vector_rsp), intent(in) :: vec_in !! Vector to be multiplied by \(\mathbf{A}\). class(abstract_vector_rsp), intent(out) :: vec_out !! Result of the matrix-vector product. end subroutine abstract_matvec_rsp end interface type, extends(abstract_linop_rsp), public :: adjoint_linop_rsp !! Utility type to define an adjoint linear operator. The definition of `matvec` and `rmatvec` !! are directly inherited from those used to define `A`. Note that this utility does not !! compute the adjoint for you. It simply provides a utility to define a new operator !! with `matvec` and `rmatvec` switched. class(abstract_linop_rsp), allocatable :: A !! Linear operator whose adjoint needs to be defined. contains private procedure, pass(self), public :: matvec => adjoint_matvec_rsp procedure, pass(self), public :: rmatvec => adjoint_rmatvec_rsp end type adjoint_linop_rsp !-------------------------------------------------------------------------------------------- !----- Definition of an abstract real(sp) exponential propagator with kind=sp ----- !-------------------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_rsp), public :: abstract_exptA_linop_rsp !! Utility type to define the exponential propagator \( \mathbf{\Phi}_\tau \) which is the linear map !! corresponding to the matrix exponential of the (possibly time-dependent) system Jacobian !! \( \mathbf{L}(t) \) over a time horizon \( \tau \) as: !! !! $$ \mathbf{\Phi}_\tau = \int_0^\tau \mathbf{L}(t) \: \text{d}t $$ !! !! Note that explicit knowledge or definition of the Jacobian is not required. This utility function !! is intended for the use in conjuction with a time-stepper algorithm that computes the integral !! directly. !! !! @warning !! While it is not necessary to use this utility operator, it is strongly recommended for operators !! that correspond to exponential propagators to extend from this abstract type to allow for more !! rigorous type checks in the application. !! @endwarning real(sp), public :: tau !! Time horizon for the temporal integration. This variable must be set when the operator is instantiated. end type abstract_exptA_linop_rsp !------------------------------------------------------------------------------ !----- Definition of an abstract real(dp) operator with kind=dp ----- !------------------------------------------------------------------------------ type, abstract, extends(abstract_linop), public :: abstract_linop_rdp !! Base type to extend in order to define a real(dp)-valued linear operator. contains private ! User defined procedures procedure(abstract_matvec_rdp), pass(self), deferred, public :: matvec !! Procedure to compute the matrix-vector product \( \mathbf{y} = \mathbf{Ax} \). procedure(abstract_matvec_rdp), pass(self), deferred, public :: rmatvec !! Procedure to compute the reversed matrix-vector product \( \mathbf{y} = \mathbf{A}^H \mathbf{x} \). ! Wrappers including counter increment procedure, pass(self), public :: apply_matvec => apply_matvec_rdp !! Wrapper for matvec including the counter increment procedure, pass(self), public :: apply_rmatvec => apply_rmatvec_rdp !! Wrapper for rmatvec including the counter increment end type abstract_linop_rdp abstract interface subroutine abstract_matvec_rdp(self, vec_in, vec_out) !! Interface for the matrix-vector product. use lightkrylov_AbstractVectors import abstract_linop_rdp implicit none(type, external) class(abstract_linop_rdp) , intent(inout) :: self !! Linear operator \(\mathbf{A}\). class(abstract_vector_rdp), intent(in) :: vec_in !! Vector to be multiplied by \(\mathbf{A}\). class(abstract_vector_rdp), intent(out) :: vec_out !! Result of the matrix-vector product. end subroutine abstract_matvec_rdp end interface type, extends(abstract_linop_rdp), public :: adjoint_linop_rdp !! Utility type to define an adjoint linear operator. The definition of `matvec` and `rmatvec` !! are directly inherited from those used to define `A`. Note that this utility does not !! compute the adjoint for you. It simply provides a utility to define a new operator !! with `matvec` and `rmatvec` switched. class(abstract_linop_rdp), allocatable :: A !! Linear operator whose adjoint needs to be defined. contains private procedure, pass(self), public :: matvec => adjoint_matvec_rdp procedure, pass(self), public :: rmatvec => adjoint_rmatvec_rdp end type adjoint_linop_rdp !-------------------------------------------------------------------------------------------- !----- Definition of an abstract real(dp) exponential propagator with kind=dp ----- !-------------------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_rdp), public :: abstract_exptA_linop_rdp !! Utility type to define the exponential propagator \( \mathbf{\Phi}_\tau \) which is the linear map !! corresponding to the matrix exponential of the (possibly time-dependent) system Jacobian !! \( \mathbf{L}(t) \) over a time horizon \( \tau \) as: !! !! $$ \mathbf{\Phi}_\tau = \int_0^\tau \mathbf{L}(t) \: \text{d}t $$ !! !! Note that explicit knowledge or definition of the Jacobian is not required. This utility function !! is intended for the use in conjuction with a time-stepper algorithm that computes the integral !! directly. !! !! @warning !! While it is not necessary to use this utility operator, it is strongly recommended for operators !! that correspond to exponential propagators to extend from this abstract type to allow for more !! rigorous type checks in the application. !! @endwarning real(dp), public :: tau !! Time horizon for the temporal integration. This variable must be set when the operator is instantiated. end type abstract_exptA_linop_rdp !------------------------------------------------------------------------------ !----- Definition of an abstract complex(sp) operator with kind=sp ----- !------------------------------------------------------------------------------ type, abstract, extends(abstract_linop), public :: abstract_linop_csp !! Base type to extend in order to define a complex(sp)-valued linear operator. contains private ! User defined procedures procedure(abstract_matvec_csp), pass(self), deferred, public :: matvec !! Procedure to compute the matrix-vector product \( \mathbf{y} = \mathbf{Ax} \). procedure(abstract_matvec_csp), pass(self), deferred, public :: rmatvec !! Procedure to compute the reversed matrix-vector product \( \mathbf{y} = \mathbf{A}^H \mathbf{x} \). ! Wrappers including counter increment procedure, pass(self), public :: apply_matvec => apply_matvec_csp !! Wrapper for matvec including the counter increment procedure, pass(self), public :: apply_rmatvec => apply_rmatvec_csp !! Wrapper for rmatvec including the counter increment end type abstract_linop_csp abstract interface subroutine abstract_matvec_csp(self, vec_in, vec_out) !! Interface for the matrix-vector product. use lightkrylov_AbstractVectors import abstract_linop_csp implicit none(type, external) class(abstract_linop_csp) , intent(inout) :: self !! Linear operator \(\mathbf{A}\). class(abstract_vector_csp), intent(in) :: vec_in !! Vector to be multiplied by \(\mathbf{A}\). class(abstract_vector_csp), intent(out) :: vec_out !! Result of the matrix-vector product. end subroutine abstract_matvec_csp end interface type, extends(abstract_linop_csp), public :: adjoint_linop_csp !! Utility type to define an adjoint linear operator. The definition of `matvec` and `rmatvec` !! are directly inherited from those used to define `A`. Note that this utility does not !! compute the adjoint for you. It simply provides a utility to define a new operator !! with `matvec` and `rmatvec` switched. class(abstract_linop_csp), allocatable :: A !! Linear operator whose adjoint needs to be defined. contains private procedure, pass(self), public :: matvec => adjoint_matvec_csp procedure, pass(self), public :: rmatvec => adjoint_rmatvec_csp end type adjoint_linop_csp !-------------------------------------------------------------------------------------------- !----- Definition of an abstract complex(sp) exponential propagator with kind=sp ----- !-------------------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_csp), public :: abstract_exptA_linop_csp !! Utility type to define the exponential propagator \( \mathbf{\Phi}_\tau \) which is the linear map !! corresponding to the matrix exponential of the (possibly time-dependent) system Jacobian !! \( \mathbf{L}(t) \) over a time horizon \( \tau \) as: !! !! $$ \mathbf{\Phi}_\tau = \int_0^\tau \mathbf{L}(t) \: \text{d}t $$ !! !! Note that explicit knowledge or definition of the Jacobian is not required. This utility function !! is intended for the use in conjuction with a time-stepper algorithm that computes the integral !! directly. !! !! @warning !! While it is not necessary to use this utility operator, it is strongly recommended for operators !! that correspond to exponential propagators to extend from this abstract type to allow for more !! rigorous type checks in the application. !! @endwarning real(sp), public :: tau !! Time horizon for the temporal integration. This variable must be set when the operator is instantiated. end type abstract_exptA_linop_csp !------------------------------------------------------------------------------ !----- Definition of an abstract complex(dp) operator with kind=dp ----- !------------------------------------------------------------------------------ type, abstract, extends(abstract_linop), public :: abstract_linop_cdp !! Base type to extend in order to define a complex(dp)-valued linear operator. contains private ! User defined procedures procedure(abstract_matvec_cdp), pass(self), deferred, public :: matvec !! Procedure to compute the matrix-vector product \( \mathbf{y} = \mathbf{Ax} \). procedure(abstract_matvec_cdp), pass(self), deferred, public :: rmatvec !! Procedure to compute the reversed matrix-vector product \( \mathbf{y} = \mathbf{A}^H \mathbf{x} \). ! Wrappers including counter increment procedure, pass(self), public :: apply_matvec => apply_matvec_cdp !! Wrapper for matvec including the counter increment procedure, pass(self), public :: apply_rmatvec => apply_rmatvec_cdp !! Wrapper for rmatvec including the counter increment end type abstract_linop_cdp abstract interface subroutine abstract_matvec_cdp(self, vec_in, vec_out) !! Interface for the matrix-vector product. use lightkrylov_AbstractVectors import abstract_linop_cdp implicit none(type, external) class(abstract_linop_cdp) , intent(inout) :: self !! Linear operator \(\mathbf{A}\). class(abstract_vector_cdp), intent(in) :: vec_in !! Vector to be multiplied by \(\mathbf{A}\). class(abstract_vector_cdp), intent(out) :: vec_out !! Result of the matrix-vector product. end subroutine abstract_matvec_cdp end interface type, extends(abstract_linop_cdp), public :: adjoint_linop_cdp !! Utility type to define an adjoint linear operator. The definition of `matvec` and `rmatvec` !! are directly inherited from those used to define `A`. Note that this utility does not !! compute the adjoint for you. It simply provides a utility to define a new operator !! with `matvec` and `rmatvec` switched. class(abstract_linop_cdp), allocatable :: A !! Linear operator whose adjoint needs to be defined. contains private procedure, pass(self), public :: matvec => adjoint_matvec_cdp procedure, pass(self), public :: rmatvec => adjoint_rmatvec_cdp end type adjoint_linop_cdp !-------------------------------------------------------------------------------------------- !----- Definition of an abstract complex(dp) exponential propagator with kind=dp ----- !-------------------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_cdp), public :: abstract_exptA_linop_cdp !! Utility type to define the exponential propagator \( \mathbf{\Phi}_\tau \) which is the linear map !! corresponding to the matrix exponential of the (possibly time-dependent) system Jacobian !! \( \mathbf{L}(t) \) over a time horizon \( \tau \) as: !! !! $$ \mathbf{\Phi}_\tau = \int_0^\tau \mathbf{L}(t) \: \text{d}t $$ !! !! Note that explicit knowledge or definition of the Jacobian is not required. This utility function !! is intended for the use in conjuction with a time-stepper algorithm that computes the integral !! directly. !! !! @warning !! While it is not necessary to use this utility operator, it is strongly recommended for operators !! that correspond to exponential propagators to extend from this abstract type to allow for more !! rigorous type checks in the application. !! @endwarning real(dp), public :: tau !! Time horizon for the temporal integration. This variable must be set when the operator is instantiated. end type abstract_exptA_linop_cdp interface adjoint module procedure initialize_adjoint_rsp module procedure initialize_adjoint_rdp module procedure initialize_adjoint_csp module procedure initialize_adjoint_cdp end interface public :: adjoint !-------------------------------------------------- !----- Definition of the Identity map ----- !-------------------------------------------------- type, extends(abstract_linop_rsp), public :: Id_rsp !! Utility type to define the Identity operator. Note that the type-bound procedures !! for `matvec` and `rmatvec` do not have to be defined by the user. contains private procedure, pass(self), public :: matvec => id_matvec_rsp procedure, pass(self), public :: rmatvec => id_matvec_rsp end type Id_rsp type, extends(abstract_linop_rdp), public :: Id_rdp !! Utility type to define the Identity operator. Note that the type-bound procedures !! for `matvec` and `rmatvec` do not have to be defined by the user. contains private procedure, pass(self), public :: matvec => id_matvec_rdp procedure, pass(self), public :: rmatvec => id_matvec_rdp end type Id_rdp type, extends(abstract_linop_csp), public :: Id_csp !! Utility type to define the Identity operator. Note that the type-bound procedures !! for `matvec` and `rmatvec` do not have to be defined by the user. contains private procedure, pass(self), public :: matvec => id_matvec_csp procedure, pass(self), public :: rmatvec => id_matvec_csp end type Id_csp type, extends(abstract_linop_cdp), public :: Id_cdp !! Utility type to define the Identity operator. Note that the type-bound procedures !! for `matvec` and `rmatvec` do not have to be defined by the user. contains private procedure, pass(self), public :: matvec => id_matvec_cdp procedure, pass(self), public :: rmatvec => id_matvec_cdp end type Id_cdp !---------------------------------------------- !----- Definition of scaled linop ----- !---------------------------------------------- type, extends(abstract_linop_rsp), public :: scaled_linop_rsp !! Defines a scaled linear operator \( \mathbf{B} = \sigma \mathbf{A} \) with \( \mathbf{A} \) a real-valued operator and !! \( \sigma \in \mathbb{R} \). The definitions of `matvec` and `rmatvec` are directly inherited from those used to define !! `A` and do not have to be defined by the user. class(abstract_linop_rsp), allocatable :: A !! Base linear operator to be scaled. real(sp) :: sigma !! Scaling factor. contains private procedure, pass(self), public :: matvec => scaled_matvec_rsp procedure, pass(self), public :: rmatvec => scaled_rmatvec_rsp end type scaled_linop_rsp type, extends(abstract_linop_rdp), public :: scaled_linop_rdp !! Defines a scaled linear operator \( \mathbf{B} = \sigma \mathbf{A} \) with \( \mathbf{A} \) a real-valued operator and !! \( \sigma \in \mathbb{R} \). The definitions of `matvec` and `rmatvec` are directly inherited from those used to define !! `A` and do not have to be defined by the user. class(abstract_linop_rdp), allocatable :: A !! Base linear operator to be scaled. real(dp) :: sigma !! Scaling factor. contains private procedure, pass(self), public :: matvec => scaled_matvec_rdp procedure, pass(self), public :: rmatvec => scaled_rmatvec_rdp end type scaled_linop_rdp type, extends(abstract_linop_csp), public :: scaled_linop_csp !! Defines a scaled linear operator \( \mathbf{B} = \sigma \mathbf{A} \) with \( \mathbf{A} \) a complex(sp)-valued operator !! and \( \sigma \in \mathbb{R}\ ) or \( \mathbb{C} \). class(abstract_linop_csp), allocatable :: A !! Base linear operator to be scaled. complex(sp) :: sigma !! Scaling factor. contains private procedure, pass(self), public :: matvec => scaled_matvec_csp procedure, pass(self), public :: rmatvec => scaled_rmatvec_csp end type scaled_linop_csp type, extends(abstract_linop_cdp), public :: scaled_linop_cdp !! Defines a scaled linear operator \( \mathbf{B} = \sigma \mathbf{A} \) with \( \mathbf{A} \) a complex(dp)-valued operator !! and \( \sigma \in \mathbb{R}\ ) or \( \mathbb{C} \). class(abstract_linop_cdp), allocatable :: A !! Base linear operator to be scaled. complex(dp) :: sigma !! Scaling factor. contains private procedure, pass(self), public :: matvec => scaled_matvec_cdp procedure, pass(self), public :: rmatvec => scaled_rmatvec_cdp end type scaled_linop_cdp !------------------------------------------------ !----- Definition of axpby operator ----- !------------------------------------------------ type, extends(abstract_linop_rsp), public :: axpby_linop_rsp !! Utility type to define a composite linear operator \( \mathbf{C} = \alpha \mathbf{A} + \beta \mathbf{B} \). !! The definitions of `matvec` and `rmatvec` are directly inherited from those used to define `A` and `B`. class(abstract_linop_rsp), allocatable :: A, B !! Underlying linear operators. real(sp) :: alpha, beta !! Scaling factors. logical :: transA = .false., transB = .false. !! Logical flag to control whether \( \mathbf{A} \) and/or \( \mathbf{B} \) need to be transposed. contains private procedure, pass(self), public :: matvec => axpby_matvec_rsp procedure, pass(self), public :: rmatvec => axpby_rmatvec_rsp end type axpby_linop_rsp type, extends(abstract_linop_rdp), public :: axpby_linop_rdp !! Utility type to define a composite linear operator \( \mathbf{C} = \alpha \mathbf{A} + \beta \mathbf{B} \). !! The definitions of `matvec` and `rmatvec` are directly inherited from those used to define `A` and `B`. class(abstract_linop_rdp), allocatable :: A, B !! Underlying linear operators. real(dp) :: alpha, beta !! Scaling factors. logical :: transA = .false., transB = .false. !! Logical flag to control whether \( \mathbf{A} \) and/or \( \mathbf{B} \) need to be transposed. contains private procedure, pass(self), public :: matvec => axpby_matvec_rdp procedure, pass(self), public :: rmatvec => axpby_rmatvec_rdp end type axpby_linop_rdp type, extends(abstract_linop_csp), public :: axpby_linop_csp !! Utility type to define a composite linear operator \( \mathbf{C} = \alpha \mathbf{A} + \beta \mathbf{B} \). !! The definitions of `matvec` and `rmatvec` are directly inherited from those used to define `A` and `B`. class(abstract_linop_csp), allocatable :: A, B !! Underlying linear operators. complex(sp) :: alpha, beta !! Scaling factors. logical :: transA = .false., transB = .false. !! Logical flag to control whether \( \mathbf{A} \) and/or \( \mathbf{B} \) need to be transposed. contains private procedure, pass(self), public :: matvec => axpby_matvec_csp procedure, pass(self), public :: rmatvec => axpby_rmatvec_csp end type axpby_linop_csp type, extends(abstract_linop_cdp), public :: axpby_linop_cdp !! Utility type to define a composite linear operator \( \mathbf{C} = \alpha \mathbf{A} + \beta \mathbf{B} \). !! The definitions of `matvec` and `rmatvec` are directly inherited from those used to define `A` and `B`. class(abstract_linop_cdp), allocatable :: A, B !! Underlying linear operators. complex(dp) :: alpha, beta !! Scaling factors. logical :: transA = .false., transB = .false. !! Logical flag to control whether \( \mathbf{A} \) and/or \( \mathbf{B} \) need to be transposed. contains private procedure, pass(self), public :: matvec => axpby_matvec_cdp procedure, pass(self), public :: rmatvec => axpby_rmatvec_cdp end type axpby_linop_cdp !---------------------------------------------------------------- !----- Definition of an abstract symmetric operator ----- !---------------------------------------------------------------- type, abstract, extends(abstract_linop_rsp), public :: abstract_sym_linop_rsp !! Abstract representation of an abstract symmetric (real valued) linear operator. contains end type abstract_sym_linop_rsp !---------------------------------------------------------------- !----- Definition of an abstract symmetric operator ----- !---------------------------------------------------------------- type, abstract, extends(abstract_linop_rdp), public :: abstract_sym_linop_rdp !! Abstract representation of an abstract symmetric (real valued) linear operator. contains end type abstract_sym_linop_rdp !---------------------------------------------------------------------------------- !----- Definition of an abstract Hermitian positive definite operator ----- !---------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_csp), public :: abstract_hermitian_linop_csp !! Abstract representation of an abstract hermitian (complex-valued) linear operator. contains end type abstract_hermitian_linop_csp !---------------------------------------------------------------------------------- !----- Definition of an abstract Hermitian positive definite operator ----- !---------------------------------------------------------------------------------- type, abstract, extends(abstract_linop_cdp), public :: abstract_hermitian_linop_cdp !! Abstract representation of an abstract hermitian (complex-valued) linear operator. contains end type abstract_hermitian_linop_cdp !------------------------------------------------ !----- Convenience dense linop type ----- !------------------------------------------------ type, extends(abstract_linop_rsp), public :: dense_linop_rsp real(sp), allocatable :: data(:, :) contains procedure, pass(self), public :: matvec => dense_matvec_rsp procedure, pass(self), public :: rmatvec => dense_rmatvec_rsp end type dense_linop_rsp type, extends(abstract_linop_rdp), public :: dense_linop_rdp real(dp), allocatable :: data(:, :) contains procedure, pass(self), public :: matvec => dense_matvec_rdp procedure, pass(self), public :: rmatvec => dense_rmatvec_rdp end type dense_linop_rdp type, extends(abstract_linop_csp), public :: dense_linop_csp complex(sp), allocatable :: data(:, :) contains procedure, pass(self), public :: matvec => dense_matvec_csp procedure, pass(self), public :: rmatvec => dense_rmatvec_csp end type dense_linop_csp type, extends(abstract_linop_cdp), public :: dense_linop_cdp complex(dp), allocatable :: data(:, :) contains procedure, pass(self), public :: matvec => dense_matvec_cdp procedure, pass(self), public :: rmatvec => dense_rmatvec_cdp end type dense_linop_cdp interface dense_linop module procedure initialize_dense_linop_from_array_rsp module procedure initialize_dense_linop_from_array_rdp module procedure initialize_dense_linop_from_array_csp module procedure initialize_dense_linop_from_array_cdp end interface public :: dense_linop contains !-------------------------------------------------------------- !----- Getter/Setter routines for abstract_linops ----- !-------------------------------------------------------------- pure integer function get_counter(self, trans) result(count) !! Getter function for the number of matvec calls implicit none(type, external) class(abstract_linop), intent(in) :: self logical, intent(in) :: trans !! matvec or rmatvec? if (trans) then count = self%rmatvec_counter else count = self%matvec_counter end if end function get_counter subroutine reset_counter(self, trans, procedure, counter, reset_timer, soft_reset, clean_timer) !! Setter routine to reset the matvec counter and reset timers implicit none(type, external) class(abstract_linop), intent(inout) :: self logical, intent(in) :: trans !! matvec or rmatvec? character(len=*), intent(in) :: procedure !! name of the caller routine integer, optional, intent(in) :: counter !! optional flag to reset to an integer other than zero. logical, optional, intent(in) :: reset_timer !! optional flag to reset also the timers logical, optional, intent(in) :: soft_reset !! optional flag to choose whether to save previous timing data (default: .true.) logical, optional, intent(in) :: clean_timer !! optional flag to choose whether to fully reset the timer (default: .false.) ! internals integer :: counter_, count_old logical :: reset_timer_ character(len=128) :: msg counter_ = optval(counter, 0) count_old = self%get_counter(trans) reset_timer_ = optval(reset_timer, .false.) if ( count_old /= 0 .or. counter_ /= 0) then if (trans) then write(msg,'(A,I0,A,I0,A)') 'Total number of rmatvecs: ', count_old, '. Resetting counter to ', counter_, '.' call log_message(msg, this_module, 'reset_counter('//trim(procedure)//')') self%rmatvec_counter = counter_ else write(msg,'(A,I0,A,I0,A)') 'Total number of matvecs: ', count_old, '. Resetting counter to ', counter_, '.' call log_message(msg, this_module, 'reset_counter('//trim(procedure)//')') self%matvec_counter = counter_ end if end if if (reset_timer_) call self%reset_timer(trans, soft_reset, clean_timer) end subroutine reset_counter subroutine print_timer_info(self, trans) !! Getter routine to print the current timing information for matvec/rmatvec !! Note: Wrapper of the corresponding routine from lightkrylov_timer implicit none(type, external) class(abstract_linop), intent(inout) :: self logical, optional, intent(in) :: trans !! matvec or rmatvec? ! internal logical :: transpose transpose = optval(trans, .false.) if (transpose) then call self%rmatvec_timer%print_info() else call self%matvec_timer%print_info() end if end subroutine print_timer_info subroutine reset_linop_timer(self, trans, soft, clean) !! Setter routine to reset the matvec/rmatvec timers !! Note: Wrapper of the corresponding routine from lightkrylov_timer implicit none(type, external) class(abstract_linop), intent(inout) :: self logical, optional, intent(in) :: trans !! matvec or rmatvec? logical, optional, intent(in) :: soft logical, optional, intent(in) :: clean ! internal logical :: transpose transpose = optval(trans, .false.) if (present(trans)) then if (transpose) then call self%rmatvec_timer%reset(soft, clean, verbose=.true.) else call self%matvec_timer%reset(soft, clean, verbose=.true.) end if else call self%rmatvec_timer%reset(soft, clean, verbose=.true.) call self%matvec_timer%reset(soft, clean, verbose=.true.) end if end subroutine reset_linop_timer subroutine finalize_linop_timer(self) !! Finalize the matvec/rmatvec timers !! Note: Wrapper of the corresponding routine from lightkrylov_timer implicit none(type, external) class(abstract_linop), intent(inout) :: self call self%matvec_timer%finalize() call self%rmatvec_timer%finalize() end subroutine finalize_linop_timer !--------------------------------------------------------------------- !----- Wrappers for matvec/rmatvec to increment counters ----- !--------------------------------------------------------------------- subroutine apply_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out ! internal character(len=128) :: msg self%matvec_counter = self%matvec_counter + 1 write(msg,'(I0,1X,A)') self%matvec_counter, 'start' call log_debug(msg, this_module, 'matvec') call self%matvec_timer%start() call self%matvec(vec_in, vec_out) call self%matvec_timer%stop() write(msg,'(I0,1X,A)') self%matvec_counter, 'end' call log_debug(msg, this_module, 'matvec') return end subroutine apply_matvec_rsp subroutine apply_rmatvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out ! internal character(len=128) :: msg self%rmatvec_counter = self%rmatvec_counter + 1 write(msg,'(I0,1X,A)') self%rmatvec_counter, 'start' call log_debug(msg, this_module, 'rmatvec') call self%rmatvec_timer%start() call self%rmatvec(vec_in, vec_out) call self%rmatvec_timer%stop() write(msg,'(I0,1X,A)') self%rmatvec_counter, 'end' call log_debug(msg, this_module, 'rmatvec') end subroutine apply_rmatvec_rsp subroutine apply_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out ! internal character(len=128) :: msg self%matvec_counter = self%matvec_counter + 1 write(msg,'(I0,1X,A)') self%matvec_counter, 'start' call log_debug(msg, this_module, 'matvec') call self%matvec_timer%start() call self%matvec(vec_in, vec_out) call self%matvec_timer%stop() write(msg,'(I0,1X,A)') self%matvec_counter, 'end' call log_debug(msg, this_module, 'matvec') return end subroutine apply_matvec_rdp subroutine apply_rmatvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out ! internal character(len=128) :: msg self%rmatvec_counter = self%rmatvec_counter + 1 write(msg,'(I0,1X,A)') self%rmatvec_counter, 'start' call log_debug(msg, this_module, 'rmatvec') call self%rmatvec_timer%start() call self%rmatvec(vec_in, vec_out) call self%rmatvec_timer%stop() write(msg,'(I0,1X,A)') self%rmatvec_counter, 'end' call log_debug(msg, this_module, 'rmatvec') end subroutine apply_rmatvec_rdp subroutine apply_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out ! internal character(len=128) :: msg self%matvec_counter = self%matvec_counter + 1 write(msg,'(I0,1X,A)') self%matvec_counter, 'start' call log_debug(msg, this_module, 'matvec') call self%matvec_timer%start() call self%matvec(vec_in, vec_out) call self%matvec_timer%stop() write(msg,'(I0,1X,A)') self%matvec_counter, 'end' call log_debug(msg, this_module, 'matvec') return end subroutine apply_matvec_csp subroutine apply_rmatvec_csp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out ! internal character(len=128) :: msg self%rmatvec_counter = self%rmatvec_counter + 1 write(msg,'(I0,1X,A)') self%rmatvec_counter, 'start' call log_debug(msg, this_module, 'rmatvec') call self%rmatvec_timer%start() call self%rmatvec(vec_in, vec_out) call self%rmatvec_timer%stop() write(msg,'(I0,1X,A)') self%rmatvec_counter, 'end' call log_debug(msg, this_module, 'rmatvec') end subroutine apply_rmatvec_csp subroutine apply_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out ! internal character(len=128) :: msg self%matvec_counter = self%matvec_counter + 1 write(msg,'(I0,1X,A)') self%matvec_counter, 'start' call log_debug(msg, this_module, 'matvec') call self%matvec_timer%start() call self%matvec(vec_in, vec_out) call self%matvec_timer%stop() write(msg,'(I0,1X,A)') self%matvec_counter, 'end' call log_debug(msg, this_module, 'matvec') return end subroutine apply_matvec_cdp subroutine apply_rmatvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(abstract_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out ! internal character(len=128) :: msg self%rmatvec_counter = self%rmatvec_counter + 1 write(msg,'(I0,1X,A)') self%rmatvec_counter, 'start' call log_debug(msg, this_module, 'rmatvec') call self%rmatvec_timer%start() call self%rmatvec(vec_in, vec_out) call self%rmatvec_timer%stop() write(msg,'(I0,1X,A)') self%rmatvec_counter, 'end' call log_debug(msg, this_module, 'rmatvec') end subroutine apply_rmatvec_cdp !------------------------------------------------------------------------------ !----- Concrete matvec/rmatvec implementations for special linops ----- !------------------------------------------------------------------------------ subroutine id_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(Id_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out call copy(vec_out, vec_in) end subroutine id_matvec_rsp subroutine id_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(Id_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out call copy(vec_out, vec_in) end subroutine id_matvec_rdp subroutine id_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(Id_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out call copy(vec_out, vec_in) end subroutine id_matvec_csp subroutine id_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(Id_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out call copy(vec_out, vec_in) end subroutine id_matvec_cdp subroutine scaled_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_matvec_rsp subroutine scaled_rmatvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_rmatvec_rsp subroutine scaled_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_matvec_rdp subroutine scaled_rmatvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_rmatvec_rdp subroutine scaled_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_matvec_csp subroutine scaled_rmatvec_csp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_rmatvec_csp subroutine scaled_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_matvec_cdp subroutine scaled_rmatvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(scaled_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) ; call vec_out%scal(self%sigma) end subroutine scaled_rmatvec_cdp subroutine axpby_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out ! Working array. class(abstract_vector_rsp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_rmatvec(vec_in, wrk) else call self%A%apply_matvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_rmatvec(vec_in, vec_out) else call self%B%apply_matvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_matvec_rsp subroutine axpby_rmatvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out ! Working array. class(abstract_vector_rsp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_matvec(vec_in, wrk) else call self%A%apply_rmatvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_matvec(vec_in, vec_out) else call self%B%apply_rmatvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_rmatvec_rsp subroutine axpby_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out ! Working array. class(abstract_vector_rdp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_rmatvec(vec_in, wrk) else call self%A%apply_matvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_rmatvec(vec_in, vec_out) else call self%B%apply_matvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_matvec_rdp subroutine axpby_rmatvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out ! Working array. class(abstract_vector_rdp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_matvec(vec_in, wrk) else call self%A%apply_rmatvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_matvec(vec_in, vec_out) else call self%B%apply_rmatvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_rmatvec_rdp subroutine axpby_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out ! Working array. class(abstract_vector_csp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_rmatvec(vec_in, wrk) else call self%A%apply_matvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_rmatvec(vec_in, vec_out) else call self%B%apply_matvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_matvec_csp subroutine axpby_rmatvec_csp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out ! Working array. class(abstract_vector_csp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_matvec(vec_in, wrk) else call self%A%apply_rmatvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_matvec(vec_in, vec_out) else call self%B%apply_rmatvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_rmatvec_csp subroutine axpby_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out ! Working array. class(abstract_vector_cdp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_rmatvec(vec_in, wrk) else call self%A%apply_matvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_rmatvec(vec_in, vec_out) else call self%B%apply_matvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_matvec_cdp subroutine axpby_rmatvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(axpby_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out ! Working array. class(abstract_vector_cdp), allocatable :: wrk ! Allocate working array. allocate(wrk, mold=vec_in) ; call wrk%zero() ! w = A @ x if (self%transA) then call self%A%apply_matvec(vec_in, wrk) else call self%A%apply_rmatvec(vec_in, wrk) endif ! y = B @ x if (self%transB) then call self%B%apply_matvec(vec_in, vec_out) else call self%B%apply_rmatvec(vec_in, vec_out) endif ! y = alpha*w + beta*y call vec_out%axpby(self%alpha, wrk, self%beta) end subroutine axpby_rmatvec_cdp !------------------------------------------------- !----- ADJOINT TYPE-BOUND PROCEDURES ----- !------------------------------------------------- function initialize_adjoint_rsp(A) result(B) implicit none(type, external) class(abstract_linop_rsp), intent(in) :: A class(adjoint_linop_rsp), allocatable :: B allocate(B) ; B%A = A end function initialize_adjoint_rsp subroutine adjoint_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) end subroutine adjoint_matvec_rsp subroutine adjoint_rmatvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) end subroutine adjoint_rmatvec_rsp function initialize_adjoint_rdp(A) result(B) implicit none(type, external) class(abstract_linop_rdp), intent(in) :: A class(adjoint_linop_rdp), allocatable :: B allocate(B) ; B%A = A end function initialize_adjoint_rdp subroutine adjoint_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) end subroutine adjoint_matvec_rdp subroutine adjoint_rmatvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) end subroutine adjoint_rmatvec_rdp function initialize_adjoint_csp(A) result(B) implicit none(type, external) class(abstract_linop_csp), intent(in) :: A class(adjoint_linop_csp), allocatable :: B allocate(B) ; B%A = A end function initialize_adjoint_csp subroutine adjoint_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) end subroutine adjoint_matvec_csp subroutine adjoint_rmatvec_csp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) end subroutine adjoint_rmatvec_csp function initialize_adjoint_cdp(A) result(B) implicit none(type, external) class(abstract_linop_cdp), intent(in) :: A class(adjoint_linop_cdp), allocatable :: B allocate(B) ; B%A = A end function initialize_adjoint_cdp subroutine adjoint_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out call self%A%apply_rmatvec(vec_in, vec_out) end subroutine adjoint_matvec_cdp subroutine adjoint_rmatvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(adjoint_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out call self%A%apply_matvec(vec_in, vec_out) end subroutine adjoint_rmatvec_cdp !-------------------------------------------------------------------------- !----- Type-bound procedures for convenience dense linop type ----- !-------------------------------------------------------------------------- subroutine dense_matvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_rsp) select type(vec_out) type is(dense_vector_rsp) block integer :: m, n real(sp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_rsp ; beta = zero_rsp vec_out = vec_in call gemv("N", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_rsp','OUT',this_module,'dense_matvec_rsp') end select class default call type_error('vec_in','dense_vector_rsp','IN',this_module,'dense_matvec_rsp') end select end subroutine dense_matvec_rsp subroutine dense_rmatvec_rsp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_rsp), intent(inout) :: self class(abstract_vector_rsp), intent(in) :: vec_in class(abstract_vector_rsp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_rsp) select type(vec_out) type is(dense_vector_rsp) block integer :: m, n real(sp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_rsp ; beta = zero_rsp vec_out = vec_in call gemv("T", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_rsp','OUT',this_module,'dense_rmatvec_rsp') end select class default call type_error('vec_in','dense_vector_rsp','IN',this_module,'dense_rmatvec_rsp') end select end subroutine dense_rmatvec_rsp subroutine dense_matvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_rdp) select type(vec_out) type is(dense_vector_rdp) block integer :: m, n real(dp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_rdp ; beta = zero_rdp vec_out = vec_in call gemv("N", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_rdp','OUT',this_module,'dense_matvec_rdp') end select class default call type_error('vec_in','dense_vector_rdp','IN',this_module,'dense_matvec_rdp') end select end subroutine dense_matvec_rdp subroutine dense_rmatvec_rdp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_rdp), intent(inout) :: self class(abstract_vector_rdp), intent(in) :: vec_in class(abstract_vector_rdp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_rdp) select type(vec_out) type is(dense_vector_rdp) block integer :: m, n real(dp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_rdp ; beta = zero_rdp vec_out = vec_in call gemv("T", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_rdp','OUT',this_module,'dense_rmatvec_rdp') end select class default call type_error('vec_in','dense_vector_rdp','IN',this_module,'dense_rmatvec_rdp') end select end subroutine dense_rmatvec_rdp subroutine dense_matvec_csp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_csp) select type(vec_out) type is(dense_vector_csp) block integer :: m, n complex(sp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_csp ; beta = zero_csp vec_out = vec_in call gemv("N", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_csp','OUT',this_module,'dense_matvec_csp') end select class default call type_error('vec_in','dense_vector_csp','IN',this_module,'dense_matvec_csp') end select end subroutine dense_matvec_csp subroutine dense_rmatvec_csp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_csp), intent(inout) :: self class(abstract_vector_csp), intent(in) :: vec_in class(abstract_vector_csp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_csp) select type(vec_out) type is(dense_vector_csp) block integer :: m, n complex(sp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_csp ; beta = zero_csp vec_out = vec_in call gemv("C", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_csp','OUT',this_module,'dense_rmatvec_csp') end select class default call type_error('vec_in','dense_vector_csp','IN',this_module,'dense_rmatvec_csp') end select end subroutine dense_rmatvec_csp subroutine dense_matvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_cdp) select type(vec_out) type is(dense_vector_cdp) block integer :: m, n complex(dp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_cdp ; beta = zero_cdp vec_out = vec_in call gemv("N", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_cdp','OUT',this_module,'dense_matvec_cdp') end select class default call type_error('vec_in','dense_vector_cdp','IN',this_module,'dense_matvec_cdp') end select end subroutine dense_matvec_cdp subroutine dense_rmatvec_cdp(self, vec_in, vec_out) implicit none(type, external) class(dense_linop_cdp), intent(inout) :: self class(abstract_vector_cdp), intent(in) :: vec_in class(abstract_vector_cdp), intent(out) :: vec_out select type(vec_in) type is(dense_vector_cdp) select type(vec_out) type is(dense_vector_cdp) block integer :: m, n complex(dp) :: alpha, beta m = size(self%data, 1) ; n = size(self%data, 2) alpha = one_cdp ; beta = zero_cdp vec_out = vec_in call gemv("C", m, n, alpha, self%data, m, vec_in%data, 1, beta, vec_out%data, 1) end block class default call type_error('vec_out','dense_vector_cdp','OUT',this_module,'dense_rmatvec_cdp') end select class default call type_error('vec_in','dense_vector_cdp','IN',this_module,'dense_rmatvec_cdp') end select end subroutine dense_rmatvec_cdp function initialize_dense_linop_from_array_rsp(A) result(linop) implicit none(type, external) real(sp), intent(in) :: A(:, :) type(dense_linop_rsp) :: linop linop%data = A end function initialize_dense_linop_from_array_rsp function initialize_dense_linop_from_array_rdp(A) result(linop) implicit none(type, external) real(dp), intent(in) :: A(:, :) type(dense_linop_rdp) :: linop linop%data = A end function initialize_dense_linop_from_array_rdp function initialize_dense_linop_from_array_csp(A) result(linop) implicit none(type, external) complex(sp), intent(in) :: A(:, :) type(dense_linop_csp) :: linop linop%data = A end function initialize_dense_linop_from_array_csp function initialize_dense_linop_from_array_cdp(A) result(linop) implicit none(type, external) complex(dp), intent(in) :: A(:, :) type(dense_linop_cdp) :: linop linop%data = A end function initialize_dense_linop_from_array_cdp end module LightKrylov_AbstractLinops