Is there a way to do this natively in Mathematica? I read in several places that this is possible with LibraryLink, but I don't have a lot of experience with C and so I'm having trouble implementing this approach.

I came across the documentation for Affine Covariant Newton Method where it states

When an iterative method like "Krylov" is specified as a method to LinearSolve, then the factorization of the Jacobian is not stored when possible, but computed anew every time it is needed.

which seems to indicate that it stores the factorization when using a direct method.

  • 2
    $\begingroup$ I'll only note that if you try to inspect the result of something like LinearSolve[mat, Method -> "Pardiso"] or LinearSolve[m, Method -> "Krylov"], the input matrix mat is stored along with a SparseArray`SparseLUDecomposition[] header in the former, and just the matrix itself in the latter, but with flags indicating which sparse method was used. $\endgroup$
    – J. M.'s torpor
    Apr 19 '20 at 6:50

No, this is not built-in and I have been waiting quite long for such a feature.


With the help of Szabolcs's package "LTemplate`" I had written such an interface for the Intel MKL Pardiso that is shipped with Mathematica a couple of years ago. I spent all day to bring it into a shape that has at least a small chance to run on a computer different from mine. You can find it on github. So far, the package supports only matrices of real numbers.

I have to say that I am a really, really bad C++ programmer, so I give no guarantees for anything. Use with caution. I tested the package only on macos 10.15.4 with Mathematica 12.0 and the clang compiler and (only very briefly) with Mathematica 11.3 on Linux. I have absolutely no idea about Windows.


I don't have any fancy installer, yet. Just download the package from github, unzip it, and drop it into the folder FileNameJoin[{$UserBaseDirectory, "Applications"}].

You might have to adjust the values of "CompileOptions" and "LinkerOptions" in the file BuildSettings.m in the package's base directory. You can get the appropriate values for your system by answering a couple of questions asked by this tool: https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/.

If you find any problems, just contact me. Unfortunately, I do not have much time to maintain the package. Actually, my hope is that somebody at Wolfram Research will implement this feature in the near future. ;)


Okay, here is a simple usage example:

n = 1000000;
m = 20;
L = SparseArray[{Band[{1, 1}] -> -2. n, Band[{2, 1}] -> 1. n, Band[{1, 2}] -> 1. n}, {n, n}, 0.];
L[[1, 1]] = 1.;
L[[1, 2]] = 0.;
L[[-2, 1]] = 0.;
L[[n, n]] = 1.;
(*The following guarantees that L is CSR conform.*)
L = SparseArray[L];
b = ConstantArray[1./n, n];
B = RandomReal[{-1, 1}, {n, m}];
b[[{-1, 1}]] = {0., 0.};

So, we have got a matrix and a couple of right hand sides. We can factorize the matrix with

P = Pardiso[L];

This results in a Pardiso object that is pretty much similar to LinearSolveFunction objects. But in constrast to almost everything in Mathematica, Pardiso is a pointer type with all its consequences. (In fact I have not implemented any deep copy routine for it. And it is almost impossible because the MKL Pardiso is pretty discrete about were it stores its data.) So use with caution!

You can solve a single equation with

x = P[b];

Or a couple of right hand sides contained in the matrix B by

X = P[B];

You can solve a linear system with the transpose of L with

y =P[b, "T"]

You can overwrite the stored matrix with another matrix with the same sparsity pattern as follows; this well keep the symbolic factorization:


You will get an error message if you use that with a matrix that has not the same number of nonzero entries as the original matrix. However, for performance reasons, the sparsity patterns are not compared. So use with caution!

user21 has recently built in a way to use alternative linear solvers in NDSolve. Here is an example how you can use Pardiso for that:

ClearAll[x, y, u];
f = NDSolveValue[{
   Laplacian[u[x, y], {x, y}] == 
    Sin[3 Pi x + 6 Pi y] + 2 Sin[7 Pi x + 5 Pi y] - 1,
   DirichletCondition[u[x, y] == 0, True]
  {x, y} \[Element] Disk[],
  Method -> {"PDEDiscretization" -> {"FiniteElement",
      "PDESolveOptions" -> {"LinearSolver" -> {Pardiso}}

It is not faster than LinearSolve[#, Method ->"Pardiso"] (it is the same backend!), but it should be able to provide quite a speed-up for for transient PDEs. In fact, the MKL Pardiso provides the feature to apply CG interations with the factorization as preconditioner (with automated refresh) for such a setup. I have not implemented it in PardisoLink, though. (And I am not going to do that in the near future, so don't ask ;) )

Per default, Pardiso performs something like an $LU$-decomposition (not really, but in the spirit). This corresponds to calling it with

P["Update"[L], "MatrixType"->11]

In principle, it should also support structurally symmetric matrices (not symmetric but with symmetric sparsity pattern, "MatrixType"->1), symmetric, but potentially indefinite matrices ("MatrixType"->-2), and symmetric, positive definite matrices ("MatrixType"->2). However, as I am writing this, I realize that the following throws an error. =O

n = 1000;
m = 20;
L = SparseArray[{Band[{1, 1}] -> -2. n, Band[{2, 1}] -> 1. n, Band[{1, 2}] -> 1. n}, {n, n}, 0.];
b = RandomReal[{-1, 1}, n];
P = Pardiso[L, "MatrixType" -> 2]

But "MatrixType"->-2 seems to work. I have to see what I can do to fix that...


Hah! Not a bug, but a feature: Pardiso told me that my assumption that the finite-difference Laplacian is positive-definite is actually wrong! It factorizes

P = Pardiso[-L, "MatrixType" -> 2]

without any problems.

Caution: "MatrixType" -> 2 and "MatrixType" -> -2 reads just the upper right triangle (including the diagonal) of the input matrix and neglects the lower left triangle (excluding the diagonal). Hence

Pardiso[-L, "MatrixType" -> 2] 


Pardiso[-UpperTriangulize[L], "MatrixType" -> 2] 

should have the same effect so that one can spare oneself the assembly of the lower left triangle.


You can find a short tutorial here:

NotebookOpen@FileNameJoin[{PardisoLink`Private`$packageDirectory, "Documentation", "Examples", "Example.nb"}]
  • 1
    $\begingroup$ @Szabolcs Your "LTemplate`" package at work. Would be great to have some feedback from you! $\endgroup$ Apr 19 '20 at 21:02
  • $\begingroup$ @Szabolcs Is there a way to return an empty vector {} (as MTensor) from a LibraryLink function? $\endgroup$ Apr 19 '20 at 21:24
  • 1
    $\begingroup$ Just noting that we're not getting notifications from the @mentions unless we've participated in the discussion thread already. I found this by accident. I will definitely take a look tomorrow (it's late now). $\endgroup$
    – Szabolcs
    Apr 19 '20 at 21:37
  • $\begingroup$ Yes, you can absolutely return a 0-size vector. Do you want to return that as a matrix, is that the problem? $\endgroup$
    – Szabolcs
    Apr 19 '20 at 21:37
  • 1
    $\begingroup$ This is really great @HenrikSchumacher! I think that there is a small mistake in the example file and in the implementation of "Update". After calling P["Update"[-L]] and computing y=P[-b] we should have y==x. The last line of the implementation of Update for SparseArrays reads P["Update"[Developer`ToPackedArray]] and so it's not passing anything to the implementation of Update for Vector inputs (I fixed this by calling P["Update"[A["NonzeroValues"]]]). I think that in the Vector implementation it should be a=Developer`ToPackedArray[N[a0]] instead of ToPackedArrayQ). $\endgroup$
    – beets
    Apr 21 '20 at 0:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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