Repository: scipy/scipy
Pull request: ENH: sparse arrays for hstack, vstack, bmat, block_diag. New block.
This PR makes `hstack`, `vstack` and `block_diags` work with sparse arrays. The return value is a sparse matrix unless any of the block inputs is a sparse array, in which case the return value is a sparse array.

It also makes `bmat` work with sparse arrays. But the name `bmat` is not very descriptive, **and** some inputs to `bmat` do not indicate whether to return sparse arrays or matrices. So, a new function `block_array` is introduced to build sparse arrays. It always returns a sparse array.  It does accept sparse matrices as input, but doesn't return sparse matrices.

In the interest of allowing partial updates from sparse matrix to sparse array, this PR also updates `bmat` to return a sparse array if any input is a sparse array. Otherwise it returns a sparse matrix. Thus it is backward compatible, yet still works with code that is using both sparse arrays and sparse matrices.  

#### Other related Issues and PRs:
Semantics of sparse array creation #18592
Add _array version of diags creation functions #18538

Fixes #18839

Problem statement:
ENH: concatenating sparse arrays should return sparse arrays
### Is your feature request related to a problem? Please describe.

Concatenating sparse arrays (via `hstack`, `vstack`, `bmat`) returns a sparse matrix.

### Describe the solution you'd like.

I would like them to maintain the input classes.

### Describe alternatives you've considered.

`hstack` and `vstack` could be replaced with something more `concat`, which is in the `array-api`, unlike `hstack`, `vstack`.

cc: @dschult

### Additional context (e.g. screenshots, GIFs)

_No response_

Hints from issue discussion:
The alternative you describe sounds like a good idea in any case. Moving towards supporting the Array API is something we might want to do in general.
The specification would be around what the outcome would be if I supply CSR, CSC and DIA arrays. Also note that, `np.hstack([a, b])` is order of magnitude more readable than `np.concat([a, b], axis=1)` so not sure that we should get rid of them that quick. 
> Also note that, np.hstack([a, b]) is order of magnitude more readable than np.concat([a, b], axis=1) so not sure that we should get rid of them that quick.

I was not suggesting a removal, just moving towards supporting the Array API. That does not necessarily imply that we would remove methods not part of the standard.

Commit message: clean up import and docs

Patch to review:
diff --git a/scipy/sparse/__init__.py b/scipy/sparse/__init__.py
index 875ce6dd8298..d27ccc4fd773 100644
--- a/scipy/sparse/__init__.py
+++ b/scipy/sparse/__init__.py
@@ -95,7 +95,8 @@
    block_diag - Build a block diagonal sparse matrix
    tril - Lower triangular portion of a matrix in sparse format
    triu - Upper triangular portion of a matrix in sparse format
-   bmat - Build a sparse matrix from sparse sub-blocks
+   block_array - Build a sparse array from sub-blocks
+   bmat - Build a sparse matrix from sub-blocks
    hstack - Stack sparse matrices horizontally (column wise)
    vstack - Stack sparse matrices vertically (row wise)
    rand - Random values in a given shape
diff --git a/scipy/sparse/_construct.py b/scipy/sparse/_construct.py
index c685c11c7741..5731c0ba2e10 100644
--- a/scipy/sparse/_construct.py
+++ b/scipy/sparse/_construct.py
@@ -5,7 +5,7 @@
 
 __all__ = ['spdiags', 'eye', 'identity', 'kron', 'kronsum',
            'hstack', 'vstack', 'bmat', 'rand', 'random', 'diags', 'block_diag',
-           'diags_array']
+           'diags_array', 'block_array']
 
 import numbers
 from functools import partial
@@ -16,12 +16,12 @@
 
 from ._sparsetools import csr_hstack
 from ._bsr import bsr_matrix
-from ._coo import coo_matrix
-from ._csc import csc_matrix
-from ._csr import csr_matrix
+from ._coo import coo_matrix, coo_array
+from ._csc import csc_matrix, csc_array
+from ._csr import csr_matrix, csr_array
 from ._dia import dia_matrix, dia_array
 
-from ._base import issparse
+from ._base import issparse, sparray
 
 
 def spdiags(data, diags, m=None, n=None, format=None):
@@ -481,9 +481,9 @@ def kronsum(A, B, format=None):
     return (L+R).asformat(format)  # since L + R is not always same format
 
 
-def _compressed_sparse_stack(blocks, axis):
+def _compressed_sparse_stack(blocks, axis, container):
     """
-    Stacking fast path for CSR/CSC matrices
+    Stacking fast path for CSR/CSC matrices or arrays
     (i) vstack for CSR, (ii) hstack for CSC.
     """
     other_axis = 1 if axis == 0 else 0
@@ -507,11 +507,20 @@ def _compressed_sparse_stack(blocks, axis):
         sum_dim += b.shape[axis]
         last_indptr += b.indptr[-1]
     indptr[-1] = last_indptr
+    # TODO remove this if-structure when sparse matrices removed
+    if container is coo_matrix:
+        if axis == 0:
+            return csr_matrix((data, indices, indptr),
+                              shape=(sum_dim, constant_dim))
+        else:
+            return csc_matrix((data, indices, indptr),
+                              shape=(constant_dim, sum_dim))
+
     if axis == 0:
