code: def rho_block_D_inv_A(A, Dinv):
      
      if not hasattr(A, 'rho_block_D_inv'):
          from scipy.sparse.linalg import LinearOperator
  
          blocksize = Dinv.shape[1]
          if Dinv.shape[1] != Dinv.shape[2]:
              raise ValueError('Dinv has incorrect dimensions')
          elif Dinv.shape[0] != int(A.shape[0]/blocksize):
              raise ValueError('Dinv and A have incompatible dimensions')
  
          Dinv = sp.sparse.bsr_matrix((Dinv,
                                       sp.arange(Dinv.shape[0]),
                                       sp.arange(Dinv.shape[0]+1)),
                                      shape=A.shape)
  
          # Don't explicitly form Dinv*A
          def matvec(x):
              return Dinv*(A*x)
          D_inv_A = LinearOperator(A.shape, matvec, dtype=A.dtype)
  
          A.rho_block_D_inv = approximate_spectral_radius(D_inv_A)
  
      return A.rho_block_D_inv
docs: Return the (approx.) spectral radius of block D^-1 * A.
  
      Parameters
      ----------
      A : sparse-matrix
          size NxN
      Dinv : array
          Inverse of diagonal blocks of A
          size (N/blocksize, blocksize, blocksize)
  
      Returns
      -------
      approximate spectral radius of (Dinv A)
  
      Examples
      --------
      >>> from pyamg.gallery import poisson
      >>> from pyamg.relaxation.smoothing import rho_block_D_inv_A
      >>> from pyamg.util.utils import get_block_diag
      >>> A = poisson((10,10), format='csr')
      >>> Dinv = get_block_diag(A, blocksize=4, inv_flag=True)