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

  • 1
    $\begingroup$ 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 $\endgroup$
    – Lukas Lang
    Commented May 28, 2018 at 22:02
  • $\begingroup$ thanks for telling me, i'll do that now $\endgroup$
    – Chalky
    Commented May 28, 2018 at 22:09
  • 1
    $\begingroup$ 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. $\endgroup$ Commented May 28, 2018 at 22:15
  • $\begingroup$ 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. $\endgroup$ Commented May 28, 2018 at 22:16
  • $\begingroup$ 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) $\endgroup$
    – Chalky
    Commented May 28, 2018 at 22:24

1 Answer 1


Turning the comments into an answer.

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


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


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]]



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]]



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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.