beyn.py

import solver as Solver
import numpy  as np
import numpy.matlib

def simple(operator, radius, origin,
           nodes=100, maxIt=10, lStart=1, lStep=1, rankTol=1e-4):
    """Solves an eigenvalue problem using Beyn's algorithm (simple version)

    Keyword arguments:
    operator -- the solver defining the operator to use
    radius -- the radius of the circular contour used to search the eigenvalues
    origin -- the origin (in the complex plane) of the above circular contour
    nodes -- the number of nodes for the trapezoidal integration rule (optional)
    lStart -- the number of columns used for A0 when algorithm starts (optional)
    lStep -- the step used for increasing the number of columns of A0 (optional)
    rankTol -- the tolerance for the rank test

    Returns the computed eigenvalues
    """

    # Initialise A0 search
    myPath = path(nodes, origin, radius)
    hasK   = False
    it     = 0
    m      = operator.size()
    l      = lStart
    k      = -1

    # Search A0
    while(not hasK and it != maxIt):
        print "Iteration: " + str(it+1)

        vHat = randomMatrix(m ,l)                       # Take a random VHat
        A0   = integrate(operator, myPath, 0, vHat)     # Compute A0
        k    = np.linalg.matrix_rank(A0, tol=rankTol)   # Rank test

        if(k == 0):                                     # Rank is zero?
            raise RuntimeError(zeroRankErrorStr())      #  -> Stop
        elif(k == m):                                   # Maximum rank reached?
            raise RuntimeError(maxRankErrorStr())       #  -> Stop
        elif(k < l):                                    # Found a null SV?
            hasK = True                                 #  -> We have A0
        else:                                           # Matrix is full rank?
            l = l + lStep                               #  -> Increase A0 size

        it += 1                                         # Keep on searching A0

    # Check if maxIt was reached
    if(it == maxIt):
        raise RuntimeError(maxItErrorStr())
    else:
        print "Constructing linear EVP..."

    # Compute V, S and Wh
    #  NB: For SVD(A) = V*S*Wh, numpy computes {v, s, w}, such that:
    #      v = V; diag(s) = S and w = Wh
    V, S, Wh = np.linalg.svd(A0, full_matrices=False, compute_uv=1)

    # Extract V0, W0 and S0Inv
    V0    = np.delete(V,  l-1, 1)
    W0    = np.delete(Wh, l-1, 0).H
    S0Inv = np.matrix(np.diag(1/np.delete(S, l-1, 0)))

    # Compute A1 and B
    A1 = integrate(operator, myPath, 1, vHat)
    B  = V0.H * A1 * W0 * S0Inv

    # Eigenvalues of B
    print "Solving linear EVP..."
    myLambda, QHat = numpy.linalg.eig(B)

    # Done
    return myLambda


## Import only simple (other functions are just helpers)
__all__ = ['simple']


## Helper functions
def path(nodes, origin, radius):
    """Returns a list with the coordinates of a circular contour

    Keyword arguments:
    nodes -- the number of nodes used to discretise the contour
    radius -- the radius of the circular contour
    origin -- the origin (in the complex plane) of the circular contour

    The returned list contains nodes+1 elements,
    such that the first and the last are identical (up the machine precision)
    """

    step         = 1.0 / nodes
    nodesPlusOne = nodes + 1

    path = list()
    for i in range(nodesPlusOne):
        path.append(origin + radius * np.exp(1j * 2 * np.pi * i * step))

    return path


def integrate(operator, path, order, vHat):
    """Computes the countour integral of Beyn's method, that is matrix A_p

    Keyword arguments:
    operator -- the solver defining the operator to use
    path -- the path to integrate on
    order -- the order of Beyn's integral (that is the 'p' defining matrix A_p)
    vHat -- the RHS matrix defing Beyn's integral
    """

    # Initialise I
    I = np.matlib.zeros(vHat.shape, dtype=complex)

    # Initialise integration loop
    F1  = multiSolve(operator, vHat, path[0])
    F1 *= np.power(path[0], order)

    # Integration loop
    pathSizeMinus = len(path) - 1;
    for i in range(pathSizeMinus):
        F2  = multiSolve(operator, vHat, path[i + 1])
        F2 *= np.power(path[i + 1], order)

        I += (F1 + F2) * (path[i + 1] - path[i])
        F1 = F2

    # Done
    return I


def multiSolve(solver, B, w):
    """Solves for multiple RHS

    Keyword arguments:
    solver -- the solver to use
    B -- the matrix of RHS
    w -- the complex frequency to use
    """

    # Number of solve
    nSolve = B.shape[1]

    # Initialise X
    size = solver.size()
    X    = np.matlib.empty((size, nSolve), dtype=complex)

    # Loop and solve
    for i in range(nSolve):
        b = B[:, i]
        solver.solve(b, w)
        X[:, i] = solver.solution()

    # Done
    return X


def randomMatrix(n, m):
    """Returns a random complex matrix of the given size (n x m)"""
    return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j


def zeroRankErrorStr():
    """Returns a string explaining the probable reason of a zero rank"""
    return ("Found a rank of zero: " +
            "the contour is probably enclosing no eigenvalues")


def maxRankErrorStr():
    """Returns a string explaining the probable reason of a maximal rank"""
    return ("Maximal rank found: " +
            "the contour is probably enclosing to many eigenvalues")


def maxItErrorStr():
    """Returns a string claiming: the maximum number of iterations is reached"""
    return "Maximum number iterations is reached!"