Python scipy.sparse 模块，kron() 实例源码

def sparse_cH(terms, ldim=2):
    """Construct a sparse cyclic nearest-neighbour Hamiltonian

    :param terms: List of nearst-neighbour terms (square array or MPO,
        see return value of :func:`cXY_local_terms`)
    :param ldim: Local dimension

    :returns: The Hamiltonian as sparse matrix

    """
    H = 0
    N = len(terms)
    for pos, term in enumerate(terms[:-1]):
        if hasattr(term, 'lt'):
            # Convert MPO to regular matrix
            term = term.to_array_global().reshape((ldim**2, ldim**2))
        left = sp.eye(ldim**pos)
        right = sp.eye(ldim**(N - pos - 2))
        H += sp.kron(left, sp.kron(term, right))
    # The last term acts on the first and last site.
    cyc = terms[-1]
    middle = sp.eye(ldim**pos)
    for i in range(cyc.ranks[0]):
        H += sp.kron(cyc.lt[0][0, ..., i], sp.kron(middle, cyc.lt[1][i, ..., 0]))
    return H

def edgeCurl(self):
        """The edgeCurl property."""
        if self.nCy > 1:
            raise NotImplementedError('Edge curl not yet implemented for '
                                      'nCy > 1')
        if getattr(self, '_edgeCurl', None) is None:
            # 1D Difference matricies
            dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1, 0],
                            self.nCx, self.nCx, format="csr")
            dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0, 1],
                            self.nCz, self.nCz+1, format="csr")

            # 2D Difference matricies
            Dr = sp.kron(sp.identity(self.nNz), dr)
            Dz = -sp.kron(dz, sp.identity(self.nCx))

            A = self.area
            E = self.edge
            # Edge curl operator
            self._edgeCurl = (utils.sdiag(1/A)*sp.vstack((Dz, Dr)) *
                              utils.sdiag(E))
        return self._edgeCurl

def aveE2CC(self):
        "Construct the averaging operator on cell edges to cell centers."
        if getattr(self, '_aveE2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if self.isSymmetric:
                avR = utils.av(n[0])[:, 1:]
                avR[0, 0] = 1.
                self._aveE2CC = sp.kron(utils.av(n[2]), avR, format="csr")
            else:
                raise NotImplementedError('wrapping in the averaging is not '
                                          'yet implemented')
                # self._aveE2CC = (1./3)*sp.hstack((utils.kron3(utils.av(n[2]),
                #                                               utils.av(n[1]),
                #                                               utils.speye(n[0])),
                #                                   utils.kron3(utils.av(n[2]),
                #                                               utils.speye(n[1]),
                #                                               utils.av(n[0])),
                #                                   utils.kron3(utils.speye(n[2]),
                #                                               utils.av(n[1]),
                #                                               utils.av(n[0]))),
                #                                  format="csr")
        return self._aveE2CC

def faceDiv(self):
        """
        Construct divergence operator (face-stg to cell-centres).
        """
        if getattr(self, '_faceDiv', None) is None:
            n = self.vnC
            # Compute faceDivergence operator on faces
            if(self.dim == 1):
                D = ddx(n[0])
            elif(self.dim == 2):
                D1 = sp.kron(speye(n[1]), ddx(n[0]))
                D2 = sp.kron(ddx(n[1]), speye(n[0]))
                D = sp.hstack((D1, D2), format="csr")
            elif(self.dim == 3):
                D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
                D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
                D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
                D = sp.hstack((D1, D2, D3), format="csr")
            # Compute areas of cell faces & volumes
            S = self.area
            V = self.vol
            self._faceDiv = sdiag(1/V)*D*sdiag(S)
        return self._faceDiv

def faceDivx(self):
        """
        Construct divergence operator in the x component (face-stg to
        cell-centres).
        """
        if getattr(self, '_faceDivx', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            # Compute faceDivergence operator on faces
            if(self.dim == 1):
                D1 = ddx(n[0])
            elif(self.dim == 2):
                D1 = sp.kron(speye(n[1]), ddx(n[0]))
            elif(self.dim == 3):
                D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
            # Compute areas of cell faces & volumes
            S = self.r(self.area, 'F', 'Fx', 'V')
            V = self.vol
            self._faceDivx = sdiag(1/V)*D1*sdiag(S)

        return self._faceDivx

def faceDivy(self):
        if(self.dim < 2):
            return None
        if getattr(self, '_faceDivy', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            # Compute faceDivergence operator on faces
            if(self.dim == 2):
                D2 = sp.kron(ddx(n[1]), speye(n[0]))
            elif(self.dim == 3):
                D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
            # Compute areas of cell faces & volumes
            S = self.r(self.area, 'F', 'Fy', 'V')
            V = self.vol
            self._faceDivy = sdiag(1/V)*D2*sdiag(S)
        return self._faceDivy

def nodalGrad(self):
        """
        Construct gradient operator (nodes to edges).
        """
        if getattr(self, '_nodalGrad', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            # Compute divergence operator on faces
            if(self.dim == 1):
                G = ddx(n[0])
            elif(self.dim == 2):
                D1 = sp.kron(speye(n[1]+1), ddx(n[0]))
                D2 = sp.kron(ddx(n[1]), speye(n[0]+1))
                G = sp.vstack((D1, D2), format="csr")
            elif(self.dim == 3):
                D1 = kron3(speye(n[2]+1), speye(n[1]+1), ddx(n[0]))
                D2 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]+1))
                D3 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]+1))
                G = sp.vstack((D1, D2, D3), format="csr")
            # Compute lengths of cell edges
            L = self.edge
            self._nodalGrad = sdiag(1/L)*G
        return self._nodalGrad

