Python linalg.splu方法代码示例

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_solve(self):
        # Test whether the lu_solve command segfaults, as reported by Nils
        # Wagner for a 64-bit machine, 02 March 2005 (EJS)
        n = 20
        np.random.seed(0)  # make tests repeatable
        A = zeros((n,n), dtype=complex)
        x = np.random.rand(n)
        y = np.random.rand(n-1)+1j*np.random.rand(n-1)
        r = np.random.rand(n)
        for i in range(len(x)):
            A[i,i] = x[i]
        for i in range(len(y)):
            A[i,i+1] = y[i]
            A[i+1,i] = conjugate(y[i])
        A = self.spmatrix(A)
        x = splu(A).solve(r)
        assert_almost_equal(A*x,r) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_expm(self):
        M = array([[1, 0, 2], [0, 0, 3], [-4, 5, 6]], float)
        sM = self.spmatrix(M, shape=(3,3), dtype=float)
        Mexp = scipy.linalg.expm(M)

        N = array([[3., 0., 1.], [0., 2., 0.], [0., 0., 0.]])
        sN = self.spmatrix(N, shape=(3,3), dtype=float)
        Nexp = scipy.linalg.expm(N)

        with suppress_warnings() as sup:
            sup.filter(SparseEfficiencyWarning, "splu requires CSC matrix format")
            sup.filter(SparseEfficiencyWarning,
                       "spsolve is more efficient when sparse b is in the CSC matrix format")
            sup.filter(SparseEfficiencyWarning,
                       "spsolve requires A be CSC or CSR matrix format")
            sMexp = expm(sM).todense()
            sNexp = expm(sN).todense()

        assert_array_almost_equal((sMexp - Mexp), zeros((3, 3)))
        assert_array_almost_equal((sNexp - Nexp), zeros((3, 3))) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_inv(self):
        def check(dtype):
            M = array([[1, 0, 2], [0, 0, 3], [-4, 5, 6]], dtype)
            with suppress_warnings() as sup:
                sup.filter(SparseEfficiencyWarning,
                           "spsolve requires A be CSC or CSR matrix format")
                sup.filter(SparseEfficiencyWarning,
                           "spsolve is more efficient when sparse b is in the CSC matrix format")
                sup.filter(SparseEfficiencyWarning,
                           "splu requires CSC matrix format")
                sM = self.spmatrix(M, shape=(3,3), dtype=dtype)
                sMinv = inv(sM)
            assert_array_almost_equal(sMinv.dot(sM).todense(), np.eye(3))
            assert_raises(TypeError, inv, M)
        for dtype in [float]:
            check(dtype) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_solve(self):
        # Test whether the lu_solve command segfaults, as reported by Nils
        # Wagner for a 64-bit machine, 02 March 2005 (EJS)
        n = 20
        np.random.seed(0)  # make tests repeatable
        A = zeros((n,n), dtype=complex)
        x = np.random.rand(n)
        y = np.random.rand(n-1)+1j*np.random.rand(n-1)
        r = np.random.rand(n)
        for i in range(len(x)):
            A[i,i] = x[i]
        for i in range(len(y)):
            A[i,i+1] = y[i]
            A[i+1,i] = conjugate(y[i])
        A = self.spmatrix(A)
        with suppress_warnings() as sup:
            sup.filter(SparseEfficiencyWarning, "splu requires CSC matrix format")
            x = splu(A).solve(r)
        assert_almost_equal(A*x,r) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I+factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float) : abbrev. for the node-to-node stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        L = splu(sp.eye(self.params.nvars, format='csc') + factor * self.A)
        me.values = L.solve(rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I+factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float) : abbrev. for the node-to-node stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver (not used here so far)
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        L = splu(sp.eye(self.params.nvars, format='csc') - factor * self.A)
        me.values = L.solve(rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-factor*A)u = rhs

        Args:
            rhs (dtype_f): right-hand side for the linear system
            factor (float): abbrev. for the local stepsize (or any other factor required)
            u0 (dtype_u): initial guess for the iterative solver
            t (float): current time (e.g. for time-dependent BCs)

        Returns:
            dtype_u: solution as mesh
        """

        me = self.dtype_u(self.init)
        L = splu(sp.eye(self.params.nvars, format='csc') - factor * self.A)
        me.values = L.solve(rhs.values)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def __init__(self, work):
        """
        Initialize structure for KKT system solution
        """
        # Construct reduced KKT matrix
        KKT = spspa.vstack([
              spspa.hstack([work.data.P + work.settings.sigma *
                            spspa.eye(work.data.n), work.data.A.T]),
              spspa.hstack([work.data.A, -spspa.diags(work.rho_inv_vec)])])

        # Initialize structure
        self.kkt_factor = spla.splu(KKT.tocsc())
        #  self.lu, self.piv = sp.linalg.lu_factor(KKT.todense()) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_preconditioner(self):
        # Check that preconditioning works
        pc = splu(Am.tocsc())
        M = LinearOperator(matvec=pc.solve, shape=A.shape, dtype=A.dtype)

        x0, count_0 = do_solve()
        x1, count_1 = do_solve(M=M)

        assert_(count_1 == 3)
        assert_(count_1 < count_0/2)
        assert_(allclose(x1, x0, rtol=1e-14)) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def solve_nonlinear(self, inputs, outputs):
        """
        Use numpy to solve Ax=b for x.

        Parameters
        ----------
        inputs : Vector
            unscaled, dimensional input variables read via inputs[key]
        outputs : Vector
            unscaled, dimensional output variables read via outputs[key]
        """
        # lu factorization for use with solve_linear
        K = self.assemble_CSC_K(inputs)
        self._lup = splu(K)
        outputs['disp_aug'] = self._lup.solve(inputs['forces']) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def sparse_solve_kkt_inverse(H_, A_, C_tilde, rx, rs, rz, ry, ns):
    nineq, nz, neq, nBatch = ns

    if neq > 0:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = torch.cat([rz, ry], 1).squeeze(0).numpy()
    else:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = rz.squeeze(0).numpy()

    full_mat = bmat([[H_, A_.transpose()],
                     [A_, C_tilde]], format='csc')
    full_res = np.concatenate([g_, h_], 0)
    sol = splu(full_mat).solve(full_res)
    # sol = spsolve(full_mat, full_res)

    dx = sol[:nz]
    ds = sol[nz:nz+nineq]
    dz = sol[nz+nineq:nz+nineq+nineq]
    dy = sol[nz+nineq+nineq:] if neq > 0 else None

    dx = torch.DoubleTensor(dx).unsqueeze(0)
    ds = torch.DoubleTensor(ds).unsqueeze(0)
    dz = torch.DoubleTensor(dz).unsqueeze(0)
    dy = torch.DoubleTensor(dy).unsqueeze(0) if neq > 0 else None

    return dx, ds, dz, dy 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def sparse_factor_solve_kkt(Q_tilde, D_tilde, A_, C_tilde, rx, rs, rz, ry, ns):
    nineq, nz, neq, nBatch = ns

    # H_ = csc_matrix((nz + nineq, nz + nineq))
    # H_[:nz, :nz] = Q_tilde
    # H_[-nineq:, -nineq:] = D_tilde
    H_ = block_diag([Q_tilde, D_tilde], format='csc')
    if neq > 0:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = torch.cat([rz, ry], 1).squeeze(0).numpy()
    else:
        g_ = torch.cat([rx, rs], 1).squeeze(0).numpy()
        h_ = rz.squeeze(0).numpy()

    H_LU = splu(H_)

    invH_A_ = csc_matrix(H_LU.solve(A_.todense().transpose()))
    invH_g_ = H_LU.solve(g_)

    S_ = A_.dot(invH_A_) + C_tilde
    S_LU = splu(S_)
    # t_ = invH_g_[np.newaxis].dot(A_.transpose()).squeeze(0) - h_
    t_ = A_.dot(invH_g_) - h_
    w_ = -S_LU.solve(t_)
    # t_ = -g_ - w_[np.newaxis].dot(A_).squeeze(0)
    t_ = -g_ - A_.transpose().dot(w_)
    v_ = H_LU.solve(t_)

    dx = v_[:nz]
    ds = v_[nz:]
    dz = w_[:nineq]
    dy = w_[nineq:] if neq > 0 else None

    dx = torch.DoubleTensor(dx).unsqueeze(0)
    ds = torch.DoubleTensor(ds).unsqueeze(0)
    dz = torch.DoubleTensor(dz).unsqueeze(0)
    dy = torch.DoubleTensor(dy).unsqueeze(0) if neq > 0 else None

    return dx, ds, dz, dy 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def test_preconditioner(self):
        # Check that preconditioning works
        pc = splu(Am.tocsc())
        M = LinearOperator(matvec=pc.solve, shape=A.shape, dtype=A.dtype)

        x0, count_0 = do_solve()
        x1, count_1 = do_solve(M=M)

        assert_equal(count_1, 3)
        assert_(count_1 < count_0/2)
        assert_(allclose(x1, x0, rtol=1e-14)) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def laplaceTria(v, t, k=10):
    """
    Compute linear finite-element method Laplace-Beltrami spectrum
    """
    from scipy.sparse.linalg import LinearOperator, eigsh, splu
    useCholmod = True
    try:
        from sksparse.cholmod import cholesky
    except ImportError:
        useCholmod = False
    if useCholmod:
        print("Solver: cholesky decomp - performance optimal ...")
    else:
        print("Package scikit-sparse not found (Cholesky decomp)")
        print("Solver: spsolve (LU decomp) - performance not optimal ...")
    # import numpy as np
    # from shapeDNA import computeABtria
    A, M = computeABtria(v, t, lump=True)
    # turns out it is much faster to use cholesky and pass operator
    sigma = -0.01
    if useCholmod:
        chol = cholesky(A - sigma * M)
        OPinv = LinearOperator(matvec=chol, shape=A.shape, dtype=A.dtype)
    else:
        lu = splu(A - sigma * M)
        OPinv = LinearOperator(matvec=lu.solve, shape=A.shape, dtype=A.dtype)
    eigenvalues, eigenvectors = eigsh(A, k, M, sigma=sigma, OPinv=OPinv)
    return eigenvalues, eigenvectors 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import splu [as 别名]
def __init__(self, A):
    LOGGER.debug(
      'computing the LU decomposition of a %s by %s sparse matrix with %s '
      'nonzeros ' % (A.shape + (A.nnz,)))
    self.factor = spla.splu(A)