Complex matrix support

[twhughes]

Hi, I'm looking to use this for complex128 csr matrices, but it is not currently supported.

Is this in the pipeline? Is it straightforward to implement?

Right now, I'm performing a workaround by constructing a twice as large, real matrix, which works but is significantly slower.

Thanks!

[ianwilliamson]

I made an attempt at getting this to work by changing the matrix type flag to 6 (complex symmetric) and commenting out the checks on the data types. It seems like Pardiso will successfully perform the solve (no error thrown) but there is still an issue with passing the data back into Python from C.

Generating random complex symmetric matrices and solving with random right hand sides gives solutions that look like:

array([[       nan+nanj],
       [       nan+nanj],
       [       nan+nanj],
       [1.42110623+nanj],
       [       nan+nanj],
       [0.         +0.j],
       [0.         +0.j],
       [0.         +0.j],
       [0.         +0.j],
       [0.         +0.j]])

Note that scipy's spsolve() is able to handle these same problems so it shouldn't be an issue with the matrix symmetry.

[ianwilliamson]

My bad... in the above comment I was actually not setting mtype correctly. Those outputs correspond to feeding complex inputs when Pardiso is expecting real inputs (still having mtype=11).

My approach to trying to get complex numbers to work was:

1. Comment out the checks in PyPardiso that see if A and b are of type float64

2. Update the pointers in _call_pardiso(self, A, b) to c_float64_p = np.ctypeslib.ndpointer(np.complex128)

3. Run the following

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
from pypardiso import *

ps = PyPardisoSolver(mtype=6)

np.random.seed(13)
n = 6
density = 0.1

A = sp.csr_matrix(sp.rand(n, n, density) + sp.eye(n)*1j)
A = A + A.transpose()
b = np.random.rand(n,1)+1j

x_pardiso = ps.solve(A,b)
print(x_pardiso)

x_scipy = spsolve(A,b)
print(x_scipy)

This gives a Segmentation fault (core dumped) from somewhere inside self._mkl_pardiso(). Note also that using c_float64_p = np.ctypeslib.ndpointer(np.float64) also works for the real case (with mtype=11), so it doesn't seem that the issue lies with using np.ctypeslib.ndpointer.

[ianwilliamson]

See my pull request #9

[victorminden]

@ianwilliamson @haasad I don't know if anyone is still interested in complex types in pypardiso now (2 years later) but I stumbled upon this as I was trying to wrap complex types for PARDISO 6 and got curious so I was playing around.

I think that this pull request #9 might actually work (more or less -- there is a test failing due to the removed conversion from half-precision that used to happen behind the scenes)!

The problem with the example code in #8 (comment) (I think) is that the matrix type is specified as a symmetric type mtype=6 which means that we should be passing a triangular matrix to PARDISO rather than passing the redundant other-half of the matrix as well. I verified on my local machine that passing either the triangular matrix or changing the mtype to nonsymmetric works.

(Possible that there were other problems I don't know about here)

[ianwilliamson]

@victorminden Feel free to take over that pull request, if you wish.

Just FYI... there is also this package: https://github.com/dwfmarchant/pyMKL that wraps (what I think is the same solver) and supports complex numbers. We have used it successfully over the past few years.

[victorminden]

@ianwilliamson Thanks for the note, I have been looking at PyMKL as well. My interest is as part of a "which of these two projects would be more straightforward for trying out PARDISO 6 with complex support via Python", so when I saw this and #10 I was trying to investigate low-hanging fruit here. Will spend more time with PyMKL as well.

[ianwilliamson]

Ah, makes sense. It would be really useful to have a working interface to pardiso version 6 from Python.

[andrejstmh]

Complex cdouble support:

https://github.com/andrejstmh/PyPardisoProject
https://github.com/dwfmarchant/pyMKL