-        return csr_matrix((data, indices, indptr),
+        return csr_array((data, indices, indptr),
                           shape=(sum_dim, constant_dim))
     else:
-        return csc_matrix((data, indices, indptr),
+        return csc_array((data, indices, indptr),
                           shape=(constant_dim, sum_dim))
 
 
@@ -565,10 +574,10 @@ def _stack_along_minor_axis(blocks, axis):
         data = np.empty(0, dtype=data_cat.dtype)
 
     if axis == 0:
-        return csc_matrix((data, indices, indptr),
+        return blocks[0]._csc_container((data, indices, indptr),
                           shape=(sum_dim, constant_dim))
     else:
-        return csr_matrix((data, indices, indptr),
+        return blocks[0]._csr_container((data, indices, indptr),
                           shape=(constant_dim, sum_dim))
 
 
@@ -602,23 +611,27 @@ def hstack(blocks, format=None, dtype=None):
            [3, 4, 6]])
 
     """
-    return bmat([blocks], format=format, dtype=dtype)
+    blocks = np.asarray(blocks, dtype='object')
+    if any(isinstance(b, sparray) for b in blocks.flat):
+        return _block_array([blocks], format, dtype, coo_array)
+    else:
+        return _block_array([blocks], format, dtype, coo_matrix)
 
 
 def vstack(blocks, format=None, dtype=None):
     """
-    Stack sparse matrices vertically (row wise)
+    Stack sparse arrays vertically (row wise)
 
     Parameters
     ----------
     blocks
-        sequence of sparse matrices with compatible shapes
+        sequence of sparse arrays with compatible shapes
     format : str, optional
         sparse format of the result (e.g., "csr")
-        by default an appropriate sparse matrix format is returned.
+        by default an appropriate sparse array format is returned.
         This choice is subject to change.
     dtype : dtype, optional
-        The data-type of the output matrix. If not given, the dtype is
+        The data-type of the output array. If not given, the dtype is
         determined from that of `blocks`.
 
     See Also
@@ -627,21 +640,29 @@ def vstack(blocks, format=None, dtype=None):
 
     Examples
     --------
-    >>> from scipy.sparse import coo_matrix, vstack
-    >>> A = coo_matrix([[1, 2], [3, 4]])
-    >>> B = coo_matrix([[5, 6]])
+    >>> from scipy.sparse import coo_array, vstack
+    >>> A = coo_array([[1, 2], [3, 4]])
+    >>> B = coo_array([[5, 6]])
     >>> vstack([A, B]).toarray()
     array([[1, 2],
            [3, 4],
            [5, 6]])
 
     """
-    return bmat([[b] for b in blocks], format=format, dtype=dtype)
+    blocks = np.asarray(blocks, dtype='object')
+    if any(isinstance(b, sparray) for b in blocks.flat):
+        return _block_array([[b] for b in blocks], format, dtype, coo_array)
+    else:
+        return _block_array([[b] for b in blocks], format, dtype, coo_matrix)
 
 
 def bmat(blocks, format=None, dtype=None):
     """
-    Build a sparse matrix from sparse sub-blocks
+    Build a sparse array or matrix from sparse sub-blocks
+
+    Note: `block_array` is preferred over `bmat`. They are the same function
+    except that `bmat` can return a deprecated sparse matrix.
+    `bmat` returns a coo_matrix if none of the inputs are a sparse array.
 
     Parameters
     ----------
@@ -658,18 +679,23 @@ def bmat(blocks, format=None, dtype=None):
 
     Returns
     -------
-    bmat : sparse matrix
+    bmat : sparse matrix or array
+        If any block in blocks is a sparse array, return a sparse array.
+        Otherwise return a sparse matrix.
+
+        If you want a sparse array built from blocks that are not sparse
+        arrays, use `block_array()`.
 
     See Also
     --------
-    block_diag, diags
+    block_array
 
     Examples
     --------
-    >>> from scipy.sparse import coo_matrix, bmat
-    >>> A = coo_matrix([[1, 2], [3, 4]])
-    >>> B = coo_matrix([[5], [6]])
-    >>> C = coo_matrix([[7]])
+    >>> from scipy.sparse import coo_array, bmat
+    >>> A = coo_array([[1, 2], [3, 4]])
+    >>> B = coo_array([[5], [6]])
+    >>> C = coo_array([[7]])
     >>> bmat([[A, B], [None, C]]).toarray()
     array([[1, 2, 5],
            [3, 4, 6],
@@ -681,7 +707,60 @@ def bmat(blocks, format=None, dtype=None):
            [0, 0, 7]])
 
     """
+    blocks = np.asarray(blocks, dtype='object')
+    if any(isinstance(b, sparray) for b in blocks.flat):
+        return _block_array(blocks, format, dtype, coo_array)
+    else:
+        return _block_array(blocks, format, dtype, coo_matrix)
 
+
+def block_array(blocks, *, format=None, dtype=None):
+    """
+    Build a sparse array from sparse sub-blocks
+
+    Parameters
+    ----------
+    blocks : array_like
+        Grid of sparse arrays with compatible shapes.
+        An entry of None implies an all-zero array.
+    format : {'bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil'}, optional
+        The sparse format of the result (e.g. "csr"). By default an
+        appropriate sparse array format is returned.
+        This choice is subject to change.
+    dtype : dtype, optional
+        The data-type of the output array. If not given, the dtype is
+        determined from that of `blocks`.
+
+    Returns
+    -------
+    block_array : sparse array
+
+    See Also
+    --------
+    block_diag : specify blocks along the main diagonals
+    diags : specify (possibly offset) diagonals
+
+    Examples
+    --------
+    >>> from scipy.sparse import coo_array, block_array
+    >>> A = coo_array([[1, 2], [3, 4]])
+    >>> B = coo_array([[5], [6]])
+    >>> C = coo_array([[7]])
+    >>> block_array([[A, B], [None, C]]).toarray()
+    array([[1, 2, 5],
+           [3, 4, 6],
+           [0, 0, 7]])
+
+    >>> block_array([[A, None], [None, C]]).toarray()
+    array([[1, 2, 0],
+           [3, 4, 0],
+           [0, 0, 7]])
+
+    """
+    return _block_array(blocks, format, dtype, coo_array)
+
+
+def _block_array(blocks, format, dtype, coo_container):
     blocks = np.asarray(blocks, dtype='object')
 
     if blocks.ndim != 2:
@@ -690,29 +769,29 @@ def bmat(blocks, format=None, dtype=None):
     M,N = blocks.shape
 
     # check for fast path cases
-    if (format in (None, 'csr') and all(isinstance(b, csr_matrix)
-                                        for b in blocks.flat)):
+    if (format in (None, 'csr') and
+        all(issparse(b) and b.format == "csr" for b in blocks.flat)
+    ):
         if N > 1:
-            # stack along columns (axis 1):
-            blocks = [[_stack_along_minor_axis(blocks[b, :], 1)]
-                      for b in range(M)]   # must have shape: (M, 1)
+            # stack along columns (axis 1): must have shape (M, 1)
+            blocks = [[_stack_along_minor_axis(blocks[b, :], 1)] for b in range(M)]
             blocks = np.asarray(blocks, dtype='object')
 
         # stack along rows (axis 0):
-        A = _compressed_sparse_stack(blocks[:, 0], 0)
+        A = _compressed_sparse_stack(blocks[:, 0], 0, coo_container)
         if dtype is not None:
             A = A.astype(dtype)
         return A
-    elif (format in (None, 'csc') and all(isinstance(b, csc_matrix)
-                                          for b in blocks.flat)):
+    elif (format in (None, 'csc') and
+          all(issparse(b) and b.format == "csc" for b in blocks.flat)
+    ):
         if M > 1:
-            # stack along rows (axis 0):
-            blocks = [[_stack_along_minor_axis(blocks[:, b], 0)
-                       for b in range(N)]]   # must have shape: (1, N)
+            # stack along rows (axis 0): must have shape (1, N)
+            blocks = [[_stack_along_minor_axis(blocks[:, b], 0) for b in range(N)]]
             blocks = np.asarray(blocks, dtype='object')
 
         # stack along columns (axis 1):
-        A = _compressed_sparse_stack(blocks[0, :], 1)
+        A = _compressed_sparse_stack(blocks[0, :], 1, coo_container)
         if dtype is not None:
             A = A.astype(dtype)
         return A
@@ -725,7 +804,7 @@ def bmat(blocks, format=None, dtype=None):
     for i in range(M):
         for j in range(N):
             if blocks[i,j] is not None:
-                A = coo_matrix(blocks[i,j])
+                A = coo_array(blocks[i,j])
                 blocks[i,j] = A
                 block_mask[i,j] = True
 
@@ -771,27 +850,28 @@ def bmat(blocks, format=None, dtype=None):
         np.add(B.col, col_offsets[j], out=col[idx], dtype=idx_dtype)
         nnz += B.nnz
 
-    return coo_matrix((data, (row, col)), shape=shape).asformat(format)
+    return coo_container((data, (row, col)), shape=shape).asformat(format)
 
 
 def block_diag(mats, format=None, dtype=None):
     """
-    Build a block diagonal sparse matrix from provided matrices.
+    Build a block diagonal sparse matrix or array from provided matrices.
 
     Parameters
     ----------
-    mats : sequence of matrices
-        Input matrices.
+    mats : sequence of matrices or arrays
+        Input matrices or arrays.
     format : str, optional
-        The sparse format of the result (e.g., "csr"). If not given, the matrix
+        The sparse format of the result (e.g., "csr"). If not given, the result
         is returned in "coo" format.
     dtype : dtype specifier, optional
-        The data-type of the output matrix. If not given, the dtype is
+        The data-type of the output. If not given, the dtype is
         determined from that of `blocks`.
 
     Returns
     -------
-    res : sparse matrix
+    res : sparse matrix or array
+        If any inputs are arrays, the output is an array
 
     Notes
     -----
@@ -800,14 +880,14 @@ def block_diag(mats, format=None, dtype=None):
 
     See Also
     --------
-    bmat, diags
+    block_array, diags
 
     Examples
     --------
-    >>> from scipy.sparse import coo_matrix, block_diag
-    >>> A = coo_matrix([[1, 2], [3, 4]])
-    >>> B = coo_matrix([[5], [6]])
-    >>> C = coo_matrix([[7]])
+    >>> from scipy.sparse import coo_array, block_diag
+    >>> A = coo_array([[1, 2], [3, 4]])
+    >>> B = coo_array([[5], [6]])
+    >>> C = coo_array([[7]])
     >>> block_diag((A, B, C)).toarray()
     array([[1, 2, 0, 0],
            [3, 4, 0, 0],
@@ -816,6 +896,11 @@ def block_diag(mats, format=None, dtype=None):
            [0, 0, 0, 7]])
 
     """
+    if any(isinstance(a, sparray) for a in mats):
+        container = coo_array
+    else:
+        container = coo_matrix
+
     row = []
     col = []
     data = []
@@ -823,7 +908,7 @@ def block_diag(mats, format=None, dtype=None):
     c_idx = 0
     for a in mats:
         if isinstance(a, (list, numbers.Number)):
-            a = coo_matrix(a)
+            a = coo_array(a)
         nrows, ncols = a.shape
         if issparse(a):
             a = a.tocoo()
@@ -840,7 +925,7 @@ def block_diag(mats, format=None, dtype=None):
     row = np.concatenate(row)
     col = np.concatenate(col)
     data = np.concatenate(data)
-    return coo_matrix((data, (row, col)),
+    return container((data, (row, col)),
                       shape=(r_idx, c_idx),
                       dtype=dtype).asformat(format)
 
diff --git a/scipy/sparse/tests/test_construct.py b/scipy/sparse/tests/test_construct.py
index 5449c28676a3..a33ba97d41c0 100644
--- a/scipy/sparse/tests/test_construct.py
+++ b/scipy/sparse/tests/test_construct.py
@@ -10,9 +10,9 @@
 from scipy._lib._util import check_random_state
 
 from scipy.sparse import (csr_matrix, coo_matrix,
+                          csr_array, coo_array,
                           _construct as construct)
 from scipy.sparse._construct import rand as sprand
-from scipy.sparse._sputils import matrix
 
 sparse_formats = ['csr','csc','coo','bsr','dia','lil','dok']
 
@@ -340,8 +340,8 @@ def test_kronsum(self):
 
     def test_vstack(self):
 
-        A = coo_matrix([[1,2],[3,4]])
-        B = coo_matrix([[5,6]])
+        A = coo_array([[1,2],[3,4]])
+        B = coo_array([[5,6]])
 
         expected = array([[1, 2],
                           [3, 4],
@@ -364,10 +364,15 @@ def test_vstack(self):
         assert_equal(result.indices.dtype, np.int32)
         assert_equal(result.indptr.dtype, np.int32)
 
+        assert construct.vstack([A, B])._is_array
+        assert not construct.vstack([coo_matrix(A), coo_matrix(B)])._is_array
+        assert construct.vstack([A, coo_matrix(B)])._is_array
+        assert construct.vstack([coo_matrix(A), B])._is_array
+
     def test_hstack(self):
 
-        A = coo_matrix([[1,2],[3,4]])
-        B = coo_matrix([[5],[6]])
+        A = coo_array([[1,2],[3,4]])
+        B = coo_array([[5],[6]])
 
         expected = array([[1, 2, 5],
                           [3, 4, 6]])
@@ -385,85 +390,136 @@ def test_hstack(self):
                                       dtype=np.float32).dtype,
                      np.float32)
 
-    def test_bmat(self):
+        assert construct.hstack([A, B])._is_array
+        assert not construct.hstack([coo_matrix(A), coo_matrix(B)])._is_array
+        assert construct.hstack([A, coo_matrix(B)])._is_array
+        assert construct.hstack([coo_matrix(A), B])._is_array
+
+    @pytest.mark.parametrize("block_array", (construct.bmat, construct.block_array))
+    def test_block_creation(self, block_array):
 
-        A = coo_matrix([[1, 2], [3, 4]])
-        B = coo_matrix([[5],[6]])
-        C = coo_matrix([[7]])
-        D = coo_matrix((0, 0))
+        A = coo_array([[1, 2], [3, 4]])
+        B = coo_array([[5],[6]])
+        C = coo_array([[7]])
+        D = coo_array((0, 0))
 
         expected = array([[1, 2, 5],
                           [3, 4, 6],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, B], [None, C]]).toarray(), expected)
-        E = csr_matrix((1, 2), dtype=np.int32)
-        assert_equal(construct.bmat([[A.tocsr(), B.tocsr()],
+        assert_equal(block_array([[A, B], [None, C]]).toarray(), expected)
+        E = csr_array((1, 2), dtype=np.int32)
+        assert_equal(block_array([[A.tocsr(), B.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), B.tocsc()],
+        assert_equal(block_array([[A.tocsc(), B.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
         expected = array([[1, 2, 0],
                           [3, 4, 0],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, None], [None, C]]).toarray(),
+        assert_equal(block_array([[A, None], [None, C]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsr(), E.T.tocsr()],
+        assert_equal(block_array([[A.tocsr(), E.T.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), E.T.tocsc()],
+        assert_equal(block_array([[A.tocsc(), E.T.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
-        Z = csr_matrix((1, 1), dtype=np.int32)
+        Z = csr_array((1, 1), dtype=np.int32)
         expected = array([[0, 5],
                           [0, 6],
                           [7, 0]])
-        assert_equal(construct.bmat([[None, B], [C, None]]).toarray(),
+        assert_equal(block_array([[None, B], [C, None]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsr(), B.tocsr()],
+        assert_equal(block_array([[E.T.tocsr(), B.tocsr()],
                                      [C.tocsr(), Z]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsc(), B.tocsc()],
+        assert_equal(block_array([[E.T.tocsc(), B.tocsc()],
                                      [C.tocsc(), Z.tocsc()]]).toarray(),
                      expected)
 
-        expected = matrix(np.empty((0, 0)))
-        assert_equal(construct.bmat([[None, None]]).toarray(), expected)
-        assert_equal(construct.bmat([[None, D], [D, None]]).toarray(),
+        expected = np.empty((0, 0))
+        assert_equal(block_array([[None, None]]).toarray(), expected)
+        assert_equal(block_array([[None, D], [D, None]]).toarray(),
                      expected)
 
         # test bug reported in gh-5976
         expected = array([[7]])
-        assert_equal(construct.bmat([[None, D], [C, None]]).toarray(),
+        assert_equal(block_array([[None, D], [C, None]]).toarray(),
                      expected)
 
         # test failure cases
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A], [B]])
+            block_array([[A], [B]])
         excinfo.match(r'Got blocks\[1,0\]\.shape\[1\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr()], [B.tocsr()]])
+            block_array([[A.tocsr()], [B.tocsr()]])
         excinfo.match(r'incompatible dimensions for axis 1')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc()], [B.tocsc()]])
+            block_array([[A.tocsc()], [B.tocsc()]])
         excinfo.match(r'Mismatching dimensions along axis 1: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A, C]])
+            block_array([[A, C]])
         excinfo.match(r'Got blocks\[0,1\]\.shape\[0\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr(), C.tocsr()]])
+            block_array([[A.tocsr(), C.tocsr()]])
         excinfo.match(r'Mismatching dimensions along axis 0: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc(), C.tocsc()]])
+            block_array([[A.tocsc(), C.tocsc()]])
         excinfo.match(r'incompatible dimensions for axis 0')
 
+    def test_block_array_return_type(self):
+        block = construct.block_array
+
+        Flist, Glist = [[1, 2],[3, 4]], [[7], [5]]
+        Fm, Gm = csr_matrix(Flist), csr_matrix(Glist)
+        assert block([[None, Flist], [Glist, None]], format="csr")._is_array
+        assert block([[None, Fm], [Gm, None]], format="csr")._is_array
+        assert block([[Fm, Gm]], format="csr")._is_array
+
+    def test_bmat_return_type(self):
+        """This can be removed after sparse matrix is removed"""
+        block = construct.bmat
+        # check return type. if any input _is_array output array, else matrix
+        Flist, Glist = [[1, 2],[3, 4]], [[7], [5]]
+        Fm, Gm = csr_matrix(Flist), csr_matrix(Glist)
+        Fa, Ga = csr_array(Flist), csr_array(Glist)
+        assert block([[Fa, Ga]], format="csr")._is_array
+        assert not block([[Fm, Gm]], format="csr")._is_array
+        assert block([[None, Fa], [Ga, None]], format="csr")._is_array
+        assert block([[None, Fm], [Ga, None]], format="csr")._is_array
+        assert not block([[None, Fm], [Gm, None]], format="csr")._is_array
+        assert not block([[None, Flist], [Glist, None]], format="csr")._is_array
+
+        # type returned by _compressed_sparse_stack (all csr)
+        assert block([[Ga, Ga]], format="csr")._is_array
+        assert block([[Gm, Ga]], format="csr")._is_array
+        assert block([[Ga, Gm]], format="csr")._is_array
+        assert not block([[Gm, Gm]], format="csr")._is_array
+        # shape is 2x2 so no _stack_along_minor_axis
+        assert block([[Fa, Fm]], format="csr")._is_array
+        assert not block([[Fm, Fm]], format="csr")._is_array
+
+        # type returned by _compressed_sparse_stack (all csc)
+        assert block([[Gm.tocsc(), Ga.tocsc()]], format="csc")._is_array
+        assert not block([[Gm.tocsc(), Gm.tocsc()]], format="csc")._is_array
+        # shape is 2x2 so no _stack_along_minor_axis
+        assert block([[Fa.tocsc(), Fm.tocsc()]], format="csr")._is_array
+        assert not block([[Fm.tocsc(), Fm.tocsc()]], format="csr")._is_array
+
+        # type returned when mixed input
+        assert block([[Glist, Ga]], format="csr")._is_array
+        assert block([[Gm.tocsc(), Ga]], format="csr")._is_array
+        assert not block([[Gm.tocsc(), Gm]], format="csr")._is_array
+        assert not block([[Gm, Gm]], format="csc")._is_array
+
     @pytest.mark.slow
     @pytest.mark.xfail_on_32bit("Can't create large array for test")
     def test_concatenate_int32_overflow(self):
@@ -471,7 +527,7 @@ def test_concatenate_int32_overflow(self):
         check_free_memory(30000)
 
         n = 33000
-        A = csr_matrix(np.ones((n, n), dtype=bool))
+        A = csr_array(np.ones((n, n), dtype=bool))
         B = A.copy()
         C = construct._compressed_sparse_stack((A,B), 0)
 
@@ -481,9 +537,9 @@ def test_concatenate_int32_overflow(self):
 
     def test_block_diag_basic(self):
         """ basic test for block_diag """
-        A = coo_matrix([[1,2],[3,4]])
-        B = coo_matrix([[5],[6]])
-        C = coo_matrix([[7]])
+        A = coo_array([[1,2],[3,4]])
+        B = coo_array([[5],[6]])
+        C = coo_array([[7]])
 
         expected = array([[1, 2, 0, 0],
                           [3, 4, 0, 0],
@@ -511,21 +567,26 @@ def test_block_diag_1(self):
         assert_equal(construct.block_diag([1]).toarray(),
                      array([[1]]))
 
-    def test_block_diag_sparse_matrices(self):
-        """ block_diag with sparse matrices """
+    def test_block_diag_sparse_arrays(self):
+        """ block_diag with sparse arrays """
 
-        sparse_col_matrices = [coo_matrix(([[1, 2, 3]]), shape=(1, 3)),
-                               coo_matrix(([[4, 5]]), shape=(1, 2))]
-        block_sparse_cols_matrices = construct.block_diag(sparse_col_matrices)
-        assert_equal(block_sparse_cols_matrices.toarray(),
+        A = coo_array([[1, 2, 3]], shape=(1, 3))
+        B = coo_array([[4, 5]], shape=(1, 2))
+        assert_equal(construct.block_diag([A, B]).toarray(),
                      array([[1, 2, 3, 0, 0], [0, 0, 0, 4, 5]]))
 
-        sparse_row_matrices = [coo_matrix(([[1], [2], [3]]), shape=(3, 1)),
-                               coo_matrix(([[4], [5]]), shape=(2, 1))]
-        block_sparse_row_matrices = construct.block_diag(sparse_row_matrices)
-        assert_equal(block_sparse_row_matrices.toarray(),
+        A = coo_array([[1], [2], [3]], shape=(3, 1))
+        B = coo_array([[4], [5]], shape=(2, 1))
+        assert_equal(construct.block_diag([A, B]).toarray(),
                      array([[1, 0], [2, 0], [3, 0], [0, 4], [0, 5]]))
 
+    def test_block_diag_return_type(self):
+        A, B = coo_array([[1, 2, 3]]), coo_matrix([[2, 3, 4]])
+        assert construct.block_diag([A, A])._is_array
+        assert construct.block_diag([A, B])._is_array
+        assert construct.block_diag([B, A])._is_array
+        assert not construct.block_diag([B, B])._is_array
+
     def test_random_sampling(self):
         # Simple sanity checks for sparse random sampling.
         for f in sprand, _sprandn:


Code review comments:
diff_hunk: @@ -681,7 +707,60 @@ def bmat(blocks, format=None, dtype=None):
            [0, 0, 7]])
 
     """
+    blocks = np.asarray(blocks, dtype='object')
+    if any(isinstance(b, sparray) for b in blocks.flat):
+        return _block_array(blocks, format, dtype, coo_array)
+    else:
+        return _block_array(blocks, format, dtype, coo_matrix)
 
+
+def block_array(blocks, *, format=None, dtype=None):
path: scipy/sparse/_construct.py
text: @user1:
I would like to push back a bit on `block_array` as the name. I think the name is a little confusing considering the classes are all named `{*}_array` and blocked arrays (and block sparse arrays) are a thing.

Additionally, numpy's `block` (as well as `dask.array.block`) has similar behavior to what we want:

```python
import numpy as np
from scipy import sparse


arrays = [
    [np.ones((3, 3)), np.ones((3, 3)) * 2],
    [np.zeros((2, 3)), np.ones((2, 3))]
]

assert np.array_equal(np.block(arrays), sparse.bmat(arrays).toarray())

sp_arrays = [[sparse.csr_matrix(a) for a in row] for row in arrays]

assert np.array_equal(np.block(arrays), sparse.bmat(sp_arrays).toarray())
```

so why not `sparse.block` now?

Looking for inspiration in other places: Julia has `hvstack` as well as `hvnstack`

@author:
I believe the idea was to rename the construction functions using e.g. `eye_array` when they return arrays. I did think about using `block` to match other libraries. But I decided it was a poor choice of name because the word block doesn't say what we are doing/returning/wanting and it means too many things. In english "block" can be a verb, noun, adjective which can be confusing. It's an adjective instead of a noun or verb in this case. `block_array` at least makes that clear. And it describes what we expect to get from this. 

@user1:
> I believe the idea was to rename the construction functions using e.g. eye_array when they return arrays.

Yes, though in all those cases there is an equivalent existing function that returns matrices. I believe we also intend to change those names later (remove the `_array`), but we are not doing this for the concatenation functions, e.g. `vstack_array`.

> In english "block" can be a verb, noun, adjective which can be confusing.

I also don't love the name `block`. Personal preference would be something like `block_concat`. But, I do think there is a lot of value to naming functions the same way as `numpy` and `dask`.

@user2:
Typical famous names for these are `spdiags`, and other `spXXXX` stuff following mathworks. And provides a better separation from the dense ones.  `sp_concat` might be an option but a very loose suggestion, please ignore freely. 

@user3:
Given that we already have `scipy.sparse.hstack` and `scipy.sparse.vstack`, I think that `scipy.sparse.block` makes sense here. We do have the `diags`/`diags_array` situation at present, but eventually the plan is to end up renaming `diags_array` back to `diags`, so if we don't need the _array suffix here then I'm happy to avoid it from the start.

I'd like to stay away from Matlab-isms if possible, especially if we can instead align more closely with the numpy-style naming conventions.

@author:
I changed the name to `block` instead of `block_array`.

@user1:
Thanks for the change!

diff_hunk: @@ -690,29 +769,29 @@ def bmat(blocks, format=None, dtype=None):
     M,N = blocks.shape
 
     # check for fast path cases
-    if (format in (None, 'csr') and all(isinstance(b, csr_matrix)
-                                        for b in blocks.flat)):
+    if (format in (None, 'csr') and
+        all(issparse(b) and b.format == "csr" for b in blocks.flat)
path: scipy/sparse/_construct.py
text: @user1:
Minor suggestion:

```suggestion
        all(isinstance(b, (csr_matrix, csr_array)) for b in blocks.flat)
```

I think it's a cleaner and more pythonic way to identify types. Similar elsewhere if you agree on this.

@author:
This question is also bigger than this issue.  I know that you prefer type checking rather than attribute checking. But that is not universal and I'm just going with what I understand is the standard way to check.  The current form will not require any change after `csr_matrix` is removed. (I don't know the other parts of that discussion.)

@user1:
I'm okay with either here. I think both versions have reasons to exist, so just preference based.

diff_hunk: @@ -340,8 +340,8 @@ def test_kronsum(self):
 
     def test_vstack(self):
path: scipy/sparse/tests/test_construct.py
text: @user2:
To maintain coverage on matrix cases, I would suggest:

```python
    @user1.mark.parametrize("coo_cls", [coo_matrix, coo_array])
    def test_vstack(self, coo_cls):
        result_cls = sparray if issubclass(coo_cls, sparray) else spmatrix
        A = coo_cls([[1,2],[3,4]])
        B = coo_cls([[5,6]])

        ...


        assert isinstance(construct.vstack([A, B]), result_cls)
        assert isinstance(construct.vstack([coo_matrix(A), coo_matrix(B)]), spmatrix)
        assert isinstance(construct.vstack([A, coo_matrix(B)]), result_cls)
        assert isinstance(construct.vstack([coo_matrix(A), B]), result_cls)
```

@author:
I've been warned that the combinatorics of the sparse subpackage makes testing expensive.  So, I'd like to avoid duplication of testing where possible.
Here, we test the base code by testing the arrays interface. And then we test the spmatrix aspects of the code. So I think the coverage is the same as your suggestion of more tests.

@user2:
> I've been warned that the combinatorics of the sparse subpackage makes testing expensive. So, I'd like to avoid duplication of testing where possible.

I think this can be true, but not sure this case will have much effect. Trying to time this test:

<details>
<summary> Full output </summary>

(scipy-dev) isaac@user3:~/github/scipy ‹block_stack›
$ python dev.py test -s sparse  -- -k test_vstack --durations=0 -vv
💻  ninja -C /Users/isaac/github/scipy/build -j8
ninja: Entering directory `/Users/isaac/github/scipy/build'
[4/4] Generating scipy/generate-version with a custom command
Build OK
💻  meson install -C build --only-changed
Installing, see meson-install.log...
Installation OK
SciPy from development installed path at: /Users/isaac/github/scipy/build-install/lib/python3.10/site-packages
Running tests for scipy version:1.12.0.dev0+1321.c969aa1, installed at:/Users/isaac/github/scipy/build-install/lib/python3.10/site-packages/scipy
============================================================================ test session starts ============================================================================
platform darwin -- Python 3.10.11, pytest-7.3.1, pluggy-1.0.0 -- /Users/isaac/miniconda3/envs/scipy-dev/bin/python3
cachedir: .pytest_cache
hypothesis profile 'default' -> database=DirectoryBasedExampleDatabase('/Users/isaac/github/scipy/build-install/lib/python3.10/site-packages/.hypothesis/examples')
metadata: {'Python': '3.10.11', 'Platform': 'macOS-13.4.1-x86_64-i386-64bit', 'Packages': {'pytest': '7.3.1', 'pluggy': '1.0.0'}, 'Plugins': {'timeout': '2.1.0', 'json-report': '1.5.0', 'xdist': '3.3.1', 'cov': '4.0.0', 'hypothesis': '6.81.2', 'metadata': '2.0.4'}}
rootdir: /Users/isaac/github/scipy
configfile: pytest.ini
plugins: timeout-2.1.0, json-report-1.5.0, xdist-3.3.1, cov-4.0.0, hypothesis-6.81.2, metadata-2.0.4
collected 10535 items / 10534 deselected / 1 selected                                                                                                                       

scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack PASSED                                                                                          [100%]

============================================================================= slowest durations =============================================================================
0.00s setup    build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
0.00s call     build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
0.00s teardown build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
==================================================================== 1 passed, 10534 deselected in 1.62s ====================================================================

</details>

```
=============================================================================
0.00s setup    build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
0.00s call     build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
0.00s teardown build-install/lib/python3.10/site-packages/scipy/sparse/tests/test_construct.py::TestConstructUtils::test_vstack
```

> Here, we test the base code by testing the arrays interface. And then we test the spmatrix aspects of the code. So I think the coverage is the same as your suggestion of more tests.

Yes, but this assumes the underlying codepath for `vstack` on matrices is the same as arrays, and that the wiring to get to that codepath is correct. The changes here would ensure correctness without having to rely on that assumption.

Broadly, my idea here is "testing the interface instead of the implementation"

@author:
I've parametrized the vstack and hstack tests to run for both matrix and array classes. I separated the tests of return matrix/array into a separate test method.

diff_hunk: @@ -385,93 +390,144 @@ def test_hstack(self):
                                       dtype=np.float32).dtype,
                      np.float32)
 
-    def test_bmat(self):
+        assert construct.hstack([A, B])._is_array
+        assert not construct.hstack([coo_matrix(A), coo_matrix(B)])._is_array
+        assert construct.hstack([A, coo_matrix(B)])._is_array
+        assert construct.hstack([coo_matrix(A), B])._is_array
+
+    @pytest.mark.parametrize("block_array", (construct.bmat, construct.block_array))
+    def test_block_creation(self, block_array):
 
-        A = coo_matrix([[1, 2], [3, 4]])
-        B = coo_matrix([[5],[6]])
-        C = coo_matrix([[7]])
-        D = coo_matrix((0, 0))
+        A = coo_array([[1, 2], [3, 4]])
+        B = coo_array([[5],[6]])
+        C = coo_array([[7]])
+        D = coo_array((0, 0))
 
         expected = array([[1, 2, 5],
                           [3, 4, 6],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, B], [None, C]]).toarray(), expected)
-        E = csr_matrix((1, 2), dtype=np.int32)
-        assert_equal(construct.bmat([[A.tocsr(), B.tocsr()],
+        assert_equal(block_array([[A, B], [None, C]]).toarray(), expected)
+        E = csr_array((1, 2), dtype=np.int32)
+        assert_equal(block_array([[A.tocsr(), B.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), B.tocsc()],
+        assert_equal(block_array([[A.tocsc(), B.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
         expected = array([[1, 2, 0],
                           [3, 4, 0],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, None], [None, C]]).toarray(),
+        assert_equal(block_array([[A, None], [None, C]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsr(), E.T.tocsr()],
+        assert_equal(block_array([[A.tocsr(), E.T.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), E.T.tocsc()],
+        assert_equal(block_array([[A.tocsc(), E.T.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
-        Z = csr_matrix((1, 1), dtype=np.int32)
+        Z = csr_array((1, 1), dtype=np.int32)
         expected = array([[0, 5],
                           [0, 6],
                           [7, 0]])
-        assert_equal(construct.bmat([[None, B], [C, None]]).toarray(),
+        assert_equal(block_array([[None, B], [C, None]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsr(), B.tocsr()],
+        assert_equal(block_array([[E.T.tocsr(), B.tocsr()],
                                      [C.tocsr(), Z]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsc(), B.tocsc()],
+        assert_equal(block_array([[E.T.tocsc(), B.tocsc()],
                                      [C.tocsc(), Z.tocsc()]]).toarray(),
                      expected)
 
-        expected = matrix(np.empty((0, 0)))
-        assert_equal(construct.bmat([[None, None]]).toarray(), expected)
-        assert_equal(construct.bmat([[None, D], [D, None]]).toarray(),
+        expected = np.empty((0, 0))
+        assert_equal(block_array([[None, None]]).toarray(), expected)
+        assert_equal(block_array([[None, D], [D, None]]).toarray(),
                      expected)
 
         # test bug reported in gh-5976
         expected = array([[7]])
-        assert_equal(construct.bmat([[None, D], [C, None]]).toarray(),
+        assert_equal(block_array([[None, D], [C, None]]).toarray(),
                      expected)
 
         # test failure cases
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A], [B]])
+            block_array([[A], [B]])
path: scipy/sparse/tests/test_construct.py
text: @user1:
This is nitpicky and probably out of scope. But while you're here...

I would suggest: change `assert_raises` to `pytest.raises` (which is a much more common usage that `from pytest import raises as assert_raises`) and calling it like:

`with pytest.raises(ValueError, match=r"..."):` so you don't need the excinfo on the next line.

You could also leave this for me.

@author:
I'd prefer to leave this suggestion for another PR.

diff_hunk: @@ -385,93 +390,144 @@ def test_hstack(self):
                                       dtype=np.float32).dtype,
                      np.float32)
 
-    def test_bmat(self):
+        assert construct.hstack([A, B])._is_array
+        assert not construct.hstack([coo_matrix(A), coo_matrix(B)])._is_array
+        assert construct.hstack([A, coo_matrix(B)])._is_array
+        assert construct.hstack([coo_matrix(A), B])._is_array
+
+    @pytest.mark.parametrize("block_array", (construct.bmat, construct.block_array))
+    def test_block_creation(self, block_array):
 
-        A = coo_matrix([[1, 2], [3, 4]])
-        B = coo_matrix([[5],[6]])
-        C = coo_matrix([[7]])
-        D = coo_matrix((0, 0))
+        A = coo_array([[1, 2], [3, 4]])
+        B = coo_array([[5],[6]])
+        C = coo_array([[7]])
+        D = coo_array((0, 0))
 
         expected = array([[1, 2, 5],
                           [3, 4, 6],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, B], [None, C]]).toarray(), expected)
-        E = csr_matrix((1, 2), dtype=np.int32)
-        assert_equal(construct.bmat([[A.tocsr(), B.tocsr()],
+        assert_equal(block_array([[A, B], [None, C]]).toarray(), expected)
+        E = csr_array((1, 2), dtype=np.int32)
+        assert_equal(block_array([[A.tocsr(), B.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), B.tocsc()],
+        assert_equal(block_array([[A.tocsc(), B.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
         expected = array([[1, 2, 0],
                           [3, 4, 0],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, None], [None, C]]).toarray(),
+        assert_equal(block_array([[A, None], [None, C]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsr(), E.T.tocsr()],
+        assert_equal(block_array([[A.tocsr(), E.T.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), E.T.tocsc()],
+        assert_equal(block_array([[A.tocsc(), E.T.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
-        Z = csr_matrix((1, 1), dtype=np.int32)
+        Z = csr_array((1, 1), dtype=np.int32)
         expected = array([[0, 5],
                           [0, 6],
                           [7, 0]])
-        assert_equal(construct.bmat([[None, B], [C, None]]).toarray(),
+        assert_equal(block_array([[None, B], [C, None]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsr(), B.tocsr()],
+        assert_equal(block_array([[E.T.tocsr(), B.tocsr()],
                                      [C.tocsr(), Z]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsc(), B.tocsc()],
+        assert_equal(block_array([[E.T.tocsc(), B.tocsc()],
                                      [C.tocsc(), Z.tocsc()]]).toarray(),
                      expected)
 
-        expected = matrix(np.empty((0, 0)))
-        assert_equal(construct.bmat([[None, None]]).toarray(), expected)
-        assert_equal(construct.bmat([[None, D], [D, None]]).toarray(),
+        expected = np.empty((0, 0))
+        assert_equal(block_array([[None, None]]).toarray(), expected)
+        assert_equal(block_array([[None, D], [D, None]]).toarray(),
                      expected)
 
         # test bug reported in gh-5976
         expected = array([[7]])
-        assert_equal(construct.bmat([[None, D], [C, None]]).toarray(),
+        assert_equal(block_array([[None, D], [C, None]]).toarray(),
                      expected)
 
         # test failure cases
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A], [B]])
+            block_array([[A], [B]])
         excinfo.match(r'Got blocks\[1,0\]\.shape\[1\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr()], [B.tocsr()]])
+            block_array([[A.tocsr()], [B.tocsr()]])
         excinfo.match(r'incompatible dimensions for axis 1')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc()], [B.tocsc()]])
+            block_array([[A.tocsc()], [B.tocsc()]])
         excinfo.match(r'Mismatching dimensions along axis 1: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A, C]])
+            block_array([[A, C]])
         excinfo.match(r'Got blocks\[0,1\]\.shape\[0\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr(), C.tocsr()]])
+            block_array([[A.tocsr(), C.tocsr()]])
         excinfo.match(r'Mismatching dimensions along axis 0: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc(), C.tocsc()]])
+            block_array([[A.tocsc(), C.tocsc()]])
         excinfo.match(r'incompatible dimensions for axis 0')
 
+    def test_block_array_return_type(self):
+        block = construct.block_array
+
+        Flist, Glist = [[1, 2],[3, 4]], [[7], [5]]
+        Fm, Gm = csr_matrix(Flist), csr_matrix(Glist)
+        assert block([[None, Flist], [Glist, None]], format="csr")._is_array
+        assert block([[None, Fm], [Gm, None]], format="csr")._is_array
+        assert block([[Fm, Gm]], format="csr")._is_array
+
+    def test_bmat_return_type(self):
+        """This can be removed after sparse matrix is removed"""
+        block = construct.bmat
path: scipy/sparse/tests/test_construct.py
text: @user1:
```suggestion
        block = construct.bmat
```

This confused me while looking through the test. Can we use the `bmat` identifier here?

diff_hunk: @@ -507,11 +507,20 @@ def _compressed_sparse_stack(blocks, axis):
         sum_dim += b.shape[axis]
         last_indptr += b.indptr[-1]
     indptr[-1] = last_indptr
+    # TODO remove this if-structure when sparse matrices removed
+    if container is coo_matrix:
path: scipy/sparse/_construct.py
text: @user1:
Should this be `if isinstance(container, spmatrix)`?

@user1:
Or I suppose maybe `container` could be changed to a boolean `return_sparse_array` parameter, given how you're calling this. That would make this part simpler, I think.

@user2:
Alternatively, could this always return an `sparray`, then have `bmat` handle the conversion to matrix?

So `bmat` becomes something like:

```python
    result = _block_array([blocks], format, dtype)
    if any(isinstance(b, sparray) for b in blocks.flat):
        return result
    else:
        return getattr(sparse, f"{result.format}_matrix")(result)
```

@author:
I changed the passed value to be boolean `return_spmatrix` with a default value of `False` so the array tests can remain the same when this testing removes spmatrix.

diff_hunk: @@ -481,9 +481,9 @@ def kronsum(A, B, format=None):
     return (L+R).asformat(format)  # since L + R is not always same format
 
 
-def _compressed_sparse_stack(blocks, axis):
+def _compressed_sparse_stack(blocks, axis, container):
path: scipy/sparse/_construct.py
text: @user1:
Could this get tested with 64 bit indices? I'm not actually sure the behavior would be correct here. I think you may also need to add the `indices` to `get_index_dtype`.

I *think* the behavior for indices here should be:

* upcast if the size requires it
* otherwise rely on type promotion
* never downcast

@author:
This PR doesn't change the `idx_dtype` behavior.
But also, I don't think we need to add `indices` to `get_index_dtype` because `indptr` will indicate the current type and the max value needed. So I think this is OK as is...


@user1:
This PR enables concatenation for the array classes, which have different index behavior that matrix classes. So there will be no existing tests for index dtype behavior.

I believe `indptr` will only indicate the correct index type after certain operations that check for this. `indptr` and `indices` sharing a `dtype` is not an invariant of these classes:

```python
csr = sparse.random(10, 10, density=0.25, format="csr")
csr.indices = csr.indices.astype(np.int64)
(csr.indices.dtype, csr.indptr.dtype)
# (dtype('int64'), dtype('int32'))
```

---

Here's an example case:

```python
from scipy import sparse
import numpy as np

csr1 = sparse.csr_array(sparse.random(10, 10, density=0.25, format="csr"))
csr2 = sparse.csr_array(sparse.random(10, 10, density=0.25, format="csr"))
csr1.indices = csr1.indices.astype(np.int64)

result = sparse.vstack([csr1, csr2])
result.indices.dtype
```

What should the resulting index dtype be? I *think* it should be `int64` following type promotion rules.

The alternative is to downcast to the smallest possible size, but this is the kind of behavior we removed.

diff_hunk: @@ -385,93 +390,144 @@ def test_hstack(self):
                                       dtype=np.float32).dtype,
                      np.float32)
 
-    def test_bmat(self):
+        assert construct.hstack([A, B])._is_array
+        assert not construct.hstack([coo_matrix(A), coo_matrix(B)])._is_array
+        assert construct.hstack([A, coo_matrix(B)])._is_array
+        assert construct.hstack([coo_matrix(A), B])._is_array
+
+    @pytest.mark.parametrize("block_array", (construct.bmat, construct.block_array))
+    def test_block_creation(self, block_array):
 
-        A = coo_matrix([[1, 2], [3, 4]])
-        B = coo_matrix([[5],[6]])
-        C = coo_matrix([[7]])
-        D = coo_matrix((0, 0))
+        A = coo_array([[1, 2], [3, 4]])
+        B = coo_array([[5],[6]])
+        C = coo_array([[7]])
+        D = coo_array((0, 0))
 
         expected = array([[1, 2, 5],
                           [3, 4, 6],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, B], [None, C]]).toarray(), expected)
-        E = csr_matrix((1, 2), dtype=np.int32)
-        assert_equal(construct.bmat([[A.tocsr(), B.tocsr()],
+        assert_equal(block_array([[A, B], [None, C]]).toarray(), expected)
+        E = csr_array((1, 2), dtype=np.int32)
+        assert_equal(block_array([[A.tocsr(), B.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), B.tocsc()],
+        assert_equal(block_array([[A.tocsc(), B.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
         expected = array([[1, 2, 0],
                           [3, 4, 0],
                           [0, 0, 7]])
-        assert_equal(construct.bmat([[A, None], [None, C]]).toarray(),
+        assert_equal(block_array([[A, None], [None, C]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsr(), E.T.tocsr()],
+        assert_equal(block_array([[A.tocsr(), E.T.tocsr()],
                                      [E, C.tocsr()]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[A.tocsc(), E.T.tocsc()],
+        assert_equal(block_array([[A.tocsc(), E.T.tocsc()],
                                      [E.tocsc(), C.tocsc()]]).toarray(),
                      expected)
 
-        Z = csr_matrix((1, 1), dtype=np.int32)
+        Z = csr_array((1, 1), dtype=np.int32)
         expected = array([[0, 5],
                           [0, 6],
                           [7, 0]])
-        assert_equal(construct.bmat([[None, B], [C, None]]).toarray(),
+        assert_equal(block_array([[None, B], [C, None]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsr(), B.tocsr()],
+        assert_equal(block_array([[E.T.tocsr(), B.tocsr()],
                                      [C.tocsr(), Z]]).toarray(),
                      expected)
-        assert_equal(construct.bmat([[E.T.tocsc(), B.tocsc()],
+        assert_equal(block_array([[E.T.tocsc(), B.tocsc()],
                                      [C.tocsc(), Z.tocsc()]]).toarray(),
                      expected)
 
-        expected = matrix(np.empty((0, 0)))
-        assert_equal(construct.bmat([[None, None]]).toarray(), expected)
-        assert_equal(construct.bmat([[None, D], [D, None]]).toarray(),
+        expected = np.empty((0, 0))
+        assert_equal(block_array([[None, None]]).toarray(), expected)
+        assert_equal(block_array([[None, D], [D, None]]).toarray(),
                      expected)
 
         # test bug reported in gh-5976
         expected = array([[7]])
-        assert_equal(construct.bmat([[None, D], [C, None]]).toarray(),
+        assert_equal(block_array([[None, D], [C, None]]).toarray(),
                      expected)
 
         # test failure cases
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A], [B]])
+            block_array([[A], [B]])
         excinfo.match(r'Got blocks\[1,0\]\.shape\[1\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr()], [B.tocsr()]])
+            block_array([[A.tocsr()], [B.tocsr()]])
         excinfo.match(r'incompatible dimensions for axis 1')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc()], [B.tocsc()]])
+            block_array([[A.tocsc()], [B.tocsc()]])
         excinfo.match(r'Mismatching dimensions along axis 1: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A, C]])
+            block_array([[A, C]])
         excinfo.match(r'Got blocks\[0,1\]\.shape\[0\] == 1, expected 2')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsr(), C.tocsr()]])
+            block_array([[A.tocsr(), C.tocsr()]])
         excinfo.match(r'Mismatching dimensions along axis 0: ({1, 2}|{2, 1})')
 
         with assert_raises(ValueError) as excinfo:
-            construct.bmat([[A.tocsc(), C.tocsc()]])
+            block_array([[A.tocsc(), C.tocsc()]])
         excinfo.match(r'incompatible dimensions for axis 0')
 
+    def test_block_array_return_type(self):
+        block = construct.block_array
+
+        Flist, Glist = [[1, 2],[3, 4]], [[7], [5]]
+        Fm, Gm = csr_matrix(Flist), csr_matrix(Glist)
+        assert block([[None, Flist], [Glist, None]], format="csr")._is_array
+        assert block([[None, Fm], [Gm, None]], format="csr")._is_array
+        assert block([[Fm, Gm]], format="csr")._is_array
path: scipy/sparse/tests/test_construct.py
text: @user1:
```suggestion
        assert isinstance(block([[None, Flist], [Glist, None]], format="csr"), sparse.csr_array)
        assert isinstance(block([[None, Fm], [Gm, None]], format="csr"), sparse.csr_array)
        assert isinstance(block([[Fm, Gm]], format="csr"), sparse.csr_array)
```

Considering you are specifiying the return type in the call, I think it would be worthwhile to test the exact result type.

It could also be worth doing a combinatorial test checking that all cases are at least `isinstance(block(...), sparray)`

@author:
I prefer using the attributes (both `.format` and `._is_array`) as I believe they are the intended mechanism to get this kind of information. But I could be wrong...

@user1:
I've opened https://github.com/scipy/scipy/issues/18921 to discuss `_is_array`, so happy to defer that for now.

I think the point about "why have the call specify `"csr"` if you don't check it" still stands. I think the check would make sense as is if the test was parameterized by `format`.

@author:
The "csr" format is used because that format allows us to test the special as well as general cases in the code. I understand your desire to test the api rather than the code, and I'm fine with approaching testing that way. But my understanding is that the combinatorial explosion is hurting our existing tests (see test_base.py). With so many formats and so many different dtypes we end up with a lot of tests.

This test is not intended to test the `format` option. It is intended to test the return object. If we decide we should test all the formats, then we could easily decide to test all the dtypes as well -- even though this test isn't intended to test handling of dtype or format. What's the right approach to testing here? 

@user1:
> The "csr" format is used because that format allows us to test the special as well as general cases in the code

Thanks for the additional context. My confusion here was that I didn't understand why the format was being specified. If the intent is to hit a certain codepath, that makes sense. I do think a comment explaining this could help so future editors know if/ when this should be changed.

I do think it would make sense to test the format somewhere (which I assume happens?), and that those tests *could* also cover the array/matrix types at the same time via instance checks. Looking locally it seems like there is coverage here so not a necessary change.