Source code for pymor.algorithms.krylov

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

"""Module for computing (rational) Krylov subspaces' bases."""

from pymor.algorithms.gram_schmidt import gram_schmidt


[docs]def rational_arnoldi(A, E, b, sigma, trans=False): r"""Rational Arnoldi algorithm. If `trans == False`, using Arnoldi process, computes a real orthonormal basis for the rational Krylov subspace .. math:: \mathrm{span}\{ (\sigma_1 E - A)^{-1} b, (\sigma_2 E - A)^{-1} b, \ldots, (\sigma_r E - A)^{-1} b \}, otherwise, computes the same for .. math:: \mathrm{span}\{ (\sigma_1 E - A)^{-T} b^T, (\sigma_2 E - A)^{-T} b^T, \ldots, (\sigma_r E - A)^{-T} b^T \}. Interpolation points in `sigma` are allowed to repeat (in any order). Then, in the above expression, .. math:: \underbrace{ (\sigma_i E - A)^{-1} b, \ldots, (\sigma_i E - A)^{-1} b }_{m \text{ times}} is replaced by .. math:: (\sigma_i E - A)^{-1} b, (\sigma_i E - A)^{-1} E (\sigma_i E - A)^{-1} b, \ldots, \left((\sigma_i E - A)^{-1} E\right)^{m - 1} (\sigma_i E - A)^{-1} b. Analogously for the `trans == True` case. Parameters ---------- A Real |Operator| A. E Real |Operator| E. b Real vector-like operator (if trans is False) or functional (if trans is True). sigma Sequence of interpolation points (closed under conjugation). trans Boolean, see above. Returns ------- V Orthonormal basis for the Krylov subspace |VectorArray|. """ assert A.source == A.range assert E.source == A.source assert E.range == A.source assert (b.range if not trans else b.source) == A.source assert not trans and b.source.dim == 1 or trans and b.range.dim == 1 r = len(sigma) V = A.source.empty(reserve=r) v = b.as_vector() v.scal(1 / v.l2_norm()[0]) for i in range(r): if sigma[i].imag < 0: continue if sigma[i].imag == 0: sEmA = sigma[i].real * E - A else: sEmA = sigma[i] * E - A if not trans: v = sEmA.apply_inverse(v if len(V) == 0 else E.apply(v)) else: v = sEmA.apply_inverse_adjoint(v if len(V) == 0 else E.apply_adjoint(v)) if sigma[i].imag == 0: V.append(v) gram_schmidt(V, atol=0, rtol=0, offset=len(V) - 1, copy=False) else: V.append(v.real) V.append(v.imag) gram_schmidt(V, atol=0, rtol=0, offset=len(V) - 2, copy=False) v = V[-1] return V
[docs]def tangential_rational_krylov(A, E, B, b, sigma, trans=False, orth=True): r"""Tangential Rational Krylov subspace. If `trans == False`, computes a real basis for the rational Krylov subspace .. math:: \mathrm{span}\{ (\sigma_1 E - A)^{-1} B b_1, (\sigma_2 E - A)^{-1} B b_2, \ldots, (\sigma_r E - A)^{-1} B b_r \}, otherwise, computes the same for .. math:: \mathrm{span}\{ (\sigma_1 E - A)^{-T} B^T b_1, (\sigma_2 E - A)^{-T} B^T b_2, \ldots, (\sigma_r E - A)^{-T} B^T b_r \}. Interpolation points in `sigma` are assumed to be pairwise distinct. Parameters ---------- A Real |Operator| A. E Real |Operator| E. B Real |Operator| B. b |VectorArray| from `B.source`, if `trans == False`, or `B.range`, if `trans == True`. sigma Sequence of interpolation points (closed under conjugation), of the same length as `b`. trans Boolean, see above. orth If `True`, orthonormalizes the basis using :meth:`pymor.algorithms.gram_schmidt.gram_schmidt`. Returns ------- V Optionally orthonormal basis for the Krylov subspace |VectorArray|. """ assert A.source == A.range assert E.source == A.source assert E.range == A.source assert (B.range if not trans else B.source) == A.source assert b in (B.source if not trans else B.range) assert len(b) == len(sigma) r = len(sigma) V = A.source.empty(reserve=r) for i in range(r): if sigma[i].imag == 0: sEmA = sigma[i].real * E - A if not trans: Bb = B.apply(b.real[i]) V.append(sEmA.apply_inverse(Bb)) else: BTb = B.apply_adjoint(b.real[i]) V.append(sEmA.apply_inverse_adjoint(BTb)) elif sigma[i].imag > 0: sEmA = sigma[i] * E - A if not trans: Bb = B.apply(b[i]) v = sEmA.apply_inverse(Bb) else: BTb = B.apply_adjoint(b[i].conj()) v = sEmA.apply_inverse_adjoint(BTb) V.append(v.real) V.append(v.imag) if orth: gram_schmidt(V, atol=0, rtol=0, copy=False) return V