def aveCC2F(self):
        "Construct the averaging operator on cell cell centers to faces."
        if getattr(self, '_aveCC2F', None) is None:
            if self.dim == 1:
                self._aveCC2F = av_extrap(self.nCx)
            elif self.dim == 2:
                self._aveCC2F = sp.vstack((
                    sp.kron(speye(self.nCy), av_extrap(self.nCx)),
                    sp.kron(av_extrap(self.nCy), speye(self.nCx))
                ), format="csr")
            elif self.dim == 3:
                self._aveCC2F = sp.vstack((
                    kron3(
                        speye(self.nCz), speye(self.nCy), av_extrap(self.nCx)
                    ),
                    kron3(
                        speye(self.nCz), av_extrap(self.nCy), speye(self.nCx)
                    ),
                    kron3(
                        av_extrap(self.nCz), speye(self.nCy), speye(self.nCx)
                    )
                ), format="csr")
        return self._aveCC2F

def aveN2E(self):
        """
        Construct the averaging operator on cell nodes to cell edges, keeping
        each dimension separate.
        """

        if getattr(self, '_aveN2E', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveN2E = av(n[0])
            elif(self.dim == 2):
                self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
                                          sp.kron(av(n[1]), speye(n[0]+1))),
                                         format="csr")
            elif(self.dim == 3):
                self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1),
                                                av(n[0])),
                                          kron3(speye(n[2]+1), av(n[1]),
                                                speye(n[0]+1)),
                                          kron3(av(n[2]), speye(n[1]+1),
                                                speye(n[0]+1))),
                                         format="csr")
        return self._aveN2E

def _calculate_jacobian_mask(nx, ny, dof):
    from scipy.sparse import diags, kron

    N = nx * ny

    j_structure = diags(
        np.ones(shape=(7,)),
        [-nx + 1, -nx, -1, 0, 1, +nx, +nx - 1],
        shape=(N, N),
        format='csr',
    )

    j_structure = kron(
        j_structure,
        np.ones(shape=(dof, dof)),
        format='csr',
    )

    return j_structure

