Added complex matrix support in cholmoddense2numpy and numpy2cholmoddense so that solve() works.

Author: yig

CHANGELOG.md

@@ -6,6 +6,10 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 
 ## [Unreleased]
 
+## [v1.5.1] - 2026-02-03
+### Added
+- Added complex matrix support in `cholmoddense2numpy` and `numpy2cholmoddense` so that `solve()` works.
+
 ## [v1.5] - 2026-02-03
 ### Added
 - Added complex matrix support in `cholmodsparse2scipy` and `scipy2cholmodsparse`.

sparseqr/__init__.py

@@ -17,7 +17,7 @@
 
 from __future__ import absolute_import
 
-__version__ = '1.5'
+__version__ = '1.5.1'
 
 # import the important things into the package's top-level namespace.
 from .sparseqr import qr, rz, solve, permutation_vector_to_matrix, qr_factorize,qmult

sparseqr/sparseqr.py

@@ -170,19 +170,41 @@ def numpy2cholmoddense( numpy_A ):
     nrow = numpy_A.shape[0]
     ncol = numpy_A.shape[1]
     lda  = nrow  # cholmod_dense is column-oriented
-    chol_A = lib.cholmod_l_allocate_dense( nrow, ncol, lda, lib.CHOLMOD_REAL, cc )
+
+    # Check if the array has complex entries and adapt cholmod type accordingly
+    is_complex = numpy.issubdtype(numpy_A.dtype, numpy.complexfloating)
+    cholmod_dtype = lib.CHOLMOD_COMPLEX if is_complex else lib.CHOLMOD_REAL
+
+    chol_A = lib.cholmod_l_allocate_dense( nrow, ncol, lda, cholmod_dtype, cc )
     if chol_A == ffi.NULL:
         raise RuntimeError("Failed to allocate chol_A")
     Adata = ffi.cast( "double*", chol_A.x )
-    for j in range(ncol):  # FIXME inefficient?
-        Adata[(j*lda):((j+1)*lda)] = numpy_A[:,j]
+
+    if is_complex:
+        # chol_A.x has size 2*nrow*ncol and real and imag parts are interleaved [real, imag, real, imag, ...]
+        array_element_size = ffi.sizeof(ffi.typeof(Adata).item)
+        Adata_view = numpy.frombuffer(ffi.buffer(Adata, 2*nrow*ncol * array_element_size), ctype2dtype["double"])
+        for j in range(ncol):
+            col_data = numpy_A[:,j]
+            Adata_view[(j*lda*2):((j+1)*lda*2):2] = col_data.real
+            Adata_view[(j*lda*2)+1:((j+1)*lda*2):2] = col_data.imag
+    else:
+        for j in range(ncol):  # FIXME inefficient?
+            Adata[(j*lda):((j+1)*lda)] = numpy_A[:,j]
     return chol_A
 
 def cholmoddense2numpy( chol_A ):
     '''Convert a CHOLMOD dense matrix to a NumPy array.'''
     Adata = ffi.cast( "double*", chol_A.x )
 
-    result = asarray( ffi, Adata, chol_A.nrow*chol_A.ncol ).copy()
+    if chol_A.xtype == lib.CHOLMOD_COMPLEX:
+        # read real and imag part from the buffer and create view as complex datatype
+        result = asarray(ffi, Adata, 2*chol_A.nrow*chol_A.ncol).copy()
+        complex_dtype = numpy.dtype(f"c{result.itemsize * 2}")
+        result = result.view(complex_dtype)
+    else:
+        result = asarray( ffi, Adata, chol_A.nrow*chol_A.ncol ).copy()
+
     result = result.reshape( (chol_A.nrow, chol_A.ncol), order='F' )
     return result
 

test/test_complex.py

@@ -2,7 +2,7 @@
 
 import numpy as np
 import pytest
-from scipy.sparse import random as sparse_random
+import scipy.sparse
 
 import sparseqr
 
@@ -19,23 +19,23 @@ class TestComplexQR:
 
     def test_complex_qr_unitary_q(self):
         """Test that Q is unitary for complex matrix (Q @ Q^H = I)."""
