v0.8.1

.gitignore

@@ -20,6 +20,7 @@ parts/
 sdist/
 var/
 wheels/
+recipe/
 pip-wheel-metadata/
 share/python-wheels/
 *.egg-info/
@@ -109,6 +110,8 @@ venv/
 ENV/
 env.bak/
 venv.bak/
+.idea/
+.vscode/
 
 # Spyder project settings
 .spyderproject

CHANGELOG.md

@@ -1,3 +1,7 @@
+### Version 0.8.1
+
+* `cast=True` will now cast to compatible complex floats if one array is complex and one is real
+
 ### Version 0.8.0
 
 * Added support for complex data types

setup.py

@@ -2,7 +2,7 @@
 from setuptools import setup, find_packages
 
 DISTNAME = 'sparse_dot_mkl'
-VERSION = '0.8.0'
+VERSION = '0.8.1'
 DESCRIPTION = "Intel MKL wrapper for sparse matrix multiplication"
 MAINTAINER = 'Chris Jackson'
 MAINTAINER_EMAIL = '[email protected]'

sparse_dot_mkl/_mkl_interface/_common.py

@@ -1,6 +1,6 @@
 from ._cfunctions import MKL, mkl_library_name
 from ._constants import *
-from ._structs import *
+from ._structs import sparse_matrix_t, MKL_Complex8, MKL_Complex16
 
 import numpy as _np
 import ctypes as _ctypes
@@ -572,58 +572,73 @@ def _sanity_check(matrix_a, matrix_b, allow_vector=False):
         raise ValueError(err_msg)
 
 
-def _cast_to_float64(matrix):
+def _cast_to(matrix, dtype):
     """ Make a copy of the array as double precision floats or return the reference if it already is"""
-    return matrix.astype(_np.float64) if matrix.dtype != _np.float64 else matrix
-
-def _cast_to_cdouble(matrix):
-    """ Make a copy of the array as double precision complex floats or return the reference if it already is"""
-    return matrix.astype(_np.cdouble) if matrix.dtype != _np.cdouble else matrix
+    return matrix.astype(dtype) if matrix.dtype != dtype else matrix
+
+def _is_valid_dtype(matrix, complex_dtype=False, all_dtype=False):
+    """ Check to see if it's a usable float dtype """
+    if all_dtype:
+        return matrix.dtype in NUMPY_FLOAT_DTYPES + NUMPY_COMPLEX_DTYPES
+    elif complex_dtype:
+        return matrix.dtype in NUMPY_COMPLEX_DTYPES
+    else:
+        return matrix.dtype in NUMPY_FLOAT_DTYPES
 
 def _type_check(matrix_a, matrix_b=None, cast=False):
     """
     Make sure that both matrices are single precision floats or both are double precision floats
     If not, convert to double precision floats if cast is True, or raise an error if cast is False
     """
 
-    all_dtypes = NUMPY_FLOAT_DTYPES + NUMPY_COMPLEX_DTYPES
     _n_complex = _np.iscomplexobj(matrix_a) + _np.iscomplexobj(matrix_b)
 
     # If there's no matrix B and matrix A is valid dtype, return it
-    if matrix_b is None and matrix_a.dtype in all_dtypes:
+    if matrix_b is None and _is_valid_dtype(matrix_a, all_dtype=True):
         return matrix_a
     # If matrix A is complex but not csingle or cdouble, and cast is True, convert it to a cdouble
     elif matrix_b is None and cast and _n_complex == 1:
-        return _cast_to_cdouble(matrix_a)
+        return _cast_to(matrix_a, _np.cdouble)
     # If matrix A is real but not float32 or float64, and cast is True, convert it to a float64
     elif matrix_b is None and cast:
-        return _cast_to_float64(matrix_a)
+        return _cast_to(matrix_a, _np.float64)
     # Raise an error - the dtype is invalid and cast is False
     elif matrix_b is None:
-        err_msg = "Matrix data type must be float32, float64, csingle, or cdouble; {a} provided".format(a=matrix_a.dtype)
-        raise ValueError(err_msg)
+        _err_msg = f"Matrix data type must be float32, float64, csingle, or cdouble; {matrix_a.dtype} provided"
+        raise ValueError(_err_msg)
 
     # If Matrix A & B have the same valid dtype, return them
-    if matrix_a.dtype in all_dtypes and matrix_a.dtype == matrix_b.dtype:
+    if _is_valid_dtype(matrix_a, all_dtype=True) and matrix_a.dtype == matrix_b.dtype:
         return matrix_a, matrix_b
+
     # If neither matrix is complex and cast is True, convert to float64s and return them
     elif cast and _n_complex == 0:
