Skip to content

Fix design problem multi_shift #147

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gnrraphi wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Fix design problem multi_shift #147

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 40 additions & 27 deletions lib/gpt/algorithms/inverter/multi_shift_fgmres.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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__()
Expand All @@ -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 = []
Expand All @@ -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]

Expand All @@ -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):
Expand Down Expand Up @@ -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 = []
Expand All @@ -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")
Expand All @@ -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())
Expand All @@ -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}"
Expand All @@ -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:
Expand All @@ -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,
Expand Down
22 changes: 13 additions & 9 deletions lib/gpt/algorithms/inverter/multi_shift_fom.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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__()
Expand All @@ -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 = []
Expand Down Expand Up @@ -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)

Expand All @@ -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")
Expand All @@ -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,
Expand Down
14 changes: 9 additions & 5 deletions tests/algorithms/solvers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
}
Expand Down