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 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.
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.
there are definitely a lot of potential areas for both research and high-performance applications!