Will this project include iterative Krylov solvers? Or an high-level framework for them?

[enricofacca]

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

[jvdp1]

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]

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]

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 AbstractOperators

It 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 DVectors

Pros of this choice:

• it is possible to "hide" complex and efficient procedures from the user.
• structures for parallellization can be hidden as well within the type.
  Cons:
• user are forced to used this type and its procedure to handle and manipulate the vectors. This reduces code readability can be partially.
  However, this issue can be partially solved declaring operators such as +,-,*,/

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:

• we only need to define efficient basic linear algebra operations (scalar product, norm). "Set" and "get" procedures are already defined, easy to use, and super easy to read.
  CONS:
  Structures handling parallelization need to put elsewhere

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 IdentityOperators

module 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 ExampleOperators

Now, 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 pcg

Here 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