Compute (perturbed) Arnoldi relation for projected operator.

    Assume that you have computed an Arnoldi relation

    .. math ::

        A V_n = V_{n+1} \\underline{H}_n

    where :math:`V_{n+1}\\in\\mathbb{C}^{N,n+1}` has orthogonal columns
    (with respect to an inner product :math:`\\langle\\cdot,\\cdot\\rangle`)
    and :math:`\\underline{H}_n\\in\\mathbb{C}^{n+1,n}` is an extended
    upper Hessenberg matrix.

    For :math:`k<n` you choose full rank matrices
    :math:`X\\in\\mathbb{C}^{n-1,k}` and :math:`Y\\in\\mathbb{C}^{n,k}` and
    define :math:`\\tilde{X}:=A V_{n_1}X = V_n \\underline{H}_{n-1} X` and
    :math:`\\tilde{Y}:=V_n Y` such that
    :math:`\\langle \\tilde{Y}, \\tilde{X} \\rangle = Y^*\\underline{H}_{n-1} X`
    is invertible. Then the projections :math:`P` and :math:`\\tilde{P}`
    characterized by

    * :math:`\\tilde{P}x = x -
      \\tilde{X} \\langle \\tilde{Y},\\tilde{X} \\rangle^{-1}
      \\langle\\tilde{Y},x\\rangle`
    * :math:`P = I - \\underline{H}_{n-1}X (Y^*\\underline{H}_{n-1}X)^{-1}Y^*`

    are well defined and :math:`\\tilde{P}V_{n+1} = [V_n P, v_{n+1}]` holds.

    This method computes for :math:`i<n-k` the Arnoldi relation

    .. math ::

        (\\tilde{P}A + E_i) W_i
        = W_{i+1} \\underline{G}_i

    where :math:`W_{i+1}=V_n U_{i+1}` has orthogonal columns with respect
    to :math:`\\langle\\cdot,\\cdot\\rangle`,
    :math:`\\underline{G}_i` is an extended upper Hessenberg matrix
    and :math:`E_i x = v_{n+1} F_i \\langle W_i,x\\rangle` with
    :math:`F_i=[f_1,\ldots,f_i]\\in\\mathbb{C}^{1,i}`.

    The perturbed Arnoldi relation can also be generated with the operator
    :math:`P_{V_n} \\tilde{P} A`:

    .. math ::

        P_{V_n} \\tilde{P} A W_i
        = W_{i+1} \\underline{G}_i.

    In a sense the perturbed Arnoldi relation is the best prediction for the
    behavior of the Krylov subspace :math:`K_i(\\tilde{P}A,\\tilde{P}v_1)`
    that can be generated only with the data from :math:`K_{n+1}(A,v_1)` and
    without carrying out further matrix-vector multiplications with A.

    :param H: the extended upper Hessenberg matrix
      :math:`\\underline{H}_n` with ``shape==(n+1,n)``.
    :param P: the projection
      :math:`P:\\mathbb{C}^n\\longrightarrow\\mathbb{C}^n` (has to be
      compatible with :py:meth:`get_linearoperator`).
    :param k: the dimension of the null space of P.
    :returns: U, G, F where

      * U is the coefficient matrix :math:`U_{i+1}` with ``shape==(n,i+1)``,
      * G is the extended upper Hessenberg matrix :math:`\\underline{G}_i`
        with ``shape==(i+1,i)``,
      * F is the error matrix :math:`F_i` with ``shape==(1,i)``.

