Skip to content

Commit

Permalink
Merge pull request #2402 from Ericgig/misc.krylov_step_computation
Browse files Browse the repository at this point in the history
Fix raise on too large min_step in krylov
  • Loading branch information
Ericgig committed Apr 23, 2024
2 parents 1a3140a + b9bc56d commit 12bf7f7
Show file tree
Hide file tree
Showing 2 changed files with 21 additions and 9 deletions.
14 changes: 8 additions & 6 deletions qutip/solver/integrator/krylov.py
Expand Up @@ -51,7 +51,7 @@ def _prepare(self):
krylov_tridiag, krylov_basis = \
self._lanczos_algorithm(rand_ket(N).data)
if (
krylov_tridiag.shape[0] < self.options["krylov_dim"]
krylov_tridiag.shape[0] < krylov_dim
or krylov_tridiag.shape[0] == N
):
self._max_step = np.inf
Expand Down Expand Up @@ -138,20 +138,22 @@ def krylov_error(t):
self._compute_psi(t, *reduced_state)
) / self.options["atol"])

dt = self.options["min_step"]
# Under 0 will cause an infinite loop in the while loop bellow.
dt = max(self.options["min_step"], 1e-14)
max_step = max(self.options["max_step"], dt)
err = krylov_error(dt)
if err > 0:
ValueError(
raise ValueError(
f"With the krylov dim of {self.options['krylov_dim']}, the "
f"error with the minimum step {dt} is {err}, higher than the "
f"desired tolerance of {self.options['atol']}."
)

while krylov_error(dt * 10) < 0 and dt < self.options["max_step"]:
while krylov_error(dt * 10) < 0 and dt < max_step:
dt *= 10

if dt > self.options["max_step"]:
return self.options["max_step"]
if dt > max_step:
return max_step

sol = root_scalar(f=krylov_error, bracket=[dt, dt * 10],
method="brentq", xtol=self.options['atol'])
Expand Down
16 changes: 13 additions & 3 deletions qutip/tests/solver/test_sesolve.py
Expand Up @@ -301,9 +301,19 @@ def test_krylovsolve(always_compute_step):
e_op.dims = H.dims
tlist = np.linspace(0, 1, 11)
ref = sesolve(H, psi0, tlist, e_ops=[e_op]).expect[0]
options = {"always_compute_step", always_compute_step}
krylov_sol = krylovsolve(H, psi0, tlist, 20, e_ops=[e_op]).expect[0]
np.testing.assert_allclose(ref, krylov_sol)
options = {"always_compute_step": always_compute_step}
krylov_sol = krylovsolve(H, psi0, tlist, 20, e_ops=[e_op], options=options)
np.testing.assert_allclose(ref, krylov_sol.expect[0])


def test_krylovsolve_error():
H = qutip.rand_herm(256, density=0.2)
psi0 = qutip.basis([256], [255])
tlist = np.linspace(0, 1, 11)
options = {"min_step": 1e10}
with pytest.raises(ValueError) as err:
krylovsolve(H, psi0, tlist, 20, options=options)
assert "error with the minimum step" in str(err.value)


def test_feedback():
Expand Down

0 comments on commit 12bf7f7

Please sign in to comment.