def get(self, Wt, full=False):
        r'''Get Arnoldi relation for a deflation subspace choice.

        :param Wt: the coefficients :math:`\tilde{W}` of the deflation vectors
          in the basis :math:`[V_n,U]` with ``Wt.shape == (n+d, k)``, i.e., the
          deflation vectors are :math:`W=[V_n,U]\tilde{W}`. Must fulfill
          :math:`\tilde{W}^*\tilde{W}=I_k`.
        :param full: (optional) should the full Arnoldi
          basis and the full perturbation be returned? Defaults to ``False``.

        :return:

          * ``Hh``: the Hessenberg matrix with ``Hh.shape == (n+d-k, n+d-k)``.
          * ``Rh``: the perturbation core matrix with
            ``Rh.shape == (l, n+d-k)``.
          * ``q_norm``: norm :math:`\|q\|_2`.
          * ``vdiff_norm``: the norm of the difference of the initial vectors
            :math:`\tilde{v}-\hat{v}`.
          * ``PWAW_norm``: norm of the projection
            :math:`P_{\mathcal{W}^\perp,A\mathcal{W}}`.
          * ``Vh``: (if ``full == True``) the Arnoldi basis with ``Vh.shape ==
            (N, n+d-k)``.
          * ``F``: (if ``full == True``) the perturbation matrix
            :math:`F=-Z\hat{R}\hat{V}_n^* - \hat{V}_n\hat{R}^*Z^*`.
        '''
        n = self.n
        n_ = self.n_
        d = self.d
        k = Wt.shape[1]

        # get orthonormal basis of Wt and Wt^\perp
        if k > 0:
            Wto, _ = scipy.linalg.qr(Wt)
            Wt = Wto[:, :k]
            Wto = Wto[:, k:]
        else:
            Wto = numpy.eye(Wt.shape[0])

        deflated_solver = self._deflated_solver

        Pt = utils.Projection(self.L.dot(Wt), self.J.T.conj().dot(Wt)) \
            .operator_complement()
        if d > 0:
            qt = Pt*(numpy.r_[[[deflated_solver.MMlr0_norm]],
                              numpy.zeros((self.n_-1, 1)),
                              numpy.linalg.solve(deflated_solver.E,
                                                 deflated_solver.UMlr)])
        else:
            qt = Pt*(numpy.r_[[[deflated_solver.MMlr0_norm]],
                              numpy.zeros((self.n_-1, 1))])
        q = Wto.T.conj().dot(self.J.dot(qt))

        # TODO: q seems to suffer from round-off errors and thus the first
        #       vector in the computed basis may differ from the exact
        #       projected one more than on the level of unit round-off.

        # rotate closest vector in [V_n,U] to first column
        Q = utils.House(q)
        q_norm = Q.xnorm

        # Arnoldify
        WtoQ = Q.apply(Wto.T.conj()).T.conj()

        from scipy.linalg import hessenberg
        Hh, T = hessenberg(
            Q.apply(Wto.T.conj().dot(self.J).dot(Pt*(self.L.dot(WtoQ)))),
            calc_q=True)

        QT = Q.apply(T)

        # construct residual
        Rh = self.N.dot(Pt*self.L.dot(Wto.dot(QT)))

        # norm of difference between initial vectors
        vdiff = self.N.dot(qt)
        vdiff_norm = 0 if vdiff.size == 0 else numpy.linalg.norm(vdiff, 2)

        # compute norm of projection P_{W^\perp,AW}
        if k > 0:
            # compute coefficients of orthonormalized AW in the basis [V,Z]
            Y = numpy.bmat([[numpy.eye(n_), deflated_solver.B_],
                            [numpy.zeros((d, n_)), deflated_solver.E],
                            [numpy.zeros((self.R12.shape[0], n_)), self.R12]])
            YL_Q, _ = scipy.linalg.qr(Y.dot(self.L.dot(Wt)), mode='economic')

            # compute <W,X> where X is an orthonormal basis of AW
            WX = Wt.T.conj().dot(numpy.r_[YL_Q[:n, :],
                                          YL_Q[n_:n_+d, :]])
            PWAW_norm = 1./numpy.min(scipy.linalg.svdvals(WX))
        else:
            PWAW_norm = 1.

        if full:
            Vh = numpy.c_[deflated_solver.V[:, :n],
                          deflated_solver.projection.U].dot(Wto.dot(QT))
            ip_Minv_B = deflated_solver.linear_system.get_ip_Minv_B()

            def _apply_F(x):
                '''Application of the perturbation.'''
                return -(self.Z.dot(Rh.dot(utils.inner(Vh, x,
                                                       ip_B=ip_Minv_B)))
                         + Vh.dot(Rh.T.conj().dot(utils.inner(self.Z, x,
                                                              ip_B=ip_Minv_B)))
                         )
            F = utils.LinearOperator((Vh.shape[0], Vh.shape[0]),
                                     dtype=deflated_solver.dtype,
                                     dot=_apply_F
                                     )
            return Hh, Rh, q_norm, vdiff_norm, PWAW_norm, Vh, F
        return Hh, Rh, q_norm, vdiff_norm, PWAW_norm
