diff --git a/dedalus/core/coords.py b/dedalus/core/coords.py
index b4d78b6b..39ffa2c9 100644
--- a/dedalus/core/coords.py
+++ b/dedalus/core/coords.py
@@ -138,7 +138,7 @@ def forward_vector_intertwiner(self, subaxis, group):
             else:
                 factors.append(np.identity(cs.dim))
             start_axis += cs.dim
-        return sparse_block_diag(factors).A
+        return sparse_block_diag(factors).toarray()
 
     def backward_vector_intertwiner(self, subaxis, group):
         factors = []
@@ -149,7 +149,7 @@ def backward_vector_intertwiner(self, subaxis, group):
             else:
                 factors.append(np.identity(cs.dim))
             start_axis += cs.dim
-        return sparse_block_diag(factors).A
+        return sparse_block_diag(factors).toarray()
 
     @CachedAttribute
     def default_nonconst_groups(self):
diff --git a/dedalus/core/solvers.py b/dedalus/core/solvers.py
index 364eea56..3de09ee3 100644
--- a/dedalus/core/solvers.py
+++ b/dedalus/core/solvers.py
@@ -172,8 +172,8 @@ def print_subproblem_ranks(self, subproblems=None, target=0):
         for i, sp in enumerate(subproblems):
             if not hasattr(sp, 'L_min'):
                 continue
-            L = sp.L_min.A
-            M = sp.M_min.A
+            L = sp.L_min.toarray()
+            M = sp.M_min.toarray()
             A = L + target * M
             print(f"MPI rank: {self.dist.comm.rank}, subproblem: {i}, group: {sp.group}, matrix rank: {np.linalg.matrix_rank(A)}/{A.shape[0]}, cond: {np.linalg.cond(A):.1e}")
 
@@ -205,15 +205,15 @@ def solve_dense(self, subproblem, rebuild_matrices=False, left=False, normalize_
         if rebuild_matrices or not hasattr(sp, 'L_min'):
             self.build_matrices([sp], ['M', 'L'])
         # Solve as dense general eigenvalue problem
-        A = sp.L_min.A
-        B = - sp.M_min.A
+        A = sp.L_min.toarray()
+        B = - sp.M_min.toarray()
         eig_output = scipy.linalg.eig(A, b=B, left=left, **kw)
         # Unpack output
         if left:
             self.eigenvalues, pre_left_evecs, pre_right_evecs = eig_output
             self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
-            self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
-            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
+            self.left_eigenvectors = sp.pre_left.conj().toarray() @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().toarray() @ pre_left_evecs
             if normalize_left:
                 norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
                 self.left_eigenvectors /= np.conj(norms)
@@ -270,8 +270,8 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
             # Note: this definition of "left eigenvectors" is consistent with the documentation for scipy.linalg.eig
             self.eigenvalues, pre_right_evecs, self.left_eigenvalues, pre_left_evecs = eig_output
             self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
-            self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
-            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
+            self.left_eigenvectors = sp.pre_left.conj().toarray() @ pre_left_evecs
+            self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).conj().toarray() @ pre_left_evecs
             # Check that eigenvalues match
             if not np.allclose(self.eigenvalues, np.conjugate(self.left_eigenvalues)):
                 if raise_on_mismatch:
@@ -363,7 +363,7 @@ def print_subproblem_ranks(self, subproblems=None):
         for i, sp in enumerate(subproblems):
             if not hasattr(sp, 'L_min'):
                 continue
-            L = sp.L_min.A
+            L = sp.L_min.toarray()
             print(f"MPI rank: {self.dist.comm.rank}, subproblem: {i}, group: {sp.group}, matrix rank: {np.linalg.matrix_rank(L)}/{L.shape[0]}, cond: {np.linalg.cond(L):.1e}")
 
     def solve(self, subproblems=None, rebuild_matrices=False):
@@ -464,7 +464,7 @@ def print_subproblem_ranks(self, subproblems=None):
         for i, sp in enumerate(subproblems):
             if not hasattr(sp, 'dF_min'):
                 continue
