From 8b7ce3bc71c1a911d8d6ba89c656f88270fcfa1d Mon Sep 17 00:00:00 2001 From: gnrraphi Date: Sat, 24 Feb 2024 20:40:30 +0100 Subject: [PATCH] fix design --- .../algorithms/inverter/multi_shift_fgmres.py | 67 +++++++++++-------- .../algorithms/inverter/multi_shift_fom.py | 22 +++--- tests/algorithms/solvers.py | 14 ++-- 3 files changed, 62 insertions(+), 41 deletions(-) diff --git a/lib/gpt/algorithms/inverter/multi_shift_fgmres.py b/lib/gpt/algorithms/inverter/multi_shift_fgmres.py index 8328c5b79..dd330596c 100644 --- a/lib/gpt/algorithms/inverter/multi_shift_fgmres.py +++ b/lib/gpt/algorithms/inverter/multi_shift_fgmres.py @@ -1,7 +1,7 @@ # # GPT - Grid Python Toolkit -# Copyright (C) 2022 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) -# 2022 Raphael Lehner (raphael.lehner@physik.uni-regensburg.de) +# Copyright (C) 2022-24 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) +# 2022-24 Raphael Lehner (raphael.lehner@physik.uni-regensburg.de) # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -17,6 +17,8 @@ # with this program; if not, write to the Free Software Foundation, Inc., # 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. # +# Reference: https://doi.org/10.1137/140979927 +# import gpt as g import numpy as np @@ -118,7 +120,7 @@ class multi_shift_fgmres(base_iterative): shifts=[], checkres=True, prec=None, - rhos=False, + saverhos=False, ) def __init__(self, params): super().__init__() @@ -129,7 +131,8 @@ def __init__(self, params): self.shifts = params["shifts"] self.checkres = params["checkres"] self.prec = params["prec"] - self.rhos = params["rhos"] + self.saverhos = params["saverhos"] + self.rhos = None def arnoldi(self, mat, V, rlen, mmp, sfgmres, prec, idx, t): H = [] @@ -139,7 +142,8 @@ def arnoldi(self, mat, V, rlen, mmp, sfgmres, prec, idx, t): Z = [sfgmres[j].Z[i] for j in idx] + [mmp] for z in Z: z[:] = 0 - rhos = prec(Z, [V[i]]) + prec(Z, V[i]) + rhos = self.prec.rhos for j, rho in zip(idx, rhos[0:-1]): sfgmres[j].gamma[i] = rho / rhos[-1] @@ -162,20 +166,20 @@ def setup_prec(self, mat): if self.prec is not None: prec_shifts = self.shifts + [0.0] self.prec.shifts = prec_shifts - self.prec.rhos = True - return self.prec(mat).mat + self.prec.saverhos = True + return self.prec(mat) return None - def restart_prec(self, mat, sfgmres): - shifts_prec = [self.shifts[0]] - idx = [0] - for j, fgmres in enumerate(sfgmres[1:]): - if fgmres.converged is False: + def restart_prec(self, mat, sfgmres, s0): + shifts_prec = [] + idx = [] + for j, fgmres in enumerate(sfgmres): + if fgmres.converged is False or j == s0: shifts_prec += [fgmres.s] - idx += [j + 1] + idx += [j] shifts_prec += [0.0] self.prec.shifts = shifts_prec - prec = self.prec(mat).mat + prec = self.prec(mat) return prec, idx def __call__(self, mat): @@ -214,7 +218,10 @@ def inv(psi, src, t): # prec prec = self.setup_prec(mat) plen = len(self.shifts) - idx = [i for i in range(plen)] + pidx = [i for i in range(plen)] + + # sorted shifts + sidx = np.argsort(self.shifts) # shifted systems sfgmres = [] @@ -225,15 +232,15 @@ def inv(psi, src, t): V = [g.copy(src) for i in range(rlen + 1)] V[0] /= r2**0.5 - # return rhos for prec fgmres - rr = self.rhos + # save rhos for prec multi_shift_fgmres + sr = self.saverhos for k in range(0, self.maxiter, rlen): # arnoldi - H = self.arnoldi(mat, V, rlen, mmp, sfgmres, prec, idx, t) + H = self.arnoldi(mat, V, rlen, mmp, sfgmres, prec, pidx, t) t("hessenberg") - fgmres = sfgmres[0] + fgmres = sfgmres[sidx[0]] Hs = fgmres.hessenberg(H, prec) t("qr") @@ -252,8 +259,9 @@ def inv(psi, src, t): t("other") self.log_convergence((k, 0), r2_new, rsq) - for j, fgmres in enumerate(sfgmres[1:]): - if fgmres.converged is False or rr: + for j, i in enumerate(sidx[1:]): + fgmres = sfgmres[i] + if fgmres.converged is False or sr is True: t("hessenberg") Hs = fgmres.hessenberg(H, prec) Hs.append(vr.copy()) @@ -268,7 +276,8 @@ def inv(psi, src, t): self.log_convergence((k, j + 1), fgmres.r2, rsq) t("other") - for fgmres in sfgmres: + for i in sidx: + fgmres = sfgmres[i] msg = fgmres.check(rsq) if msg: msg += f" at iteration {k+rlen}" @@ -281,23 +290,26 @@ def inv(psi, src, t): if all([fgmres.converged for fgmres in sfgmres]): self.log(f"converged in {k+rlen} iterations") - return [fgmres.rho for fgmres in sfgmres] if rr else None + if sr is True: + self.rhos = [fgmres.rho for fgmres in sfgmres] + return if self.maxiter != rlen: t("restart") r2 = g.norm2(r) V[0] @= r / r2**0.5 - if prec is not None and rr is False: + if prec is not None and sr is False: t("restart_prec") plen_new = sum([not fgmres.converged for fgmres in sfgmres[1:]]) + 1 if plen_new != plen: plen = plen_new - prec, idx = self.restart_prec(mat, sfgmres) + prec, pidx = self.restart_prec(mat, sfgmres, sidx[0]) self.debug("performed restart") t("other") - for fgmres in sfgmres: + for i in sidx: + fgmres = sfgmres[i] if fgmres.converged is False: msg = f"shift {fgmres.s} NOT converged in {k+rlen} iterations" if self.maxiter != rlen: @@ -309,7 +321,8 @@ def inv(psi, src, t): cs = sum([fgmres.converged for fgmres in sfgmres if True]) ns = len(self.shifts) self.log(f"NOT converged in {k+rlen} iterations; {cs} / {ns} converged shifts") - return [fgmres.rho for fgmres in sfgmres] if rr else None + if sr is True: + self.rhos = [fgmres.rho for fgmres in sfgmres] return g.matrix_operator( mat=inv, diff --git a/lib/gpt/algorithms/inverter/multi_shift_fom.py b/lib/gpt/algorithms/inverter/multi_shift_fom.py index 5246eea8b..0efd2d6c5 100644 --- a/lib/gpt/algorithms/inverter/multi_shift_fom.py +++ b/lib/gpt/algorithms/inverter/multi_shift_fom.py @@ -1,7 +1,7 @@ # # GPT - Grid Python Toolkit -# Copyright (C) 2022 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) -# 2022 Raphael Lehner (raphael.lehner@physik.uni-regensburg.de) +# Copyright (C) 2022-24 Christoph Lehner (christoph.lehner@ur.de, https://github.com/lehner/gpt) +# 2022-24 Raphael Lehner (raphael.lehner@physik.uni-regensburg.de) # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -79,7 +79,7 @@ class multi_shift_fom(base_iterative): restartlen=20, shifts=[], checkres=True, - rhos=False, + saverhos=False, ) def __init__(self, params): super().__init__() @@ -89,7 +89,8 @@ def __init__(self, params): self.restartlen = params["restartlen"] self.shifts = params["shifts"] self.checkres = params["checkres"] - self.rhos = params["rhos"] + self.saverhos = params["saverhos"] + self.rhos = None def arnoldi(self, mat, V, rlen): H = [] @@ -155,15 +156,15 @@ def inv(psi, src, t): V = [g.copy(src) for i in range(rlen + 1)] V[0] /= r2**0.5 - # return rhos for prec fgmres - rr = self.rhos + # save rhos for prec multi_shift_fgmres + sr = self.saverhos for k in range(0, self.maxiter, rlen): t("arnoldi") H = self.arnoldi(mat, V, rlen) for j, fom in enumerate(sfoms): - if fom.converged is False or rr: + if fom.converged is False or sr is True: t("solve_hessenberg") fom.solve_hessenberg(H, r2) @@ -187,7 +188,9 @@ def inv(psi, src, t): if all([fom.converged for fom in sfoms]): self.log(f"converged in {k+rlen} iterations") - return [fom.rho for fom in sfoms] if rr else None + if sr is True: + self.rhos = [fom.rho for fom in sfoms] + return if self.maxiter != rlen: t("restart") @@ -207,7 +210,8 @@ def inv(psi, src, t): cs = sum([fom.converged for fom in sfoms if True]) ns = len(self.shifts) self.log(f"NOT converged in {k+rlen} iterations; {cs} / {ns} converged shifts") - return [fom.rho for fom in sfoms] if rr else None + if sr is True: + self.rhos = [fom.rho for fom in sfoms] return g.matrix_operator( mat=inv, diff --git a/tests/algorithms/solvers.py b/tests/algorithms/solvers.py index 811b96859..14f5aadf7 100755 --- a/tests/algorithms/solvers.py +++ b/tests/algorithms/solvers.py @@ -179,7 +179,7 @@ def test(slv, name): # ordering of dst fields. cg = inv.cg({"eps": 1e-8, "maxiter": 500}) -shifts = [0.5, 1.0, 1.7] +shifts = [0.5, 1.7, 1.0] mat = eo2_odd(w).Mpc # also test with multiple sources @@ -197,24 +197,28 @@ def test(slv, name): mscg = inv.multi_shift_cg({"eps": 1e-8, "maxiter": 1024, "shifts": shifts}) g.default.set_verbose("multi_shift_fom") -msfom = inv.multi_shift_fom({"eps": 1e-8, "maxiter": 1024, "restartlen": 10, "shifts": shifts}) +msfom = inv.multi_shift_fom({"eps": 1e-8, "maxiter": 1024, "restartlen": 4, "shifts": shifts}) g.default.set_verbose("multi_shift_fgmres") msfgmres = inv.multi_shift_fgmres( - {"eps": 1e-8, "maxiter": 1024, "restartlen": 10, "shifts": shifts} + {"eps": 1e-8, "maxiter": 1024, "restartlen": 4, "shifts": shifts} ) +g.default.set_verbose("multi_shift_fom", False) prec_fom = inv.multi_shift_fom({"maxiter": 4, "restartlen": 2}) +g.default.set_verbose("multi_shift_fgmres") msfgmres_fom = inv.multi_shift_fgmres( - {"eps": 1e-8, "maxiter": 512, "restartlen": 5, "shifts": shifts, "prec": prec_fom} + {"eps": 1e-8, "maxiter": 512, "restartlen": 2, "shifts": shifts, "prec": prec_fom} ) +g.default.set_verbose("multi_shift_fgmres", False) prec_fgmres = inv.multi_shift_fgmres({"maxiter": 4, "restartlen": 2}) +g.default.set_verbose("multi_shift_fgmres") msfgmres_fgmres = inv.multi_shift_fgmres( { "eps": 1e-8, "maxiter": 512, - "restartlen": 5, + "restartlen": 2, "shifts": shifts, "prec": prec_fgmres, }