-        debug_print("Recasting matrix data types {a} and {b} to _np.float64".format(a=matrix_a.dtype, b=matrix_b.dtype))
-        return _cast_to_float64(matrix_a), _cast_to_float64(matrix_b)
+        debug_print(f"Recasting matrix data types {matrix_a.dtype} and {matrix_b.dtype} to _np.float64")
+        return _cast_to(matrix_a, _np.float64), _cast_to(matrix_b, _np.float64)
+
     # If both matrices are complex and cast is True, convert to cdoubles and return them
     elif cast and _n_complex == 2:
-        debug_print("Recasting matrix data types {a} and {b} to _np.cdouble".format(a=matrix_a.dtype, b=matrix_b.dtype))
-        return _cast_to_cdouble(matrix_a), _cast_to_cdouble(matrix_b)
-    # Can't cast reals and complex matrices together
+        debug_print(f"Recasting matrix data types {matrix_a.dtype} and {matrix_b.dtype} to _np.cdouble")
+        return _cast_to(matrix_a, _np.cdouble), _cast_to(matrix_b, _np.cdouble)
+
+    # Cast reals and complex matrices together
+    elif cast and _n_complex == 1 and _is_valid_dtype(matrix_a, complex_dtype=True):
+        debug_print(f"Recasting matrix data type {matrix_b.dtype} to {matrix_a.dtype}")
+        return matrix_a, _cast_to(matrix_b, matrix_a.dtype)
+    elif cast and _n_complex == 1 and _is_valid_dtype(matrix_b, complex_dtype=True):
+        debug_print(f"Recasting matrix data type {matrix_a.dtype} to {matrix_b.dtype}")
+        return _cast_to(matrix_a, matrix_b.dtype), matrix_b
     elif cast and _n_complex == 1:
-        err_msg = "Cannot cast reals and complexes together; {a} and {b} provided".format(a=matrix_a.dtype,
-                                                                                          b=matrix_b.dtype)
-        raise ValueError(err_msg)
+        debug_print(f"Recasting matrix data type {matrix_a.dtype} and {matrix_b.dtype} to _np.cdouble")
+        return _cast_to(matrix_a, _np.cdouble), _cast_to(matrix_b, _np.cdouble)
+
     # If cast is False, can't cast anything together
     elif not cast:
-        err_msg = "Matrix data types must be in concordance; {a} and {b} provided".format(a=matrix_a.dtype,
-                                                                                          b=matrix_b.dtype)
-        raise ValueError(err_msg)
+        raise ValueError(
+            "Matrix data type must be float32, float64, csingle, or cdouble, " +
+            "and must be the same if cast=False; " +
+            f"{matrix_a.dtype} & {matrix_b.dtype} provided"
+        )
 
 def _mkl_scalar(scalar, complex_type, double_precision):
     """Turn complex scalars into appropriate precision MKL scalars or leave floats as floats"""
@@ -748,4 +763,4 @@ def _empty_output_check(matrix_a, matrix_b):
 
     # Neither trivial condition
     else:
-        return False
+        return False

sparse_dot_mkl/tests/test_mkl.py

@@ -14,7 +14,7 @@
 import scipy.sparse as _spsparse
 from sparse_dot_mkl import dot_product_mkl
 from sparse_dot_mkl._mkl_interface import (_create_mkl_sparse, _export_mkl, sparse_matrix_t, _destroy_mkl_handle,
-                                           _convert_to_csr, _order_mkl_handle, MKL)
+                                           _convert_to_csr, _order_mkl_handle, MKL, _type_check)
 
 SEED = 86
 
@@ -227,6 +227,123 @@ def test_create_convert_bsr(self):
         self.is_sparse_identical_A(self.mat1.astype(np.float32), cycle_3)
 
 
