(ctx, n0, n1, A, Q, shift)
This subroutine executes a single implicitly shifted QR step applied to an
    upper Hessenberg matrix A. Given A and shift as input, first an QR
    decomposition is calculated:

      Q R = A - shift * 1 .

    The output is then following matrix:

      R Q + shift * 1

    parameters:
      n0, n1    (input) Two integers which specify the submatrix A[n0:n1,n0:n1]
                on which this subroutine operators. The subdiagonal elements
                to the left and below this submatrix must be deflated (i.e. zero).
                following restriction is imposed: n1>=n0+2
      A         (input/output) On input, A is an upper Hessenberg matrix.
                On output, A is replaced by "R Q + shift * 1"
      Q         (input/output) The parameter Q is multiplied by the unitary matrix
                Q arising from the QR decomposition. Q can also be false, in which
                case the unitary matrix Q is not computated.
      shift     (input) a complex number specifying the shift. idealy close to an
                eigenvalue of the bottemmost part of the submatrix A[n0:n1,n0:n1].

    references:
      Stoer, Bulirsch - Introduction to Numerical Analysis.
      Kresser : Numerical Methods for General and Structured Eigenvalue Problems
This subroutine executes a single implicitly shifted QR step applied to an
    upper Hessenberg matrix A. Given A and shift as input, first an QR
    decomposition is calculated:
This subroutine executes a single implicitly shifted QR step applied to an upper Hessenberg matrix A . Given A and shift as input first an QR decomposition is calculated :
def qr_step(ctx, n0, n1, A, Q, shift):
    """
    This subroutine executes a single implicitly shifted QR step applied to an
    upper Hessenberg matrix A. Given A and shift as input, first an QR
    decomposition is calculated:

      Q R = A - shift * 1 .

    The output is then following matrix:

      R Q + shift * 1

    parameters:
      n0, n1    (input) Two integers which specify the submatrix A[n0:n1,n0:n1]
                on which this subroutine operators. The subdiagonal elements
                to the left and below this submatrix must be deflated (i.e. zero).
                following restriction is imposed: n1>=n0+2
      A         (input/output) On input, A is an upper Hessenberg matrix.
                On output, A is replaced by "R Q + shift * 1"
      Q         (input/output) The parameter Q is multiplied by the unitary matrix
                Q arising from the QR decomposition. Q can also be false, in which
                case the unitary matrix Q is not computated.
      shift     (input) a complex number specifying the shift. idealy close to an
                eigenvalue of the bottemmost part of the submatrix A[n0:n1,n0:n1].

    references:
      Stoer, Bulirsch - Introduction to Numerical Analysis.
      Kresser : Numerical Methods for General and Structured Eigenvalue Problems
    """

    # implicitly shifted and bulge chasing is explained at p.398/399 in "Stoer, Bulirsch - Introduction to Numerical Analysis"
    # for bulge chasing see also "Watkins - The Matrix Eigenvalue Problem" sec.4.5,p.173

    # the Givens rotation we used is determined as follows: let c,s be two complex
    # numbers. then we have following relation:
    #
    #     v = sqrt(|c|^2 + |s|^2)
    #
    #     1/v [ c~  s~]  [c] = [v]
    #         [-s   c ]  [s]   [0]
    #
    # the matrix on the left is our Givens rotation.

    n = A.rows

    # first step

    # calculate givens rotation
    c = A[n0  ,n0] - shift
    s = A[n0+1,n0]

    v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s)))

    if v == 0:
        v = 1
        c = 1
        s = 0
    else:
        c /= v
        s /= v

    for k in xrange(n0, n):
        # apply givens rotation from the left
        x = A[n0  ,k]
        y = A[n0+1,k]
        A[n0  ,k] = ctx.conj(c) * x + ctx.conj(s) * y
        A[n0+1,k] =         -s  * x +          c  * y

    for k in xrange(min(n1, n0+3)):
        # apply givens rotation from the right
        x = A[k,n0  ]
        y = A[k,n0+1]
        A[k,n0  ] =           c  * x +          s  * y
        A[k,n0+1] = -ctx.conj(s) * x + ctx.conj(c) * y

    if not isinstance(Q, bool):
        for k in xrange(n):
            # eigenvectors
            x = Q[k,n0  ]
            y = Q[k,n0+1]
            Q[k,n0  ] =           c  * x +          s  * y
            Q[k,n0+1] = -ctx.conj(s) * x + ctx.conj(c) * y

    # chase the bulge

    for j in xrange(n0, n1 - 2):
        # calculate givens rotation

        c = A[j+1,j]
        s = A[j+2,j]

        v = ctx.hypot(ctx.hypot(ctx.re(c), ctx.im(c)), ctx.hypot(ctx.re(s), ctx.im(s)))

        if v == 0:
            A[j+1,j] = 0
            v = 1
            c = 1
            s = 0
        else:
            A[j+1,j] = v
            c /= v
            s /= v

        A[j+2,j] = 0

        for k in xrange(j+1, n):
            # apply givens rotation from the left
            x = A[j+1,k]
            y = A[j+2,k]
            A[j+1,k] = ctx.conj(c) * x + ctx.conj(s) * y
            A[j+2,k] =         -s  * x +          c  * y

        for k in xrange(0, min(n1, j+4)):
            # apply givens rotation from the right
            x = A[k,j+1]
            y = A[k,j+2]
            A[k,j+1] =           c  * x +          s  * y
            A[k,j+2] = -ctx.conj(s) * x + ctx.conj(c) * y

        if not isinstance(Q, bool):
            for k in xrange(0, n):
                # eigenvectors
                x = Q[k,j+1]
                y = Q[k,j+2]
                Q[k,j+1] =           c  * x +          s  * y
                Q[k,j+2] = -ctx.conj(s) * x + ctx.conj(c) * y
