Python linalg.cg方法代码示例

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def test():
    print("\nTesting benchmark functions...")
    A, b, x0 = setup_input(n=50, sparsity=7)  # dense A
    Asp = to_sparse(A)
    x1, _ = calc_arrayfire(A, b, x0)
    x2, _ = calc_arrayfire(Asp, b, x0)
    if af.sum(af.abs(x1 - x2)/x2 > 1e-5):
        raise ValueError("arrayfire test failed")
    if np:
        An = to_numpy(A)
        bn = to_numpy(b)
        x0n = to_numpy(x0)
        x3, _ = calc_numpy(An, bn, x0n)
        if not np.allclose(x3, x1.to_list()):
            raise ValueError("numpy test failed")
    if sp:
        Asc = to_scipy_sparse(Asp)
        x4, _ = calc_scipy_sparse(Asc, bn, x0n)
        if not np.allclose(x4, x1.to_list()):
            raise ValueError("scipy.sparse test failed")
        x5, _ = calc_scipy_sparse_linalg_cg(Asc, bn, x0n)
        if not np.allclose(x5, x1.to_list()):
            raise ValueError("scipy.sparse.linalg.cg test failed")
    print("    all tests passed...") 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [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)
        me.values = cg(sp.eye(self.params.nvars[0] * self.params.nvars[1], format='csc') - factor * self.A,
                       rhs.values.flatten(), x0=u0.values.flatten(), tol=1E-12)[0]
        me.values = me.values.reshape(self.params.nvars)
        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def cg_diffusion(qsims, Wn, alpha = 0.99, maxiter = 10, tol = 1e-3):
    Wnn = eye(Wn.shape[0]) - alpha * Wn
    out_sims = []
    for i in range(qsims.shape[0]):
        #f,inf = s_linalg.cg(Wnn, qsims[i,:], tol=tol, maxiter=maxiter)
        f,inf = s_linalg.minres(Wnn, qsims[i,:], tol=tol, maxiter=maxiter)
        out_sims.append(f.reshape(-1,1))
    out_sims = np.concatenate(out_sims, axis = 1)
    ranks = np.argsort(-out_sims, axis = 0)
    return ranks 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def SetSolver(self,linear_solver="direct", linear_solver_type="umfpack",
        apply_preconditioner=False, preconditioner="amg_smoothed_aggregation",
        iterative_solver_tolerance=1.0e-12, reduce_matrix_bandwidth=False,
        geometric_discretisation=None):
        """

            input:
                linear_solver:          [str] type of solver either "direct",
                                        "iterative", "petsc" or "amg"

                linear_solver_type      [str] type of direct or linear solver to
                                        use, for instance "umfpack", "superlu" or
                                        "mumps" for direct solvers, or "cg", "gmres"
                                        etc for iterative solvers or "amg" for algebraic
                                        multigrid solver. See WhichSolvers method for
                                        the complete set of available linear solvers

                preconditioner:         [str] either "smoothed_aggregation",
                                        or "ruge_stuben" or "rootnode" for
                                        a preconditioner based on algebraic multigrid
                                        or "ilu" for scipy's spilu linear
                                        operator

                geometric_discretisation:
                                        [str] type of geometric discretisation used, for
                                        instance for FEM discretisations this would correspond
                                        to "tri", "quad", "tet", "hex" etc

        """

        self.solver_type = linear_solver
        self.solver_subtype = "umfpack"
        self.iterative_solver_tolerance = iterative_solver_tolerance
        self.apply_preconditioner = apply_preconditioner
        self.requires_cuthill_mckee = reduce_matrix_bandwidth
        self.geometric_discretisation = geometric_discretisation 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def WhichLinearSolvers(self):
        return {"direct":["superlu", "umfpack", "mumps", "pardiso"],
                "iterative":["cg", "bicg", "cgstab", "bicgstab", "gmres", "lgmres"],
                "amg":["cg", "bicg", "cgstab", "bicgstab", "gmres", "lgmres"],
                "petsc":["cg", "bicgstab", "gmres"]} 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def calc_scipy_sparse_linalg_cg(A, b, x0, maxiter=10):
    x = np.zeros(len(b), dtype=np.float32)
    x, _ = linalg.cg(A, b, x, tol=0., maxiter=maxiter)
    res = x0 - x
    return x, np.dot(res, res) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def bench(n=4*1024, sparsity=7, maxiter=10, iters=10):

    # generate data
    print("\nGenerating benchmark data for n = %i ..." %n)
    A, b, x0 = setup_input(n, sparsity)  # dense A
    Asp = to_sparse(A)  # sparse A
    input_info(A, Asp)

    # make benchmarks
    print("Benchmarking CG solver for n = %i ..." %n)
    t1 = timeit(calc_arrayfire, iters, args=(A, b, x0, maxiter))
    print("    arrayfire - dense:            %f ms" %t1)
    t2 = timeit(calc_arrayfire, iters, args=(Asp, b, x0, maxiter))
    print("    arrayfire - sparse:           %f ms" %t2)
    if np:
        An = to_numpy(A)
        bn = to_numpy(b)
        x0n = to_numpy(x0)
        t3 = timeit(calc_numpy, iters, args=(An, bn, x0n, maxiter))
        print("    numpy     - dense:            %f ms" %t3)
    if sp:
        Asc = to_scipy_sparse(Asp)
        t4 = timeit(calc_scipy_sparse, iters, args=(Asc, bn, x0n, maxiter))
        print("    scipy     - sparse:           %f ms" %t4)
        t5 = timeit(calc_scipy_sparse_linalg_cg, iters, args=(Asc, bn, x0n, maxiter))
        print("    scipy     - sparse.linalg.cg: %f ms" %t5) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [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
        """

        class context:
            num_iter = 0

        def callback(xk):
            context.num_iter += 1
            return context.num_iter

        me = self.dtype_u(self.init)

        Id = sp.eye(self.params.nvars[0] * self.params.nvars[1])

        me.values = cg(Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(), tol=self.params.lin_tol,
                       maxiter=self.params.lin_maxiter, callback=callback)[0]
        me.values = me.values.reshape(self.params.nvars)

        self.lin_ncalls += 1
        self.lin_itercount += context.num_iter

        return me


# noinspection PyUnusedLocal 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def solve_system_1(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
        """

        class context:
            num_iter = 0

        def callback(xk):
            context.num_iter += 1
            return context.num_iter

        me = self.dtype_u(self.init)

        Id = sp.eye(self.params.nvars[0] * self.params.nvars[1])

        me.values = cg(Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(), tol=self.params.lin_tol,
                       maxiter=self.params.lin_maxiter, callback=callback)[0]
        me.values = me.values.reshape(self.params.nvars)

        self.lin_ncalls += 1
        self.lin_itercount += context.num_iter

        return me 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def get_offline_result(i):
    ids = trunc_ids[i]
    trunc_lap = lap_alpha[ids][:, ids]
    scores, _ = linalg.cg(trunc_lap, trunc_init, tol=1e-6, maxiter=20)
    ranks = np.argsort(-scores)
    scores = scores[ranks]
    ranks = ids[ranks]
    return scores, ranks 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def ngstep(x0, obj0, objgrad0, obj_and_kl_func, hvpx0_func, max_kl, damping, max_cg_iter, enable_bt):
    '''
    Natural gradient step using hessian-vector products

    Args:
        x0: current point
        obj0: objective value at x0
        objgrad0: grad of objective value at x0
        obj_and_kl_func: function mapping a point x to the objective and kl values
        hvpx0_func: function mapping a vector v to the KL Hessian-vector product H(x0)v
        max_kl: max kl divergence limit. Triggers a line search.
        damping: multiple of I to mix with Hessians for Hessian-vector products
        max_cg_iter: max conjugate gradient iterations for solving for natural gradient step
    '''

    assert x0.ndim == 1 and x0.shape == objgrad0.shape

    # Solve for step direction
    damped_hvp_func = lambda v: hvpx0_func(v) + damping*v
    hvpop = ssl.LinearOperator(shape=(x0.shape[0], x0.shape[0]), matvec=damped_hvp_func)
    step, _ = ssl.cg(hvpop, -objgrad0, maxiter=max_cg_iter)
    fullstep = step / np.sqrt(.5 * step.dot(damped_hvp_func(step)) / max_kl + 1e-8)

    # Line search on objective with a hard KL wall
    if not enable_bt:
        return x0+fullstep, 0

    def barrierobj(p):
        obj, kl = obj_and_kl_func(p)
        return np.inf if kl > 2*max_kl else obj
    xnew, num_bt_steps = btlinesearch(
        f=barrierobj,
        x0=x0,
        fx0=obj0,
        g=objgrad0,
        dx=fullstep,
        accept_ratio=.1, shrink_factor=.5, max_steps=10)
    return xnew, num_bt_steps 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def _cg_wrapper(A, b, x0=None, tol=1e-5, maxiter=None):
        return cg(A, b, x0=x0, tol=tol, maxiter=maxiter, atol=0.0) 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def flow_matrix_row(G, weight='weight', dtype=float, solver='lu'):
    # Generate a row of the current-flow matrix
    import numpy as np
    from scipy import sparse
    from scipy.sparse import linalg
    solvername={"full" :FullInverseLaplacian,
                "lu": SuperLUInverseLaplacian,
                "cg": CGInverseLaplacian}
    n = G.number_of_nodes()
    L = laplacian_sparse_matrix(G, nodelist=range(n), weight=weight, 
                                dtype=dtype, format='csc')
    C = solvername[solver](L, dtype=dtype) # initialize solver
    w = C.w # w is the Laplacian matrix width
    # row-by-row flow matrix
    for u,v,d in G.edges_iter(data=True):
        B = np.zeros(w, dtype=dtype)
        c = d.get(weight,1.0)
        B[u%w] = c
        B[v%w] = -c
        # get only the rows needed in the inverse laplacian 
        # and multiply to get the flow matrix row
        row = np.dot(B, C.get_rows(u,v))  
        yield row,(u,v) 


# Class to compute the inverse laplacian only for specified rows
# Allows computation of the current-flow matrix without storing entire
# inverse laplacian matrix 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def solve(self,rhs):
        s = np.zeros(rhs.shape, dtype=self.dtype)
        s[1:]=linalg.cg(self.L1, rhs[1:], M=self.M)[0]
        return s 

# 需要导入模块: from scipy.sparse import linalg [as 别名]
# 或者: from scipy.sparse.linalg import cg [as 别名]
def solve_inverse(self,r):
        rhs = np.zeros(self.n, self.dtype)
        rhs[r] = 1
        return linalg.cg(self.L1, rhs[1:], M=self.M)[0]


# graph laplacian, sparse version, will move to linalg/laplacianmatrix.py