Open enricofacca opened 2 years ago
jvdp1
commented
2 years ago Thank you @enricofacca for the proposition. A while ago, when I was testing Coarrays, I did something similar to your proposition (see here) The idea was that the coefficient matrix and preconditioner (for PCG) would be a DT with a method to multiply it by a vector.
Would you like to propose a first API, e.g., for PCG?
enricofacca
commented
2 years ago Dear @jvdp1 You wrote something really close to what I was thinking. I would also include a further abstract DT called "vector" having methods that will handle the basis linear algebra operation used in iterative solver. Moreover, I was thinking to gather matrix and preconditioner in one class with an "apply" method.
I will be happy to put these ideas into code and write the API.
enricofacca
commented
2 years ago I will present my proposal focusing on the Krylov solver for double precision arrays. I apologe for the length of this comments. The source code with CmakeList file is available here. I would build Krylov solvers on the following abstract type.
module AbstractOperators
!! Module containing the abstract class parent of
!! all operator classes defined on Dvector derived type
use DVectors ! see below for its definition
implicit none
public :: abstract_operator, application
type, abstract :: abstract_operator
!! Abstract type defining base general operator
!! Members shared by any operator:
!! Number of columns (the domain of the application)
integer :: ncol=0
!! Number of rows (the co-domain of the application)
integer :: nrow=0
!! Logical flag for symmetry
logical :: is_symmetric=.false.
!! Real estimating the nonlinearity of the application, defined
!! as: nonlinearity = |OP(x)+OP(y)-OP(x+y)|/(max(|x|,|y|,|x+y|)
!! It is zero (the default value) for linear operators (up to
!! machine precision). It may be greater that zero for quasi
!! linear operator like inverse operator.
double precision :: nonlinearity = 0.0d0
contains
!! Procedure for the application of the operator to a vector.
procedure(application), deferred :: apply
end type abstract_operator
abstract interface
!! Abstract procedure defining the interface for a general
!! application actions
!!
!! y = M(x) = ( M*x if M is linear)
subroutine application(this,vec_in,vec_out,ierr)
import abstract_operator
import dvector
implicit none
!! The class with all quantities describing the operator
class(abstract_operator), intent(inout) :: this
!! The input" vector
class(dvector), intent(in ) :: vec_in
!! The output vector y= "operator apllied on" x
class(dvector), 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
end interface
end module AbstractOperatorsIt is aimed to be be as general as possible to include linear and nonlinear operators. The latter may include inexact solvers (like Krylov or stationary solvers) that may be used as a preconditioner in the FlexiblePCG and FlaxibleGMRES algorithms.
For Krylov methods such as PCG, BICGSTAB, etc. , we will basically need to define two vector-vector operations: 1 - y=ax+by were a,b are real and x,y are vectors. 2 - Scalar product of two vectors.
Now, the most important decision becomes how to define a structure describing vectors, those appearing in the definition of the abstract procedure "apply". The design of this structure is crucial to define concrete operators that will extend the abstract class abstract_operator. When we will need to define the action of any operator, it is fundamental to have access to the vector entries, either reading and writing them. Once the structure is defined, it should be easy to code the basic linear algebra operations mentioned above can be easily built.
My first proposal of structure describing real vectors is a derived type containing the basic array info and vector coefficients. I introduced a constructor, a destructor, and the linear algebra procedure mentioned above. I also introduced two procedures, "set" and "get", to define the interface of reading and writing the vector coefficient.
module DVectors
private
public :: dvector
public :: axpby_dvector
public :: dot_dvector
public :: norm_dvector
!! simple real vector
type :: dvector
!! vector length
integer :: size
!! vector coefficients
double precision, allocatable :: coefficients(:)
contains
!! constructor
procedure, public, pass :: init => init_dvector
!! destructor
procedure, public, pass :: kill => kill_dvector
!! axpby procedure
procedure, public, nopass :: axpby => axpby_dvector
!! get procedure
procedure, public, pass :: get => get_dvector
!! set procedure
procedure, public, pass :: set => set_dvector
end type dvector
!procedure, public :: scalar_product => scalar_product_dvector
contains
!! constructor
subroutine init_dvector(this,size)
implicit none
class(dvector), intent(inout) :: this
integer, intent(in ) :: size
this%size=size
allocate(this%coefficients(size))
end subroutine init_dvector
!! destructor
subroutine kill_dvector(this)
implicit none
class(dvector), intent(inout) :: this
this%size=0
deallocate(this%coefficients)
end subroutine kill_dvector
subroutine axpby_dvector(alpha,x,beta,y)
use Precision
implicit none
double precision, intent(in ) :: alpha
class(dvector), intent(in ) :: x
double precision, intent(in ) :: beta
class(dvector), intent(inout) :: y
!! not optimized
y%coefficients=alpha*x%coefficients+beta*y%coefficients
end subroutine axpby_dvector
!! get one entry procedure
function get_dvector(x, index) result(res)
implicit none
!!
class(dvector), intent(in ) :: x
integer, intent(in ) :: index
double precision :: res
res = x%coefficients(index)
end function get_dvector
!! set one entry procedure
subroutine set_dvector(x, index, value)
implicit none
!!
class(dvector), intent(inout) :: x
integer, intent(in ) :: index
double precision, intent(in ) :: value
x%coefficients(index)=value
end subroutine set_dvector
!! scalar product of two vectors
function dot_dvector(x,y) result(res)
implicit none
class(dvector), intent(in ) :: x
class(dvector), intent(in ) :: y
double precision :: res
! local
integer :: i
res=0.0d0
do i=1,x%size
res=res+x%coefficients(i)*y%coefficients(i)
end do
end function dot_dvector
!! vector norm
function norm_dvector(x) result(res)
implicit none
class(dvector), intent(in ) :: x
double precision :: res
res=sqrt(dot_dvector(x,x))
end function norm_dvector
end module DVectorsPros of this choice:
This is more or less the choice done by the PSBLAS library and the Trilinos Library
An alternative could be to use get standard Fortran arrays or coarrays as @ did. PROS:
It is possible to see positive and negative aspects of both choices by looking at following operators.
module IdentityOperators
use DVectors
use AbstractOperators
implicit none
!! identity operator
type, public, extends(abstract_operator) :: eye
contains
!! constructor
procedure, public, pass :: init => init_eye
!! procedure to be overrided
procedure, public, pass :: apply => apply_eye
end type eye
contains
!! set identity dimensions
subroutine init_eye(this,size)
implicit none
class(eye), intent(inout) :: this
integer, intent(in ) :: size
!! the set domain/codomain operator dimension
this%nrow = size
this%ncol = size
this%is_symmetric=.False.
this%nonlinearity=0.d0
end subroutine init_eye
!! simply copy x into y
subroutine apply_eye(this,vec_in,vec_out,ierr)
use DVectors
implicit none
class(eye), intent(inout) :: this
class(dvector), intent(in ) :: vec_in
class(dvector), intent(inout) :: vec_out
integer, intent(inout) :: ierr
! using basic linear algebra operation defined in Dvectors module
call axpby_dvector(1.0d0,vec_in,0.0d0,vec_out)
! or we can simply use a procedure assocaite with standard Fortran vector
!vec_out%coefficinets=vec_in%coefficients
ierr=0
end subroutine apply_eye
end module IdentityOperatorsmodule ExampleOperators
use Precision
use AbstractOperators
use DVectors
private
!! Derived type describing the 1d-Laplacian
!! with periodic boundary conditions
type, public, extends(abstract_operator) :: Laplacian
contains
procedure, public, pass :: init => init_LP
!! procedure overrided
procedure, public, pass :: apply => apply_LP
end type Laplacian
!! Derived type describing the 1d-Laplacian
!! with periodic boundary conditions
type, public, extends(abstract_operator) :: DiagonalMatrix
double precision, allocatable :: coefficients(:)
contains
procedure, public, pass :: init => init_DM
procedure, public, pass :: kill => kill_DM
!! procedure overrided
procedure, public, pass :: apply => apply_DM
end type DiagonalMatrix
contains
!! constructor
subroutine init_LP(this, nnode)
implicit none
class(Laplacian), intent(inout) :: this
integer, intent(in ) :: nnode
!! the set domain/codomain operator dimension
this%nrow = nnode
this%ncol = nnode
this%is_symmetric=.True.
end subroutine init_LP
!! apply tridiagonal operator
subroutine apply_LP(this,vec_in,vec_out,ierr)
implicit none
class(Laplacian), intent(inout) :: this
class(dvector), intent(in ) :: vec_in
class(dvector), intent(inout) :: vec_out
integer, intent(inout) :: ierr
!local
double precision :: sign=1.0d0, val
integer :: i_row
!! first row
i_row=1
val = 2*vec_in%get(i_row) + sign* vec_in%get(i_row+1)
call vec_out%set(i_row,val)
!! middle rowsuse DVectors
do i_row=2,this%nrow-1
val = sign * vec_in%get(i_row-1) + 2*vec_in%get(i_row) + sign* vec_in%get(i_row+1)
call vec_out%set(i_row,val)
end do
!! last row
i_row=this%nrow
val = sign* vec_in%get(i_row-1) + 2*vec_in%get(i_row)
call vec_out%set(i_row,val)
ierr=0
end subroutine apply_LP
!! constructor
subroutine init_DM(this, size)
class(DiagonalMatrix), intent(inout) :: this
integer, intent(in ) :: size
this%nrow = size
this%ncol = size
this%is_symmetric=.True.
allocate(this%coefficients(size))
end subroutine init_DM
!! destructor
subroutine kill_DM(this)
class(DiagonalMatrix), intent(inout) :: this
this%nrow = 0
this%ncol = 0
this%is_symmetric=.False.
deallocate(this%coefficients)
end subroutine kill_DM
!! apply diagonal matrix
subroutine apply_DM(this,vec_in,vec_out,ierr)
implicit none
class(DiagonalMatrix), intent(inout) :: this
class(dvector), intent(in ) :: vec_in
class(dvector), intent(inout) :: vec_out
integer, intent(inout) :: ierr
!local
integer :: i_row
double precision :: val
do i_row=1,this%nrow
val=vec_in%get(i_row)*this%coefficients(i_row)
call vec_out%set(i_row,val)
end do
ierr=0
end subroutine apply_DM
end module ExampleOperatorsNow, the interface for the PCG solver is given by
subroutine pcg(matrix,rhs,sol,ierr,&
tolerance,max_iterations,preconditioner)
use DVectors
use AbstractOperators
use IdentityOperators
implicit none
class(abstract_operator), intent(inout) :: matrix
class(dvector), intent(in ) :: rhs
class(dvector), intent(inout) :: sol
integer, intent(inout) :: ierr
double precision, intent(in ) :: tolerance
integer, intent(in ) :: max_iterations
class(abstract_operator), target,optional, intent(inout) :: preconditioner
!
! local vars
!
! logical
logical :: exit_test
! string
character(len=256) :: msg
! integer
integer :: i
integer :: lun_err, lun_out
integer :: iprt,max_iter
integer :: nequ,dim_ker,ndir=0
integer :: iort
integer :: info_prec
integer :: iter
! real
class(dvector), pointer :: aux(:)
type(dvector), target, allocatable :: aux_local(:)
class(dvector), pointer :: axp, pres, resid, pk,scr
double precision :: alpha, beta, presnorm,resnorm,bnorm
double precision :: ptap
double precision :: normres
double precision :: tol,rort
double precision :: inverse_residum_weight
type(eye),target :: identity
class(abstract_operator), pointer :: prec
!! checks
if (.not. matrix%is_symmetric) then
ierr = 998
return
end if
if (matrix%nrow .ne. matrix%ncol ) then
ierr = 999
return
end if
nequ=matrix%nrow
tol=tolerance
max_iter=max_iterations
!! handle preconditing
if (present(preconditioner)) then
if (.not. preconditioner%is_symmetric) then
ierr = 999
return
end if
if ( preconditioner%nrow .ne. preconditioner%ncol ) then
ierr = 999
return
end if
prec => preconditioner
else
call identity%init(nequ)
prec => identity
end if
! allocate memory if required and set pointers for scratch arrays
allocate(aux_local(5))
do i=1,5
call aux_local(i)%init(nequ)
end do
aux=>aux_local
axp => aux(1)
pres => aux(2)
resid => aux(3)
pk => aux(4)
scr => aux(5)
! set tolerance
tol = tolerance
! set max iterations number
max_iter = max_iterations
exit_test = .false.
! compute rhs norm and return zero solution
bnorm = norm_dvector(rhs)
if (bnorm<1.0d-16) then
! set solution and flag
ierr=0
do i=1,nequ
call sol%set(i,0.0d0)
end do
! free memory
aux=>null()
if (allocated(aux_local))then
do i=1,6
call aux_local(i)%kill()
end do
deallocate(aux_local)
end if
return
end if
! calculate initial residual (res = rhs-M*sol)
call matrix%apply(sol,resid,ierr)
call axpby_dvector(1.0d0,rhs,-1.0d0,resid)
resnorm = norm_dvector(resid)/bnorm
!
! cycle
!
iter=0
do while (.not. exit_test)
iter = iter + 1
! compute pres = PREC (r_{k+1})
call prec%apply(resid,pres,ierr)
! calculates beta_k
if (iter.eq.1) then
beta = 0.0d0
else
beta = -dot_dvector(pres,axp)/ptap
end if
! calculate p_{k+1}:=pres_{k}+beta_{k}*p_{k}
call axpby_dvector(1.0d0,pres,beta,pk)
! calculates axp_{k+1}:= matrix * p_{k+1}
call matrix%apply(pk,axp,ierr)
! calculates \alpha_k
ptap = dot_dvector(pk,axp)
alpha = dot_dvector(pk,resid)/ptap
! calculate x_k+1 and r_k+1
call axpby_dvector(alpha,pk,1.0d0,sol)
call axpby_dvector(-alpha,axp,1.0d0,resid)
! compute residum
resnorm = norm_dvector(resid)/bnorm
! checks
exit_test = (&
iter .gt. max_iter .or.&
resnorm .le. tol)
if (iter.ge. max_iter) ierr = 1
end do
! free memory
aux=>null()
if (allocated(aux_local))then
do i=1,5
call aux_local(i)%kill()
end do
deallocate(aux_local)
end if
end subroutine pcgHere below an example program.
program main
use AbstractOperators
use DVectors
use ExampleOperators
use KrylovSolver
implicit none
! problem vars
type(Laplacian) :: laplacian_1d
type(DiagonalMatrix) :: inverse_diagonal_laplacian
type(dvector) :: sol
type(dvector) :: rhs
integer :: n
integer :: max_iterations=400
double precision :: tolerance=1.0d-5,h,sign
! local vars
integer :: ierr,i,m,j
type(dvector) :: res
n=100
h=1.0d0/(n-1)
!! init matrix
call laplacian_1d%init(n)
!! allocate solution and rhs and set the latter
call sol%init(n)
call rhs%init(n)
call res%init(n)
do i=1,n
call rhs%set(i,0.d0)
end do
call rhs%set(1,1.0d0)
call rhs%set(n,1.0d0)
do i=1,n
call sol%set(i,0.d0)
end do
!! init diagonal and set values with inverse of
!! diagonal of the Laplacian
call inverse_diagonal_laplacian%init(n)
inverse_diagonal_laplacian%coefficients=1.0d0/(2.0d0)
write(*,*) norm_dvector(rhs)
call pcg(laplacian_1d,rhs,sol,ierr,&
tolerance=1.0d-3,&
max_iterations=400)!,&
!preconditi1.0d0r=inverse_diagonal_laplacian)
!! check results
call laplacian_1d%apply(sol,res,ierr)
call axpby_dvector(-1.0d0,rhs,1.0d0,res)
write(*,*) 'Res=',norm_dvector(res)
end program main
Motivation
I could not find in the Issues or Linear Algebra section any reference to include standard iterative Krylov solvers (e.g., PCG, BICGSTAB, GMRES) in the stdlib library. I believe that it would be great being able to call this type of algorithms like in Matlab, Scipy, (something like the Matlab Preconditioned Conjugate Gradient
x = pcg(A,b,tol,maxit,M1,M2,x0)passing any type of linear operator and/or preconditioner.Prior Art
A great source of inspiration on this topic is the PSBLAS library . Nevertheless but I believe that, using more abstract classes and procedure overriding, it would be nice to build a flexible high level framework. This high level structure would not depend on the actual implementation of matrix vector product and preconditioner action called by krylov solvers. This would close the gap, in terms of ease of use and flexibility, with other implementations of these solver in Matlab, Python or C++.
Additional Information
No response