-            dF = sp.dF_min.A
+            dF = sp.dF_min.toarray()
             print(f"MPI rank: {self.dist.comm.rank}, subproblem: {i}, group: {sp.group}, matrix rank: {np.linalg.matrix_rank(dF)}/{dF.shape[0]}, cond: {np.linalg.cond(dF):.1e}")
 
     def newton_iteration(self, damping=1):
@@ -739,7 +739,7 @@ def print_subproblem_ranks(self, subproblems=None, dt=1):
         for i, sp in enumerate(subproblems):
             M = sp.M_min
             L = sp.L_min
-            A = (M + dt*L).A
+            A = (M + dt*L).toarray()
             print(f"MPI rank: {self.dist.comm.rank}, subproblem: {i}, group: {sp.group}, matrix rank: {np.linalg.matrix_rank(A)}/{A.shape[0]}, cond: {np.linalg.cond(A):.1e}")
 
     def evaluate_handlers_now(self, dt, handlers=None):
diff --git a/dedalus/libraries/dedalus_sphere/clenshaw.py b/dedalus/libraries/dedalus_sphere/clenshaw.py
index 01f9efaa..20e427f5 100644
--- a/dedalus/libraries/dedalus_sphere/clenshaw.py
+++ b/dedalus/libraries/dedalus_sphere/clenshaw.py
@@ -32,7 +32,7 @@ def jacobi_recursion(N, a, b, X):
     """
     # Jacobi matrix
     J = jacobi_matrix(N, a, b)
-    JA = J.A
+    JA = J.toarray()
     # Identity element
     if np.isscalar(X):
         I = 1
diff --git a/dedalus/libraries/dedalus_sphere/jacobi.py b/dedalus/libraries/dedalus_sphere/jacobi.py
index 635fd816..bba2d233 100644
--- a/dedalus/libraries/dedalus_sphere/jacobi.py
+++ b/dedalus/libraries/dedalus_sphere/jacobi.py
@@ -132,7 +132,7 @@ def grid_guess(n,a,b,dtype='longdouble',quick=False):
         quick = True
 
     if n == 1:
-        return operator('Z')(n,a,b).A[0]
+        return operator('Z')(n,a,b).toarray()[0]
 
     if quick:
         return np.cos(np.pi*(np.arange(4*n-1,2,-4,dtype=dtype)+2*a)/(4*n+2*(a+b+1)))
diff --git a/dedalus/libraries/dedalus_sphere/tests/test_jacobi.py b/dedalus/libraries/dedalus_sphere/tests/test_jacobi.py
index f5d994f4..7973480b 100644
--- a/dedalus/libraries/dedalus_sphere/tests/test_jacobi.py
+++ b/dedalus/libraries/dedalus_sphere/tests/test_jacobi.py
@@ -87,7 +87,7 @@ def test_Ap_Bp_commutation(N, a, b):
     Bp00 = jacobi128.operator('B+', N, a, b)
     Ap01 = jacobi128.operator('A+', N, a, b+1)
     path2 = Ap01 @ Bp00
-    assert np.allclose(path1.A, path2.A)
+    assert np.allclose(path1.toarray(), path2.toarray())
 
 
 @pytest.mark.parametrize('N', N_range)
@@ -105,5 +105,5 @@ def test_App_Bpp_commutation(N, a, b):
     Ap02 = jacobi128.operator('A+', N, a, b+2)
     Ap12 = jacobi128.operator('A+', N, a+1, b+2)
     path2 = Ap12 @ Ap02 @ Bp01 @ Bp00
-    assert np.allclose(path1.A, path2.A)
+    assert np.allclose(path1.toarray(), path2.toarray())
 
diff --git a/dedalus/libraries/matsolvers.py b/dedalus/libraries/matsolvers.py
index 082e7e34..811b4542 100644
--- a/dedalus/libraries/matsolvers.py
+++ b/dedalus/libraries/matsolvers.py
@@ -137,7 +137,7 @@ def __init__(self, matrix, solver=None):
         if self.trans == "T":
             matrix = matrix.T
         elif self.trans == "H":
-            matrix = matrix.H
+            matrix = matrix.conj().toarray()
         self.LU = spla.splu(matrix.tocsc(),
                             permc_spec=self.permc_spec,
                             diag_pivot_thresh=self.diag_pivot_thresh,
@@ -235,7 +235,7 @@ class DenseInverse(DenseSolver):
     """Dense inversion solve."""
 
     def __init__(self, matrix, solver=None):
-        self.matrix_inverse = sla.inv(matrix.A)
+        self.matrix_inverse = sla.inv(matrix.toarray())
 
     def solve(self, vector):
         return self.matrix_inverse @ vector
@@ -276,7 +276,7 @@ class ScipyDenseLU(DenseSolver):
     """Scipy dense LU factorized solve."""
 
     def __init__(self, matrix, solver=None):
-        self.LU = sla.lu_factor(matrix.A, check_finite=False)
+        self.LU = sla.lu_factor(matrix.toarray(), check_finite=False)
 
     def solve(self, vector):
         return sla.lu_solve(self.LU, vector, check_finite=False)
@@ -296,9 +296,9 @@ def __init__(self, matrix, subproblem, matsolver):
         self.V = V = np.zeros((2*R, matrix.shape[1]), dtype=matrix.dtype)
         # Remove top border, leaving upper left subblock
         U[:R, :R] = np.identity(R)
-        V[:R, R:] = matrix[:R, R:].A
+        V[:R, R:] = matrix[:R, R:].toarray()
         # Remove right border, leaving upper right and lower right subblocks
-        U[R:-R, R:] = matrix[R:-R, -R:].A
+        U[R:-R, R:] = matrix[R:-R, -R:].toarray()
         V[-R:, -R:] = np.identity(R)
         self.A = matrix - sp.csr_matrix(U) @ sp.csr_matrix(V)
         # Solve A using specified matsolver
diff --git a/dedalus/tests/test_nlbvp.py b/dedalus/tests/test_nlbvp.py
index 6d0446f3..a512ef22 100644
--- a/dedalus/tests/test_nlbvp.py
+++ b/dedalus/tests/test_nlbvp.py
@@ -324,7 +324,7 @@ def error(perts):
         R = -φcen  / (n+1)**(3/2)
         print(solver.iteration, φcen, R, err)
         dH = solver.subproblems[0].dH_min
-        print('%.2e' %np.linalg.cond(dH.A))
+        print('%.2e' %np.linalg.cond(dH.toarray()))
     if err > tolerance:
         raise ValueError("Did not converge")
     # Compare to reference solutions from Boyd
diff --git a/dedalus/tools/array.py b/dedalus/tools/array.py
index 0a196b7e..bf28cd73 100644
--- a/dedalus/tools/array.py
+++ b/dedalus/tools/array.py
@@ -433,7 +433,7 @@ def matvec(x):
         # Build sparse linear operator representing (A^H - conj(σ)B^H)^I B^H = C^-H B^H = left_D
         # Note: left_D is not equal to D^H
         def matvec_left(x):
-            return solver.solve_H(B.H.dot(x))
+            return solver.solve_H(B.conj().toarray().dot(x))
         left_D = spla.LinearOperator(dtype=A.dtype, shape=A.shape, matvec=matvec_left)
         # Solve using scipy sparse algorithm
         left_evals, left_evecs = spla.eigs(left_D, k=N, which='LM', sigma=None, **kw)
diff --git a/dedalus/tools/clenshaw.py b/dedalus/tools/clenshaw.py
index d1a40d70..af120157 100644
--- a/dedalus/tools/clenshaw.py
+++ b/dedalus/tools/clenshaw.py
@@ -79,7 +79,7 @@ def jacobi_recursion(N, a, b, X):
     """
     # Jacobi matrix
     J = jacobi.jacobi_matrix(N, a, b)
-    JA = J.A
+    JA = J.toarray()
     # Identity element
     if np.isscalar(X):
         I = 1