Python scipy.sparse.linalg.gmres() Examples

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) 

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 

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 

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 

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
# 

def gmres_no_fft(A, b):
    LOG.debug(f"Solve with GMRES for {A} without using FFT.")

    x, info = ssl.gmres(A.no_toeplitz() if isinstance(A, BlockMatrix) else A, b, atol=1e-6)

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

    return x 

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 

def gmres_linsolve(A, b):
    """

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

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 

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 

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 

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 

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 

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 

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
# 

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 

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 

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) 

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"]} 

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 

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'):
            dtype = M.dtype
        else:
            x = np.zeros(M.shape[1])
            dtype = (M * x).dtype
        LinearOperator.__init__(self, M.shape, self._matvec, dtype=dtype) 

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