-        A = sparse_random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
+        A = scipy.sparse.random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
         Q, R, P, rank = sparseqr.qr(A)
 
         QQH = (Q @ adj(Q)).toarray()
         np.testing.assert_allclose(QQH, np.eye(Q.shape[0]), atol=self.ATOL)
 
     def test_complex_qr_upper_triangular_r(self):
         """Test that R is upper triangular for complex matrix."""
-        A = sparse_random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
+        A = scipy.sparse.random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
         Q, R, P, rank = sparseqr.qr(A)
 
         lower_triangle = np.tril(R.toarray(), k=-1)
         np.testing.assert_allclose(lower_triangle, 0, atol=self.ATOL)
 
     def test_complex_qr_reconstruction(self):
         """Test that A can be reconstructed from QR decomposition: A = (QR) @ P^T."""
-        A = sparse_random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
+        A = scipy.sparse.random(m=120, n=100, density=0.1, dtype=np.complex128, random_state=42)
         Q, R, P, rank = sparseqr.qr(A)
 
         P_matrix = sparseqr.permutation_vector_to_matrix(P)
@@ -46,13 +46,60 @@ def test_complex_qr_reconstruction(self):
     def test_complex_qr_various_shapes(self, shape):
         """Test complex QR on square, tall, and wide matrices."""
         m, n = shape
-        A = sparse_random(m=m, n=n, density=0.1, dtype=np.complex128, random_state=42)
+        A = scipy.sparse.random(m=m, n=n, density=0.1, dtype=np.complex128, random_state=42)
         Q, R, P, rank = sparseqr.qr(A)
 
         P_matrix = sparseqr.permutation_vector_to_matrix(P)
         A_reconstructed = (Q @ R) @ P_matrix.T
         np.testing.assert_allclose(A_reconstructed.toarray(), A.toarray(), atol=self.ATOL)
 
 
+class TestComplexSolve:
+    """Tests for solve function with complex matrices."""
+
+    ATOL = 1e-10
+
+    def test_complex_solve_exact_system(self):
+        """Test solving an exact complex system where Ax=b has exact solution."""
+        A = scipy.sparse.random(m=50, n=50, density=0.2, dtype=np.complex128, random_state=42)
+        # Add diagonal to make it well-conditioned
+        A = A + 5 * scipy.sparse.eye(50, dtype=np.complex128)
+        x_true = np.random.RandomState(42).random(50) + 1j * np.random.RandomState(43).random(50)
+        b = A @ x_true
+
+        x = sparseqr.solve(A, b, tolerance=0)
+
+        assert x is not None, "solve() returned None"
+        np.testing.assert_allclose(x, x_true, atol=self.ATOL)
+
+    def test_complex_solve_overdetermined_least_squares(self):
+        """Test least-squares solution for overdetermined complex system."""
+        # 100 equations, 50 unknowns
+        A = scipy.sparse.random(m=100, n=50, density=0.2, dtype=np.complex128, random_state=42)
+        x_true = np.random.RandomState(42).random(50) + 1j * np.random.RandomState(43).random(50)
+        b = A @ x_true
+
+        x = sparseqr.solve(A, b, tolerance=0)
+
+        assert x is not None, "solve() returned None"
+        assert x.shape == (50,), f"Solution shape wrong: {x.shape}"
+        # For an exact system (b in range of A), solution should match
+        np.testing.assert_allclose(x, x_true, atol=self.ATOL)
+
+    def test_complex_solve_multiple_rhs(self):
+        """Test solving AX=B with multiple complex RHS vectors."""
+        A = scipy.sparse.random(m=50, n=50, density=0.2, dtype=np.complex128, random_state=42)
+        A = A + 5 * scipy.sparse.eye(50, dtype=np.complex128)
+        X_true = (np.random.RandomState(42).random((50, 3)) +
+                  1j * np.random.RandomState(43).random((50, 3)))
+        B = A @ X_true
+
+        X = sparseqr.solve(A, B, tolerance=0)
+
+        assert X is not None, "solve() returned None"
+        assert X.shape == (50, 3), f"Solution shape wrong: {X.shape}"
+        np.testing.assert_allclose(X, X_true, atol=self.ATOL)
+
+
 if __name__ == '__main__':
     pytest.main([__file__, '-v'])