Source code for gemseo.algos.linear_solvers.lib_scipy_linalg

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

# Contributors:
#    INITIAL AUTHORS - API and implementation and/or documentation
#        :author: Francois Gallard
#    OTHER AUTHORS   - MACROSCOPIC CHANGES

"""Wrappers for scipy's linear solvers."""

from __future__ import division, unicode_literals

import logging
from typing import Any, Dict, List, Mapping, Optional, Tuple, Union

from numpy import find_common_type, ndarray
from scipy.sparse import spmatrix
from scipy.sparse.base import issparse
from scipy.sparse.linalg import LinearOperator, bicg, bicgstab, gmres, lgmres, qmr, splu

from gemseo.algos.linear_solvers.linear_solver_lib import LinearSolverLib

LOGGER = logging.getLogger(__name__)


[docs]class ScipyLinalgAlgos(LinearSolverLib): """Wrapper for scipy linalg sparse linear solvers. Attributes: save_fpath (str): The path to the file to saved the problem when it is not converged and the option save_when_fail is active. methods_map (Dict[str, Callable]): The mapping between the solver names and the solvers methods in scipy.sparse. """ BASE_INFO_MSG = "scipy linear solver algorithm stop info: " OPTIONS_MAP = { "max_iter": "maxiter", "preconditioner": "M", "store_outer_av": "store_outer_Av", } LGMRES_SPEC_OPTS = ( "inner_m", "outer_k", "outer_v", "store_outer_av", "prepend_outer_v", ) __WEBSITE = "https://docs.scipy.org/doc/scipy/reference/generated/{}.html" __WEBPAGE = "scipy.sparse.linalg.{}" __WEBPAGES = { "BICG": __WEBPAGE.format("bicg"), "GMRES": __WEBPAGE.format("gmres"), "LGMRES": __WEBPAGE.format("lgmres"), "QMR": __WEBPAGE.format("qmr"), "BICGSTAB": __WEBPAGE.format("bicgstab"), "DEFAULT": __WEBPAGE.format("splu"), } def __init__(self): # type: (...) -> None super(ScipyLinalgAlgos, self).__init__() self.methods_map = { "LGMRES": lgmres, "GMRES": gmres, "BICG": bicg, "QMR": qmr, "BICGSTAB": bicgstab, "DEFAULT": self._run_default_solver, } self.lib_dict = { name: self.get_default_properties(name) for name in self.methods_map.keys() } self.lib_dict["DEFAULT"]["description"] = ( "This starts by LGMRES, but if it fails, " "switches to GMRES, then direct method super LU factorization." )
[docs] @classmethod def get_default_properties( cls, algo_name # type: str ): # type: (...) -> Dict[str, Union[bool,str]] """Return the properties of the algorithm. It states if it requires symmetric, or positive definite matrices for instance. Args: algo_name: The algorithm name. Returns: The properties of the solver. """ return { cls.LHS_MUST_BE_POSITIVE_DEFINITE: False, cls.LHS_MUST_BE_SYMMETRIC: False, cls.LHS_CAN_BE_LINEAR_OPERATOR: True, cls.INTERNAL_NAME: algo_name, "description": "Linear solver implemented in the SciPy library.", "website": cls.__WEBSITE.format(algo_name), }
def _get_options( self, max_iter=1000, # type: int preconditioner=None, # type: Optional[Union[ndarray, LinearOperator]] tol=1e-12, # type: float atol=None, # type: Optional[float] x0=None, # type: Optional[ndarray] use_ilu_precond=True, # type: bool inner_m=30, # type: int outer_k=3, # type: int outer_v=None, # type: Optional[List[Tuple]] store_outer_av=True, # type: bool prepend_outer_v=False, # type: bool save_when_fail=False, # type: bool store_residuals=False, # type: bool ): # type: (...) -> Dict[str, Any] """Check the options and sets the default values. See https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html Args: max_iter: The maximum number of iterations. preconditioner: The preconditionner, approximation of RHS^-1. If None, no preconditioner is used. tol: The relative tolerance for convergence, norm(RHS.dot(sol)) <= max(tol*norm(LHS), atol). atol: The absolute tolerance for convergence, norm(RHS.dot(sol)) <= max(tol*norm(LHS), atol). x0: The initial guess for the solution. M{sparse matrix, dense matrix, LinearOperator}. If None, solvers usually start from the null vector. inner_m int: The number of inner GMRES iterations per outer iteration. outer_k: The number of vectors to carry between inner GMRES iterations. outer_v: The data used to augment the Krylov subspace. store_outer_A: Whether LGMRES should store also A*v in addition to the vectors v in outer_v. prepend_outer_v: Whether to put outer_v augmentation vectors before the Krylov iterates. use_ilu_precond: Whether to use superLU to compute an incomplete LU factorization as preconditioner. save_when_fail: Whether to save the linear system to the disk when the solver failed to converge. store_residuals: Whether to store the residuals convergence history. Returns: The options. Raises: ValueError: When the LHR and RHS shapes are inconsistent, or when the preconditioner options are inconsistent. """ if preconditioner is not None and ( preconditioner.shape != self.problem.lhs.shape ): msg = "Inconsistent Preconditioner shape: %s != %s" raise ValueError(msg.format(preconditioner.shape, self.problem.lhs.shape)) if x0 is not None and len(x0) != len(self.problem.rhs): msg = "Inconsistent initial guess shape: %s != %s" raise ValueError(msg.format(x0.shape, self.problem.rhs.shape)) if use_ilu_precond and preconditioner is not None: raise ValueError( "Use either 'use_ilu_precond' or provide 'preconditioner'," " but not both." ) return self._process_options( max_iter=max_iter, preconditioner=preconditioner, tol=tol, atol=atol, x0=x0, inner_m=inner_m, outer_k=outer_k, outer_v=outer_v, store_outer_av=store_outer_av, prepend_outer_v=prepend_outer_v, use_ilu_precond=use_ilu_precond, save_when_fail=save_when_fail, store_residuals=store_residuals, ) def _run( self, **options # type: Union[None, bool, int, float, ndarray] ): # type: (...) -> ndarray """Run the algorithm. Args: **options: The options for the algorithm. Returns: The solution of the problem. """ if issparse(self.problem.rhs): self.problem.rhs = self.problem.rhs.toarray() rhs = self.problem.rhs lhs = self.problem.lhs opts_solver = options.copy() c_dtype = None if rhs.dtype != lhs.dtype and not isinstance(lhs, LinearOperator): c_dtype = find_common_type([rhs.dtype, lhs.dtype], []) if lhs.dtype != c_dtype: lhs = lhs.astype(c_dtype) if rhs.dtype != c_dtype: rhs = rhs.astype(c_dtype) if opts_solver["use_ilu_precond"] and not isinstance(lhs, LinearOperator): opts_solver["M"] = self._build_ilu_preconditioner(lhs, c_dtype) del opts_solver["use_ilu_precond"] del opts_solver["store_residuals"] method = self.methods_map[self.algo_name] if options.get("store_residuals", False): opts_solver["callback"] = self.__store_residuals opts_solver.pop("save_when_fail") self.problem.solution, info = method(lhs, rhs, **opts_solver) self._check_solver_info(info, options) return self.problem.solution def __store_residuals( self, current_x # type: ndarray ): # type: (...) -> ndarray """Store the current iteration residuals. Args: current_x: The current solution. """ self.problem.solution = current_x self.problem.compute_residuals(True, True) def _check_solver_info( self, info, # type: int options, # type: Mapping[str, Any] ): # type: (...) -> bool """Check the info returned by the solver. Args: info: The info value, negative, 0 or positive depending on status. options: The options passed to the solver. Returns: Whether the solver converged. Raises: RuntimeError: If the inputs are illegal for the solver. """ self.problem.is_converged = info == 0 if info > 0: if self.problem.solution is not None: LOGGER.warning( "%s, residual = %s", self.BASE_INFO_MSG, self.problem.compute_residuals(True), ) LOGGER.warning("info = %s", info) else: LOGGER.warning("Failed to converge") return False # check the dimensions if info < 0: raise RuntimeError( self.BASE_INFO_MSG + "illegal input or breakdown" ", options = %s", options, ) return True def _run_default_solver( self, lhs, # type: Union[ndarray, spmatrix, LinearOperator] rhs, # type: ndarray **options # type: Any ): # type: (...) -> Tuple[ndarray, int] """Run the default solver. This starts by LGMRES, but if it fails, switches to GMRES, then direct method super LU factorization. Args: lhs: The left hand side of the equation (matrix). rhs: The right hand side of the equation. **options: The user options. Returns: The last solution found and the info. """ # Try LGMRES first best_sol, info = lgmres(A=lhs, b=rhs, **options) min_res = self.problem.compute_residuals(True, current_x=best_sol) if self._check_solver_info(info, options): # If converged, stop return best_sol, info else: # Otherwise try GMRES min_res = self.problem.compute_residuals(True, current_x=best_sol) for k in self.LGMRES_SPEC_OPTS: # Adapt options if k in options: del options[k] if min_res < 1.0: options["x0"] = best_sol sol, info = gmres(A=lhs, b=rhs, **options) res = self.problem.compute_residuals(True, current_x=sol) if res < min_res: best_sol = sol options["x0"] = best_sol if self._check_solver_info(info, options): # pragma: no cover return best_sol, info # In this case previous runs failed, trying direct method # based on super LU a_fact = splu(lhs) sol = a_fact.solve(rhs) res = self.problem.compute_residuals(True, current_x=best_sol) if res < options["tol"]: # pragma: no cover best_sol = sol info = 0 self.problem.is_converged = True else: info = 1 return best_sol, info