Factorize and solve the Laplacian on the GPU

[lukem12345]

PR #75 added GPU support to CombinatorialSpaces.jl via the CUDA.jl API.

Currently, we directly offer GPU support for the classic DEC primitive operators: $d, \star$, and $\wedge^{pp}$. And Decapodes.jl is able to automatically create bindings for higher-order operators on the GPU - such as the Laplacian - by expanding their definitions - e.g. $d\star d\star$ - and pre-multiplying those matrices. The branch llm/cuda-wedge-music is adding GPU support for $\wedge_{10}^{p d}, \wedge^{dp}_{01}$, and the interpolating musical operator $\flat\sharp$. These operators are implemented as matrix-multiplications are relatively straight-forward to port to the GPU.

However, some operators - such as the geometric Hodge star for 1-forms - are instead implemented by solving a matrix implementing that operator. $\star^{-1}$ was implemented by using gmres from Krylov.jl, and performantly solves this system, since it is "mostly" diagonal.

However, when we need to solve the heat equation by solving the sparse $\Delta_0$ matrix, gmres is not able to provide adequate performance. Further, even the factorization of this matrix on the GPU is prohibitively-expensive.

So, we need to provide a canonical way of solving this problem quickly in the CombinatorialSpaces library. The current prototype factors the matrix on the CPU, sends its components to the GPU, and in-houses an LU-solve.

We are also experimenting with approximations to the Laplacian on well-structured meshes that can be exploited for superior factorizations.

[GeorgeR227]

Looking at the CUDA.jl source I found that they actually do implement in-place sparse solvers on GPU with QR and Cholesky. They also have a eigenvalue solver. These are called csrlsvqr!, csrlsvchol! and csreigvsi respectively.

They seem to have sparse factorization support but I'm not sure how to use those to solve.

Link here: https://github.com/JuliaGPU/CUDA.jl/blob/master/lib/cusolver/sparse.jl#L101