krypy icon indicating copy to clipboard operation
krypy copied to clipboard

[Bug] find_common_dtype results casting complex to real

Open homocomputeris opened this issue 7 years ago • 0 comments

Hi, It seems that find_common_dtype returns wrong dtype when passed an operator (maybe a matrix too).

Here is the minimal "not working" example:

import numpy as np
import scipy.sparse.linalg as sla
from krypy.utils import Arnoldi

dim = 8
shp = (dim, dim)
# a complex Hermitian matrix
R = np.random.rand(dim, dim) + 1j*np.random.rand(dim, dim)
H = R + R.conj().T
# a real vector
psi = np.random.rand(dim, 1)


def mv(psi):
    return H@psi


# define an operator
HLO = sla.LinearOperator(shp, matvec=mv)

# run Arnoldi iteration
arn = Arnoldi(HLO, psi)
for k in range(dim//2):
    arn.advance()
    v, h = arn.get()

print("Operator's type is", HLO.dtype)
print("Vector's types is", psi.dtype)
print("but the common type in Arnoldi is", arn.dtype)
print('\n')

Output:

Operator's type is complex128
Vector's types is float64
but the common type in Arnoldi is float64

/usr/lib/python3.6/site-packages/krypy/utils.py:984: ComplexWarning: Casting complex values to real discards the imaginary part
  self.H[j, k] += alpha
/usr/lib/python3.6/site-packages/krypy/utils.py:1002: ComplexWarning: Casting complex values to real discards the imaginary part
  self.V[:, [k+1]] = Av / self.H[k+1, k]

Maybe NumPy's default function would be easier?

if M:
    self.dtype = numpy.common_type(A, v, M)
else:
    self.dtype = numpy.common_type(A, v)

homocomputeris avatar Feb 17 '18 14:02 homocomputeris