From 0bacdd014de4db85018ca4bc301eb67df221f9be Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Wed, 6 Dec 2023 13:39:19 +0100 Subject: [PATCH 01/34] Changed the handling of preconditioners in inverse LinearOperators. The option `pc` must now take as value a `LinearOperator`, which is the preconditioner. If `pc=None`, an `IdentityOperator` of the domain space is created automatically. The class `InverseLinearOperators` could be deprecated, since the Jacobi preconditiner should be implemented as a `LinearOperator`. The new solver `PBiConjugateGradientStabilized` (preconditioned) has been added from Struphy. `test_solvers.py` has been updated. The Jacobi preconditioner is not tested at the moment. --- psydac/linalg/solvers.py | 309 +++++++++++++++++++++++++--- psydac/linalg/tests/test_solvers.py | 16 +- 2 files changed, 289 insertions(+), 36 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 93817ac51..8519bc79f 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -44,7 +44,7 @@ def inverse(A, solver, **kwargs): """ # Check solver input - solvers = ('cg', 'pcg', 'bicg', 'bicgstab', 'minres', 'lsmr', 'gmres') + solvers = ('cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres') if solver not in solvers: raise ValueError(f"Required solver '{solver}' not understood.") @@ -65,6 +65,8 @@ def inverse(A, solver, **kwargs): obj = BiConjugateGradient(A, **kwargs) elif solver == 'bicgstab': obj = BiConjugateGradientStabilized(A, **kwargs) + elif solver == 'pbicgstab': + obj = PBiConjugateGradientStabilized(A, **kwargs) elif solver == 'minres': obj = MinimumResidual(A, **kwargs) elif solver == 'lsmr': @@ -268,7 +270,7 @@ class PConjugateGradient(InverseLinearOperator): The .dot (and also the .solve) function are based on a preconditioned conjugate gradient method. Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte - system Ax = b. It assumes that pc(r) returns the solution to Ps = r, + system Ax = b. It assumes that pc.dot(r) returns the solution to Ps = r, where P is positive definite. Parameters @@ -276,14 +278,8 @@ class PConjugateGradient(InverseLinearOperator): A : psydac.linalg.stencil.StencilMatrix Left-hand-side matrix A of linear system - pc: str - Preconditioner for A, it should approximate the inverse of A. - Can currently only be: - * The strings 'jacobi' or 'weighted_jacobi'. (rather obsolete, supply a callable instead, if possible)(14.02.: test weighted_jacobi) - The following should also be possible - * None, i.e. not pre-conditioning (this calls the standard `cg` method) - * A LinearSolver object (in which case the out parameter is used) - * A callable with two parameters (A, r), where A is the LinearOperator from above, and r is the residual. + pc: psydac.linalg.basic.LinearOperator + Preconditioner for A, it should approximate the inverse of A (can be None). x0 : psydac.linalg.basic.Vector First guess of solution for iterative solver (optional). @@ -301,7 +297,7 @@ class PConjugateGradient(InverseLinearOperator): Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems """ - def __init__(self, A, *, pc='jacobi', x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): + def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): assert isinstance(A, LinearOperator) assert A.domain.dimension == A.codomain.dimension @@ -318,6 +314,8 @@ def __init__(self, A, *, pc='jacobi', x0=None, tol=1e-6, maxiter=1000, verbose=F self._domain = domain self._codomain = codomain self._solver = 'pcg' + if pc is None: + pc = IdentityOperator(self._domain) self._options = {"x0":x0, "pc":pc, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} self._check_options(**self._options) tmps_codomain = {key: codomain.zeros() for key in ("p", "s")} @@ -329,8 +327,7 @@ def _check_options(self, **kwargs): for key, value in kwargs.items(): if key == 'pc': - assert value is not None, "pc may not be None" - assert value == 'jacobi', "unsupported preconditioner" + assert isinstance(value, LinearOperator) elif key == 'x0': if value is not None: assert isinstance(value, Vector), "x0 must be a Vector or None" @@ -357,7 +354,7 @@ def transpose(self, conjugate=False): def solve(self, b, out=None): """ Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte - system Ax = b. It assumes that pc(r) returns the solution to Ps = r, + system Ax = b. It assumes that pc.dot(r) returns the solution to Ps = r, where P is positive definite. Info can be accessed using get_info(), see :func:~`basic.InverseLinearOperator.get_info`. @@ -390,6 +387,8 @@ def solve(self, b, out=None): assert isinstance(b, Vector) assert b.space is domain + + assert isinstance(pc, LinearOperator) # First guess of solution if out is not None: @@ -398,22 +397,6 @@ def solve(self, b, out=None): x = x0.copy(out=out) - # Preconditioner - assert pc is not None - if pc == 'jacobi': - psolve = lambda r, out: InverseLinearOperator.jacobi(A, r, out) - #elif pc == 'weighted_jacobi': - # psolve = lambda r, out: InverseLinearOperator.weighted_jacobi(A, r, out) # allows for further specification not callable like this! - #elif isinstance(pc, str): - # pcfun = getattr(InverseLinearOperator, str) - # #pcfun = globals()[pc] - # psolve = lambda r: pcfun(A, r) - #elif isinstance(pc, LinearSolver): - # s = b.space.zeros() - # psolve = lambda r: pc.solve(r, out=s) - #elif hasattr(pc, '__call__'): - # psolve = lambda r: pc(A, r) - # Extract local storage v = self._tmps["v"] r = self._tmps["r"] @@ -425,7 +408,7 @@ def solve(self, b, out=None): b.copy(out=r) r -= v nrmr_sqr = r.dot(r).real - psolve(r, out=s) + pc.dot(r, out=s) am = s.dot(r) s.copy(out=p) @@ -453,7 +436,7 @@ def solve(self, b, out=None): r.mul_iadd(-l, v) # this is r -= l*v nrmr_sqr = r.dot(r).real - psolve(r, out=s) + pc.dot(r, out=s) am1 = s.dot(r) @@ -708,6 +691,268 @@ def solve(self, b, out=None): def dot(self, b, out=None): return self.solve(b, out=out) +#=============================================================================== +class PBiConjugateGradientStabilized(InverseLinearOperator): + """ + A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. + The .dot (and also the .solve) function are based on the + preconditioned Biconjugate gradient Stabilized (PBCGSTAB) algorithm for solving linear system Ax=b. + Implementation from [1], page 251. + Parameters + ---------- + A : psydac.linalg.basic.LinearOperator + Left-hand-side matrix A of linear system; individual entries A[i,j] + can't be accessed, but A has 'shape' attribute and provides 'dot(p)' + function (i.e. matrix-vector product A*p). + pc: psydac.linalg.basic.LinearOperator + Preconditioner for A, it should approximate the inverse of A (can be None). + x0 : psydac.linalg.basic.Vector + First guess of solution for iterative solver (optional). + tol : float + Absolute tolerance for 2-norm of residual r = A*x - b. + maxiter: int + Maximum number of iterations. + verbose : bool + If True, 2-norm of residual r is printed at each iteration. + References + ---------- + [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015. + """ + def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False): + + assert isinstance(A, LinearOperator) + assert A.domain.dimension == A.codomain.dimension + domain = A.codomain + codomain = A.domain + + if x0 is not None: + assert isinstance(x0, Vector) + assert x0.space is codomain + else: + x0 = codomain.zeros() + + self._A = A + self._domain = domain + self._codomain = codomain + self._solver = 'pbicgstab' + if pc is None: + pc = IdentityOperator(self._domain) + self._options = { "pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose} + self._check_options(**self._options) + self._tmps = {key: domain.zeros() for key in ("v", "r", "s", "t", + "vp", "rp", "sp", "tp", + "pp", "av", "app", "osp", + "rp0")} + self._info = None + + def _check_options(self, **kwargs): + keys = ('pc', 'x0', 'tol', 'maxiter', 'verbose') + for key, value in kwargs.items(): + idx = [key == keys[i] for i in range(len(keys))] + assert any(idx), "key not supported, check options" + true_idx = idx.index(True) + if true_idx == 0: + assert isinstance(value, LinearOperator) + elif true_idx == 1: + if value is not None: + assert isinstance(value, Vector), "x0 must be a Vector or None" + assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" + elif true_idx == 2: + assert is_real(value), "tol must be a real number" + assert value > 0, "tol must be positive" + elif true_idx == 3: + assert isinstance(value, int), "maxiter must be an int" + assert value > 0, "maxiter must be positive" + elif true_idx == 4: + assert isinstance(value, bool), "verbose must be a bool" + + def _update_options( self ): + self._options = {"pc": self._pc, "x0":self._x0, "tol":self._tol, "maxiter": self._maxiter, "verbose": self._verbose} + + def transpose(self, conjugate=False): + At = self._A.transpose(conjugate=conjugate) + solver = self._solver + options = self._options + return inverse(At, solver, **options) + + def solve(self, b, out=None): + """ + Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. + Implementation from [1], page 251. + Parameters + ---------- + b : psydac.linalg.basic.Vector + Right-hand-side vector of linear system. Individual entries b[i] need + not be accessed, but b has 'shape' attribute and provides 'copy()' and + 'dot(p)' functions (dot(p) is the vector inner product b*p ); moreover, + scalar multiplication and sum operations are available. + out : psydac.linalg.basic.Vector | NoneType + The output vector, or None (optional). + Returns + ------- + x : psydac.linalg.basic.Vector + Numerical solution of linear system. To check the convergence of the solver, + use the method InverseLinearOperator.get_info(). + + info : dict + Dictionary containing convergence information: + - 'niter' = (int) number of iterations + - 'success' = (boolean) whether convergence criteria have been met + - 'res_norm' = (float) 2-norm of residual vector r = A*x - b. + References + ---------- + [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015. + """ + + A = self._A + domain = self._domain + codomain = self._codomain + options = self._options + pc = options["pc"] + x0 = options["x0"] + tol = options["tol"] + maxiter = options["maxiter"] + verbose = options["verbose"] + + assert isinstance(b, Vector) + assert b.space is domain + + assert isinstance(pc, LinearOperator) + + # first guess of solution + if out is not None: + assert isinstance(out, Vector) + assert out.space == codomain + out *= 0 + if x0 is None: + x = out + else: + assert x0.shape == (A.shape[0],) + out += x0 + x = out + else: + if x0 is None: + x = b.copy() + x *= 0.0 + else: + assert x0.shape == (A.shape[0],) + x = x0.copy() + + # preconditioner (must have a .solve method) + assert isinstance(pc, LinearOperator) + + # extract temporary vectors + v = self._tmps['v'] + r = self._tmps['r'] + s = self._tmps['s'] + t = self._tmps['t'] + + vp = self._tmps['vp'] + rp = self._tmps['rp'] + pp = self._tmps['pp'] + sp = self._tmps['sp'] + tp = self._tmps['tp'] + + av = self._tmps['av'] + + app = self._tmps['app'] + osp = self._tmps['osp'] + + # first values: r = b - A @ x, rp = pp = PC @ r, rhop = |rp|^2 + A.dot(x, out=v) + b.copy(out=r) + r -= v + + pc.dot(r, out=rp) + rp.copy(out=pp) + + rhop = rp.dot(rp) + + # save initial residual vector rp0 + rp0 = self._tmps['rp0'] + rp.copy(out=rp0) + + # squared residual norm and squared tolerance + res_sqr = r.dot(r) + tol_sqr = tol**2 + + if verbose: + print("Pre-conditioned BICGSTAB solver:") + print("+---------+---------------------+") + print("+ Iter. # | L2-norm of residual |") + print("+---------+---------------------+") + template = "| {:7d} | {:19.2e} |" + + # iterate to convergence or maximum number of iterations + niter = 0 + + while res_sqr > tol_sqr and niter < maxiter: + + # v = A @ pp, vp = PC @ v, alphap = rhop/(vp.rp0) + A.dot(pp, out=v) + pc.dot(v, out=vp) + alphap = rhop / vp.dot(rp0) + + # s = r - alphap*v, sp = PC @ s + r.copy(out=s) + v.copy(out=av) + av *= alphap + s -= av + pc.dot(s, out=sp) + + # t = A @ sp, tp = PC @ t, omegap = (tp.sp)/(tp.tp) + A.dot(sp, out=t) + pc.dot(t, out=tp) + omegap = tp.dot(sp) / tp.dot(tp) + + # x = x + alphap*pp + omegap*sp + pp.copy(out=app) + sp.copy(out=osp) + app *= alphap + osp *= omegap + x += app + x += osp + + # r = s - omegap*t, rp = sp - omegap*tp + s.copy(out=r) + t *= omegap + r -= t + + sp.copy(out=rp) + tp *= omegap + rp -= tp + + # rhop_new = rp.rp0, betap = (alphap*rhop_new)/(omegap*rhop) + rhop_new = rp.dot(rp0) + betap = (alphap*rhop_new) / (omegap*rhop) + rhop = 1*rhop_new + + # pp = rp + betap*(pp - omegap*vp) + vp *= omegap + pp -= vp + pp *= betap + pp += rp + + # new residual norm + res_sqr = r.dot(r) + + niter += 1 + + if verbose: + print(template.format(niter, sqrt(res_sqr))) + + if verbose: + print("+---------+---------------------+") + + # convergence information + self._info = {'niter': niter, 'success': res_sqr < + tol_sqr, 'res_norm': sqrt(res_sqr)} + + return x + + def dot(self, b, out=None): + return self.solve(b, out=out) + #=============================================================================== class BiConjugateGradientStabilized(InverseLinearOperator): """ diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 803877257..793b9e89c 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -59,19 +59,24 @@ def define_data(n, p, matrix_data, dtype=float): @pytest.mark.parametrize( 'n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @pytest.mark.parametrize('dtype', [float, complex]) -@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'minres', 'lsmr', 'gmres']) +@pytest.mark.parametrize('pc', [None]) +@pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) -def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): +def test_solver_tridiagonal(n, p, dtype, pc, solver, verbose=False): #--------------------------------------------------------------------------- # PARAMETERS #--------------------------------------------------------------------------- - if solver in ['bicg', 'bicgstab', 'lsmr']: + if solver in ['bicg', 'bicgstab', 'pbicgstab', 'lsmr']: if dtype==complex: diagonals = [1-10j,6+9j,3+5j] else: diagonals = [1,6,3] + + if solver == 'pbicgstab' and dtype == complex: + # pbicgstab only works for real matrices + return elif solver == 'gmres': if dtype==complex: diagonals = [-7-2j,-6-2j,-1-10j] @@ -102,7 +107,10 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): print() #Create the solvers - solv = inverse(A, solver, tol=tol, verbose=False, recycle=True) + if solver in ['pcg', 'pbicgstab']: + solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=True) + else: + solv = inverse(A, solver, tol=1e-13, verbose=True) solvt = solv.transpose() solvh = solv.H solv2 = inverse(A@A, solver, tol=1e-13, verbose=True, recycle=True) # Test solver of composition of operators From 11965fae552fa8ebf13f181f190cdd80ba7a815b Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Wed, 6 Dec 2023 14:19:32 +0100 Subject: [PATCH 02/34] Moved _check_options() and transpose() from the solver classes to the base class `InverseLinearOperator`. All solvers call now `super().__init__` which takes the base arguments for any iterative solver. Additional arguments must be checked in the individual solvers. The base class now has the setter method for `linop`, which allows to change tha matrix A of the solver. Moreover, there is the abstract property `solver` which is the name string of the solver. --- psydac/linalg/basic.py | 80 +++- psydac/linalg/solvers.py | 818 +++++++++++++-------------------------- 2 files changed, 342 insertions(+), 556 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index dc86957b4..f8a137053 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -879,11 +879,39 @@ class InverseLinearOperator(LinearOperator): Iterative solver for square linear system Ax=b, where x and b belong to (normed) vector space V. + Parameters + ---------- + A : psydac.linalg.basic.LinearOperator + Left-hand-side matrix A of linear system. + + x0 : psydac.linalg.basic.Vector + First guess of solution for iterative solver (optional). + + tol : float + Absolute tolerance for L2-norm of residual r = A*x - b. + + maxiter: int + Maximum number of iterations. + + verbose : bool + If True, L2-norm of residual r is printed at each iteration. """ - @property - def space(self): - return self._space + def __init__(self, A, **kwargs): + + assert isinstance(A, LinearOperator) + assert A.domain.dimension == A.codomain.dimension + domain = A.codomain + codomain = A.domain + + if kwargs['x0'] is None: + kwargs['x0'] = codomain.zeros() + + self._A = A + self._domain = domain + self._codomain = codomain + self._options = kwargs + self._check_options(**self._options) @property def domain(self): @@ -900,10 +928,39 @@ def dtype(self): @property def linop(self): return self._A + + @linop.setter + def linop(self, a): + assert isinstance(a, LinearOperator) + assert a.domain is self.domain + assert a.codomain is self.codomain + self._A = a @property def options(self): return self._options + + @property + @abstractmethod + def solver(self): + "String that identifies the solver." + pass + + def _check_options(self, **kwargs): + for key, value in kwargs.items(): + + if key == 'x0': + if value is not None: + assert isinstance(value, Vector), "x0 must be a Vector or None" + assert value.space == self.codomain, "x0 belongs to the wrong VectorSpace" + elif key == 'tol': + assert is_real(value), "tol must be a real number" + assert value > 0, "tol must be positive" + elif key == 'maxiter': + assert isinstance(value, int), "maxiter must be an int" + assert value > 0, "maxiter must be positive" + elif key == 'verbose': + assert isinstance(value, bool), "verbose must be a bool" def toarray(self): raise NotImplementedError('toarray() is not defined for InverseLinearOperators.') @@ -921,13 +978,13 @@ def set_options(self, **kwargs): self._check_options(**kwargs) self._options.update(kwargs) - @abstractmethod - def _check_options(self, **kwargs): - pass - - @abstractmethod def transpose(self, conjugate=False): - pass + from psydac.linalg.solvers import inverse + + At = self.linop.transpose(conjugate=conjugate) + solver = self.solver + options = self._options + return inverse(At, solver, **options) @staticmethod def jacobi(A, b, out=None): @@ -1016,3 +1073,8 @@ def solve(self, rhs, out=None): @property def T(self): return self.transpose() + + +def is_real(x): + from numbers import Number + return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 8519bc79f..6f21761d4 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -113,51 +113,16 @@ class ConjugateGradient(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() - - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'cg' self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("v", "r", "p")} + + super().__init__(A, **self._options) + + self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p")} self._info = None - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) + @property + def solver(self): + return 'cg' def solve(self, b, out=None): """ @@ -299,57 +264,23 @@ class PConjugateGradient(InverseLinearOperator): """ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() - - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'pcg' - if pc is None: - pc = IdentityOperator(self._domain) self._options = {"x0":x0, "pc":pc, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) - tmps_codomain = {key: codomain.zeros() for key in ("p", "s")} - tmps_domain = {key: domain.zeros() for key in ("v", "r")} + + super().__init__(A, **self._options) + + if pc is None: + self._options['pc'] = IdentityOperator(self.domain) + else: + assert isinstance(pc, LinearOperator) + + tmps_codomain = {key: self.codomain.zeros() for key in ("p", "s")} + tmps_domain = {key: self.domain.zeros() for key in ("v", "r")} self._tmps = {**tmps_codomain, **tmps_domain} self._info = None - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'pc': - assert isinstance(value, LinearOperator) - elif key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) + @property + def solver(self): + return 'pcg' def solve(self, b, out=None): """ @@ -501,52 +432,17 @@ class BiConjugateGradient(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() - - self._A = A - self._Ah = A.H - self._domain = domain - self._codomain = codomain - self._solver = 'bicg' self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")} + + super().__init__(A, **self._options) + + self._Ah = A.H + self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")} self._info = None - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) + @property + def solver(self): + return 'bicg' def solve(self, b, out=None): """ @@ -691,6 +587,207 @@ def solve(self, b, out=None): def dot(self, b, out=None): return self.solve(b, out=out) +#=============================================================================== +class BiConjugateGradientStabilized(InverseLinearOperator): + """ + A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. + + The .dot (and also the .solve) function are based on the + Biconjugate gradient Stabilized (BCGSTAB) algorithm for solving linear system Ax=b. + Implementation from [1], page 175. + + Parameters + ---------- + A : psydac.linalg.basic.LinearOperator + Left-hand-side matrix A of linear system; individual entries A[i,j] + can't be accessed, but A has 'shape' attribute and provides 'dot(p)' + function (i.e. matrix-vector product A*p). + + x0 : psydac.linalg.basic.Vector + First guess of solution for iterative solver (optional). + + tol : float + Absolute tolerance for 2-norm of residual r = A*x - b. + + maxiter: int + Maximum number of iterations. + + verbose : bool + If True, 2-norm of residual r is printed at each iteration. + + recycle : bool + Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems + + References + ---------- + [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015. + + """ + def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): + + self._options = {"x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose, "recycle":recycle} + + super().__init__(A, **self._options) + + self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vr", "r0")} + self._info = None + + @property + def solver(self): + return 'bicgstab' + + def solve(self, b, out=None): + """ + Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b. + Implementation from [1], page 175. + ToDo: Add optional preconditioner + + Parameters + ---------- + b : psydac.linalg.basic.Vector + Right-hand-side vector of linear system. Individual entries b[i] need + not be accessed, but b has 'shape' attribute and provides 'copy()' and + 'dot(p)' functions (dot(p) is the vector inner product b*p ); moreover, + scalar multiplication and sum operations are available. + out : psydac.linalg.basic.Vector | NoneType + The output vector, or None (optional). + + Returns + ------- + x : psydac.linalg.basic.Vector + Numerical solution of linear system. To check the convergence of the solver, + use the method InverseLinearOperator.get_info(). + + info : dict + Dictionary containing convergence information: + - 'niter' = (int) number of iterations + - 'success' = (boolean) whether convergence criteria have been met + - 'res_norm' = (float) 2-norm of residual vector r = A*x - b. + + References + ---------- + [1] H. A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the + solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comp., 13(2):631–644, 1992. + """ + + A = self._A + domain = self._domain + codomain = self._codomain + options = self._options + x0 = options["x0"] + tol = options["tol"] + maxiter = options["maxiter"] + verbose = options["verbose"] + recycle = options["recycle"] + + assert isinstance(b, Vector) + assert b.space is domain + + # First guess of solution + if out is not None: + assert isinstance(out, Vector) + assert out.space is codomain + + x = x0.copy(out=out) + + # Extract local storage + v = self._tmps["v"] + r = self._tmps["r"] + p = self._tmps["p"] + vr = self._tmps["vr"] + r0 = self._tmps["r0"] + + # First values + A.dot(x, out=v) + b.copy(out=r) + r -= v + #r = b - A.dot(x) + r.copy(out=p) + v *= 0.0 + vr *= 0.0 + + r.copy(out=r0) + + res_sqr = r.dot(r).real + tol_sqr = tol ** 2 + + if verbose: + print("BiCGSTAB solver:") + print("+---------+---------------------+") + print("+ Iter. # | L2-norm of residual |") + print("+---------+---------------------+") + template = "| {:7d} | {:19.2e} |" + + # Iterate to convergence + for m in range(1, maxiter + 1): + + if res_sqr < tol_sqr: + m -= 1 + break + + # ----------------------- + # MATRIX-VECTOR PRODUCTS + # ----------------------- + v = A.dot(p, out=v) + # ----------------------- + + # c := (r0, r) + c = r0.dot(r) + + # a := (r0, r) / (r0, v) + a = c / (r0.dot(v)) + + # r := r - a*v + r.mul_iadd(-a, v) + + # vr := A*r + vr = A.dot(r, out=vr) + + # w := (r, A*r) / (A*r, A*r) + w = r.dot(vr) / vr.dot(vr) + + # ----------------------- + # SOLUTION UPDATE + # ----------------------- + # x := x + a*p +w*r + x.mul_iadd(a, p) + x.mul_iadd(w, r) + # ----------------------- + + # r := r - w*A*r + r.mul_iadd(-w, vr) + + # ||r||_2 := (r, r) + res_sqr = r.dot(r).real + + if res_sqr < tol_sqr: + break + + # b := a / w * (r0, r)_{m+1} / (r0, r)_m + b = r0.dot(r) * a / (c * w) + + # p := r + b*p- b*w*v + p *= b + p += r + p.mul_iadd(-b * w, v) + + if verbose: + print(template.format(m, sqrt(res_sqr))) + + if verbose: + print("+---------+---------------------+") + + # Convergence information + self._info = {'niter': m, 'success': res_sqr < tol_sqr, 'res_norm': sqrt(res_sqr)} + + if recycle: + x.copy(out=self._options["x0"]) + + return x + + def dot(self, b, out=None): + return self.solve(b, out=out) + #=============================================================================== class PBiConjugateGradientStabilized(InverseLinearOperator): """ @@ -720,61 +817,28 @@ class PBiConjugateGradientStabilized(InverseLinearOperator): """ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain + self._options = { "pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose} - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() + super().__init__(A, **self._options) - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'pbicgstab' if pc is None: - pc = IdentityOperator(self._domain) - self._options = { "pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("v", "r", "s", "t", + self._options['pc'] = IdentityOperator(self.domain) + else: + assert isinstance(pc, LinearOperator) + + self._tmps = {key: self.domain.zeros() for key in ("v", "r", "s", "t", "vp", "rp", "sp", "tp", "pp", "av", "app", "osp", "rp0")} self._info = None - def _check_options(self, **kwargs): - keys = ('pc', 'x0', 'tol', 'maxiter', 'verbose') - for key, value in kwargs.items(): - idx = [key == keys[i] for i in range(len(keys))] - assert any(idx), "key not supported, check options" - true_idx = idx.index(True) - if true_idx == 0: - assert isinstance(value, LinearOperator) - elif true_idx == 1: - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif true_idx == 2: - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif true_idx == 3: - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif true_idx == 4: - assert isinstance(value, bool), "verbose must be a bool" + @property + def solver(self): + return 'pbicgstab' def _update_options( self ): self._options = {"pc": self._pc, "x0":self._x0, "tol":self._tol, "maxiter": self._maxiter, "verbose": self._verbose} - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) - def solve(self, b, out=None): """ Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. @@ -953,242 +1017,6 @@ def solve(self, b, out=None): def dot(self, b, out=None): return self.solve(b, out=out) -#=============================================================================== -class BiConjugateGradientStabilized(InverseLinearOperator): - """ - A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. - - The .dot (and also the .solve) function are based on the - Biconjugate gradient Stabilized (BCGSTAB) algorithm for solving linear system Ax=b. - Implementation from [1], page 175. - - Parameters - ---------- - A : psydac.linalg.basic.LinearOperator - Left-hand-side matrix A of linear system; individual entries A[i,j] - can't be accessed, but A has 'shape' attribute and provides 'dot(p)' - function (i.e. matrix-vector product A*p). - - x0 : psydac.linalg.basic.Vector - First guess of solution for iterative solver (optional). - - tol : float - Absolute tolerance for 2-norm of residual r = A*x - b. - - maxiter: int - Maximum number of iterations. - - verbose : bool - If True, 2-norm of residual r is printed at each iteration. - - recycle : bool - Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems - - References - ---------- - [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015. - - """ - def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() - - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'bicgstab' - self._options = {"x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose, "recycle":recycle} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("v", "r", "p", "vr", "r0")} - self._info = None - - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) - - def solve(self, b, out=None): - """ - Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b. - Implementation from [1], page 175. - ToDo: Add optional preconditioner - - Parameters - ---------- - b : psydac.linalg.basic.Vector - Right-hand-side vector of linear system. Individual entries b[i] need - not be accessed, but b has 'shape' attribute and provides 'copy()' and - 'dot(p)' functions (dot(p) is the vector inner product b*p ); moreover, - scalar multiplication and sum operations are available. - out : psydac.linalg.basic.Vector | NoneType - The output vector, or None (optional). - - Returns - ------- - x : psydac.linalg.basic.Vector - Numerical solution of linear system. To check the convergence of the solver, - use the method InverseLinearOperator.get_info(). - - info : dict - Dictionary containing convergence information: - - 'niter' = (int) number of iterations - - 'success' = (boolean) whether convergence criteria have been met - - 'res_norm' = (float) 2-norm of residual vector r = A*x - b. - - References - ---------- - [1] H. A. van der Vorst. Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the - solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comp., 13(2):631–644, 1992. - """ - - A = self._A - domain = self._domain - codomain = self._codomain - options = self._options - x0 = options["x0"] - tol = options["tol"] - maxiter = options["maxiter"] - verbose = options["verbose"] - recycle = options["recycle"] - - assert isinstance(b, Vector) - assert b.space is domain - - # First guess of solution - if out is not None: - assert isinstance(out, Vector) - assert out.space is codomain - - x = x0.copy(out=out) - - # Extract local storage - v = self._tmps["v"] - r = self._tmps["r"] - p = self._tmps["p"] - vr = self._tmps["vr"] - r0 = self._tmps["r0"] - - # First values - A.dot(x, out=v) - b.copy(out=r) - r -= v - #r = b - A.dot(x) - r.copy(out=p) - v *= 0.0 - vr *= 0.0 - - r.copy(out=r0) - - res_sqr = r.dot(r).real - tol_sqr = tol ** 2 - - if verbose: - print("BiCGSTAB solver:") - print("+---------+---------------------+") - print("+ Iter. # | L2-norm of residual |") - print("+---------+---------------------+") - template = "| {:7d} | {:19.2e} |" - - # Iterate to convergence - for m in range(1, maxiter + 1): - - if res_sqr < tol_sqr: - m -= 1 - break - - # ----------------------- - # MATRIX-VECTOR PRODUCTS - # ----------------------- - v = A.dot(p, out=v) - # ----------------------- - - # c := (r0, r) - c = r0.dot(r) - - # a := (r0, r) / (r0, v) - a = c / (r0.dot(v)) - - # r := r - a*v - r.mul_iadd(-a, v) - - # vr := A*r - vr = A.dot(r, out=vr) - - # w := (r, A*r) / (A*r, A*r) - w = r.dot(vr) / vr.dot(vr) - - # ----------------------- - # SOLUTION UPDATE - # ----------------------- - # x := x + a*p +w*r - x.mul_iadd(a, p) - x.mul_iadd(w, r) - # ----------------------- - - # r := r - w*A*r - r.mul_iadd(-w, vr) - - # ||r||_2 := (r, r) - res_sqr = r.dot(r).real - - if res_sqr < tol_sqr: - break - - # b := a / w * (r0, r)_{m+1} / (r0, r)_m - b = r0.dot(r) * a / (c * w) - - # p := r + b*p- b*w*v - p *= b - p += r - p.mul_iadd(-b * w, v) - - if verbose: - print(template.format(m, sqrt(res_sqr))) - - if verbose: - print("+---------+---------------------+") - - # Convergence information - self._info = {'niter': m, 'success': res_sqr < tol_sqr, 'res_norm': sqrt(res_sqr)} - - if recycle: - x.copy(out=self._options["x0"]) - - return x - - def dot(self, b, out=None): - return self.solve(b, out=out) - #=============================================================================== class MinimumResidual(InverseLinearOperator): """ @@ -1237,52 +1065,16 @@ class MinimumResidual(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - assert A.domain.dtype == float - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() + self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose} - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'minres' - self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")} + super().__init__(A, **self._options) + + self._tmps = {key: self.domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")} self._info = None - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) + @property + def solver(self): + return 'minres' def solve(self, b, out=None): """ @@ -1569,67 +1361,34 @@ class LSMR(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, conlim=1e8, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain - - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() - - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'lsmr' self._options = {"x0":x0, "tol":tol, "atol":atol, "btol":btol, "maxiter":maxiter, "conlim":conlim, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) + + super().__init__(A, **self._options) + + # check additional options + if atol is not None: + assert is_real(atol), "atol must be a real number" + assert atol >= 0, "atol must not be negative" + if btol is not None: + assert is_real(btol), "btol must be a real number" + assert btol >= 0, "btol must not be negative" + assert is_real(conlim), "conlim must be a real number" # actually an integer? + assert conlim > 0, "conlim must be positive" # supposedly + self._info = None self._successful = None - tmps_domain = {key: domain.zeros() for key in ("u", "u_work")} - tmps_codomain = {key: codomain.zeros() for key in ("v", "v_work", "h", "hbar")} + tmps_domain = {key: self.domain.zeros() for key in ("u", "u_work")} + tmps_codomain = {key: self.codomain.zeros() for key in ("v", "v_work", "h", "hbar")} self._tmps = {**tmps_codomain, **tmps_domain} + @property + def solver(self): + return 'lsmr' + def get_success(self): return self._successful - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - if value is not None: - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" # suppose atol/btol must also be positive numbers - elif key == 'atol' or key == 'btol': - if value is not None: - assert is_real(value), "atol/btol must be a real number" - assert value >= 0, "atol/btol must not be negative" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'conlim': - assert is_real(value), "conlim must be a real number" # actually an integer? - assert value > 0, "conlim must be positive" # supposedly - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) - def solve(self, b, out=None): """Iterative solver for least-squares problems. lsmr solves the system of linear equations ``Ax = b``. If the system @@ -1956,55 +1715,20 @@ class GMRES(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=100, verbose=False, recycle=False): - assert isinstance(A, LinearOperator) - assert A.domain.dimension == A.codomain.dimension - domain = A.codomain - codomain = A.domain + self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - if x0 is not None: - assert isinstance(x0, Vector) - assert x0.space is codomain - else: - x0 = codomain.zeros() + super().__init__(A, **self._options) - self._A = A - self._domain = domain - self._codomain = codomain - self._solver = 'gmres' - self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} - self._check_options(**self._options) - self._tmps = {key: domain.zeros() for key in ("r", "p")} + self._tmps = {key: self.domain.zeros() for key in ("r", "p")} # Initialize upper Hessenberg matrix self._H = np.zeros((self._options["maxiter"] + 1, self._options["maxiter"]), dtype=A.domain.dtype) self._Q = [] self._info = None - def _check_options(self, **kwargs): - for key, value in kwargs.items(): - - if key == 'x0': - if value is not None: - assert isinstance(value, Vector), "x0 must be a Vector or None" - assert value.space == self._codomain, "x0 belongs to the wrong VectorSpace" - elif key == 'tol': - assert is_real(value), "tol must be a real number" - assert value > 0, "tol must be positive" - elif key == 'maxiter': - assert isinstance(value, int), "maxiter must be an int" - assert value > 0, "maxiter must be positive" - elif key == 'verbose': - assert isinstance(value, bool), "verbose must be a bool" - elif key == 'recycle': - assert isinstance(value, bool), "recycle must be a bool" - else: - raise ValueError(f"Key '{key}' not understood. See self._options for allowed keys.") - - def transpose(self, conjugate=False): - At = self._A.transpose(conjugate=conjugate) - solver = self._solver - options = self._options - return inverse(At, solver, **options) + @property + def solver(self): + return 'gmres' def solve(self, b, out=None): """ From 81de0d7d274b086c5f43c3c66517676563371430 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Dec 2023 15:52:11 +0100 Subject: [PATCH 03/34] Write new diagonal method for class StencilMatrix --- psydac/linalg/stencil.py | 76 +++++++++++++++++++++++++++++----------- 1 file changed, 56 insertions(+), 20 deletions(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index cf4fcddc1..9499f4bf3 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1368,30 +1368,66 @@ def _exchange_assembly_data_serial(self): # Copy data from left to right idx_to = tuple( idx_front + [slice( m*p, m*p+p)] + idx_back ) - idx_from = tuple( idx_front + [ slice(-m*p,-m*p+p) if (-m*p+p)!=0 else slice(-m*p,None)] + idx_back ) + idx_from = tuple( idx_front + [slice(-m*p,-m*p+p) if (-m*p+p)!=0 else slice(-m*p,None)] + idx_back ) self._data[idx_to] += self._data[idx_from] - def diagonal(self): - if self._diag_indices is None: - cm = self.codomain.shifts - dm = self.domain.shifts - ss = self.codomain.starts - pp = [compute_diag_len(p,mj,mi)-(p+1) for p,mi,mj in zip(self._pads, cm, dm)] - nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) - indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2*len(nrows))] - l = 0 - for xx in np.ndindex(*nrows): - ii = [m*p+x for m,p,x in zip(self.domain.shifts, self.domain.pads, xx)] - jj = [p+x+s-((x+s)//mi)*mj for (x,mi,mj,p,s) in zip(xx,cm,dm,pp,ss)] - for k in range(len(nrows)): - indices[k][l] = ii[k] - indices[k+len(nrows)][l] = jj[k] - l += 1 - self._diag_indices = tuple(indices) + # ... + def diagonal(self, *, inverse = False, out = None): + """Get the coefficients on the main diagonal as another StencilMatrix object. + + Parameters + ---------- + inverse : bool + If True, get the inverse of the diagonal. (Default: False). + + out : StencilMatrix + If provided, write the diagonal entries into this matrix. (Default: None). + + Returns + ------- + StencilMatrix + The matrix which contains the main diagonal of self (or its inverse). + + """ + + V, W = self.domain, self.codomain + if inverse: + V, W = W, V + + if out is not None: + assert isinstance(out, StencilMatrix) + assert out.domain is V + assert out.codomain is W + assert all(p == 0 for p in out.pads) # is this really needed? else: - nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) + out = StencilMatrix(V, W, pads = [0] * V.ndim) - return self._data[self._diag_indices].reshape(nrows) + index = tuple(slice(s, e) for s, e in zip(W.starts, W.ends)) + (0,) * V.ndim + out[index] = 1.0 / M[index] if inverse else M[index] + return out + + +# def diagonal(self): +# if self._diag_indices is None: +# cm = self.codomain.shifts +# dm = self.domain.shifts +# ss = self.codomain.starts +# pp = [compute_diag_len(p,mj,mi)-(p+1) for p,mi,mj in zip(self._pads, cm, dm)] +# nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) +# indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2*len(nrows))] +# l = 0 +# for xx in np.ndindex(*nrows): +# ii = [m*p+x for m,p,x in zip(self.domain.shifts, self.domain.pads, xx)] +# jj = [p+x+s-((x+s)//mi)*mj for (x,mi,mj,p,s) in zip(xx,cm,dm,pp,ss)] +# for k in range(len(nrows)): +# indices[k][l] = ii[k] +# indices[k+len(nrows)][l] = jj[k] +# l += 1 +# self._diag_indices = tuple(indices) +# else: +# nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) +# +# return self._data[self._diag_indices].reshape(nrows) # ... def topetsc(self): From 3e36e33639c4387264882428837c67f4e5c5ff9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Dec 2023 15:53:13 +0100 Subject: [PATCH 04/34] Add diagonal method to class BlockLinearOperator --- psydac/linalg/block.py | 41 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index eed6c7906..94640d719 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -777,6 +777,42 @@ def __sub__(self, M): #-------------------------------------- # New properties/methods #-------------------------------------- + def diagonal(self, *, inverse = False, out = None): + """Get the coefficients on the main diagonal as another BlockLinearOperator object. + + Parameters + ---------- + inverse : bool + If True, get the inverse of the diagonal. (Default: False). + + out : BlockLinearOperator + If provided, write the diagonal entries into this matrix. (Default: None). + + Returns + ------- + BlockLinearOperator + The matrix which contains the main diagonal of self (or its inverse). + + """ + V, W = self.domain, self.codomain + if inverse: + V, W = W, V + + if out is not None: + assert isinstance(out, BlockLinearOperator) + assert out.domain is V + assert out.codomain is W + assert all(i==j for i, j in out._blocks.keys()) # is this really needed? + for i in range(min(self.n_block_rows, self.n_block_cols)): + self[i, i].diagonal(inverse = inverse, out = out[i, i]) + else: + out = BlockLinearOperator(V, W) + for i in range(min(self.n_block_rows, self.n_block_cols)): + out[i, i] = self[i, i].diagonal(inverse = inverse) + + return out + + # ... @property def blocks(self): """ Immutable 2D view (tuple of tuples) of the linear operator, @@ -868,7 +904,8 @@ def __setitem__(self, key, value): assert value.codomain is self.codomain[i] self._blocks[i,j] = value - + + # ... def transform(self, operation): """ Applies an operation on each block in this BlockLinearOperator. @@ -881,6 +918,7 @@ def transform(self, operation): blocks = {ij: operation(Bij) for ij, Bij in self._blocks.items()} return BlockLinearOperator(self.domain, self.codomain, blocks=blocks) + # ... def backend(self): return self._backend @@ -909,7 +947,6 @@ def copy(self, out=None): out.set_backend(self._backend) - # ... def __imul__(self, a): for Bij in self._blocks.values(): From f52a291fc9e904f2a845a82e80aa73c463b7e4f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Dec 2023 15:57:54 +0100 Subject: [PATCH 05/34] Use inverse diagonal as preconditioner for P-CG and P-BiCGStab in tests --- psydac/linalg/tests/test_solvers.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 793b9e89c..78863c70e 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -59,10 +59,9 @@ def define_data(n, p, matrix_data, dtype=float): @pytest.mark.parametrize( 'n', [5, 10, 13] ) @pytest.mark.parametrize('p', [2, 3]) @pytest.mark.parametrize('dtype', [float, complex]) -@pytest.mark.parametrize('pc', [None]) @pytest.mark.parametrize('solver', ['cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres']) -def test_solver_tridiagonal(n, p, dtype, pc, solver, verbose=False): +def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): #--------------------------------------------------------------------------- # PARAMETERS @@ -108,9 +107,10 @@ def test_solver_tridiagonal(n, p, dtype, pc, solver, verbose=False): #Create the solvers if solver in ['pcg', 'pbicgstab']: - solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=True) + pc = A.diagonal(inverse=True) + solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=True) else: - solv = inverse(A, solver, tol=1e-13, verbose=True) + solv = inverse(A, solver, tol=1e-13, verbose=True) solvt = solv.transpose() solvh = solv.H solv2 = inverse(A@A, solver, tol=1e-13, verbose=True, recycle=True) # Test solver of composition of operators From bf198ef839d611e481881faf13b2ce832bbb54b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 15 Dec 2023 07:50:05 +0100 Subject: [PATCH 06/34] Fix typos --- psydac/linalg/stencil.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 9499f4bf3..b0332e691 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1402,8 +1402,8 @@ def diagonal(self, *, inverse = False, out = None): else: out = StencilMatrix(V, W, pads = [0] * V.ndim) - index = tuple(slice(s, e) for s, e in zip(W.starts, W.ends)) + (0,) * V.ndim - out[index] = 1.0 / M[index] if inverse else M[index] + index = tuple(slice(s, e+1) for s, e in zip(W.starts, W.ends)) + (0,) * V.ndim + out[index] = 1.0 / self[index] if inverse else self[index] return out From 1fafaabf8e031cf1cff7d6f8ff71c67e4d0f96ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 15 Dec 2023 07:57:21 +0100 Subject: [PATCH 07/34] Improve efficiency --- psydac/linalg/stencil.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index b0332e691..d3c3afef1 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1403,7 +1403,11 @@ def diagonal(self, *, inverse = False, out = None): out = StencilMatrix(V, W, pads = [0] * V.ndim) index = tuple(slice(s, e+1) for s, e in zip(W.starts, W.ends)) + (0,) * V.ndim - out[index] = 1.0 / self[index] if inverse else self[index] + if inverse: + np.divide(1, self[index], out=out[index]) + else: + out[index] = self[index] + return out From c55d6b6267991f2c1d35d0103f9927a4723638a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 15 Dec 2023 08:06:18 +0100 Subject: [PATCH 08/34] Fix some tests --- psydac/api/tests/test_api_feec_1d.py | 3 ++- psydac/api/tests/test_api_feec_2d.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index 260bb0381..3c06204c7 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -336,7 +336,8 @@ def discrete_energies(self, e, b): step_ampere_1d = dt * ( M0_dir_inv @ D0_T_dir @ M1 ) elif bc_mode == 'penalization': - M0_M0_bc_inv = inverse(M0+M0_bc, 'pcg', pc='jacobi', **kwargs) + M0_M0_bc = M0 + M0_bc + M0_M0_bc_inv = inverse(M0_M0_bc, 'pcg', pc = M0_M0_bc.diagonal(inverse=True), **kwargs) step_ampere_1d = dt * ( M0_M0_bc_inv @ D0_T @ M1 ) half_step_faraday_1d = (dt/2) * D0 diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 5c27b678e..1453211f1 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -519,7 +519,8 @@ def discrete_energies(self, e, b): M1_inv = inverse(M1, 'cg', **kwargs) step_ampere_2d = dt * (M1_inv @ D1_T @ M2) else: - M1_M1_bc_inv = inverse(M1 + M1_bc, 'pcg', pc='jacobi', **kwargs) + M1_M1_bc = M1 + M1_bc + M1_M1_bc_inv = inverse(M1_M1_bc, 'pcg', pc = M1_M1_bc.diagonal(inverse=True), **kwargs) step_ampere_2d = dt * (M1_M1_bc_inv @ D1_T @ M2) half_step_faraday_2d = (dt/2) * D1 From 68b51ac8f6e583d23e30251a5eb479ae3baf2727 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 15 Dec 2023 12:20:08 +0100 Subject: [PATCH 09/34] Fix some tests --- psydac/linalg/tests/test_linalg.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index d929e42cf..b3b9614a1 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -623,6 +623,7 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): tol = 1e-5 C = inverse(B, 'cg', tol=tol) + pc = B.T.diagonal(inverse=True) # -1,T & T,-1 -> equal assert isinstance(C.T, ConjugateGradient) @@ -634,8 +635,8 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): assert isinstance(inverse(C.T, 'bicg'), BlockLinearOperator) assert np.array_equal(inverse(C.T, 'bicg').dot(u).toarray(), B.T.dot(u).toarray()) # T,-1,-1 -> equal T - assert isinstance(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc='jacobi'), BlockLinearOperator) - assert np.array_equal(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc='jacobi').dot(u).toarray(), B.T.dot(u).toarray()) + assert isinstance(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc=pc), BlockLinearOperator) + assert np.array_equal(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc=pc).dot(u).toarray(), B.T.dot(u).toarray()) # T,-1,T -> equal -1 assert isinstance(inverse(B.T, 'cg', tol=tol).T, ConjugateGradient) assert np.sqrt(sum(((inverse(B.T, 'cg', tol=tol).dot(u) - C.dot(u)).toarray())**2)) < tol @@ -647,6 +648,7 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): tol = 1e-5 C = inverse(S, 'cg', tol=tol) + pc = S.T.diagonal(inverse=True) # -1,T & T,-1 -> equal assert isinstance(C.T, ConjugateGradient) @@ -658,8 +660,8 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): assert isinstance(inverse(C.T, 'bicg'), StencilMatrix) assert np.array_equal(inverse(C.T, 'bicg').dot(v).toarray(), S.T.dot(v).toarray()) # T,-1,-1 -> equal T - assert isinstance(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc='jacobi'), StencilMatrix) - assert np.array_equal(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc='jacobi').dot(v).toarray(), S.T.dot(v).toarray()) + assert isinstance(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc=pc), StencilMatrix) + assert np.array_equal(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc=pc).dot(v).toarray(), S.T.dot(v).toarray()) # T,-1,T -> equal -1 assert isinstance(inverse(S.T, 'cg', tol=tol).T, ConjugateGradient) assert np.sqrt(sum(((inverse(S.T, 'cg', tol=tol).dot(v) - C.dot(v)).toarray())**2)) < tol @@ -833,8 +835,8 @@ def test_operator_evaluation(n1, n2, p1, p2): S_cg = inverse(S, 'cg', tol=tol) B_cg = inverse(B, 'cg', tol=tol) - S_pcg = inverse(S, 'pcg', pc='jacobi', tol=tol) - B_pcg = inverse(B, 'pcg', pc='jacobi', tol=tol) + S_pcg = inverse(S, 'pcg', pc=S.diagonal(inverse=True), tol=tol) + B_pcg = inverse(B, 'pcg', pc=B.diagonal(inverse=True), tol=tol) S_bicg = inverse(S, 'bicg', tol=tol) B_bicg = inverse(B, 'bicg', tol=tol) S_lsmr = inverse(S, 'lsmr', tol=tol) @@ -944,7 +946,7 @@ def test_x0update(solver): # Create Inverse tol = 1e-6 if solver == 'pcg': - A_inv = inverse(A, solver, pc='jacobi', tol=tol) + A_inv = inverse(A, solver, pc=A.diagonal(inverse=True), tol=tol) else: A_inv = inverse(A, solver, tol=tol) options = A_inv.options From 885f0b3edb71404e81df618d40a105521bed4597 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 8 Jan 2024 12:45:42 +0100 Subject: [PATCH 10/34] Pass correct preconditioner in test_2d_multipatch_mapping_maxwell.py --- psydac/api/tests/test_2d_multipatch_mapping_maxwell.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py index 90dbc6b0a..61b3fa5d5 100644 --- a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py +++ b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py @@ -93,7 +93,8 @@ def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename= equation_h = discretize(equation, domain_h, [Vh, Vh], backend=PSYDAC_BACKEND_GPYCCEL) l2norm_h = discretize(l2norm, domain_h, Vh, backend=PSYDAC_BACKEND_GPYCCEL) - equation_h.set_solver('pcg', pc='jacobi', tol=1e-8) + jacobi_pc = equation_h.lhs.diagonal(inverse=True) + equation_h.set_solver('pcg', pc=jacobi_pc, tol=1e-8) uh = equation_h.solve() From 10fb4c764b75bf5bdd23eb4cd866ecc67705377a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 8 Jan 2024 13:05:01 +0100 Subject: [PATCH 11/34] Fix typo in test_2d_multipatch_mapping_maxwell.py --- psydac/api/tests/test_2d_multipatch_mapping_maxwell.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py index 61b3fa5d5..8301d0a4b 100644 --- a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py +++ b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py @@ -93,7 +93,7 @@ def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename= equation_h = discretize(equation, domain_h, [Vh, Vh], backend=PSYDAC_BACKEND_GPYCCEL) l2norm_h = discretize(l2norm, domain_h, Vh, backend=PSYDAC_BACKEND_GPYCCEL) - jacobi_pc = equation_h.lhs.diagonal(inverse=True) + jacobi_pc = equation_h.linear_system.lhs.diagonal(inverse=True) equation_h.set_solver('pcg', pc=jacobi_pc, tol=1e-8) uh = equation_h.solve() From 9089039e2758457b976d4d628b4bbac69e2d6ceb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 8 Jan 2024 19:37:46 +0100 Subject: [PATCH 12/34] Assemble a DiscreteEquation before accessing the linear_system attribute --- psydac/api/tests/test_2d_multipatch_mapping_maxwell.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py index 8301d0a4b..a9759ba86 100644 --- a/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py +++ b/psydac/api/tests/test_2d_multipatch_mapping_maxwell.py @@ -93,11 +93,19 @@ def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename= equation_h = discretize(equation, domain_h, [Vh, Vh], backend=PSYDAC_BACKEND_GPYCCEL) l2norm_h = discretize(l2norm, domain_h, Vh, backend=PSYDAC_BACKEND_GPYCCEL) + # Explicitly assemble the linear system + equation_h.assemble() + + # Extract the matrix' diagonal to build a Jacobi preconditioner jacobi_pc = equation_h.linear_system.lhs.diagonal(inverse=True) + + # Choose a linear solver and pass any flags to it equation_h.set_solver('pcg', pc=jacobi_pc, tol=1e-8) + # Solve the linear system and obtain the solution as a FEM field uh = equation_h.solve() + # Compute the L2 norm of the error l2_error = l2norm_h.assemble(F=uh) return l2_error, uh From 11aba5e32da4394a8924fa40b9ef1d65de2abfe8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 8 Jan 2024 19:39:14 +0100 Subject: [PATCH 13/34] Pass correct preconditioner in test_api_2d_compatible_spaces --- psydac/api/tests/test_api_2d_compatible_spaces.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psydac/api/tests/test_api_2d_compatible_spaces.py b/psydac/api/tests/test_api_2d_compatible_spaces.py index 1638af024..e84e3782e 100644 --- a/psydac/api/tests/test_api_2d_compatible_spaces.py +++ b/psydac/api/tests/test_api_2d_compatible_spaces.py @@ -399,7 +399,8 @@ def run_maxwell_time_harmonic_2d_dir(uex, f, alpha, ncells, degree): #+++++++++++++++++++++++++++++++ # Solve linear system - M_inv = inverse(M, 'pcg', pc='jacobi', tol=1e-8) + jacobi_pc = M.diagonal(inverse=True) + M_inv = inverse(M, 'pcg', pc=jacobi_pc, tol=1e-8) sol = M_inv @ b uh = FemField( Vh, sol ) From 15c028d4e8d61543ccb41ffca105d5e368256c64 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 10 Jan 2024 15:00:58 +0100 Subject: [PATCH 14/34] Fix wrong parameters to inverse() in test_solvers.py --- psydac/linalg/tests/test_solvers.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/psydac/linalg/tests/test_solvers.py b/psydac/linalg/tests/test_solvers.py index 78863c70e..6d4b37e73 100644 --- a/psydac/linalg/tests/test_solvers.py +++ b/psydac/linalg/tests/test_solvers.py @@ -108,12 +108,12 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): #Create the solvers if solver in ['pcg', 'pbicgstab']: pc = A.diagonal(inverse=True) - solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=True) + solv = inverse(A, solver, pc=pc, tol=1e-13, verbose=verbose, recycle=True) else: - solv = inverse(A, solver, tol=1e-13, verbose=True) + solv = inverse(A, solver, tol=1e-13, verbose=verbose, recycle=True) solvt = solv.transpose() solvh = solv.H - solv2 = inverse(A@A, solver, tol=1e-13, verbose=True, recycle=True) # Test solver of composition of operators + solv2 = inverse(A@A, solver, tol=1e-13, verbose=verbose, recycle=True) # Test solver of composition of operators # Manufacture right-hand-side vector from exact solution be = A @ xe @@ -160,7 +160,6 @@ def test_solver_tridiagonal(n, p, dtype, solver, verbose=False): if solver != 'pcg': bc = A @ A @ xc - err = b - be err_norm = np.linalg.norm( err.toarray() ) err2 = b2 - be2 From 9cf15012257e8e391c233c7cbd33582b06a6ced5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 11 Jan 2024 16:25:15 +0100 Subject: [PATCH 15/34] Remove obsolete static method jacobi from InverseLinearOperator --- psydac/linalg/basic.py | 56 ------------------------------------------ 1 file changed, 56 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index f8a137053..e1ebb724f 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -986,62 +986,6 @@ def transpose(self, conjugate=False): options = self._options return inverse(At, solver, **options) - @staticmethod - def jacobi(A, b, out=None): - """ - Jacobi preconditioner. - - A : psydac.linalg.stencil.StencilMatrix | psydac.linalg.block.BlockLinearOperator - Left-hand-side matrix A of linear system. - - b : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector - Right-hand-side vector of linear system. - - Returns - ------- - x : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector - Preconditioner solution - - """ - from psydac.linalg.block import BlockLinearOperator, BlockVector - from psydac.linalg.stencil import StencilMatrix, StencilVector - - # In case A is None we return a zero vector - if A is None: - return b.space.zeros() - - # Sanity checks - assert isinstance(A, (StencilMatrix, BlockLinearOperator)) - assert isinstance(b, (StencilVector, BlockVector)) - assert A.codomain.dimension == A.domain.dimension - assert A.codomain == b.space - - #------------------------------------------------------------- - # Handle the case of a block linear system - if isinstance(A, BlockLinearOperator): - if out is not None: - for i, bi in enumerate(b.blocks): - InverseLinearOperator.jacobi(A[i,i], bi, out=out[i]) - return out - else: - x = [InverseLinearOperator.jacobi(A[i, i], bi) for i, bi in enumerate(b.blocks)] - y = BlockVector(b.space, blocks=x) - return y - #------------------------------------------------------------- - - V = b.space - i = tuple(slice(s, e + 1) for s, e in zip(V.starts, V.ends)) - - if out is not None: - b.copy(out=out) - out[i] /= A.diagonal() - out.update_ghost_regions() - else: - out = b.copy() - out[i] /= A.diagonal() - out.update_ghost_regions() - return out - #=============================================================================== class LinearSolver(ABC): """ From f3e8807549643ae4946f92b4b8a57939f27a08a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 11 Jan 2024 16:47:15 +0100 Subject: [PATCH 16/34] Clean up inverse() factory function --- psydac/linalg/solvers.py | 58 ++++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 23 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 6f21761d4..47404b96b 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -9,12 +9,23 @@ from psydac.linalg.basic import Vector, LinearOperator, InverseLinearOperator, IdentityOperator, ScaledLinearOperator from psydac.linalg.utilities import _sym_ortho -__all__ = ('ConjugateGradient', 'PConjugateGradient', 'BiConjugateGradient', 'BiConjugateGradientStabilized', 'MinimumResidual', 'LSMR', 'GMRES') +__all__ = ( + 'ConjugateGradient', + 'PConjugateGradient', + 'BiConjugateGradient', + 'BiConjugateGradientStabilized', + 'PBiConjugateGradientStabilized', + 'MinimumResidual', + 'LSMR', + 'GMRES' +) +#=============================================================================== def is_real(x): from numbers import Number return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool) +#=============================================================================== def inverse(A, solver, **kwargs): """ A function to create objects of all InverseLinearOperator subclasses. @@ -34,21 +45,35 @@ def inverse(A, solver, **kwargs): solver : str Preferred iterative solver. Options are: 'cg', 'pcg', 'bicg', - 'bicgstab', 'minres', 'lsmr', 'gmres'. + 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres'. Returns ------- obj : psydac.linalg.basic.InverseLinearOperator - More specifically: Returns the chosen subclass, for example psydac.linalg.solvers.ConjugateGradient - A linear operator acting as the inverse of A. + A linear operator acting as the inverse of A, of the chosen subclass + (for example psydac.linalg.solvers.ConjugateGradient). """ + + # Map each possible value of the `solver` string with a specific + # `InverseLinearOperator` subclass in this module: + solvers_dict = { + 'cg' : ConjugateGradient, + 'pcg' : PConjugateGradient, + 'bicg' : BiConjugateGradient, + 'bicgstab' : BiConjugateGradientStabilized, + 'pbicgstab': PBiConjugateGradientStabilized, + 'minres' : MinimumResidual, + 'lsmr' : LSMR, + 'gmres' : GMRES, + } + # Check solver input - solvers = ('cg', 'pcg', 'bicg', 'bicgstab', 'pbicgstab', 'minres', 'lsmr', 'gmres') - if solver not in solvers: + if solver not in solvers_dict: raise ValueError(f"Required solver '{solver}' not understood.") assert isinstance(A, LinearOperator) + if isinstance(A, IdentityOperator): return A elif isinstance(A, ScaledLinearOperator): @@ -57,22 +82,9 @@ def inverse(A, solver, **kwargs): return A.linop # Instantiate object of correct solver class - if solver == 'cg': - obj = ConjugateGradient(A, **kwargs) - elif solver == 'pcg': - obj = PConjugateGradient(A, **kwargs) - elif solver == 'bicg': - obj = BiConjugateGradient(A, **kwargs) - elif solver == 'bicgstab': - obj = BiConjugateGradientStabilized(A, **kwargs) - elif solver == 'pbicgstab': - obj = PBiConjugateGradientStabilized(A, **kwargs) - elif solver == 'minres': - obj = MinimumResidual(A, **kwargs) - elif solver == 'lsmr': - obj = LSMR(A, **kwargs) - elif solver == 'gmres': - obj = GMRES(A, **kwargs) + cls = solvers_dict[solver] + obj = cls(A, **kwargs) + return obj #=============================================================================== @@ -1897,4 +1909,4 @@ def apply_givens_rotation(self, k, sn, cn): def dot(self, b, out=None): return self.solve(b, out=out) - \ No newline at end of file + From 9b2e24027eb72c16524adcdf06bfcb73a3587204 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 11 Jan 2024 16:55:38 +0100 Subject: [PATCH 17/34] Add recycle option to PBiConjugateGradientStabilized --- psydac/linalg/solvers.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 47404b96b..847c454b9 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -827,9 +827,9 @@ class PBiConjugateGradientStabilized(InverseLinearOperator): ---------- [1] A. Maister, Numerik linearer Gleichungssysteme, Springer ed. 2015. """ - def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False): + def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - self._options = { "pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose} + self._options = {"pc": pc, "x0": x0, "tol": tol, "maxiter": maxiter, "verbose": verbose, "recycle": recycle} super().__init__(A, **self._options) @@ -848,9 +848,6 @@ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False def solver(self): return 'pbicgstab' - def _update_options( self ): - self._options = {"pc": self._pc, "x0":self._x0, "tol":self._tol, "maxiter": self._maxiter, "verbose": self._verbose} - def solve(self, b, out=None): """ Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. @@ -889,6 +886,7 @@ def solve(self, b, out=None): tol = options["tol"] maxiter = options["maxiter"] verbose = options["verbose"] + recycle = options["recycle"] assert isinstance(b, Vector) assert b.space is domain @@ -1024,6 +1022,9 @@ def solve(self, b, out=None): self._info = {'niter': niter, 'success': res_sqr < tol_sqr, 'res_norm': sqrt(res_sqr)} + if recycle: + x.copy(out=self._options["x0"]) + return x def dot(self, b, out=None): From 7b947d9230f66f86f494c1d48eb5b4b7a6c36c5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 11 Jan 2024 21:45:26 +0100 Subject: [PATCH 18/34] Add recycle option to MinimumResidual --- psydac/linalg/solvers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 847c454b9..1f64f8877 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -1078,7 +1078,7 @@ class MinimumResidual(InverseLinearOperator): """ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): - self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose} + self._options = {"x0":x0, "tol":tol, "maxiter":maxiter, "verbose":verbose, "recycle":recycle} super().__init__(A, **self._options) From 1867ba176443763fca7f6dfb6855a2407e36f701 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 22 Jan 2024 18:59:29 +0100 Subject: [PATCH 19/34] Introduce class StencilDiagonalMatrix --- psydac/linalg/stencil.py | 182 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 181 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index d3c3afef1..43c1c0d93 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -8,7 +8,7 @@ import numpy as np from types import MappingProxyType -from scipy.sparse import coo_matrix +from scipy.sparse import coo_matrix, diags as sp_diags from mpi4py import MPI from psydac.linalg.basic import VectorSpace, Vector, LinearOperator @@ -1853,6 +1853,186 @@ def set_backend(self, backend): self._args.pop('ndiags') self._func = dot.func +#=============================================================================== +class StencilDiagonalMatrix(LinearOperator): + """ + Linear operator which operates between stencil vector spaces, and which can + be represented by a matrix with non-zero entries only on its main diagonal. + As such this operator is completely local and requires no data communication. + + We assume that the vectors in the domain and the codomain have the same + shape and are distributed in the same way. + + Parameters + ---------- + V : psydac.linalg.stencil.StencilVectorSpace + Domain of the new linear operator. + + W : psydac.linalg.stencil.StencilVectorSpace + Codomain of the new linear operator. + + """ + def __init__(self, V, W, data): + + # Check domain and codomain + assert isinstance(V, StencilVectorSpace) + assert isinstance(W, StencilVectorSpace) + assert V.starts == W.starts + assert V.ends == W.ends + + data = np.asarray(data) + + # Check shape of provided data + shape = tuple(e - s + 1 for s, e in zip(V.starts, V.ends)) + assert data.shape == shape + + # Store info in object + self._domain = V + self._codomain = W + self._data = data + + #-------------------------------------- + # Abstract interface + #-------------------------------------- + @property + def domain(self): + return self._domain + + @property + def codomain(self): + return self._codomain + + @property + def dtype(self): + return self._data.dtype + + def tosparse(self): + return sp_diags(self._data.ravel()) + + def toarray(self): + return self._data.copy() + + def dot(self, v, out=None): + + assert isinstance(v, StencilVector) + assert v.space is self.domain + + if out is not None: + assert isinstance(out, StencilVector) + assert out.space is self.codomain + else: + out = self.codomain.zeros() + + V = self.domain + i = tuple(slice(s, e + 1) for s, e in zip(V.starts, V.ends)) + np.multiply(self._data, v[i], out=out[i]) + + out.ghost_regions_in_sync = False + + return out + + # ... + # TODO [YG 22.01.2024]: idot function will require a dedicated kernel + # ... + + def transpose(self, *, conjugate=False, out=None): + + assert isinstance(conjugate, bool) + + if out is not None: + assert isinstance(out, StencilDiagonalMatrix) + assert out.domain is self.codomain + assert out.codomain is self.domain + + if not (conjugate and self.dtype is complex): + + if out is None: + data = self._data.copy() + else: + np.copyto(out._data, self._data, casting='no') + + else: + + if out is None: + data = np.conjugate(self._data, casting='no') + else: + np.conjugate(self._data, out=out._data, casting='no') + + if out is None: + out = StencilDiagonalMatrix(self.codomain, self.domain, data) + + return out + + #-------------------------------------- + # Other properties/methods + #-------------------------------------- + def copy(self, *, out=None): + + if out is self: + return self + + if out is None: + data = self._data.copy() + out = StencilDiagonalMatrix(self.domain, self.codomain, data) + else: + assert isinstance(out, StencilDiagonalMatrix) + assert out.domain is self.domain + assert out.codomain is self.codomain + np.copyto(out._data, self._data, casting='no') + + return out + + def diagonal(self, *, inverse = False, out = None): + """ + Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. + + In the default case (inverse=False, out=None) self is returned. + + Parameters + ---------- + inverse : bool + If True, get the inverse of the diagonal. (Default: False). + + out : StencilDiagonalMatrix + If provided, write the diagonal entries into this matrix. (Default: None). + + Returns + ------- + StencilDiagonalMatrix + Either self, or another StencilDiagonalMatrix with the diagonal inverse. + + """ + # Check `inverse` argument + assert isinstance(inverse, bool) + + # Determine domain and codomain of the `out` matrix + V, W = self.domain, self.codomain + if inverse: + V, W = W, V + + # Check `out` argument and identify `data` array of output vector + if out is None: + data = None + else: + assert isinstance(out, StencilDiagonalMatrix) + assert out.domain is V + assert out.codomain is W + data = out._data + + # Calculate entries, or set `out=self` in default case + if inverse: + data = np.divide(1, diag, out=data) + elif out: + np.copyto(data, diag) + else: + out = self + + # If needed create a new StencilDiagonalMatrix object + if out is None: + out = StencilDiagonalMatrix(V, W, data) + + return out + #=============================================================================== # TODO [YG, 28.01.2021]: # - Check if StencilMatrix should be subclassed From 7bac5cc8cbee8240ae151f10b6cb046a7a1029ca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 22 Jan 2024 19:01:08 +0100 Subject: [PATCH 20/34] Make StencilMatrix.diagonal() return a StencilDiagonalMatrix object --- psydac/linalg/stencil.py | 74 ++++++++++++++++++++++++++++++++++------ 1 file changed, 63 insertions(+), 11 deletions(-) diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 43c1c0d93..8c3cfd824 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1373,44 +1373,58 @@ def _exchange_assembly_data_serial(self): # ... def diagonal(self, *, inverse = False, out = None): - """Get the coefficients on the main diagonal as another StencilMatrix object. + """ + Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. Parameters ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). - out : StencilMatrix + out : StencilDiagonalMatrix If provided, write the diagonal entries into this matrix. (Default: None). Returns ------- - StencilMatrix + StencilDiagonalMatrix The matrix which contains the main diagonal of self (or its inverse). """ + # Check `inverse` argument + assert isinstance(inverse, bool) + # Determine domain and codomain of the StencilDiagonalMatrix V, W = self.domain, self.codomain if inverse: V, W = W, V + # Check `out` argument if out is not None: - assert isinstance(out, StencilMatrix) + assert isinstance(out, StencilDiagonalMatrix) assert out.domain is V assert out.codomain is W - assert all(p == 0 for p in out.pads) # is this really needed? - else: - out = StencilMatrix(V, W, pads = [0] * V.ndim) - index = tuple(slice(s, e+1) for s, e in zip(W.starts, W.ends)) + (0,) * V.ndim + + # Extract diagonal data from self and identify output array + diagonal_indices = self._get_diagonal_indices() + diag = self._data[diagonal_indices] + data = out._data if out else None + + # Calculate entries of StencilDiagonalMatrix if inverse: - np.divide(1, self[index], out=out[index]) + data = np.divide(1, diag, out=data) + elif out: + np.copyto(data, diag) else: - out[index] = self[index] + data = diag.copy() - return out + # If needed create a new StencilDiagonalMatrix object + if out is None: + out = StencilDiagonalMatrix(V, W, data) + return out +# # ... # def diagonal(self): # if self._diag_indices is None: # cm = self.codomain.shifts @@ -1853,6 +1867,44 @@ def set_backend(self, backend): self._args.pop('ndiags') self._func = dot.func + # ... + def _get_diagonal_indices(self): + """ + Compute the indices which should be applied to self._data in order to + get the matrix entries on the main diagonal. The result is also stored + in self._diag_indices, and retrieved from there on successive calls. + + Returns + ------- + tuple[numpy.ndarray, ndim] + The diagonal indices as a tuple of NumPy arrays of identical shape + (n1, n2, n3, ...). + + """ + + if self._diag_indices is None: + + dp = self.domain.pads + dm = self.domain.shifts + cm = self.codomain.shifts + ss = self.codomain.starts + pp = [compute_diag_len(p, mj, mi) - p - 1 for p, mi, mj in zip(self._pads, cm, dm)] + nrows = [e - s + 1 for s, e in zip(self.codomain.starts, self.codomain.ends)] + ndim = self.domain.ndim + + indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2 * ndim)] + + for l, xx in enumerate(np.ndindex(*nrows)): + ii = [m * p + x for m, p, x in zip(dm, dp, xx)] + jj = [p + x + s - ((x+s) // mi) * mj for x, mi, mj, p, s in zip(xx, cm, dm, pp, ss)] + for k in range(ndim): + indices[k][l] = ii[k] + indices[k + ndim][l] = jj[k] + + self._diag_indices = tuple(idx.reshape(nrows) for idx in indices) + + return self._diag_indices + #=============================================================================== class StencilDiagonalMatrix(LinearOperator): """ From e90fb503395dda4fe72d4ff90a12d55b69c6a72b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 23 Jan 2024 09:43:17 +0100 Subject: [PATCH 21/34] Do not call diagonal() on empty block --- psydac/linalg/block.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 94640d719..02552e4cc 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -794,21 +794,24 @@ def diagonal(self, *, inverse = False, out = None): The matrix which contains the main diagonal of self (or its inverse). """ + # Determine domain and codomain of result V, W = self.domain, self.codomain if inverse: V, W = W, V + # Check the `out` argument, if `None` create a new BlockLinearOperator if out is not None: assert isinstance(out, BlockLinearOperator) assert out.domain is V assert out.codomain is W assert all(i==j for i, j in out._blocks.keys()) # is this really needed? - for i in range(min(self.n_block_rows, self.n_block_cols)): - self[i, i].diagonal(inverse = inverse, out = out[i, i]) else: out = BlockLinearOperator(V, W) - for i in range(min(self.n_block_rows, self.n_block_cols)): - out[i, i] = self[i, i].diagonal(inverse = inverse) + + # Store the diagonal (or its inverse) into `out` + for i, j in self.nonzero_block_indices: + if i == j: + out[i, i] = self[i, i].diagonal(inverse = inverse, out = out[i, i]) return out From f0da5b949bd766677d72b01d7cd8929c9b0a5112 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 29 Jan 2024 09:11:48 +0100 Subject: [PATCH 22/34] Minor cleanup --- psydac/linalg/basic.py | 1 - psydac/linalg/stencil.py | 23 ----------------------- 2 files changed, 24 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 3ccaa408a..2f2537e47 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -448,7 +448,6 @@ def __new__(cls, domain, codomain=None): return BlockLinearOperator(domain, domain, blocks) else: return super().__new__(cls) - def __init__(self, domain, codomain=None): diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 8c3cfd824..2c8c6fa96 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1424,29 +1424,6 @@ def diagonal(self, *, inverse = False, out = None): return out -# # ... -# def diagonal(self): -# if self._diag_indices is None: -# cm = self.codomain.shifts -# dm = self.domain.shifts -# ss = self.codomain.starts -# pp = [compute_diag_len(p,mj,mi)-(p+1) for p,mi,mj in zip(self._pads, cm, dm)] -# nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) -# indices = [np.zeros(np.product(nrows), dtype=int) for _ in range(2*len(nrows))] -# l = 0 -# for xx in np.ndindex(*nrows): -# ii = [m*p+x for m,p,x in zip(self.domain.shifts, self.domain.pads, xx)] -# jj = [p+x+s-((x+s)//mi)*mj for (x,mi,mj,p,s) in zip(xx,cm,dm,pp,ss)] -# for k in range(len(nrows)): -# indices[k][l] = ii[k] -# indices[k+len(nrows)][l] = jj[k] -# l += 1 -# self._diag_indices = tuple(indices) -# else: -# nrows = tuple(e-s+1 for s,e in zip(self.codomain.starts, self.codomain.ends)) -# -# return self._data[self._diag_indices].reshape(nrows) - # ... def topetsc(self): """ Convert to PETSc data structure. From 982663671ff4417d4719b86b0f75a2657af1ff1e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 29 Jan 2024 09:17:54 +0100 Subject: [PATCH 23/34] Remove solver property from InverseLinearOperator base class --- psydac/linalg/basic.py | 14 +++----------- psydac/linalg/solvers.py | 32 -------------------------------- 2 files changed, 3 insertions(+), 43 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 2f2537e47..72b3c5e29 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -939,12 +939,6 @@ def linop(self, a): def options(self): return self._options - @property - @abstractmethod - def solver(self): - "String that identifies the solver." - pass - def _check_options(self, **kwargs): for key, value in kwargs.items(): @@ -978,12 +972,10 @@ def set_options(self, **kwargs): self._options.update(kwargs) def transpose(self, conjugate=False): - from psydac.linalg.solvers import inverse - - At = self.linop.transpose(conjugate=conjugate) - solver = self.solver + cls = type(self) + At = self.linop.transpose(conjugate=conjugate) options = self._options - return inverse(At, solver, **options) + return cls(At, **options) #=============================================================================== class LinearSolver(ABC): diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index 1f64f8877..c4c571851 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -132,10 +132,6 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p")} self._info = None - @property - def solver(self): - return 'cg' - def solve(self, b, out=None): """ Conjugate gradient algorithm for solving linear system Ax=b. @@ -290,10 +286,6 @@ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False self._tmps = {**tmps_codomain, **tmps_domain} self._info = None - @property - def solver(self): - return 'pcg' - def solve(self, b, out=None): """ Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte @@ -452,10 +444,6 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vs", "rs", "ps")} self._info = None - @property - def solver(self): - return 'bicg' - def solve(self, b, out=None): """ Biconjugate gradient (BCG) algorithm for solving linear system Ax=b. @@ -644,10 +632,6 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("v", "r", "p", "vr", "r0")} self._info = None - @property - def solver(self): - return 'bicgstab' - def solve(self, b, out=None): """ Biconjugate gradient stabilized method (BCGSTAB) algorithm for solving linear system Ax=b. @@ -844,10 +828,6 @@ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False "rp0")} self._info = None - @property - def solver(self): - return 'pbicgstab' - def solve(self, b, out=None): """ Preconditioned biconjugate gradient stabilized method (PBCGSTAB) algorithm for solving linear system Ax=b. @@ -1085,10 +1065,6 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle self._tmps = {key: self.domain.zeros() for key in ("res_old", "res_new", "w_new", "w_work", "w_old", "v", "y")} self._info = None - @property - def solver(self): - return 'minres' - def solve(self, b, out=None): """ Use MINimum RESidual iteration to solve Ax=b @@ -1395,10 +1371,6 @@ def __init__(self, A, *, x0=None, tol=None, atol=None, btol=None, maxiter=1000, tmps_codomain = {key: self.codomain.zeros() for key in ("v", "v_work", "h", "hbar")} self._tmps = {**tmps_codomain, **tmps_domain} - @property - def solver(self): - return 'lsmr' - def get_success(self): return self._successful @@ -1739,10 +1711,6 @@ def __init__(self, A, *, x0=None, tol=1e-6, maxiter=100, verbose=False, recycle= self._Q = [] self._info = None - @property - def solver(self): - return 'gmres' - def solve(self, b, out=None): """ Generalized minimal residual algorithm for solving linear system Ax=b. From 54f467123afa06f34a93616e1643d9588937aa69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 30 Jan 2024 14:03:34 +0100 Subject: [PATCH 24/34] Add docstring to linop property of InverseLinearOperator: The text clarifies that replacing the original linear operator through the setter is a dangerous operation which should only be made in extreme cases. --- psydac/linalg/basic.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 72b3c5e29..b23417068 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -926,6 +926,16 @@ def dtype(self): @property def linop(self): + """ + The linear operator L of which this object is the inverse L^{-1}. + + The linear operator L can be modified in place, or replaced entirely + through the setter. A substitution should only be made in cases where + no other options are viable, as it breaks the one-to-one map between + the original linear operator L (passed to the constructor) and the + current `InverseLinearOperator` object L^{-1}. Use with extreme care! + + """ return self._A @linop.setter From 24b9be29def76b05600c29c4faddfbf9565f1341 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 30 Jan 2024 19:02:50 +0100 Subject: [PATCH 25/34] Set off-diagonal blocks to zero in 'out' argument of diagonal() --- psydac/linalg/block.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 02552e4cc..e8a9c44d6 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -804,7 +804,11 @@ def diagonal(self, *, inverse = False, out = None): assert isinstance(out, BlockLinearOperator) assert out.domain is V assert out.codomain is W - assert all(i==j for i, j in out._blocks.keys()) # is this really needed? + + # Set any off-diagonal blocks to zero + for i, j in out.nonzero_block_indices: + if i != j: + out[i, j] = None else: out = BlockLinearOperator(V, W) From f514ef43e152e1c6b6e56d9b503f5a550ad4749d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 30 Jan 2024 19:25:40 +0100 Subject: [PATCH 26/34] Improve docstring of PConjugateGradient class --- psydac/linalg/solvers.py | 29 +++++++++++++++++++---------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index c4c571851..efcb1b620 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -242,32 +242,41 @@ class PConjugateGradient(InverseLinearOperator): A LinearOperator subclass. Objects of this class are meant to be created using :func:~`solvers.inverse`. The .dot (and also the .solve) function are based on a preconditioned conjugate gradient method. - Preconditioned Conjugate Gradient (PCG) solves the symetric positive definte - system Ax = b. It assumes that pc.dot(r) returns the solution to Ps = r, - where P is positive definite. + The Preconditioned Conjugate Gradient (PCG) algorithm solves the linear + system A x = b where A is a symmetric and positive-definite matrix, i.e. + A = A^T and y A y > 0 for any vector y. The preconditioner P is a matrix + which approximates the inverse of A. The algorithm assumes that P is also + symmetric and positive definite. + + Since this is a matrix-free iterative method, both A and P are provided as + `LinearOperator` objects which must implement the `dot` method. Parameters ---------- - A : psydac.linalg.stencil.StencilMatrix - Left-hand-side matrix A of linear system + A : psydac.linalg.basic.LinearOperator + Left-hand-side matrix A of the linear system. This should be symmetric + and positive definite. pc: psydac.linalg.basic.LinearOperator - Preconditioner for A, it should approximate the inverse of A (can be None). + Preconditioner which should approximate the inverse of A (optional). + Like A, the preconditioner should be symmetric and positive definite. x0 : psydac.linalg.basic.Vector First guess of solution for iterative solver (optional). tol : float - Absolute tolerance for L2-norm of residual r = A*x - b. + Absolute tolerance for L2-norm of residual r = A x - b. (Default: 1e-6) maxiter: int - Maximum number of iterations. + Maximum number of iterations. (Default: 1000) verbose : bool - If True, L2-norm of residual r is printed at each iteration. + If True, the L2-norm of the residual r is printed at each iteration. + (Default: False) recycle : bool - Stores a copy of the output in x0 to speed up consecutive calculations of slightly altered linear systems + If True, a copy of the output is stored in x0 to speed up consecutive + calculations of slightly altered linear systems. (Default: False) """ def __init__(self, A, *, pc=None, x0=None, tol=1e-6, maxiter=1000, verbose=False, recycle=False): From ca6593bd0937cc761a025d539bd93fcf6d63deec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 10:37:34 +0100 Subject: [PATCH 27/34] Avoid unnecessary transpose operations in test_inverse_transpose_interaction --- psydac/linalg/tests/test_linalg.py | 50 +++++++++++++++++------------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index b3b9614a1..a011e8f83 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -623,23 +623,26 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): tol = 1e-5 C = inverse(B, 'cg', tol=tol) - pc = B.T.diagonal(inverse=True) + P = B.diagonal(inverse=True) + + B_T = B.T + C_T = C.T # -1,T & T,-1 -> equal - assert isinstance(C.T, ConjugateGradient) - assert isinstance(inverse(B.T, 'cg', tol=tol), ConjugateGradient) - assert np.sqrt(sum(((C.T.dot(u) - inverse(B.T, 'cg', tol=tol).dot(u)).toarray())**2)) < 2*tol + assert isinstance(C_T, ConjugateGradient) + assert isinstance(inverse(B_T, 'cg', tol=tol), ConjugateGradient) + assert np.sqrt(sum(((C_T.dot(u) - inverse(B_T, 'cg', tol=tol).dot(u)).toarray())**2)) < 2*tol # -1,T,T -> equal -1 - assert np.sqrt(sum(((C.T.T.dot(u) - C.dot(u)).toarray())**2)) < 2*tol + assert np.sqrt(sum(((C_T.T.dot(u) - C.dot(u)).toarray())**2)) < 2*tol # -1,T,-1 -> equal T - assert isinstance(inverse(C.T, 'bicg'), BlockLinearOperator) - assert np.array_equal(inverse(C.T, 'bicg').dot(u).toarray(), B.T.dot(u).toarray()) + assert isinstance(inverse(C_T, 'bicg'), BlockLinearOperator) + assert np.array_equal(inverse(C_T, 'bicg').dot(u).toarray(), B_T.dot(u).toarray()) # T,-1,-1 -> equal T - assert isinstance(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc=pc), BlockLinearOperator) - assert np.array_equal(inverse(inverse(B.T, 'cg', tol=tol), 'pcg', pc=pc).dot(u).toarray(), B.T.dot(u).toarray()) + assert isinstance(inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P), BlockLinearOperator) + assert np.array_equal(inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P).dot(u).toarray(), B_T.dot(u).toarray()) # T,-1,T -> equal -1 - assert isinstance(inverse(B.T, 'cg', tol=tol).T, ConjugateGradient) - assert np.sqrt(sum(((inverse(B.T, 'cg', tol=tol).dot(u) - C.dot(u)).toarray())**2)) < tol + assert isinstance(inverse(B_T, 'cg', tol=tol).T, ConjugateGradient) + assert np.sqrt(sum(((inverse(B_T, 'cg', tol=tol).dot(u) - C.dot(u)).toarray())**2)) < tol ### ### StencilMatrix Transpose - Inverse Tests @@ -648,23 +651,26 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): tol = 1e-5 C = inverse(S, 'cg', tol=tol) - pc = S.T.diagonal(inverse=True) + P = S.diagonal(inverse=True) + + S_T = S.T + C_T = C.T # -1,T & T,-1 -> equal - assert isinstance(C.T, ConjugateGradient) - assert isinstance(inverse(S.T, 'cg', tol=tol), ConjugateGradient) - assert np.sqrt(sum(((C.T.dot(v) - inverse(S.T, 'cg', tol=tol).dot(v)).toarray())**2)) < 2*tol + assert isinstance(C_T, ConjugateGradient) + assert isinstance(inverse(S_T, 'cg', tol=tol), ConjugateGradient) + assert np.sqrt(sum(((C_T.dot(v) - inverse(S_T, 'cg', tol=tol).dot(v)).toarray())**2)) < 2*tol # -1,T,T -> equal -1 - assert np.sqrt(sum(((C.T.T.dot(v) - C.dot(v)).toarray())**2)) < 2*tol + assert np.sqrt(sum(((C_T.T.dot(v) - C.dot(v)).toarray())**2)) < 2*tol # -1,T,-1 -> equal T - assert isinstance(inverse(C.T, 'bicg'), StencilMatrix) - assert np.array_equal(inverse(C.T, 'bicg').dot(v).toarray(), S.T.dot(v).toarray()) + assert isinstance(inverse(C_T, 'bicg'), StencilMatrix) + assert np.array_equal(inverse(C_T, 'bicg').dot(v).toarray(), S_T.dot(v).toarray()) # T,-1,-1 -> equal T - assert isinstance(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc=pc), StencilMatrix) - assert np.array_equal(inverse(inverse(S.T, 'cg', tol=tol), 'pcg', pc=pc).dot(v).toarray(), S.T.dot(v).toarray()) + assert isinstance(inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P), StencilMatrix) + assert np.array_equal(inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P).dot(v).toarray(), S_T.dot(v).toarray()) # T,-1,T -> equal -1 - assert isinstance(inverse(S.T, 'cg', tol=tol).T, ConjugateGradient) - assert np.sqrt(sum(((inverse(S.T, 'cg', tol=tol).dot(v) - C.dot(v)).toarray())**2)) < tol + assert isinstance(inverse(S_T, 'cg', tol=tol).T, ConjugateGradient) + assert np.sqrt(sum(((inverse(S_T, 'cg', tol=tol).dot(v) - C.dot(v)).toarray())**2)) < tol #=============================================================================== @pytest.mark.parametrize('n1', [3, 5]) From fd61f5f1e8e841f7b69c9e30e9ff0c816ea9e13f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 11:15:08 +0100 Subject: [PATCH 28/34] Clean up test_inverse_transpose_interaction --- psydac/linalg/tests/test_linalg.py | 38 ++++++++++++++++++++++-------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index a011e8f83..da9945488 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -631,18 +631,27 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): # -1,T & T,-1 -> equal assert isinstance(C_T, ConjugateGradient) assert isinstance(inverse(B_T, 'cg', tol=tol), ConjugateGradient) - assert np.sqrt(sum(((C_T.dot(u) - inverse(B_T, 'cg', tol=tol).dot(u)).toarray())**2)) < 2*tol + diff = C_T @ u - inverse(B_T, 'cg', tol=tol) @ u + assert diff.dot(diff) == 0 + # -1,T,T -> equal -1 - assert np.sqrt(sum(((C_T.T.dot(u) - C.dot(u)).toarray())**2)) < 2*tol + diff = C_T.T @ u - C @ u + assert diff.dot(diff) == 0 + # -1,T,-1 -> equal T assert isinstance(inverse(C_T, 'bicg'), BlockLinearOperator) - assert np.array_equal(inverse(C_T, 'bicg').dot(u).toarray(), B_T.dot(u).toarray()) + diff = inverse(C_T, 'bicg') @ u - B_T @ u + assert diff.dot(diff) == 0 + # T,-1,-1 -> equal T assert isinstance(inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P), BlockLinearOperator) - assert np.array_equal(inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P).dot(u).toarray(), B_T.dot(u).toarray()) + diff = inverse(inverse(B_T, 'cg', tol=tol), 'pcg', pc=P) @ u - B_T @ u + assert diff.dot(diff) == 0 + # T,-1,T -> equal -1 assert isinstance(inverse(B_T, 'cg', tol=tol).T, ConjugateGradient) - assert np.sqrt(sum(((inverse(B_T, 'cg', tol=tol).dot(u) - C.dot(u)).toarray())**2)) < tol + diff = inverse(B_T, 'cg', tol=tol) @ u - C @ u + assert diff.dot(diff) == 0 ### ### StencilMatrix Transpose - Inverse Tests @@ -659,18 +668,27 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): # -1,T & T,-1 -> equal assert isinstance(C_T, ConjugateGradient) assert isinstance(inverse(S_T, 'cg', tol=tol), ConjugateGradient) - assert np.sqrt(sum(((C_T.dot(v) - inverse(S_T, 'cg', tol=tol).dot(v)).toarray())**2)) < 2*tol + diff = C_T @ v - inverse(S_T, 'cg', tol=tol) @ v + assert diff.dot(diff) == 0 + # -1,T,T -> equal -1 - assert np.sqrt(sum(((C_T.T.dot(v) - C.dot(v)).toarray())**2)) < 2*tol + diff = C_T.T @ v - C @ v + assert diff.dot(diff) == 0 + # -1,T,-1 -> equal T assert isinstance(inverse(C_T, 'bicg'), StencilMatrix) - assert np.array_equal(inverse(C_T, 'bicg').dot(v).toarray(), S_T.dot(v).toarray()) + diff = inverse(C_T, 'bicg') @ v - S_T @ v + assert diff.dot(diff) == 0 + # T,-1,-1 -> equal T assert isinstance(inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P), StencilMatrix) - assert np.array_equal(inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P).dot(v).toarray(), S_T.dot(v).toarray()) + diff = inverse(inverse(S_T, 'cg', tol=tol), 'pcg', pc=P) @ v - S_T @ v + assert diff.dot(diff) == 0 + # T,-1,T -> equal -1 assert isinstance(inverse(S_T, 'cg', tol=tol).T, ConjugateGradient) - assert np.sqrt(sum(((inverse(S_T, 'cg', tol=tol).dot(v) - C.dot(v)).toarray())**2)) < tol + diff = inverse(S_T, 'cg', tol=tol) @ v - C @ v + assert diff.dot(diff) == 0 #=============================================================================== @pytest.mark.parametrize('n1', [3, 5]) From ee5070f64195eaeee94a9a29bd77f821e3bb8b70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 11:33:12 +0100 Subject: [PATCH 29/34] Minor cleanup in psydac.linalg.utilities --- psydac/linalg/utilities.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py index b486bb8e6..515cef8c9 100644 --- a/psydac/linalg/utilities.py +++ b/psydac/linalg/utilities.py @@ -1,12 +1,18 @@ # coding: utf-8 import numpy as np -from math import sqrt +from math import sqrt + from psydac.linalg.stencil import StencilVectorSpace, StencilVector from psydac.linalg.block import BlockVector, BlockVectorSpace -__all__ = ('array_to_psydac', 'petsc_to_psydac', '_sym_ortho') +__all__ = ( + 'array_to_psydac', + 'petsc_to_psydac', + '_sym_ortho' +) +#============================================================================== def array_to_psydac(x, Xh): """ converts a numpy array to StencilVector or BlockVector format""" @@ -47,6 +53,7 @@ def array_to_psydac(x, Xh): u.update_ghost_regions() return u +#============================================================================== def petsc_to_psydac(vec, Xh): """ converts a petsc Vec object to a StencilVector or a BlockVector format. We gather the petsc global vector in all the processes and extract the chunk owned by the Psydac Vector. @@ -129,6 +136,7 @@ def petsc_to_psydac(vec, Xh): u.update_ghost_regions() return u +#============================================================================== def _sym_ortho(a, b): """ Stable implementation of Givens rotation. From 4b0f821f217d6ee3a987cc0c3c35eabe381ff808 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 11:33:48 +0100 Subject: [PATCH 30/34] Move is_real function to psydac.utilities.utils --- psydac/linalg/basic.py | 27 ++++++++++++++++++--------- psydac/linalg/solvers.py | 13 +++++-------- psydac/utilities/utils.py | 25 ++++++++++++++++++++++++- 3 files changed, 47 insertions(+), 18 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index b23417068..6acd79c74 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -3,12 +3,26 @@ # Copyright 2018 Yaman Güçlü, Jalal Lakhlili # Copyright 2022 Yaman Güçlü, Said Hadjout, Julian Owezarek -from abc import ABC, abstractmethod -from scipy.sparse import coo_matrix +from abc import ABC, abstractmethod + import numpy as np +from scipy.sparse import coo_matrix -__all__ = ('VectorSpace', 'Vector', 'LinearOperator', 'ZeroOperator', 'IdentityOperator', 'ScaledLinearOperator', - 'SumLinearOperator', 'ComposedLinearOperator', 'PowerLinearOperator', 'InverseLinearOperator', 'LinearSolver') +from psydac.utilities.utils import is_real + +__all__ = ( + 'VectorSpace', + 'Vector', + 'LinearOperator', + 'ZeroOperator', + 'IdentityOperator', + 'ScaledLinearOperator', + 'SumLinearOperator', + 'ComposedLinearOperator', + 'PowerLinearOperator', + 'InverseLinearOperator', + 'LinearSolver' +) #=============================================================================== class VectorSpace(ABC): @@ -1018,8 +1032,3 @@ def solve(self, rhs, out=None): @property def T(self): return self.transpose() - - -def is_real(x): - from numbers import Number - return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool) diff --git a/psydac/linalg/solvers.py b/psydac/linalg/solvers.py index efcb1b620..aa0b5e21b 100644 --- a/psydac/linalg/solvers.py +++ b/psydac/linalg/solvers.py @@ -1,13 +1,15 @@ # coding: utf-8 """ -This module provides iterative solvers and precondionners. +This module provides iterative solvers and preconditioners. """ -from math import sqrt import numpy as np +from math import sqrt -from psydac.linalg.basic import Vector, LinearOperator, InverseLinearOperator, IdentityOperator, ScaledLinearOperator +from psydac.utilities.utils import is_real from psydac.linalg.utilities import _sym_ortho +from psydac.linalg.basic import (Vector, LinearOperator, + InverseLinearOperator, IdentityOperator, ScaledLinearOperator) __all__ = ( 'ConjugateGradient', @@ -20,11 +22,6 @@ 'GMRES' ) -#=============================================================================== -def is_real(x): - from numbers import Number - return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool) - #=============================================================================== def inverse(A, solver, **kwargs): """ diff --git a/psydac/utilities/utils.py b/psydac/utilities/utils.py index 623120569..279b1bd60 100644 --- a/psydac/utilities/utils.py +++ b/psydac/utilities/utils.py @@ -3,8 +3,31 @@ # Copyright 2018 Yaman Güçlü import numpy as np +from numbers import Number -__all__ = ('refine_array_1d', 'unroll_edges', 'split_space', 'split_field', 'animate_field') +__all__ = ( + 'refine_array_1d', + 'unroll_edges', + 'split_space', + 'split_field', + 'animate_field' +) + +#============================================================================== +def is_real(x): + """Determine whether the given input represents a real number. + + Parameters + ---------- + x : Any + + Returns + ------- + bool + True if x is real, False otherwise. + + """ + return isinstance(x, Number) and np.isrealobj(x) and not isinstance(x, bool) #=============================================================================== def refine_array_1d(x, n, remove_duplicates=True): From 0fa1958d6eddb2038395940c1bb4248431f4ce51 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 12:10:57 +0100 Subject: [PATCH 31/34] Improve docstring of InverseLinearOperator base class --- psydac/linalg/basic.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 6acd79c74..c2e95be0d 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -889,9 +889,17 @@ def dot(self, v, out=None): #=============================================================================== class InverseLinearOperator(LinearOperator): """ - Iterative solver for square linear system Ax=b, where x and b belong to (normed) + Abstract base class for the (approximate) inverse A_inv := A^{-1} of a + square matrix A. The result of A_inv.dot(b) is the (approximate) solution x + of the linear system A x = b, where x and b belong to the same (normed) vector space V. + We assume that the linear system is solved by an iterative method, which + needs a first guess `x0` and an exit condition based on `tol` and `maxiter`. + + Concrete subclasses of this class must implement the `dot` method and take + care of any internal storage which might be necessary. + Parameters ---------- A : psydac.linalg.basic.LinearOperator From 9fe92f11caf46f0d757f12c2e17d93879bb931ec Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 12:14:05 +0100 Subject: [PATCH 32/34] Remove useless 'options' property from InverseLinearOperator --- psydac/linalg/basic.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index c2e95be0d..0a64a9bed 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -931,8 +931,9 @@ def __init__(self, A, **kwargs): self._A = A self._domain = domain self._codomain = codomain + + self._check_options(**kwargs) self._options = kwargs - self._check_options(**self._options) @property def domain(self): @@ -967,10 +968,6 @@ def linop(self, a): assert a.codomain is self.codomain self._A = a - @property - def options(self): - return self._options - def _check_options(self, **kwargs): for key, value in kwargs.items(): @@ -997,9 +994,13 @@ def get_info(self): return self._info def get_options(self): + """Get a copy of the solver options. + """ return self._options.copy() def set_options(self, **kwargs): + """Set the solver options by passing keyword arguments. + """ self._check_options(**kwargs) self._options.update(kwargs) From 7def82786ce804fbdc323078914c5b26b0aeb0c1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 13:23:03 +0100 Subject: [PATCH 33/34] Improve get_options method of InverseLinearOperator, fix bug --- psydac/linalg/basic.py | 22 +++++++++++++++++++--- psydac/linalg/tests/test_linalg.py | 24 +++++++++++------------- 2 files changed, 30 insertions(+), 16 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 0a64a9bed..57d7adcd8 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -993,10 +993,26 @@ def tosparse(self): def get_info(self): return self._info - def get_options(self): - """Get a copy of the solver options. + def get_options(self, key=None): + """Get a copy of all the solver options, or a specific value of interest. + + Parameters + ---------- + x0 : str | None + Name of the specific option of interest (default: None). + + Returns + ------- + dict | type(self._options['key'] | None + If `key` is given, get the specific option of interest. If there is + no such option, `None` is returned instead. If `key` is not given, + get a copy of all the solver options in a dictionary. + """ - return self._options.copy() + if key is None: + return self._options.copy() + else: + return self._options.get(key) def set_options(self, **kwargs): """Set the solver options by passing keyword arguments. diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index da9945488..8b30802b8 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -973,28 +973,26 @@ def test_x0update(solver): A_inv = inverse(A, solver, pc=A.diagonal(inverse=True), tol=tol) else: A_inv = inverse(A, solver, tol=tol) - options = A_inv.options - x0_init = options["x0"] + # Check whether x0 is not None + x0_init = A_inv.get_options("x0") assert x0_init is not None + # Apply inverse and check x0 x = A_inv @ b - options = A_inv.options - x0_new1 = options["x0"] + x0_new1 = A_inv.get_options("x0") assert x0_new1 is x0_init + # Change x0, apply A_inv and check for x0 A_inv.set_options(x0 = b) - options = A_inv.options - assert options["x0"] is b + assert A_inv.get_options("x0") is b + x = A_inv @ b - options = A_inv.options - x0_new2 = options["x0"] - assert x0_new2 is b + assert A_inv.get_options("x0") is b + # Apply inverse using out=x0 and check for updated x0 - x = A_inv.dot(b, out=x0_new2) - options = A_inv.options - x0_new3 = options["x0"] - assert x0_new3 is x + x = A_inv.dot(b, out=b) + assert A_inv.get_options('x0') is x #=============================================================================== # SCRIPT FUNCTIONALITY From 1e67613cceec7b4f892ba31df098b4c0f52753a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 31 Jan 2024 14:24:43 +0100 Subject: [PATCH 34/34] Fix docstring of get_options --- psydac/linalg/basic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 57d7adcd8..352b4e380 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -998,12 +998,12 @@ def get_options(self, key=None): Parameters ---------- - x0 : str | None + key : str | None Name of the specific option of interest (default: None). Returns ------- - dict | type(self._options['key'] | None + dict | type(self._options['key']) | None If `key` is given, get the specific option of interest. If there is no such option, `None` is returned instead. If `key` is not given, get a copy of all the solver options in a dictionary.