Python linalg.gmres方法代碼示例

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [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

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

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def solve_system(self,rhs,factor,u0,t):
        """
        Simple linear solver for (I-dtA)u = rhs

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

        Returns:
            solution as mesh
        """

        b = rhs.values.flatten()
        # NOTE: A = -M, therefore solve Id + factor*M here
        sol, info =  LA.gmres( self.Id + factor*self.c_s*self.M, b, x0=u0.values.flatten(), tol=1e-13, restart=10, maxiter=20)
        me = mesh(self.nvars)
        me.values = unflatten(sol, 3, self.N[0], self.N[1])

        return me 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [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)

        if self.params.direct_solver:
            me.values = spsolve(self.Id - factor * self.A, rhs.values.flatten())
        else:
            me.values = gmres(self.Id - factor * self.A, rhs.values.flatten(), x0=u0.values.flatten(),
                              tol=self.params.lintol, maxiter=self.params.liniter)[0]
        me.values = me.values.reshape(self.params.nvars)
        return me 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def __init__(self, M, 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(M.dtype).eps
        self.M = M
        self.ifunc = ifunc
        self.tol = tol
        if hasattr(M, 'dtype'):
            self.dtype = M.dtype
        else:
            x = np.zeros(M.shape[1])
            self.dtype = (M * x).dtype
        self.shape = M.shape 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [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 gmres [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 gmres [as 別名]
def gmres_loose(A, b, tol):
    """
    gmres with looser termination condition.
    """
    b = np.asarray(b)
    min_tol = 1000 * np.sqrt(b.size) * np.finfo(b.dtype).eps
    return gmres(A, b, tol=max(tol, min_tol), atol=0) 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def solve_system(self, rhs, factor, u0, t):
        """
        Simple linear solver for (I-dtA)u = rhs using GMRES

        Args:
            rhs (dtype_f): right-hand side for the nonlinear 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
        """

        b = rhs.values.flatten()
        cb = Callback()

        sol, info = gmres(self.Id - factor * self.M, b, x0=u0.values.flatten(), tol=self.params.gmres_tol_limit,
                          restart=self.params.gmres_restart, maxiter=self.params.gmres_maxiter, callback=cb)
        # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
        if factor != 0.0:
            self.gmres_logger.add(cb.getcounter())
        me = self.dtype_u(self.init)
        me.values = unflatten(sol, 4, self.N[0], self.N[1])

        return me 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def f_fast_solve(self, rhs, alpha, u0):
        cb = Callback()
        sol, info = gmres(self.problem.Id - alpha * self.problem.M, rhs, x0=u0,
                          tol=self.problem.params.gmres_tol_limit, restart=self.problem.params.gmres_restart,
                          maxiter=self.problem.params.gmres_maxiter, callback=cb)
        if alpha != 0.0:
            self.logger.add(cb.getcounter())
        return sol


#
# Trapezoidal rule
# 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def f_solve(self, b, alpha, u0):
        cb = Callback()
        sol, info = gmres(self.problem.Id - alpha * (self.problem.D_upwind + self.problem.M), b, x0=u0,
                          tol=self.problem.params.gmres_tol_limit, restart=self.problem.params.gmres_restart,
                          maxiter=self.problem.params.gmres_maxiter, callback=cb)
        if alpha != 0.0:
            self.logger.add(cb.getcounter())
        return sol


#
#  Split-Explicit method
# 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def __init__(self,
               A,
               drop_tol=0.005,
               fill_factor=2.0,
               normalize_inplace=False):
    # the spilu and gmres functions are most efficient with csc sparse. If the
    # matrix is already csc then this will do nothing
    A = sp.csc_matrix(A)
    n = row_norms(A)
    if normalize_inplace:
      divide_rows(A, n, inplace=True)
    else:
      A = divide_rows(A, n, inplace=False).tocsc()

    LOGGER.debug(
      'computing the ILU decomposition of a %s by %s sparse matrix with %s '
      'nonzeros ' % (A.shape + (A.nnz,)))
    ilu = spla.spilu(
      A,
      drop_rule='basic',
      drop_tol=drop_tol,
      fill_factor=fill_factor)
    LOGGER.debug('done')
    M = spla.LinearOperator(A.shape, ilu.solve)
    self.A = A
    self.M = M
    self.n = n 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def solve(self, b, tol=1.0e-10):
    '''
    Solve `Ax = b` for `x`

    Parameters
    ----------
    b : (n,) array

    tol : float, optional

    Returns
    -------
    (n,) array

    '''
    # solve the system using GMRES and define the callback function to
    # print info for each iteration
    def callback(res, _itr=[0]):
      l2 = np.linalg.norm(res)
      LOGGER.debug('GMRES error on iteration %s: %s' % (_itr[0], l2))
      _itr[0] += 1

    LOGGER.debug('solving the system with GMRES')
    x, info = spla.gmres(
      self.A,
      b/self.n,
      tol=tol,
      M=self.M,
      callback=callback)
    LOGGER.debug('finished GMRES with info %s' % info)
    return x 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def krylovMethod(self,tol=1e-8): 
        """
        We obtain ``pi`` by using the :func:``gmres`` solver for the system of linear equations. 
        It searches in Krylov subspace for a vector with minimal residual. The result is stored in the class attribute ``pi``.   

        Example
        -------
        >>> P = np.array([[0.5,0.5],[0.6,0.4]])
        >>> mc = markovChain(P)
        >>> mc.krylovMethod()
        >>> print(mc.pi) 
        [ 0.54545455  0.45454545]
        
        Parameters
        ----------
        tol : float, optional(default=1e-8)
            Tolerance level for the precision of the end result. A lower tolerance leads to more accurate estimate of ``pi``.        
        
        Remarks
        -------
        For large state spaces, this method may not always give a solution. 
        Code due to http://stackoverflow.com/questions/21308848/
        """            
        P       = self.getIrreducibleTransitionMatrix()
        
        #if P consists of one element, then set self.pi = 1.0
        if P.shape == (1, 1):
            self.pi = np.array([1.0]) 
            return
            
        size    = P.shape[0]
        dP      = P - eye(size)
        #Replace the first equation by the normalizing condition.
        A       = vstack([np.ones(size), dP.T[1:,:]]).tocsr()
        rhs     = np.zeros((size,))
        rhs[0]  = 1
                
        pi, info = gmres(A, rhs, tol=tol)
        if info != 0:
            raise RuntimeError("gmres did not converge")
        self.pi = pi 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def gmres_linsolve(A, b):
    """

    :param A:
    :param b:
    :return:
    """
    x, info = gmres(A, b)
    return x 

# 需要導入模塊: from scipy.sparse import linalg [as 別名]
# 或者: from scipy.sparse.linalg import gmres [as 別名]
def solve_gmres(A, b):
    LOG.debug(f"Solve with GMRES for {A}.")

    if LOG.isEnabledFor(logging.DEBUG):
        counter = Counter()
        x, info = ssl.gmres(A, b, atol=1e-6, callback=counter)
        LOG.debug(f"End of GMRES after {counter.nb_iter} iterations.")

    else:
        x, info = ssl.gmres(A, b, atol=1e-6)

    if info != 0:
        LOG.warning(f"No convergence of the GMRES. Error code: {info}")

    return x