# Generic or Kind-Agnostic linear system solvers

**URL:** <https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574>\
**Category:** Help\
**Created:** [March 8, 2024, 2:24pm UTC](https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574 "2024-03-08T14:24:45Z")\
**Posts on this page:** 1\
**Showing post:** 29

<div class="post-metadata">

**Author:** ![kimala](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/kimala/32/2707_2.png) [@kimala](https://fortran-lang.discourse.group/u/kimala)\
**Post date:** [March 3, 2025, 10:34pm UTC](https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574/29 "2025-03-03T22:34:26Z")

</div>

> [@stdlib\_LinAlgOperators: Linear Algebra Operators](https://fortran-lang.discourse.group/t/stdlib-linalgoperators-linear-algebra-operators/9286/5):
>
> > [@stdlib\_LinAlgOperators: Linear Algebra Operators](https://fortran-lang.discourse.group/t/stdlib-linalgoperators-linear-algebra-operators/9286/4):
> >
> > What would be the other paths one could go in modern fortran to accomplish what you’re doing ?
> 
> Besides type extension, the methods available for injecting behavior are:
> 
> - **Procedure arguments** (as in SLATEC SLAP, e.g. [`dcg.f`](https://www.netlib.org/slatec/lin/dcg.f))
> - **Procedure pointers** (packed in a derived type for tidyness)
> - **Reverse communication** (just like ARPACK, the venerable [templates](https://www.netlib.org/templates/), or oneMKL [RCI ISS](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2024-0/iterative-sp-solvers-reverse-comm-iface-rci-iss.html))
> 
> By means of internal procedures, it’s always possible to convert an interface using one approach into a different one. Internal procedures have some weakness (indirection, executable stack) as explained here: [Is creating nested subroutines/functions considered good practice in Fortran? - #25 by szakharin](https://fortran-lang.discourse.group/t/is-creating-nested-subroutines-functions-considered-good-practice-in-fortran/6545/25). For realistic problems I’d expect the linear operator `apply` method to dominate the runtime.

> [@stdlib\_LinAlgOperators: Linear Algebra Operators](https://fortran-lang.discourse.group/t/stdlib-linalgoperators-linear-algebra-operators/9286/8):
>
> This is also _the_ spot where you can, and want to, introduce optimizations, which IMHO advocates for simplicity of interfacing with the libraries, but I have to admit that the [monkey in my brain plays the cymbals](https://youtu.be/H3KpR6EyZX0?t=47) when I try to understand reverse communication.
> 
> When extending a derived data type (if I’m not mistaken) you can introduce additional dependencies. This means that I could extend my base `absOperator_dp` and add the module mpi\_f08 and create an `mpiOperator_dp`. In the same way I could extend `Operator_dp` with openMP, or openACC and so on. Moreover, these extension are standalone from the library and thus easily testable outside. Could reverse communication be as simple/user friendly as this?

> [@stdlib\_LinAlgOperators: Linear Algebra Operators](https://fortran-lang.discourse.group/t/stdlib-linalgoperators-linear-algebra-operators/9286/9):
>
> IMO it is just a matter of choosing your own poison. Some people avoid type extension and deferred bindings, others don’t like reverse communication. I think @victotronics [has described](https://web.archive.org/web/20220121212952/http://computer-programming-forum.com/49-fortran/cd933b17b5b4592e.htm) reverse communication fairly when he says:
> 
> > The idea is to make the function data-independent by keeping  
> > the data entirely outside its scope. Thus for every operation on that data  
> > you have to return, with a request parameter indicating what  
> > operation has to be performed outside, and you have to have a way  
> > to indicate where you resume after the operation is done.
> > 
> > I’ve written code like this, and **it’s kinda neat, but I’m not**  
> > **entirely unconvinced that it’s not totally evil. Or not.**  
> > [emphasis added]

@ivanpribec I’m moving this part of the conversation here because it’s more relevant here. If I understood correctly from 30 minutes of research and [this page](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-fortran/2025-0/cg-interface-description.html) reverse communication wouldn’t be a big deal to implement on the library side. An first test is the following:

```fortran
module stdlib_krylov_PCG
  use stdlib_kinds
  
  implicit none
  private
  !! Public methods
  public :: stdlib_linalg_cg
  !!
  interface stdlib_linalg_cg
    !!
    module procedure stdlib_krylov_CG
    module procedure stdlib_krylov_CG_reverse
    !!
  end interface
  !!
contains !!
  subroutine stdlib_krylov_CG_reverse( n, wrk, x, info, itnum )
    !!
    integer, intent(in ) :: n
    real(dp), intent(inout) :: wrk(n,3) ! or 4 for preconditioning
    real(dp), intent(inout) :: x(n)
    integer, intent( out) :: info
    integer, intent(inout) :: itnum
    !!
    real(dp) :: alpha, beta
    real(dp) :: p_dot_Ap
    real(dp), save :: r_dot_r
    real(dp), save :: r_dot_r_old
    !!
    ! Mapping (maybe associate?)
    ! p = wrk(:,1)
    ! Ap = wrk(:,2)
    ! r = wrk(:,3)
    ! Kr = wrk(:,4)

    !cg_iterations: do itnum = 1, options%itmax
      !
      if (info == 22) then
        info = 22 !request Matrix-Vector multiplication
        return
      end if 
      ! Ap = matmul( A, p ) ! Performed outside
      !
      r_dot_r = dot_product( wrk(:,3), wrk(:,3) )  
      
      if (info == 11) then ! the matvec product has done
        p_dot_Ap = dot_product( wrk(:,1), wrk(:,2) )
        alpha = r_dot_r / p_dot_Ap
        x = x + (alpha) * wrk(:,1) 
        wrk(:,3) = wrk(:,3) + (-alpha) * wrk(:,2) 
        r_dot_r_old = r_dot_r
        ! Time for test
        info = 21
        return
      end if 

      !if ( sqrt(r_dot_r) < options%tol_dp ) then ! Performed outside
      ! info = 0
      ! exit cg_iterations
      !end if

      if (info == 12) then
        r_dot_r = dot_product( wrk(:,3), wrk(:,3) )
        beta = r_dot_r/r_dot_r_old
        wrk(:,1) = beta * wrk(:,1) + 1._dp * wrk(:,3)
        itnum = itnum + 1
        ! Time to go back to Ap
        info = 22
      end if

    !end do cg_iterations

    return

  end subroutine stdlib_krylov_CG_reverse
  !!
  subroutine stdlib_krylov_CG( A, rhs, x, xzero, info )
    !!
    real(dp), intent(in ) :: A(:,:) 
    real(dp), intent(in ) :: rhs(:)
    real(dp), intent(inout) :: x(:)
    logical, intent(in ), optional :: xzero
    integer, intent( out) :: info
    !!
    !! Default options
    real(dp), allocatable :: r(:)
    real(dp), allocatable :: p(:)
    real(dp), allocatable :: Ap(:)
    !! Residual, search direction and other vectors.
    real(dp) :: alpha, beta
    real(dp) :: p_dot_Ap
    real(dp) :: r_dot_r
    real(dp) :: r_dot_r_old
    integer :: itnum
    !!
    allocate (Ap, source = rhs)
    Ap = 0._dp
    allocate (r, source = rhs)

    if ( present(xzero) .and. xzero ) then 
      Ap = matmul( A, x )
      r = 1._dp * r + (-1._dp) * Ap 
    else
      x = 0._dp
    end if
 
    allocate (p, source = r)
    
    r_dot_r = dot_product( r, r )  
    info = 1
    
    cg_iterations: do itnum = 1, 10
      
      Ap = matmul( A, p )
      p_dot_Ap = dot_product( p, Ap )
      alpha = r_dot_r / p_dot_Ap
      x = x + (alpha) * p 
      r = r + (-alpha) * Ap 
      r_dot_r_old = r_dot_r
      r_dot_r = dot_product( r, r )

      if ( sqrt(r_dot_r) < 1.e-10 ) then
        info = 0
        exit cg_iterations
      end if

      beta = r_dot_r/r_dot_r_old
      p = beta * p + 1._dp * r 

    end do cg_iterations

    deallocate (r, p, Ap)

    return

  end subroutine stdlib_krylov_CG
  !!
end module stdlib_krylov_PCG

program main
  use stdlib_kinds
  use stdlib_krylov_PCG
  implicit none

  integer :: i,j, n
  real(dp), allocatable :: A(:,:)
  real(dp), allocatable :: L(:,:)
  real(dp), allocatable :: U(:,:)
  real(dp), allocatable :: b(:), x(:)
  real(dp), allocatable :: work(:,:), ones(:)
  integer :: info
  integer :: itnum

  n = 10
  allocate( A(n,n), L(n,n), U(n,n), b(n), ones(n), x(n), work(n,3) )
  L = 0._dp
  do j = 1, n
    do i = j, n
      call random_number(L(i,j))
    end do
    L(j,j) = 3._dp + 3._dp*L(j,j)! sqrt(real(n))
  end do
  U = transpose(L)
  A = matmul( L, transpose(L) )
  ones = 1._dp
  b = matmul( A, ones)

  call stdlib_linalg_cg(A, b, x, .false., info )
  print *, x
  
  x = 0._dp
  ! This in a subroutine init
  work(:,1) = b
  work(:,2) = 0._dp
  work(:,3) = b
  info = 22

  cg_iteration: do i = 1, 3*10
  
    call stdlib_linalg_cg( n, work, x, info, itnum )
    if (info == 22) then 
       work(:,2) = matmul( A, work(:,1) )
       info = 11
    end if

    if (info == 21) then 
       if ( sqrt(dot_product(work(:,3),work(:,3))) < 1.e-10 ) then
         info = 0
         exit cg_iteration
       else 
         info = 12
       end if
    end if

  end do cg_iteration
 
  print *, x

end program

```

I see the value in this kind of implementation (there is the standard cg there for reference), but I don’t see any benefit in terms of code simplicity to be honest. I will test performance though.

---

_[View the full topic](https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574)._