def area(self):
        """Face areas"""
        if getattr(self, '_area', None) is None:
            if self.nCy > 1:
                raise NotImplementedError('area not yet implemented for 3D '
                                          'cyl mesh')
            areaR = np.kron(self.hz, 2*pi*self.vectorNx)
            areaZ = np.kron(np.ones_like(self.vectorNz), pi*(self.vectorNx**2 -
                            np.r_[0, self.vectorNx[:-1]]**2))
            self._area = np.r_[areaR, areaZ]
        return self._area

def vol(self):
        """Volume of each cell"""
        if getattr(self, '_vol', None) is None:
            if self.nCy > 1:
                raise NotImplementedError('vol not yet implemented for 3D '
                                          'cyl mesh')
            az = pi*(self.vectorNx**2 - np.r_[0, self.vectorNx[:-1]]**2)
            self._vol = np.kron(self.hz, az)
        return self._vol

    ####################################################
    # Operators
    ####################################################

def aveF2CC(self):
        "Construct the averaging operator on cell faces to cell centers."
        if getattr(self, '_aveF2CC', None) is None:
            n = self.vnC
            if self.isSymmetric:
                avR = utils.av(n[0])[:, 1:]
                avR[0, 0] = 1.
                self._aveF2CC = ((0.5)*sp.hstack((sp.kron(utils.speye(n[2]),
                                                          avR),
                                                  sp.kron(utils.av(n[2]),
                                                          utils.speye(n[0]))),
                                                 format="csr"))
            else:
                raise NotImplementedError('wrapping in the averaging is not '
                                          'yet implemented')
                # self._aveF2CC = (1./3.)*sp.hstack((utils.kron3(utils.speye(n[2]),
                #                                                utils.speye(n[1]),
                #                                                utils.av(n[0])),
                #                                    utils.kron3(utils.speye(n[2]),
                #                                                utils.av(n[1]),
                #                                                utils.speye(n[0])),
                #                                    utils.kron3(utils.av(n[2]),
                #                                                utils.speye(n[1]),
                #                                                utils.speye(n[0]))),
                #                                   format="csr")
        return self._aveF2CC

def nodalLaplacian(self):
        """
        Construct laplacian operator (nodes to edges).
        """
        if getattr(self, '_nodalLaplacian', None) is None:
            print('Warning: Laplacian has not been tested rigorously.')
            # The number of cell centers in each direction
            n = self.vnC
            # Compute divergence operator on faces
            if(self.dim == 1):
                D1 = sdiag(1./self.hx) * ddx(self.nCx)
                L = - D1.T*D1
            elif(self.dim == 2):
                D1 = sdiag(1./self.hx) * ddx(n[0])
                D2 = sdiag(1./self.hy) * ddx(n[1])
                L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
                L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
                L = L1 + L2
            elif(self.dim == 3):
                D1 = sdiag(1./self.hx) * ddx(n[0])
                D2 = sdiag(1./self.hy) * ddx(n[1])
                D3 = sdiag(1./self.hz) * ddx(n[2])
                L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
                L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
                L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
                L = L1 + L2 + L3
            self._nodalLaplacian = L
        return self._nodalLaplacian

def _cellGradxStencil(self):
        BC = ['neumann', 'neumann']
        n = self.vnC
        if(self.dim == 1):
            G1 = ddxCellGrad(n[0], BC)
        elif(self.dim == 2):
            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
        elif(self.dim == 3):
            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
        return G1

def _cellGradyStencil(self):
        if self.dim < 2: return None
        BC = ['neumann', 'neumann']
        n = self.vnC
        if(self.dim == 2):
            G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
        elif(self.dim == 3):
            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
        return G2

