Reusing PARDISO symbolic factorization

Question

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.

Answer 1

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

PardisoLink

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.

Installation

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. ;)

Example

Okay, here is a simple usage example:

Needs["PardisoLink`"];
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:

P["Update"[-L]]

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]
   },
  u,
  {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...

Edit

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] 

and

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

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

Tutorial

You can find a short tutorial here:

NotebookOpen@FileNameJoin[{PardisoLink`Private`$packageDirectory, "Documentation", "Examples", "Example.nb"}]