CUDA for linear equations

TLDR; Is there a way to solve linear equations of a sparse matrix (discretized laplace operator) efficiently using CUDALink in Mathematica? I didn't find a CUDALinearSolve or CUDAMatrixInverse or something in the documentation. This is my first time working with CUDA link.

The long version: Our prof gave us a project to solve a Poisson equation (laplace Phi = 4pi * rho) in 2D. I discretized the Laplacian to an n² x n² matrix given by ({-4 for i=j, 1 for |i-j|=1 or |i-j|=n}). (*)

I get solutions in reasonable time for up to n=100, using LinearSolve[N[laplace],rho], but our programs are judged by performance and this isn't really too well performing....

Unfortunately we were specifically tasked to only use Java, Mathematica or Python, I don't know python and well, Java is Java, so Mathematica is really my only shot. But since I am basically throwing around matrices all day I though my a GPU might be better equipped to handle this and Mathematica has CUDA integration...

I looked through the CUDALink documentation and didn't find anything helpful, which is super surprising to me, since this should be a fairly standard use-case of CUDA. Does anyone have an idea what I could use?

I found cusolver in the CUDA documentation and know that CUDAFunctionImport is a thing, but I have no clue how to bring these components together

(*) - There is something weird with this discretization though, as my results seem to represent the surface of a cylinder not a flat plane, but this is besides the point

• If you have something to add to your question, please do so by editing the question instead of adding a comment - the question should give anyone reading all necessary information without having to read through the comments – Lukas Lang May 28 '18 at 22:02
• thanks for telling me, i'll do that now – user3271145 May 28 '18 at 22:09
• 10000 x 10000 is a tiny matrix is you use SparseArray. Try also LinearSolve with the option 'Method->"Pardiso". On my machine, "Pardiso" factorizes the 5-stencil Laplacian of a $1120 \times 1120$ in about 4.2 seconds while the standard solver for sparse matrices ('Method->"Pardiso") needs about 14.8 seconds for the same task. – Henrik Schumacher May 28 '18 at 22:15
• Making that work on a GPU is highly nontrivial and as far as I can tell, no reasonably working linear solvers for the GPU get shipped with Mathematica. – Henrik Schumacher May 28 '18 at 22:16
• Thanks for the tip, I will try... n=100 worked with the naivest possible code for this problem, I hope to get it up to n=1000, which seems like a tall order (and be it just for RAM, even building the laplace matrix makes me run out of RAM) – user3271145 May 28 '18 at 22:24

Making a linear solver work on a GPU is highly nontrivial (it's a task to be assigned to someone who poisened his mother and father), and as far as I can tell, no reasonably working linear solvers for the GPU are shipped with Mathematica. But one can work reasonably well with the built-in tools for SparseArrays.

The 5-stencil Laplacian is the graph Laplacian of the underlying grid graph. We can obtain it as follows.

n = 1000;
A = 1. AdjacencyMatrix[GridGraph[{n, n}]] -
4. IdentityMatrix[n^2, SparseArray]; // AbsoluteTiming // First

0.190387

This uses the Pardiso solver from the Intel MKL to factorize the matrix and to create a LinearSolveFunction object that we can use to perform multiple solves

S = LinearSolve[A, "Method" -> "Pardiso"]; // AbsoluteTiming // First

1.90579

Creating a right-hand side for the equation and applying the LinearSolveFunction S

b = RandomReal[{-1, 1}, n^2];
x = S[b]; // AbsoluteTiming // First
Max[Abs[A.x - b]]

0.285522

1.03473*10^-13

Creating many right-hand sides and solving the equations for all of them at once

B = RandomReal[{-1, 1}, {n^2, 100}];
X = S[B]; // AbsoluteTiming // First
Max[Abs[A.X - B]]

8.24018

1.03473*10^-13

Now you might wonder why the latter needs considerably less than 100 times the time as the former. The reason for that is that S[b] and S[B] involve so-called forward- and backward-substitutions which are not parallelizable. So S[b] runs at single-core speed, while S[B] can utilize my CPU's 4 compute cores in parallel.