def edgeCurl(self):
        """
        Construct the 3D curl operator.
        """
        if getattr(self, '_edgeCurl', None) is None:
            assert self.dim > 1, "Edge Curl only programed for 2 or 3D."

            n = self.vnC  # The number of cell centers in each direction
            L = self.edge  # Compute lengths of cell edges
            S = self.area # Compute areas of cell faces

            # Compute divergence operator on faces
            if self.dim == 2:

                D21 = sp.kron(ddx(n[1]), speye(n[0]))
                D12 = sp.kron(speye(n[1]), ddx(n[0]))
                C = sp.hstack((-D21, D12), format="csr")
                self._edgeCurl = C*sdiag(1/S)

            elif self.dim == 3:

                D32 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]+1))
                D23 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]+1))
                D31 = kron3(ddx(n[2]), speye(n[1]+1), speye(n[0]))
                D13 = kron3(speye(n[2]), speye(n[1]+1), ddx(n[0]))
                D21 = kron3(speye(n[2]+1), ddx(n[1]), speye(n[0]))
                D12 = kron3(speye(n[2]+1), speye(n[1]), ddx(n[0]))

                O1 = spzeros(np.shape(D32)[0], np.shape(D31)[1])
                O2 = spzeros(np.shape(D31)[0], np.shape(D32)[1])
                O3 = spzeros(np.shape(D21)[0], np.shape(D13)[1])

                C = sp.vstack((sp.hstack((O1, -D32, D23)),
                               sp.hstack((D31, O2, -D13)),
                               sp.hstack((-D21, D12, O3))), format="csr")

                self._edgeCurl = sdiag(1/S)*(C*sdiag(L))
        return self._edgeCurl

def aveFx2CC(self):
        """
        Construct the averaging operator on cell faces in the x direction to
        cell centers.
        """

        if getattr(self, '_aveFx2CC', None) is None:
            n = self.vnC
            if self.dim == 1:
                self._aveFx2CC = av(n[0])
            elif self.dim == 2:
                self._aveFx2CC = sp.kron(speye(n[1]), av(n[0]))
            elif self.dim == 3:
                self._aveFx2CC = kron3(speye(n[2]), speye(n[1]), av(n[0]))
        return self._aveFx2CC

