search

scipy / scipy

SciPy library main repository
https://scipy.org
BSD 3-Clause "New" or "Revised" License
13.07k stars 5.18k forks source link

Should iterative linear solvers force the right-hand-side array to be an ndarray? #15657

Open michaelbynum opened 2 years ago

michaelbynum commented 2 years ago

At lease some of the iterative solvers in SciPy force the right-hand-side (RHS) array to an ndarray, even if the original RHS array is a subclass of ndarray. Should this be the case?

I am developing some custom matrix and array classes that enable parallel computing on distributed memory machines. The matrix class utilizes SciPy sparse matrices under the hood, and the array class inherits from ndarray. I thought I would be able to use SciPy's iterative linear solvers by wrapping the custom matrix class in a LinearOperator (my custom matrix can perform a parallel matvec with my custom array), but I don't have a way to convert my custom array class to an ndarray without a lot of unnecessary communication with MPI, which ruins parallel scalability.

Take the following code snippet for example: https://github.com/scipy/scipy/blob/2f1fc78b3d76fb974be57e1948ac828f7734d6f8/scipy/sparse/linalg/_isolve/utils.py#L73-L91

Line 73 uses asanarray which does allow subclasses, but then line 91 forces b to an ndarray. The comment on line 91 implies that the dtype just needs to be changed, so it might be okay to use astype here.

I would be happy to create a pull request to modify this behavior, but I wanted to know if there were reasons to force the RHS array to an ndarray that I don't know about.

zhaog6 commented 2 years ago

The comment on line 91 implies that the dtype just needs to be changed, so it might be okay to use astype here.

@michaelbynum No, the asarray conversion is to avoid the case where b is a deprecated subclass type numpy.matrix, it'll fail to be directly used as parameter of Krylov methods such as cg/cgs/bicg/bicgstab/... (Need to be converted).

michaelbynum commented 2 years ago

Thanks, @zhaog6. I appreciate the clarification. Do you see any way to enable the use of subclasses of ndarray?

zhaog6 commented 2 years ago

I have one year of experience for Python math library, and I have less experences (previously used C++). I haven't seen it for the time being in this case. But if you want to add new subclasses (a new array), I think you can try.

j-bowhay commented 1 month ago

@michaelbynum would overriding asdtype for your subclass fix this issue?

EDIT: oops sorry I misread