def qr_step ( ctx , n0 , n1 , A , Q , shift ) : # implicitly shifted and bulge chasing is explained at p.398/399 in "Stoer, Bulirsch - Introduction to Numerical Analysis" # for bulge chasing see also "Watkins - The Matrix Eigenvalue Problem" sec.4.5,p.173 # the Givens rotation we used is determined as follows: let c,s be two complex # numbers. then we have following relation: # #     v = sqrt(|c|^2 + |s|^2) # #     1/v [ c~  s~]  [c] = [v] #         [-s   c ]  [s]   [0] # # the matrix on the left is our Givens rotation. n = A . rows # first step # calculate givens rotation c = A [ n0 , n0 ] - shift s = A [ n0 + 1 , n0 ] v = ctx . hypot ( ctx . hypot ( ctx . re ( c ) , ctx . im ( c ) ) , ctx . hypot ( ctx . re ( s ) , ctx . im ( s ) ) ) if v == 0 : v = 1 c = 1 s = 0 else : c /= v s /= v for k in xrange ( n0 , n ) : # apply givens rotation from the left x = A [ n0 , k ] y = A [ n0 + 1 , k ] A [ n0 , k ] = ctx . conj ( c ) * x + ctx . conj ( s ) * y A [ n0 + 1 , k ] = - s * x + c * y for k in xrange ( min ( n1 , n0 + 3 ) ) : # apply givens rotation from the right x = A [ k , n0 ] y = A [ k , n0 + 1 ] A [ k , n0 ] = c * x + s * y A [ k , n0 + 1 ] = - ctx . conj ( s ) * x + ctx . conj ( c ) * y if not isinstance ( Q , bool ) : for k in xrange ( n ) : # eigenvectors x = Q [ k , n0 ] y = Q [ k , n0 + 1 ] Q [ k , n0 ] = c * x + s * y Q [ k , n0 + 1 ] = - ctx . conj ( s ) * x + ctx . conj ( c ) * y # chase the bulge for j in xrange ( n0 , n1 - 2 ) : # calculate givens rotation c = A [ j + 1 , j ] s = A [ j + 2 , j ] v = ctx . hypot ( ctx . hypot ( ctx . re ( c ) , ctx . im ( c ) ) , ctx . hypot ( ctx . re ( s ) , ctx . im ( s ) ) ) if v == 0 : A [ j + 1 , j ] = 0 v = 1 c = 1 s = 0 else : A [ j + 1 , j ] = v c /= v s /= v A [ j + 2 , j ] = 0 for k in xrange ( j + 1 , n ) : # apply givens rotation from the left x = A [ j + 1 , k ] y = A [ j + 2 , k ] A [ j + 1 , k ] = ctx . conj ( c ) * x + ctx . conj ( s ) * y A [ j + 2 , k ] = - s * x + c * y for k in xrange ( 0 , min ( n1 , j + 4 ) ) : # apply givens rotation from the right x = A [ k , j + 1 ] y = A [ k , j + 2 ] A [ k , j + 1 ] = c * x + s * y A [ k , j + 2 ] = - ctx . conj ( s ) * x + ctx . conj ( c ) * y if not isinstance ( Q , bool ) : for k in xrange ( 0 , n ) : # eigenvectors x = Q [ k , j + 1 ] y = Q [ k , j + 2 ] Q [ k , j + 1 ] = c * x + s * y Q [ k , j + 2 ] = - ctx . conj ( s ) * x + ctx . conj ( c ) * y