def aveEx2CC(self):
        """
        Construct the averaging operator on cell edges in the x direction to
        cell centers.
        """
        if getattr(self, '_aveEx2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveEx2CC = speye(n[0])
            elif(self.dim == 2):
                self._aveEx2CC = sp.kron(av(n[1]), speye(n[0]))
            elif(self.dim == 3):
                self._aveEx2CC = kron3(av(n[2]), av(n[1]), speye(n[0]))
        return self._aveEx2CC

def aveEy2CC(self):
        """
        Construct the averaging operator on cell edges in the y direction to
        cell centers.
        """
        if self.dim < 2:
            return None
        if getattr(self, '_aveEy2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 2):
                self._aveEy2CC = sp.kron(speye(n[1]), av(n[0]))
            elif(self.dim == 3):
                self._aveEy2CC = kron3(av(n[2]), speye(n[1]), av(n[0]))
        return self._aveEy2CC

def aveN2CC(self):
        "Construct the averaging operator on cell nodes to cell centers."
        if getattr(self, '_aveN2CC', None) is None:
            # The number of cell centers in each direction
            n = self.vnC
            if(self.dim == 1):
                self._aveN2CC = av(n[0])
            elif(self.dim == 2):
                self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
            elif(self.dim == 3):
                self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
        return self._aveN2CC

def kron3(A, B, C):
    """Three kron prods"""
    return sp.kron(sp.kron(A, B), C, format="csr")

def makePropertyTensor(M, tensor):
    if tensor is None:  # default is ones
        tensor = np.ones(M.nC)

    if isScalar(tensor):
        tensor = tensor * np.ones(M.nC)

    propType = TensorType(M, tensor)
    if propType == 1:  # Isotropic!
        Sigma = sp.kron(sp.identity(M.dim), sdiag(mkvc(tensor)))
    elif propType == 2:  # Diagonal tensor
        Sigma = sdiag(mkvc(tensor))
    elif M.dim == 2 and tensor.size == M.nC*3:  # Fully anisotropic, 2D
        tensor = tensor.reshape((M.nC, 3), order='F')
        row1 = sp.hstack((sdiag(tensor[:, 0]), sdiag(tensor[:, 2])))
        row2 = sp.hstack((sdiag(tensor[:, 2]), sdiag(tensor[:, 1])))
        Sigma = sp.vstack((row1, row2))
    elif M.dim == 3 and tensor.size == M.nC*6:  # Fully anisotropic, 3D
        tensor = tensor.reshape((M.nC, 6), order='F')
        row1 = sp.hstack(
            (sdiag(tensor[:, 0]), sdiag(tensor[:, 3]), sdiag(tensor[:, 4]))
        )
        row2 = sp.hstack(
            (sdiag(tensor[:, 3]), sdiag(tensor[:, 1]), sdiag(tensor[:, 5]))
        )
        row3 = sp.hstack(
            (sdiag(tensor[:, 4]), sdiag(tensor[:, 5]), sdiag(tensor[:, 2]))
        )
        Sigma = sp.vstack((row1, row2, row3))
    else:
        raise Exception('Unexpected shape of tensor')

    return Sigma

def _getInnerProductProjectionMatrices(self, projType, tensorType):
        """
        :param str projType: 'F' for faces 'E' for edges
        :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma)
        """
        assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.'
        assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"

        d = self.dim
        # We will multiply by sqrt on each side to keep symmetry
        V = sp.kron(sp.identity(d), sdiag(np.sqrt((2**(-d))*self.vol)))

        nodes = ['000', '100', '010', '110', '001', '101', '011',  '111'][:2**d]

        if projType == 'F':
            locs = {
                    '000': [('fXm',), ('fXm', 'fYm'), ('fXm', 'fYm', 'fZm')],
                    '100': [('fXp',), ('fXp', 'fYm'), ('fXp', 'fYm', 'fZm')],
                    '010': [  None  , ('fXm', 'fYp'), ('fXm', 'fYp', 'fZm')],
                    '110': [  None  , ('fXp', 'fYp'), ('fXp', 'fYp', 'fZm')],
                    '001': [  None  ,      None     , ('fXm', 'fYm', 'fZp')],
                    '101': [  None  ,      None     , ('fXp', 'fYm', 'fZp')],
                    '011': [  None  ,      None     , ('fXm', 'fYp', 'fZp')],
                    '111': [  None  ,      None     , ('fXp', 'fYp', 'fZp')]
                   }
            proj = getattr(self, '_getFaceP' + ('x'*d))()

        elif projType == 'E':
            locs = {
                    '000': [('eX0',), ('eX0', 'eY0'), ('eX0', 'eY0', 'eZ0')],
                    '100': [('eX0',), ('eX0', 'eY1'), ('eX0', 'eY1', 'eZ1')],
                    '010': [  None  , ('eX1', 'eY0'), ('eX1', 'eY0', 'eZ2')],
                    '110': [  None  , ('eX1', 'eY1'), ('eX1', 'eY1', 'eZ3')],
                    '001': [  None  ,      None     , ('eX2', 'eY2', 'eZ0')],
                    '101': [  None  ,      None     , ('eX2', 'eY3', 'eZ1')],
                    '011': [  None  ,      None     , ('eX3', 'eY2', 'eZ2')],
                    '111': [  None  ,      None     , ('eX3', 'eY3', 'eZ3')]
                   }
            proj = getattr(self, '_getEdgeP' + ('x'*d))()

        return [V*proj(*locs[node][d-1]) for node in nodes]

def visit_Kron(self, node):
        """ Kron(I, SpMatrix) => SpMatrix """
        node = self.generic_visit(node)
        L, R = node.children
        if isinstance(L, Eye):
            L = L.realize()
        if isinstance(L, SpMatrix) and isinstance(R, SpMatrix):
            name = "({}(x){})".format(L._name, R._name)
            log.debug('realizing kron %s x %s', L._name, R._name)
            K = spp.kron(L._matrix, R._matrix)
            return SpMatrix( node._backend, K, name=name )
        else:
            return node

def test_KronI(backend, M, N, K, density, alpha, beta):
    b = backend()
    mat_h = indigo.util.randM(M,N,density)
    A_h = spp.kron( spp.eye(K), mat_h )

    mat_d = b.SpMatrix(mat_h)
    A = b.KronI(K, mat_d)

    # forward
    x = b.rand_array((A.shape[1],K))
    y = b.rand_array((A.shape[0],K))
    y_exp = beta * y.to_host() + alpha * A_h @ x.to_host()
    A.eval(y, x, alpha=alpha, beta=beta)
    npt.assert_allclose(y.to_host(), y_exp, rtol=1e-5)

    # adjoint
    x = b.rand_array((A.shape[0],K))
    y = b.rand_array((A.shape[1],K))
    y_exp = beta * y.to_host() + alpha * np.conj(A_h.T) @ x.to_host()
    A.H.eval(y, x, alpha=alpha, beta=beta)
    npt.assert_allclose(y.to_host(), y_exp, rtol=1e-5)

    # shape
    assert A.shape == (M*K,N*K)
    assert A.H.shape == (N*K,M*K)

    # dtype
    assert A.dtype == np.dtype('complex64')

def test_nested_kroni(backend, batch, M, N, c1, c2):
    A00_h = indigo.util.randM(M, N, 0.9)
    A0_h = spp.kron( spp.eye(c1), A00_h )
    A_h = spp.kron( spp.eye(c2), A0_h )

    b = backend()
    A00 = b.SpMatrix(A00_h)
    A0 = b.KronI(c1, A00)
    A = b.KronI(c2, A0)

    V, U = A.shape
    u = np.random.rand(U) + 1j*np.random.rand(U)
    u = np.require(u, dtype=np.dtype('complex64'), requirements='F')
    v = np.random.rand(V) + 1j*np.random.rand(V)
    v = np.require(v, dtype=np.dtype('complex64'), requirements='F')

    # forward
    u_d = b.copy_array(u)
    v_d = b.copy_array(v)
    A.eval(v_d, u_d)
    v_act = v_d.to_host()
    v_exp = A_h @ u
    np.testing.assert_allclose( v_act, v_exp, rtol=1e-6 )

    # adjoint
    u_d = b.copy_array(u)
    v_d = b.copy_array(v)
    A.H.eval(u_d, v_d)
    u_act = u_d.to_host()
    u_exp = A_h.H @ v
    np.testing.assert_allclose( u_act, u_exp, rtol=1e-6 )

def reorder(array,axes=None,permutation=None):
    '''
    Reorder the axes of an array from the ordinary numpy.kron order to the correct quantum number collection order.

    Parameters
    ----------
    array : ndarray-like
        The original array in the ordinary numpy.kron order.
    axes : list of integer, optional
        The axes of the array to be reordered.
    permutation : 1d ndarray of integers, optional
        The permutation array applied to the required axes.

    Returns
    -------
    ndarray-like
        The axes-reordered array.
    '''
    result=array
    if permutation is not None:
        axes=xrange(array.ndim) if axes is None else axes
        for axis in axes:
            temp=[slice(None,None,None)]*array.ndim
            temp[axis]=permutation
            result=result[tuple(temp)]
    return result

def cellGradBC(self):
        """
        The cell centered Gradient boundary condition matrix
        """
        if getattr(self, '_cellGradBC', None) is None:
            BC = self.setCellGradBC(self._cellGradBC_list)
            n = self.vnC
            if(self.dim == 1):
                G = ddxCellGradBC(n[0], BC[0])
            elif(self.dim == 2):
                G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G = sp.block_diag((G1, G2), format="csr")
            elif(self.dim == 3):
                G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
                G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
                G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
                G = sp.block_diag((G1, G2, G3), format="csr")
            # Compute areas of cell faces & volumes
            S = self.area
            V = self.aveCC2F*self.vol  # Average volume between adjacent cells
            self._cellGradBC = sdiag(S/V)*G
        return self._cellGradBC

    # def cellGradBC():
    #     doc = "The cell centered Gradient boundary condition matrix"

    #     def fget(self):
    #         if(self._cellGradBC is None):
    #             BC = self.setCellGradBC(self._cellGradBC_list)
    #             n = self.vnC
    #             if(self.dim == 1):
    #                 G = ddxCellGradBC(n[0], BC[0])
    #             elif(self.dim == 2):
    #                 G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
    #                 G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
    #                 G = sp.block_diag((G1, G2), format="csr")
    #             elif(self.dim == 3):
    #                 G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
    #                 G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
    #                 G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
    #                 G = sp.block_diag((G1, G2, G3), format="csr")
    #             # Compute areas of cell faces & volumes
    #             S = self.area
    #             V = self.aveCC2F*self.vol  # Average volume between adjacent cells
    #             self._cellGradBC = sdiag(S/V)*G
    #         return self._cellGradBC
    #     return locals()
    # _cellGradBC = None
    # cellGradBC = property(**cellGradBC())

def random_model(shape, seed=None, anisotropy=None, its=100, bounds=None):
    """
        Create a random model by convolving a kernel with a
        uniformly distributed model.

        :param tuple shape: shape of the model.
        :param int seed: pick which model to produce, prints the seed if you don't choose.
        :param numpy.ndarray anisotropy: this is the (3 x n) blurring kernel that is used.
        :param int its: number of smoothing iterations
        :param list bounds: bounds on the model, len(list) == 2
        :rtype: numpy.ndarray
        :return: M, the model


        .. plot::

            import matplotlib.pyplot as plt
            import discretize
            plt.colorbar(plt.imshow(discretize.utils.random_model((50, 50), bounds=[-4, 0])))
            plt.title('A very cool, yet completely random model.')
            plt.show()


    """
    if bounds is None:
        bounds = [0, 1]

    if seed is None:
        seed = np.random.randint(1e3)
        print('Using a seed of: ', seed)

    if type(shape) in num_types:
        shape = (shape, ) # make it a tuple for consistency

    np.random.seed(seed)
    mr = np.random.rand(*shape)
    if anisotropy is None:
        if len(shape) is 1:
            smth = np.array([1, 10., 1], dtype=float)
        elif len(shape) is 2:
            smth = np.array([[1, 7, 1], [2, 10, 2], [1, 7, 1]], dtype=float)
        elif len(shape) is 3:
            kernal = np.array([1, 4, 1], dtype=float).reshape((1, 3))
            smth = np.array(sp.kron(sp.kron(kernal, kernal.T).todense()[:], kernal).todense()).reshape((3, 3, 3))
    else:
        assert len(anisotropy.shape) is len(shape), 'Anisotropy must be the same shape.'
        smth = np.array(anisotropy, dtype=float)

    smth = smth/smth.sum() # normalize
    mi = mr
    for i in range(its):
        mi = ndi.convolve(mi, smth)

    # scale the model to live between the bounds.
    mi = (mi - mi.min())/(mi.max()-mi.min()) # scaled between 0 and 1
    mi = mi*(bounds[1]-bounds[0])+bounds[0]

    return mi

def test_Kron_general(backend, L, Q, alpha, beta, eyeL, eyeR):
    b = backend()

    if eyeL:
        A = np.eye(L, dtype=np.complex64)
        A_d = b.Eye(L)
        if not eyeR: A_d._eval = lambda unreachable: ()
    else:
        A = indigo.util.rand64c(L,L)
        A_d = b.DenseMatrix(A)

    if eyeR:
        B = np.eye(Q, dtype=np.complex64)
        B_d = b.Eye(Q)
        if not eyeL: B_d._eval = lambda unreachable: ()
    else:
        B = indigo.util.rand64c(Q,Q)
        B = B + B.T
        B.imag = 0
        B_d = b.DenseMatrix(B)

    K = b.Kron( B_d, A_d )
    vec = lambda arr: arr.reshape((-1,1), order='F')

    # check
    X = indigo.util.rand64c(L,Q)
    Y = indigo.util.rand64c(L,Q)
    y_exp  = alpha * np.kron(B.T, A) @ vec(X) + beta * vec(Y)
    y_exp2 = alpha * vec(A @ X @ B) + beta * vec(Y)
    npt.assert_allclose(y_exp2, y_exp, rtol=1e-5)

    # forward
    x_d = b.copy_array(X)
    y_d = b.copy_array(Y)
    K.eval(y_d, x_d, alpha=alpha, beta=beta)
    y_act = y_d.to_host()
    npt.assert_allclose(vec(y_act), y_exp, rtol=1e-5)

    # adjoint
    X = indigo.util.rand64c(Q,L)
    Y = indigo.util.rand64c(Q,L)
    y_exp = alpha * np.conj(np.kron(B.T, A).T) @ vec(X) + beta * vec(Y)
    x_d = b.copy_array(X)
    y_d = b.copy_array(Y)
    K.eval(y_d, x_d, alpha=alpha, beta=beta, forward=False)
    y_act = y_d.to_host()
    npt.assert_allclose(vec(y_act), y_exp, rtol=1e-5)

def kron(m1,m2,rcs=None,timers=None):
    '''
    Kronecker product of two matrices.

    Parameters
    ----------
    m1,m2 : 2d ndarray
        The matrices.
    rcs : 1d ndarray or 3-tuple of 1d ndarray
        * When 1d ndarray
            The selected rows and columns of the kronecker product
        * When 3-tuple of 1d ndarray
            * tuple[0]: the selected rows and columns of the first matrix `m1`
            * tuple[1]: the selected rows and columns of the second matrix `m2`
            * tuple[2]: the map between the indices before and after the selection of the rows and columns of the kronecker product.
    timers : Timers, optional
        The timers to record certain procedures of this function.

    Returns
    -------
    csr_matrix
        The product.
    '''
    if rcs is None:
        result=sp.kron(m1,m2,format='csr')
    else:
        assert m1.dtype==m2.dtype and m1.shape[0]==m1.shape[1] and m2.shape[0]==m2.shape[1]
        if isinstance(rcs,np.ndarray):
            rcs1,rcs2=np.divide(rcs,m2.shape[1]),np.mod(rcs,m2.shape[1])
            slices=np.zeros(m1.shape[1]*m2.shape[1],dtype=np.int64)
            slices[rcs]=xrange(len(rcs))
        else:
            rcs1,rcs2,slices=rcs
        def csr(m1,m2):
            return sp.csr_matrix(m1),sp.csr_matrix(m2)
        def fkron(m1,m2):
            nnz=(m1.indptr[rcs1+1]-m1.indptr[rcs1]).dot(m2.indptr[rcs2+1]-m2.indptr[rcs2])
            if nnz>0:
                if m1.dtype==np.float32:
                    data,indices,indptr=fkron_r4(m1.data,m1.indices,m1.indptr,rcs1,m2.data,m2.indices,m2.indptr,rcs2,nnz,slices)
                elif m1.dtype==np.float64:
                    data,indices,indptr=fkron_r8(m1.data,m1.indices,m1.indptr,rcs1,m2.data,m2.indices,m2.indptr,rcs2,nnz,slices)
                elif m1.dtype==np.complex64:
                    data,indices,indptr=fkron_c4(m1.data,m1.indices,m1.indptr,rcs1,m2.data,m2.indices,m2.indptr,rcs2,nnz,slices)
                elif m1.dtype==np.complex128:
                    data,indices,indptr=fkron_c8(m1.data,m1.indices,m1.indptr,rcs1,m2.data,m2.indices,m2.indptr,rcs2,nnz,slices)
                else:
                    raise ValueError("_fkron_ error: only matrices with dtype being float32, float64, complex64 or complex128 are supported.")
                result=sp.csr_matrix((data,indices,indptr),shape=(len(rcs1),len(rcs1)))
            else:
                result=sp.csr_matrix((len(rcs1),len(rcs1)),dtype=m1.dtype)
            return result
        if timers is None:
            m1,m2=csr(m1,m2)
            result=fkron(m1,m2)
        else:
            with timers.get('csr'):
                m1,m2=csr(m1,m2)
            with timers.get('fkron'):
                result=fkron(m1,m2)
    return result