+class TestTypeConversions(unittest.TestCase):
+
+    dtype = np.float32
+    cast_dtype = np.float64
+
+    final_dtype=None
+    always_cast=False
+
+    def setUp(self):
+        self.mat1 = MATRIX_1.copy()
+        self.mat2 = MATRIX_2.copy()
+
+    def test_valid_pairs(self):
+
+        a, b = self.mat1.astype(self.dtype), self.mat2.astype(self.dtype)
+        c, d = _type_check(a, b, cast=self.always_cast)
+
+        if self.always_cast:
+            self.assertNotEqual(id(a), id(c))
+            self.assertNotEqual(id(b), id(d))
+        else:
+            self.assertEqual(id(a), id(c))
+            self.assertEqual(id(b), id(d))
+
+        fd = self.final_dtype if self.final_dtype is not None else self.dtype
+
+        self.assertEqual(c.dtype, fd)
+        self.assertEqual(d.dtype, fd)
+
+    def test_cast_pairs_right(self):
+
+        a, b = self.mat1.astype(self.dtype), self.mat2.astype(self.cast_dtype)
+
+        with self.assertRaises(ValueError):
+            c, d = _type_check(a, b)
+
+        c, d = _type_check(a, b, cast=True)
+
+        self.assertNotEqual(id(a), id(c))
+
+        if self.always_cast:
+            self.assertNotEqual(id(b), id(d))
+        else:
+            self.assertEqual(id(b), id(d))
+
+        fd = self.final_dtype if self.final_dtype is not None else self.cast_dtype
+
+        self.assertEqual(c.dtype, fd)
+        self.assertEqual(d.dtype, fd)
+
+    def test_cast_pairs_left(self):
+
+        a, b = self.mat1.astype(self.cast_dtype), self.mat2.astype(self.dtype)
+
+        with self.assertRaises(ValueError):
+            c, d = _type_check(a, b)
+
+        c, d = _type_check(a, b, cast=True)
+
+        if self.always_cast:
+            self.assertNotEqual(id(a), id(c))
+        else:
+            self.assertEqual(id(a), id(c))
+
+        self.assertNotEqual(id(b), id(d))
+
+        fd = self.final_dtype if self.final_dtype is not None else self.cast_dtype
+
+        self.assertEqual(c.dtype, fd)
+        self.assertEqual(d.dtype, fd)
+
+
+class TestTypeConversions2(TestTypeConversions):
+
+    dtype = np.csingle
+    cast_dtype = np.cdouble
+
+
+class TestTypeConversions3(TestTypeConversions):
+
+    dtype = np.float32
+    cast_dtype = np.cdouble
+
+
+class TestTypeConversions4(TestTypeConversions):
+
+    dtype = np.float64
+    cast_dtype = np.cdouble
+
+
+class TestTypeConversions5(TestTypeConversions):
+
+    dtype = np.int32
+    cast_dtype = np.float32
+
+    final_dtype = np.float64
+    always_cast = True
+
+
+class TestTypeConversions6(TestTypeConversions):
+
+    dtype = np.int32
+    cast_dtype = np.int64
+
+    final_dtype = np.float64
+    always_cast = True
+
+
+class TestTypeConversions7(TestTypeConversions):
+
+    dtype = np.clongdouble
+    cast_dtype = np.clongdouble
+
+    final_dtype = np.cdouble
+    always_cast = True
+
+
 def run():
     unittest.main(module='sparse_dot_mkl.tests.test_mkl')
     unittest.main(module='sparse_dot_mkl.tests.test_gram_matrix')

sparse_dot_mkl/tests/test_qr_solver.py

@@ -20,8 +20,9 @@ def setUp(self):
         self.mat3 = self.X.copy()
 
     def test_sparse_solver(self):
-        mat3 = sparse_qr_solve_mkl(self.mat1, self.mat2, debug=True)
-        npt.assert_array_almost_equal(self.mat3, mat3)
+        with self.assertWarns(DeprecationWarning):
+            mat3 = sparse_qr_solve_mkl(self.mat1, self.mat2, debug=True)
+            npt.assert_array_almost_equal(self.mat3, mat3)
 
     def test_sparse_solver_single(self):
         mat3 = sparse_qr_solve_mkl(self.mat1.astype(np.float32), self.mat2.astype(np.float32))

sparse_dot_mkl/tests/test_sparse_dense.py

@@ -1,8 +1,7 @@
 import unittest
-from unittest.case import skip
 import numpy as np
 import numpy.testing as npt
-from scipy.sparse.csr import csr_matrix
+from scipy.sparse import csr_matrix
 from sparse_dot_mkl import dot_product_mkl
 from sparse_dot_mkl.tests.test_mkl import MATRIX_1, MATRIX_2, make_matrixes
 
@@ -27,15 +26,16 @@ def setUp(self):
     def test_float32_b_sparse(self):
         d1, d2 = self.mat1_d.astype(self.single_dtype), self.mat2.astype(self.single_dtype)
 
-        mat3 = dot_product_mkl(d1, d2, debug=True)
-        mat3_np = np.dot(d1, d2.A)
+        with self.assertWarns(DeprecationWarning):
+            mat3 = dot_product_mkl(d1, d2, debug=True)
+            mat3_np = np.dot(d1, d2.A)
 
         npt.assert_array_almost_equal(mat3_np, mat3)
 
         mat3_np += 3.
         out = np.ones(mat3_np.shape, dtype=self.single_dtype, order=self.order)
 
-        mat3 = dot_product_mkl(d1, d2, out=out, out_scalar=3., debug=True)
+        mat3 = dot_product_mkl(d1, d2, out=out, out_scalar=3.)
         npt.assert_array_almost_equal(mat3_np, mat3)
         npt.assert_array_almost_equal(mat3_np, out)
 
@@ -47,7 +47,7 @@ def test_float32_b_sparse(self):
     def test_float64_b_sparse(self):
         d1, d2 = self.mat1_d.astype(self.double_dtype), self.mat2.astype(self.double_dtype)
 
-        mat3 = dot_product_mkl(d1, d2, debug=True)
+        mat3 = dot_product_mkl(d1, d2)
         mat3_np = np.dot(d1, d2.A)
 
         npt.assert_array_almost_equal(mat3_np, mat3)