Python interface.LinearOperator方法代码示例

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def __init__(self, ops, dtype=None):
        self.ops = ops
        mops = np.zeros(len(ops), dtype=np.int)
        nops = np.zeros(len(ops), dtype=np.int)
        for iop, oper in enumerate(ops):
            if not isinstance(oper, (LinearOperator, spLinearOperator)):
                self.ops[iop] = MatrixMult(oper, dtype=oper.dtype)
            nops[iop] = self.ops[iop].shape[0]
            mops[iop] = self.ops[iop].shape[1]
        self.nops = nops.sum()
        self.mops = mops.sum()
        self.nnops = np.insert(np.cumsum(nops), 0, 0)
        self.mmops = np.insert(np.cumsum(mops), 0, 0)
        self.shape = (self.nops, self.mops)
        if dtype is None:
            self.dtype = _get_dtype(ops)
        else:
            self.dtype = np.dtype(dtype)
        self.explicit = False 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def __init__(self, ops, dtype=None):
        self.ops = ops
        nops = np.zeros(len(self.ops), dtype=np.int)
        for iop, oper in enumerate(ops):
            if not isinstance(oper, (LinearOperator, spLinearOperator)):
                self.ops[iop] = MatrixMult(oper, dtype=oper.dtype)
            nops[iop] = self.ops[iop].shape[0]
        self.nops = nops.sum()
        mops = [oper.shape[1] for oper in self.ops]
        if len(set(mops)) > 1:
            raise ValueError('operators have different number of columns')
        self.mops = mops[0]
        self.nnops = np.insert(np.cumsum(nops), 0, 0)
        self.shape = (self.nops, self.mops)
        if dtype is None:
            self.dtype = _get_dtype(self.ops)
        else:
            self.dtype = np.dtype(dtype)
        self.explicit = False 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def __init__(self, ops, dtype='float64'):
        self.ops = ops
        mops = np.zeros(len(ops), dtype=np.int)
        for iop, oper in enumerate(ops):
            if not isinstance(oper, (LinearOperator, spLinearOperator)):
                self.ops[iop] = MatrixMult(oper, dtype=oper.dtype)
            mops[iop] = self.ops[iop].shape[1]
        self.mops = mops.sum()
        nops = [oper.shape[0] for oper in self.ops]
        if len(set(nops)) > 1:
            raise ValueError('operators have different number of rows')
        self.nops = nops[0]
        self.mmops = np.insert(np.cumsum(mops), 0, 0)
        self.shape = (self.nops, self.mops)
        if dtype is None:
            self.dtype = _get_dtype(self.ops)
        else:
            self.dtype = np.dtype(dtype)
        self.explicit = False 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def _check_reentrancy(solver, is_reentrant):
    def matvec(x):
        A = np.array([[1.0, 0, 0], [0, 2.0, 0], [0, 0, 3.0]])
        y, info = solver(A, x)
        assert_equal(info, 0)
        return y
    b = np.array([1, 1./2, 1./3])
    op = LinearOperator((3, 3), matvec=matvec, rmatvec=matvec,
                        dtype=b.dtype)

    if not is_reentrant:
        assert_raises(RuntimeError, solver, op, b)
    else:
        y, info = solver(op, b)
        assert_equal(info, 0)
        assert_allclose(y, [1, 1, 1])


