# 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:** 1

<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 8, 2024, 2:24pm UTC](https://fortran-lang.discourse.group/t/generic-or-kind-agnostic-linear-system-solvers/7574/1 "2024-03-08T14:24:46Z")

</div>

Dear All, today I would like to have some clarifications about handling generic procedures in Fortran without re-coding the same procedure for each set of new data types. I want to apply this to a very specific problem, that of solving a linear system with a Krylov method such as the pcg. The following code is taken from [here](https://github.com/fortran-lang/stdlib/issues/628) where Dr. Enrico Facca kindly provided a template for a Krylov API. As you can see, the input field `matrix` is of derived type `AbstractOperator` and the vectors are of type `Vectors`.

```fortran
 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(abstract_operator), target,optional, intent(inout) :: preconditioner
    [...]
    !
    ! 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 ) )
       [...]
    end do

    [...]

  end subroutine pcg

```

Since this is just a combination of matrix-vector products, scalar products and vector updates, it would be interesting if there was a way to make it “kind agnostic”. Considering only for simplicity different kind of matrices, it would be amazing if the pcg could be applied to a tabular matrix, a CSR matrix, a COO matrix, a CSC matrix or even a user defined data structure, without redefining the pcg for each type of matrix. This last possibility is, IMHO, the key for user-defined parallelization, because you would need to define the data structure and its communication procedures and simply provide it to the pcg. As an example, imagine to solve the same linear system twice: the first time you provide the full matrix to the pcg as it is, the second time you have a domain decomposition procedure that imposes a block structure to the matrix. You might want to define the block structure as a derived data type, define a matvet function handling correctly the communication and use type-bound procedures to apply it, give this data structure to the pcg. This idea could be taken even further, because the concept could be extended to the data type `DVectors` to provide, as an example, an openMP `DVector` and an OpenACC `DVector` and so on.

Let’s consider the following case: there is one module implementing the pcg; two modules containing predefined matrix representation, let’s say tabular and CSR; the last module, defines a user-defined matrix representation. I would like the pcg to accept any of the three matrix representations, so it should be in some way kind-agnostic. The tabular matrix and CSR, as well as for the pcg, could be coming from some library (hopefully one day stdlib) that the user should not be recoding or modifying directly. Finally, the last module containing the user defined data type should be the only one directly accessible to the user, and yet perfectly interchangeable with the other two.

What would be your preferred approach to do this? Any suggestion for further reading is well accepted.  
This could also be something interesting to develop in stdlib, especially if you could really have a set of solvers that only need to know the type of the matrices and vectors.

Feel free to ask for clarifications if you find this post not clear enough. I hope to come back during the weekend with a MWE.

---

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