Solving sparse linear systems in CSR format with complex values

[mhassell]

Hello, I'm working on using some of scipy's sparse solvers for both real and complex problems, and I seem to have found an issue when solving linear systems Ax=b where A is a real -valued matrix in csc format, and b is a complex numpy array. When I try to solve such a system, I get the following warning:

/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/core/numeric.py:460:
ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)

and the vector x that is returned only contains the real part of the solution. Now, if we explicitly cast A to be a complex data type by way of the dtype=np.cfloat flag in the constructor, then everything turns out fine. Now the rub is that I often use the same matrix A to solve linear systems with both real and complex right-hand-sides, and it would make sense for the output to be the same as the input, i.e. x should be real when A and b are real, and it should be complex when either A or b is complex. Casting A to a complex data type returns x as a complex array, which it may not be when b is real. Here's a minimal working example to show the warning being thrown:

import numpy as np
import scipy.sparse.linalg
import scipy.sparse

A = scipy.sparse.csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])

b = np.array([1,1,1]) +1j*np.array([1,1,1])

#throws warning and returns erroneous result
x_complex = scipy.sparse.linalg.spsolve(A,b)

[perimosocordiae]

Seems like a pretty easy fix: after making sure A is a floating-point type here, we can get the result dtype using np.promote_types and cast A and/or b to it as needed.