Calling cuSolverSp from CUDA FORTRAN

Greetings to all people here!

Recently, I tried to search/find a way to solve large sparse linear Ax=b systems using CUDA library routines.

The sparse linear system is part of a Computational Fluid Dynamics solver and can be quite large.
For a modest computational mesh of 5million nodes, and by having 5 Navier-Stokes equations for each node, a 25mil x 25mil A matrix is formed.
Thus, direct and dense solutions are out of the question.
Matrix A is always square and the diagonal is always non-zero.

The current solution process for the linear system is a good old iterative Jacobi kernel (corresponding to computational nodes or cells in the CFD solver), which is easy to implement and debug, but has kinda bad convergence properties.
I was searching for something better.

I stumbled upon the cuSolverSp library, but I cannot seem to find any CUDA-Fortran examples to help me implement/adjust it for my case and setup.
Moreover, when I include -cudalib=cusolver in my Makefile, the compilation complains that it cannot find the relevant module.
For context, I already utilize -cudalib=cublas for parallel reduction and everything compiles and runs fine.

Is cuSolver available for CUDA-Fortran?
Is there a sparse linear system example I can use?

I am familiar with CSR format, since I also utilize it to solve the linear system using the PETSc library.
However, I am not entirely familiar with the CUDA handles and the corresponding API and I was hoping for a simple example that can guide me.

The example folders that come bundled with the HPC-SDK do not have any cuSolver examples for CUDA-Fortran.

What options do I have?

Best regards,
Konstantinos

According to the docs, we don’t have Fortran interfaces for cusolverSp. We can work to fix that but it would be helpful to know exactly the functions you need so we can prioritize.

Math Libs examples can be found here, although they are primarily C/C++.

Well, I was hoping that the Fortran interfaces were there and that I had somehow missed them.
Based on the docs (which I had already read), it indeed seemed like interfaces for Fortran were not available.

I am interested in any routine that can solve a Ax=b linear system with A being square and non singular.

From the available interfaces, the following seem suitable:

  • cusolverSpScsrlsvlu
  • cusolverSpScsrlsvqr

since they only require A to be square and positive definite.

Thank you for the feedback. I don’t have an ETA, but in the meantime you could this guide to write your own interface.