def arnoldi_projected(H, P, k, ortho='mgs'):
    """Compute (perturbed) Arnoldi relation for projected operator.

    Assume that you have computed an Arnoldi relation

    .. math ::

        A V_n = V_{n+1} \\underline{H}_n

    where :math:`V_{n+1}\\in\\mathbb{C}^{N,n+1}` has orthogonal columns
    (with respect to an inner product :math:`\\langle\\cdot,\\cdot\\rangle`)
    and :math:`\\underline{H}_n\\in\\mathbb{C}^{n+1,n}` is an extended
    upper Hessenberg matrix.

    For :math:`k<n` you choose full rank matrices
    :math:`X\\in\\mathbb{C}^{n-1,k}` and :math:`Y\\in\\mathbb{C}^{n,k}` and
    define :math:`\\tilde{X}:=A V_{n_1}X = V_n \\underline{H}_{n-1} X` and
    :math:`\\tilde{Y}:=V_n Y` such that
    :math:`\\langle \\tilde{Y}, \\tilde{X} \\rangle = Y^*\\underline{H}_{n-1} X`
    is invertible. Then the projections :math:`P` and :math:`\\tilde{P}`
    characterized by

    * :math:`\\tilde{P}x = x -
      \\tilde{X} \\langle \\tilde{Y},\\tilde{X} \\rangle^{-1}
      \\langle\\tilde{Y},x\\rangle`
    * :math:`P = I - \\underline{H}_{n-1}X (Y^*\\underline{H}_{n-1}X)^{-1}Y^*`

    are well defined and :math:`\\tilde{P}V_{n+1} = [V_n P, v_{n+1}]` holds.

    This method computes for :math:`i<n-k` the Arnoldi relation

    .. math ::

        (\\tilde{P}A + E_i) W_i
        = W_{i+1} \\underline{G}_i

    where :math:`W_{i+1}=V_n U_{i+1}` has orthogonal columns with respect
    to :math:`\\langle\\cdot,\\cdot\\rangle`,
    :math:`\\underline{G}_i` is an extended upper Hessenberg matrix
    and :math:`E_i x = v_{n+1} F_i \\langle W_i,x\\rangle` with
    :math:`F_i=[f_1,\ldots,f_i]\\in\\mathbb{C}^{1,i}`.

    The perturbed Arnoldi relation can also be generated with the operator
    :math:`P_{V_n} \\tilde{P} A`:

    .. math ::

        P_{V_n} \\tilde{P} A W_i
        = W_{i+1} \\underline{G}_i.

    In a sense the perturbed Arnoldi relation is the best prediction for the
    behavior of the Krylov subspace :math:`K_i(\\tilde{P}A,\\tilde{P}v_1)`
    that can be generated only with the data from :math:`K_{n+1}(A,v_1)` and
    without carrying out further matrix-vector multiplications with A.

    :param H: the extended upper Hessenberg matrix
      :math:`\\underline{H}_n` with ``shape==(n+1,n)``.
    :param P: the projection
      :math:`P:\\mathbb{C}^n\\longrightarrow\\mathbb{C}^n` (has to be
      compatible with :py:meth:`get_linearoperator`).
    :param k: the dimension of the null space of P.
    :returns: U, G, F where

      * U is the coefficient matrix :math:`U_{i+1}` with ``shape==(n,i+1)``,
      * G is the extended upper Hessenberg matrix :math:`\\underline{G}_i`
        with ``shape==(i+1,i)``,
      * F is the error matrix :math:`F_i` with ``shape==(1,i)``.
    """
    n = H.shape[1]
    dtype = find_common_dtype(H, P)
    invariant = H.shape[0] == n
    hlast = 0 if invariant else H[-1, -1]
    H = get_linearoperator((n, n), H if invariant else H[:-1, :])
    P = get_linearoperator((n, n), P)
    v = P * numpy.eye(n, 1)
    maxiter = n-k+1
    F = numpy.zeros((1, maxiter), dtype=dtype)
    PH = lambda x: P*(H*x)
    PH = LinearOperator((n, n), dtype, PH)
    _arnoldi = Arnoldi(PH, v, maxiter=maxiter, ortho=ortho)
    while _arnoldi.iter < _arnoldi.maxiter and not _arnoldi.invariant:
        u, _ = _arnoldi.get_last()
        F[0, _arnoldi.iter] = hlast * u[-1, 0]
        _arnoldi.advance()
    U, G = _arnoldi.get()
    return U, G, F[[0], :_arnoldi.iter]