def get ( self , Wt , full = False ) : n = self . n n_ = self . n_ d = self . d k = Wt . shape [ 1 ] if k > 0 : Wto , _ = scipy . linalg . qr ( Wt ) Wt = Wto [ : , : k ] Wto = Wto [ : , k : ] else : Wto = numpy . eye ( Wt . shape [ 0 ] ) deflated_solver = self . _deflated_solver Pt = utils . Projection ( self . L . dot ( Wt ) , self . J . T . conj ( ) . dot ( Wt ) ) . operator_complement ( ) if d > 0 : qt = Pt * ( numpy . r_ [ [ [ deflated_solver . MMlr0_norm ] ] , numpy . zeros ( ( self . n_ - 1 , 1 ) ) , numpy . linalg . solve ( deflated_solver . E , deflated_solver . UMlr ) ] ) else : qt = Pt * ( numpy . r_ [ [ [ deflated_solver . MMlr0_norm ] ] , numpy . zeros ( ( self . n_ - 1 , 1 ) ) ] ) q = Wto . T . conj ( ) . dot ( self . J . dot ( qt ) ) Q = utils . House ( q ) q_norm = Q . xnorm WtoQ = Q . apply ( Wto . T . conj ( ) ) . T . conj ( ) from scipy . linalg import hessenberg Hh , T = hessenberg ( Q . apply ( Wto . T . conj ( ) . dot ( self . J ) . dot ( Pt * ( self . L . dot ( WtoQ ) ) ) ) , calc_q = True ) QT = Q . apply ( T ) Rh = self . N . dot ( Pt * self . L . dot ( Wto . dot ( QT ) ) ) vdiff = self . N . dot ( qt ) vdiff_norm = 0 if vdiff . size == 0 else numpy . linalg . norm ( vdiff , 2 ) if k > 0 : Y = numpy . bmat ( [ [ numpy . eye ( n_ ) , deflated_solver . B_ ] , [ numpy . zeros ( ( d , n_ ) ) , deflated_solver . E ] , [ numpy . zeros ( ( self . R12 . shape [ 0 ] , n_ ) ) , self . R12 ] ] ) YL_Q , _ = scipy . linalg . qr ( Y . dot ( self . L . dot ( Wt ) ) , mode = 'economic' ) WX = Wt . T . conj ( ) . dot ( numpy . r_ [ YL_Q [ : n , : ] , YL_Q [ n_ : n_ + d , : ] ] ) PWAW_norm = 1. / numpy . min ( scipy . linalg . svdvals ( WX ) ) else : PWAW_norm = 1. if full : Vh = numpy . c_ [ deflated_solver . V [ : , : n ] , deflated_solver . projection . U ] . dot ( Wto . dot ( QT ) ) ip_Minv_B = deflated_solver . linear_system . get_ip_Minv_B ( ) def _apply_F ( x ) : return - ( self . Z . dot ( Rh . dot ( utils . inner ( Vh , x , ip_B = ip_Minv_B ) ) ) + Vh . dot ( Rh . T . conj ( ) . dot ( utils . inner ( self . Z , x , ip_B = ip_Minv_B ) ) ) ) F = utils . LinearOperator ( ( Vh . shape [ 0 ] , Vh . shape [ 0 ] ) , dtype = deflated_solver . dtype , dot = _apply_F ) return Hh , Rh , q_norm , vdiff_norm , PWAW_norm , Vh , F return Hh , Rh , q_norm , vdiff_norm , PWAW_norm
r'''Get Arnoldi relation for a deflation subspace choice.

        :param Wt: the coefficients :math:`\tilde{W}` of the deflation vectors
          in the basis :math:`[V_n,U]` with ``Wt.shape == (n+d, k)``, i.e., the
          deflation vectors are :math:`W=[V_n,U]\tilde{W}`. Must fulfill
          :math:`\tilde{W}^*\tilde{W}=I_k`.
        :param full: (optional) should the full Arnoldi
          basis and the full perturbation be returned? Defaults to ``False``.

        :return:

          * ``Hh``: the Hessenberg matrix with ``Hh.shape == (n+d-k, n+d-k)``.
          * ``Rh``: the perturbation core matrix with
            ``Rh.shape == (l, n+d-k)``.
          * ``q_norm``: norm :math:`\|q\|_2`.
          * ``vdiff_norm``: the norm of the difference of the initial vectors
            :math:`\tilde{v}-\hat{v}`.
          * ``PWAW_norm``: norm of the projection
            :math:`P_{\mathcal{W}^\perp,A\mathcal{W}}`.
          * ``Vh``: (if ``full == True``) the Arnoldi basis with ``Vh.shape ==
            (N, n+d-k)``.
          * ``F``: (if ``full == True``) the perturbation matrix
            :math:`F=-Z\hat{R}\hat{V}_n^* - \hat{V}_n\hat{R}^*Z^*`.