Before starting, I have to thank @loiseaujc for the incredible support and for all the discussions we had. Moreover, this proposal is very close to his implementation of LightKrylov. I need to thank @ivanpribec, @FedericoPerini, @hkvzjal and @fluidnumerics_joe for the discussions. Also, this builds upon Enrico Facca’s PR, and it is being tested in the Krylov suite I’m working on (among other DDTs).
The following aims at laying the foundations for some abstract data type containing the notion of an operator and of the space(s) on which it acts. Nothing about this implementation is off the discussion (naming, operations to implement, implementations…), please share your comments, suggestions or doubts. Contributors are also welcome.
stdlib_LinAlgOperators: Linear Algebra Operators
From the Wikipedia page Operator:
In mathematics, an operator is generally a mapping or function that acts on elements of a space to produce elements of another space (possibly and sometimes required to be the same space).
We can write the operator \mathcal{A} as \mathcal{A} \colon U \rightarrow V where U is the domain (or set of inputs) and V is the codomain (or set of destination). Based on this definition, it is easy to envision something similar to be included in the Standard Library, for the Linear Algebra package.
Operators
To begin, lets consider Intrinsic Operators: those are the operators that take as an input a intrinsic array vec_in(:) and that provide a processed intrisic array vec_out(:). The current proposal is the same as Enrico Facca’s PR plus two extra deferred procedures and two extra default implementations.
The two extra procedures are get_domsize() and get_cdmsize() to get respectively the domain size and the codomain size: basically, the first returns the number of columns, the second one the number of rows. Notice that in some cases (like in distributed matrices/vectors, the local size of the operator does not coincide with the size of the domain and codomain, that’s why I’m proposing domsize and cdmsize and still have ncol and nrows defined somewhere else, so to have local and global.)
The two extra default operations are Axpy and Axpby, that can be defined easily through the deferred apply and some BLAS routines.
module stdlib_Operators
use stdlib_kinds
use stdlib_linalg_blas
use stdlib_linalg_lapack
!!
!! Module containing the abstract class parent for an
!! operator derived type, for different precisions,
!! uses AbstractVectors to match the precisions
!!
implicit none
private
!!
!! Public types of operator
!!
public :: Operator_dp
!!
!! real(dp) abstract operator
!!
type, abstract :: Operator_dp
!! Abstract type defining base general operator
!! Members shared by any operator:
!! Number of columns
integer :: ncol=0
!! Number of rows
integer :: nrow=0
!! Logical flag for symmetry
logical :: is_symmetric=.false.
contains
!! Procedure for the application of the operator to a vector.
procedure(application_dp), deferred :: apply
!! Procedure to inqire the domain size
procedure(get_domSize_interface_dp), deferred :: get_domsize
!! Procedure to inqire the domain size
procedure(get_cdmSize_interface_dp), deferred :: get_cdmsize
!
! Procedures defined upon apply
!!
procedure, nopass :: axpby => Axpby_dp
!!
procedure, nopass :: axpy => Axpy_dp
!
end type Operator_dp
abstract interface
!! Abstract procedure defining the interface for a general
!! application actions
!!
!! y = M(x) = ( M*x if M is linear)
subroutine application_dp(this,vec_in,vec_out,ierr)
import dp
import Operator_dp
implicit none
!! The class with all quantities describing the operator
class(Operator_dp), intent(in ) :: this
!! The input" vector
real(dp), intent(in ) :: vec_in(:)
!! The output vector y= "operator apllied on" x
real(dp), intent(inout) :: vec_out(:)
!! An integer flag to inform about error
!! It is zero if no error occured
integer, intent(inout) :: ierr
end subroutine application_dp
!! Abstract procedure defining the interface for a general
!! get_domainSize procedure
!!
!! y = M(x) = ( M*x if M is linear)
function get_domSize_interface_dp(this) result(size)
import Operator_dp
implicit none
!! The class with all quantities describing the operator
class(Operator_dp), intent(in ) :: this
!! The size of the domain
integer :: size
end function get_domSize_interface_dp
!! Abstract procedure defining the interface for a general
!! get_codomainSize procedure
!!
!! y = M(x) = ( M*x if M is linear)
function get_cdmSize_interface_dp(this) result(size)
import Operator_dp
implicit none
!! The class with all quantities describing the operator
class(Operator_dp), intent(in ) :: this
!! The size of the domain
integer :: size
end function get_cdmSize_interface_dp
end interface
contains
!! y = M(x) = ( M*x if M is linear)
subroutine Axpy_dp(this,x,y,ierr)
implicit none
!! The class with all quantities describing the operator
class(Operator_dp), intent(in ) :: this
!! The input" vector
real(dp), intent(in ) :: x(:)
!! The output vector y= "operator applied on" x
real(dp), intent(inout) :: y(:)
!! An integer flag to inform about error
!! It is zero if no error occured
integer, intent(inout) :: ierr
!
real(dp), allocatable :: vec_aux(:)
allocate( vec_aux, mold=x )
call this%apply(x, vec_aux, ierr)
call axpy( size(y), 1._dp, vec_aux, 1, y, 1 )
end subroutine Axpy_dp
!! y = M(x) = ( M*x if M is linear)
subroutine Axpby_dp(this,x,beta,y,ierr)
implicit none
!! The class with all quantities describing the operator
class(Operator_dp), intent(in ) :: this
!! The input" vector
real(dp), intent(in ) :: x(:)
real(dp), intent(in ) :: beta
!! The output vector y= "operator applied on" x
real(dp), intent(inout) :: y(:)
!! An integer flag to inform about error
!! It is zero if no error occured
integer, intent(inout) :: ierr
!
real(dp), allocatable :: vec_aux(:)
allocate( vec_aux, mold=x )
call this%apply(x, vec_aux, ierr)
call scal( size(y), beta, y, 1 )
call axpy( size(y), 1._dp, vec_aux, 1, y, 1 )
end subroutine Axpby_dp
end module stdlib_Operators
Why do we need an Operator acting on intrinsic arrays?
- In SELF (if I remember correctly) the matrices are never explicitly computed but rather only their action is known. A standard way of making a container for such actions could be helpful.
- On shared memory architectures, often you need to define the parallelization of the matrix-vector product, but the input and output vectors are completely in memory.
- In the Krylov suite we will need preconditioners, and since the sparse suite acts on intrinsic arrays, so should the preconditioner. The preconditioner could be defined as an extension of the operator type or not. As for now, it is not.
Abstract Operators and Abstract Vectors
Abstract Operators
The Abstract operator is absolutely the same as the Intrinsic Operator, it just accepts abstract vectors as an input. Clearly, an abstract preconditioner is defined in the same manner as the intrinsic preconditioner, but dealing with abstract vectors.
Abstract Vectors
This module defines abstract vectors and bases of abstract vector. The fundamental procedures needed to be implemented by the user are init and kill, axpby and axpy, dot, getval and setval and rand. Upon these, a set of default implementations is provided.
About axpy: I know that it can be implemented easily upon axpby, but you would require an extra support vector to do the trick and to be fair, the two routines are extremely similar, so I don’t think it’s going to be a burden to code the both (i.e. copy paste rename an entire routine and get rid of the bs).
An important thing is that once you have set these operations, you have almost everything you need to work with bases as well. Some operations are already implemented.
There are also two routines called localizer and globalizer: these are for now just an identity mappings, but those are thought to be handles to access local indexes starting from global indexes and vice versa (e.g. for distributed vectors).
module stdlib_AbstractVectors
use stdlib_kinds
!!
!! Module containing the abstract class parent for a
!! vector derived type, for different precisions
!!
implicit none
private
!!
!! Public types of vector
!!
public :: absVector_dp
public :: linearCombination
public :: projectionOntoBasis
public :: copy
interface linearCombination
module procedure :: linearCombination_Vector_default_dp
module procedure :: linearCombination_Basis_default_dp
end interface
interface projectionOntoBasis
module procedure :: ProjectionOntoBasis_Vector_default_dp
module procedure :: ProjectionOntoBasis_Basis_default_dp
end interface
interface copy
module procedure :: Copy_Basis_default_dp
end interface
!!
!! real(dp) abstract vector
!!
type, abstract :: absVector_dp
!! Members shared by any operator:
!! Number of elements
integer :: size = 0
!!
contains
!! Procedure for the initialization of the vector
procedure(init_vector_interface_dp), pass(self), deferred :: init
!! Procedure for the destruction of the vector
procedure(kill_vector_interface_dp), pass(self), deferred :: kill
!! Procedure to scale and add two vectors
procedure(axpby_vector_interface_dp), pass(y), deferred :: axpby
!! Procedure to scale one vector and add a constant
procedure(apby_vector_interface_dp), pass(y), deferred :: apby
!! Procedure to set the values
procedure(dot_vector_interface_dp), pass(self), deferred :: dot
!! get procedure
procedure(get_value_interface_dp), public, pass(self), deferred :: getval
!! set procedure
procedure(set_value_interface_dp), public, pass(self), deferred :: setval
!! rand procedure
procedure(randvalue_interface_dp), public, pass(self), deferred :: rand
!
! Procedures defined upon axpby
!! Procedure to scale first vector and add a second
procedure, pass(y) :: axpy => axpy_default_dp
!! Procedure to perform addition
procedure, nopass :: add => xpy_default_dp
!! Procedure to perform subtraction
procedure, nopass :: sub => xmy_default_dp
!! Procedure to perform copy from x to y
procedure, pass(to) :: copy => copy_default_dp
!
! Procedures defined upon axpb
!! Procedure to set all the values to one given scalar
procedure, pass(self) :: setall => set_scalar_default_dp
!! Procedure to scale a vector
procedure, nopass :: scal => scale_default_dp
!
! Procedures defined upon dot
!! Procedure to compute the norm
procedure, pass :: norm => norm_default_dp
!
! Localizer and Globalizer
procedure, nopass :: iglob => globalizer_default_dp
procedure, nopass :: iloc => localizer_default_dp
!
end type absVector_dp
abstract interface
!!
!! Abstract procedure defining the interface for a general
!! initialization scheme
!!
subroutine init_vector_interface_dp(self,size,source)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: self
integer, intent(in ) :: size
real(dp), optional, intent(in ) :: source(:)
!
end subroutine init_vector_interface_dp
!!
!! Abstract procedure defining the interface for a general
!! destructor
!!
subroutine kill_vector_interface_dp(self)
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: self
!
end subroutine kill_vector_interface_dp
!!
!! Procedure : axpby_vector_interface_dp(y,alpha,x,beta)
!!
!! Operation : y <- alpha*x + beta*y
!!
!! Rationale : This routine can be used for all basic operations
!! with two vectors
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! y, intent(inout) : absVector
!! alpha, intent(in ) : real(dp)
!! x, intent(in ) : absVector
!! beta, intent(in ) : real(dp)
!!
!!##### Example
!!
!! ...
!! y = [4._dp, 5._dp, 6._dp]
!! call y % axpy( alpha=3._dp, x=[1._dp, 2._dp, 3._dp], beta=1._dp )
!! print *, y ! 7, 11, 15
!! ...
!!
subroutine axpby_vector_interface_dp(y,alpha,x,beta)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: y
real(dp), intent(in ) :: alpha
class(absVector_dp), intent(in ) :: x
real(dp), intent(in ) :: beta
!
end subroutine axpby_vector_interface_dp
!!
!! Procedure : apby_vector_interface_dp(y,alpha,beta)
!!
!! Operation : y <- alpha + beta*y
!!
!! Rationale : This routine can be used to initialize an
!! abstract vector with a single value, or to
!! scale the whole abstract vector with a given
!! scalar, avoiding the allocation of a second
!! abstract vector (that would be needed with axpby)
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! y, intent(inout) : absVector
!! alpha, intent(in ) : real(dp)
!! beta, intent(in ) : real(dp)
!!
!!##### Example
!!
!! ...
!! !y = ...
!! call y % axpy( alpha=3._dp, 3._dp] )
!! print *, y ! 7, 8, 9
!! ...
!!
subroutine apby_vector_interface_dp(y,alpha,beta)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: y
real(dp), intent(in ) :: alpha
real(dp), intent(in ) :: beta
!
end subroutine apby_vector_interface_dp
!!
!! Procedure : get_value_interface_dp(self,index) result(res)
!!
!! Operation : res <- self(index)
!!
!! Rationale : Access a specific index of the (possibly local)
!! abstract vector and return its value
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!! index, intent(in ) : integer
!! res, intent( out) : real(dp)
!!
!!##### Example
!!
!! ...
!! y = [4._dp, 5._dp, 6._dp]
!! res = y%getval(2)
!! print *, res ! 5
!! ...
!!
function get_value_interface_dp(self,index) result(res)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(in ) :: self
integer, intent(in ) :: index
real(dp) :: res
!
end function get_value_interface_dp
!!
!! Procedure : set_value_interface_dp(self,index,val)
!!
!! Operation : self(index) <- val
!!
!! Rationale : Access a specific index of the (possibly local)
!! abstract vector and set the value
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!! index, intent(in ) : integer
!! val, intent(in ) : real(dp)
!!
!!##### Example
!!
!! ...
!! y = [4._dp, 5._dp, 6._dp]
!! call y%setval(2, 11._dp)
!! print *, y%getval(2) ! 11
!! ...
!!
subroutine set_value_interface_dp(self,index,val)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: self
integer, intent(in ) :: index
real(dp), intent(in ) :: val
!
end subroutine set_value_interface_dp
!!
!! Procedure : dot_vector_interface_dp(self,y) result(res)
!!
!! Operation : res <- self^T*y
!!
!! Rationale : Define the dot product between two abstract vectors
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! self, intent(in ) : absVector
!! y, intent(in ) : absVector
!! res, intent( out) : real(dp)
!!
!!##### Example
!!
!! ...
!! y = [1._dp, 2._dp, 3._dp]
!! x = [4._dp, 5._dp, 6._dp]
!! print *, y%dot(x) ! 32
!! ...
!!
function dot_vector_interface_dp(self,y) result(res)
import dp
import absVector_dp
implicit none
class(absVector_dp), intent(in ) :: self
class(absVector_dp), intent(in ) :: y
real(dp) :: res
!
end function dot_vector_interface_dp
!!
!! Procedure : randvector_interface_dp(self,y) result(res)
!!
!! Operation : self <- random_number()
!!
!! Rationale : Define a procedure to randomize the input vector
!!
!! Implementation : deferred
!!
!! Input-Output :
!!
!! self, intent(in ) : absVector
!! seed, intent(in ) : integer
!!
!!##### Example
!!
!! ...
!! [TODO]
!! ...
!!
subroutine randvalue_interface_dp(self,seed)
import absVector_dp
implicit none
class(absVector_dp), intent(inout) :: self
integer, intent(in ), optional :: seed(:)
!
end subroutine randvalue_interface_dp
end interface
contains
!!
!! Procedure : axpy_default_dp(self,alpha,x)
!!
!! Operation : y <- alpha*x + y
!!
!! Rationale : [TODO] forgot why I cared here
!!
!! Implementation : Overridable, built upon AXPBY
!! y <- alpha*x + 1*y = axpby(y,alpha,x)
!!
!! Input-Output :
!!
!! y, intent(inout) : absVector
!! alpha, intent(in ) : real(dp)
!! x, intent(in ) : absVector
!!
!!##### Example
!!
!! ...
!! y = [4._dp, 5._dp, 6._dp]
!! call y % axpy( alpha=3._dp, x=[1._dp, 2._dp, 3._dp] )
!! print *, y ! 7, 11, 15
!! ...
!!
subroutine axpy_default_dp(y,alpha,x)
class(absVector_dp), intent(inout) :: y
real(dp), intent(in ) :: alpha
class(absVector_dp), intent(in ) :: x
!
call y%axpby(alpha,x,1._dp)
!
end subroutine axpy_default_dp
!!
!! Procedure : xpy_default_dp(self,x)
!!
!! Operation : y <- x + y
!!
!! Rationale : This procedure can be used to override the (+) operator
!!
!! Implementation : Overridable, built upon AXPBY
!! y <- 1*x + 1*y = axpby(1._dp,y,1_dp,x)
!!
!! Input-Output :
!!
!! y, intent(inout) : absVector
!! x, intent(in ) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! y = [4._dp, 5._dp, 6._dp]
!! call y % xpy( x )
!! print *, y ! [5._dp, 7._dp, 9._dp]
!! ...
!!
subroutine xpy_default_dp(self,x)
class(absVector_dp), intent(inout) :: self
class(absVector_dp), intent(in ) :: x
!
call self%axpby(1._dp,x,1._dp)
!
end subroutine xpy_default_dp
!!
!! Procedure : xpy_default_dp(self,x)
!!
!! Operation : y <- x - y
!!
!! Rationale : This procedure can be used to override the (-) operator
!!
!! Implementation : Overridable, built upon AXPBY
!! y <- (-1)*x + 1*y = axpby(1._dp,y,-1_dp,x)
!!
!! Input-Output :
!!
!! y, intent(inout) : absVector
!! x, intent(in ) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! y = [4._dp, 5._dp, 6._dp]
!! call y % xmy( x )
!! print *, y ! [-3._dp, -3._dp, -3._dp]
!! ...
!!
subroutine xmy_default_dp(self,x)
class(absVector_dp), intent(inout) :: self
class(absVector_dp), intent(in ) :: x
!
call self%axpby(1._dp,x,-1._dp)
!
end subroutine xmy_default_dp
!!
!! Procedure : copy_default_dp(to,from)
!!
!! Operation : to <- from
!!
!! Rationale : This procedure can be used to override the (=) operator
!!
!! Implementation : Overridable, built upon AXPBY
!! y <- 1*x + 0*y = axpby(0._dp,y,1_dp,x)
!!
!! Input-Output :
!!
!! to, intent( out) : absVector
!! from, intent(in ) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! y = [4._dp, 5._dp, 6._dp]
!! call y % copy( x )
!! print *, y ! [1._dp, 2._dp, 3._dp]
!! ...
!!
subroutine copy_default_dp(to,from)
class(absVector_dp), intent( out) :: to
class(absVector_dp), intent(in ) :: from
!
call to%axpby(1._dp,from,0._dp)
!
end subroutine copy_default_dp
!!
!! Procedure : scale_default_dp(self, beta)
!!
!! Operation : self <- beta*self
!!
!! Rationale : This procedure can be used to initialise and to scale
!!
!! Implementation : Overridable, built upon AXPB
!! self <- 0 + beta*self = axpb(self,beta)
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! call x % scal( 3._dp )
!! print *, x ! [3._dp, 6._dp, 9._dp]
!! ...
!!
subroutine scale_default_dp(self, beta)
class(absVector_dp), intent(inout) :: self
real(dp), intent(in ) :: beta
!
call self%apby(0._dp,beta)
!
end subroutine scale_default_dp
!!
!! Procedure : set_scalar_default_dp(self, beta)
!!
!! Operation : self <- beta
!!
!! Rationale : This procedure can be used to initialise and to scale
!!
!! Implementation : Overridable, built upon AXPB
!! self <- 0*self + beta = axpb(self,0._dp,beta)
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! call x % setall( 3._dp )
!! print *, x ! [3._dp, 3._dp, 3._dp]
!! ...
!!
subroutine set_scalar_default_dp(self, beta)
class(absVector_dp), intent(inout) :: self
real(dp), intent(in ) :: beta
!
call self%apby(0._dp,beta)
!
end subroutine set_scalar_default_dp
!!
!! Procedure : norm_default_dp(self) result(res)
!!
!! Operation : res <- || self ||
!!
!! Rationale : This procedure implements the norm of a vector
!!
!! Implementation : Overridable, built upon DOT
!! self <- sqrt( self^T * self ) = sqrt(dot(self,self))
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!!
!!##### Example
!!
!! ...
!! x = [1._dp, 2._dp, 3._dp]
!! print *, x % norm( ) ! sqrt(14)§
!! ...
!!
function norm_default_dp(self) result(res)
class(absVector_dp), intent(in ) :: self
real(dp) :: res
!
res = sqrt(self%dot(self))
!
end function norm_default_dp
!!
!! Procedure : globalizer_default_dp(iloc) result(iglob)
!!
!! Operation : global_index <- local_index
!!
!! Rationale : This routine converts the local index into a global
!! index. It can be useful (for Householder and Givens
!! algorithms) in distributed memory environments. On
!! shared memory is the identity map.
!!
!! Implementation : Overridable, identity mapping
!! iglob <- iloc
!!
!! Input-Output :
!!
!! iloc, intent(in ) : integer
!! iglob, intent( out) : integer
!!
integer function globalizer_default_dp(iloc) result(iglob)
implicit none
integer, intent(in ) :: iloc
!
iglob=iloc
!
end function globalizer_default_dp
!!
!! Procedure : localizer_default_dp(iglob) result(iloc)
!!
!! Operation : local_index <- global_index
!!
!! Rationale : This routine converts the global index into a local
!! index. It can be useful (for Householder and Givens
!! algorithms) in distributed memory environments. On
!! shared memory is the identity map.
!!
!! Implementation : Overridable, identity mapping
!! iloc <- iglob
!!
!! Input-Output :
!!
!! iglob, intent(in ) : integer
!! iloc, intent( out) : integer
!!
integer function localizer_default_dp(iglob) result(iloc)
implicit none
integer, intent(in ) :: iglob
!
iloc=iglob
!
end function localizer_default_dp
!!
!! Procedure : linearCombination_Basis_default_dp(self, coefs, basis)
!!
!! Operation : self <- \sum_{i} basis(i)*coefs(i)
!!
!! Rationale : Procedure implementing linear combinations.
!!
!! Implementation : Overridable, built upon DOT
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector(:)
!! coefs, intent( out) : real(dp)(:,:)
!! basis, intent(in ) : absVector(:)
!!
!!
subroutine linearCombination_Vector_default_dp(self, coefs, basis)
class(absVector_dp), intent(inout) :: self
real(dp), intent(in ) :: coefs(:)
class(absVector_dp), intent(in ) :: basis(:)
!
integer :: i
!
if ( basis(1)%size /= self%size ) STOP "Incompatible dimensions"
call self%setall(0._dp)
do i = 1, size(basis)
call self%axpby(coefs(i), basis(i), 1._dp)
end do
!
end subroutine linearCombination_Vector_default_dp
!!
!! Procedure : linearCombination_Basis_default_dp(self, coefs, basis)
!!
!! Operation : self(j) <- \sum_{i} basis(i)*coefs(i,j) for all j
!!
!! Rationale : Procedure implementing linear combinations. Notice that
!! the coefficients are stored column-wise (i.e. each column
!! contains one set of coefficients for a linear combination)
!!
!! Implementation : Overridable, built upon DOT
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector(:)
!! coefs, intent( out) : real(dp)(:,:)
!! basis, intent(in ) : absVector(:)
!!
!!
subroutine linearCombination_Basis_default_dp(self, coefs, basis)
class(absVector_dp), intent(inout) :: self(:)
real(dp), intent(in ) :: coefs(:,:)
class(absVector_dp), intent(in ) :: basis(:)
!
integer :: i, j
!
if ( basis(1)%size /= self(1)%size ) STOP "Incompatible dimensions"
do j = 1, size(self)
call self(j)%setall(0._dp)
do i = 1, size(basis)
call self(j)%axpby(coefs(i,j), basis(i), 1._dp)
end do
end do
!
end subroutine linearCombination_Basis_default_dp
!!
!! Procedure : ProjectionOntoBasis_Vector_default_dp(self, coefs, basis)
!!
!! Operation : coefs(i) <- self * basis(i) for all i
!!
!! Rationale : Project one abstract vector onto a basis of abstract vectors.
!!
!! Implementation : Overridable, built upon DOT
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector
!! coefs, intent( out) : real(dp)(:)
!! basis, intent(in ) : absVector(:)
!!
!!
subroutine ProjectionOntoBasis_Vector_default_dp(self, coefs, basis)
class(absVector_dp), intent(inout) :: self
real(dp), intent( out) :: coefs(:)
class(absVector_dp), intent(in ) :: basis(:)
!
integer :: i
!
coefs = 0._dp
if ( basis(1)%size /= self%size ) STOP "Incompatible dimensions"
do i = 1, size(basis)
coefs(i) = self%dot( basis(i) )
end do
!
end subroutine ProjectionOntoBasis_Vector_default_dp
!!
!! Procedure : ProjectionOntoBasis_Basis_default_dp(self, coefs, basis)
!!
!! Operation : coefs(i,j) <- self(j) * basis(i) for all i and j
!!
!! Rationale : Project one array of abstract vectors onto a basis of abstract
!! vectors. Notice that the coefficients are stored column-wise
!! (i.e. each column contains one set of coefficients for a
!! linear combination)
!!
!! Implementation : Overridable, built upon DOT
!!
!! Input-Output :
!!
!! self, intent(inout) : absVector(:)
!! coefs, intent( out) : real(dp)(:,:)
!! basis, intent(in ) : absVector(:)
!!
!!
subroutine ProjectionOntoBasis_Basis_default_dp(self, coefs, basis)
class(absVector_dp), intent(inout) :: self(:)
real(dp), intent( out) :: coefs(:,:)
class(absVector_dp), intent(in ) :: basis(:)
!
integer :: i, j
!
coefs = 0._dp
if ( basis(1)%size /= self(1)%size ) STOP "Incompatible dimensions"
do j = 1, size(self)
call self(j)%setall(0._dp)
do i = 1, size(basis)
coefs(i,j) = self(j)%dot( basis(i) )
end do
end do
!
end subroutine ProjectionOntoBasis_Basis_default_dp
!!
!! Procedure : Copy_Basis_default_dp(to,from)
!!
!! Operation : to(i) <- from(i) for all i
!!
!! Rationale : Copy one array of abstract vectors onto
!! another array
!!
!! Implementation : Overridable, built upon COPY
!!
!! Input-Output :
!!
!! to, intent( out) : absVector(:)
!! from, intent(in ) : absVector(:)
!!
subroutine Copy_Basis_default_dp(to,from)
class(absVector_dp), intent( out) :: to(:)
class(absVector_dp), intent(in ) :: from(:)
!
integer :: i
!
if ( from(1)%size /= to(1)%size ) STOP "Incompatible dimensions"
do i = 1, size(to)
call to(i)%copy(from(i) )
end do
!
end subroutine Copy_Basis_default_dp
end module stdlib_AbstractVectors