#------------------------------------------------------------------------------ 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def _check_reentrancy(solver, is_reentrant):
    def matvec(x):
        A = np.array([[1.0, 0, 0], [0, 2.0, 0], [0, 0, 3.0]])
        y, info = solver(A, x)
        assert_equal(info, 0)
        return y
    b = np.array([1, 1./2, 1./3])
    op = LinearOperator((3, 3), matvec=matvec, rmatvec=matvec,
                        dtype=b.dtype)

    if not is_reentrant:
        assert_raises(RuntimeError, solver, op, b)
    else:
        y, info = solver(op, b)
        assert_equal(info, 0)
        assert_allclose(y, [1, 1, 1]) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def test_matmul(self):
        if not TEST_MATMUL:
            pytest.skip("matmul is only tested in Python 3.5+")

        D = {'shape': self.A.shape,
             'matvec': lambda x: np.dot(self.A, x).reshape(self.A.shape[0]),
             'rmatvec': lambda x: np.dot(self.A.T.conj(),
                                         x).reshape(self.A.shape[1]),
             'matmat': lambda x: np.dot(self.A, x)}
        A = interface.LinearOperator(**D)
        B = np.array([[1, 2, 3],
                      [4, 5, 6],
                      [7, 8, 9]])
        b = B[0]

        assert_equal(operator.matmul(A, b), A * b)
        assert_equal(operator.matmul(A, B), A * B)
        assert_raises(ValueError, operator.matmul, A, 2)
        assert_raises(ValueError, operator.matmul, 2, A) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def compute_dr_wrt(self, wrt):
        if wrt is self.a:
            if False:
                from scipy.sparse.linalg.interface import LinearOperator
                return LinearOperator((self.size, wrt.size), lambda x : self.reorder(x.reshape(self.a.shape)).ravel())
            else:
                a = self.a
                asz = a.size
                ashape = a.shape
                key = self.unique_reorder_id()
                if key not in self.dr_lookup or key is None:
                    JS = self.reorder(np.arange(asz).reshape(ashape))
                    IS = np.arange(JS.size)
                    data = np.ones_like(IS)
                    shape = JS.shape
                    self.dr_lookup[key] = sp.csc_matrix((data, (IS, JS.ravel())), shape=(self.r.size, wrt.r.size))
                return self.dr_lookup[key] 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def __init__(self, A, M, sigma, ifunc=gmres, tol=0):
        if tol <= 0:
            # when tol=0, ARPACK uses machine tolerance as calculated
            # by LAPACK's _LAMCH function.  We should match this
            tol = 2 * np.finfo(A.dtype).eps
        self.A = A
        self.M = M
        self.sigma = sigma
        self.ifunc = ifunc
        self.tol = tol

        def mult_func(x):
            return A.matvec(x) - sigma * M.matvec(x)

        def mult_func_M_None(x):
            return A.matvec(x) - sigma * x

        x = np.zeros(A.shape[1])
        if M is None:
            dtype = mult_func_M_None(x).dtype
            self.OP = LinearOperator(self.A.shape,
                                     mult_func_M_None,
                                     dtype=dtype)
        else:
            dtype = mult_func(x).dtype
            self.OP = LinearOperator(self.A.shape,
                                     mult_func,
                                     dtype=dtype)
        self.shape = A.shape 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [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.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def check_precond_dummy(solver, case):
    tol = 1e-8

    def identity(b,which=None):
        """trivial preconditioner"""
        return b

    A = case.A

    M,N = A.shape
    D = spdiags([1.0/A.diagonal()], [0], M, N)

    b = arange(A.shape[0], dtype=float)
    x0 = 0*b

    precond = LinearOperator(A.shape, identity, rmatvec=identity)

    if solver is qmr:
        x, info = solver(A, b, M1=precond, M2=precond, x0=x0, tol=tol)
    else:
        x, info = solver(A, b, M=precond, x0=x0, tol=tol)
    assert_equal(info,0)
    assert_normclose(A.dot(x), b, tol)

    A = aslinearoperator(A)
    A.psolve = identity
    A.rpsolve = identity

    x, info = solver(A, b, x0=x0, tol=tol)
    assert_equal(info,0)
    assert_normclose(A*x, b, tol=tol) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def test_leftright_precond(self):
        """Check that QMR works with left and right preconditioners"""

        with warnings.catch_warnings():
            warnings.simplefilter("ignore", category=SparseEfficiencyWarning)
            from scipy.sparse.linalg.dsolve import splu
            from scipy.sparse.linalg.interface import LinearOperator

            n = 100

            dat = ones(n)
            A = spdiags([-2*dat, 4*dat, -dat], [-1,0,1],n,n)
            b = arange(n,dtype='d')

            L = spdiags([-dat/2, dat], [-1,0], n, n)
            U = spdiags([4*dat, -dat], [0,1], n, n)

            L_solver = splu(L)
            U_solver = splu(U)

            def L_solve(b):
                return L_solver.solve(b)

            def U_solve(b):
                return U_solver.solve(b)

            def LT_solve(b):
                return L_solver.solve(b,'T')

            def UT_solve(b):
                return U_solver.solve(b,'T')

            M1 = LinearOperator((n,n), matvec=L_solve, rmatvec=LT_solve)
            M2 = LinearOperator((n,n), matvec=U_solve, rmatvec=UT_solve)

            x,info = qmr(A, b, tol=1e-8, maxiter=15, M1=M1, M2=M2)

            assert_equal(info,0)
            assert_normclose(A*x, b, tol=1e-8) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [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.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def test_repr():
    A = interface.LinearOperator(shape=(1, 1), matvec=lambda x: 1)
    repr_A = repr(A)
    assert_('unspecified dtype' not in repr_A, repr_A) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def test_attributes():
    A = interface.aslinearoperator(np.arange(16).reshape(4, 4))

    def always_four_ones(x):
        x = np.asarray(x)
        assert_(x.shape == (3,) or x.shape == (3, 1))
        return np.ones(4)

    B = interface.LinearOperator(shape=(4, 3), matvec=always_four_ones)

    for op in [A, B, A * B, A.H, A + A, B + B, A ** 4]:
        assert_(hasattr(op, "dtype"))
        assert_(hasattr(op, "shape"))
        assert_(hasattr(op, "_matvec")) 

# 需要导入模块: from scipy.sparse.linalg import interface [as 别名]
# 或者: from scipy.sparse.linalg.interface import LinearOperator [as 别名]
def test_pickle():
    import pickle

    for protocol in range(pickle.HIGHEST_PROTOCOL + 1):
        A = interface.LinearOperator((3, 3), matvec)
        s = pickle.dumps(A, protocol=protocol)
        B = pickle.loads(s)

        for k in A.__dict__:
            assert_equal(getattr(A, k), getattr(B, k))