compute A0 and A1 in a single run

bin/beyn.py

@@ -13,14 +13,14 @@ import numpy  as np
 import numpy.matlib
 
 def simple(operator, origin, radius,
-           nodes=100, maxIt=10, lStart=1, lStep=1, rankTol=1e-4, verbose=True):
+           nNode=100, maxIt=10, lStart=1, lStep=1, rankTol=1e-4, verbose=True):
     """Solves an eigenvalue problem using Beyn's algorithm (simple version)
 
     Keyword arguments:
     operator -- the solver defining the operator to use
     origin -- the origin (in the complex plane) of the above circular contour
     radius -- the radius of the circular contour used to search the eigenvalues
-    nodes -- the number of nodes for the trapezoidal integration rule (optional)
+    nNode -- the number of nodes for the trapezoidal integration rule (optional)
     maxIt -- the maximal number of iteration for constructing A0 (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)
@@ -32,15 +32,14 @@ def simple(operator, origin, radius,
     """
 
     # Display the parameter used
-    if(verbose): display(nodes, maxIt, lStart, lStep, rankTol, origin, radius)
+    if(verbose): display(nNode, maxIt, lStart, lStep, rankTol, origin, radius)
 
     # Initialise A0 search
-    myPath = path(nodes, origin, radius)
-    hasK   = False
-    it     = 0
-    m      = operator.size()
-    l      = lStart
-    k      = -1
+    hasK = False
+    it   = 0
+    m    = operator.size()
+    l    = lStart
+    k    = -1
 
     # Search A0
     if(verbose): print "Searching A0..."
@@ -48,8 +47,9 @@ def simple(operator, origin, radius,
         if(verbose): print " # Iteration " + str(it+1) + ": ",
         if(verbose): sys.stdout.flush()
 
-        vHat = randomMatrix(m, l)                       # Take a random VHat
-        A0   = integrate(operator, myPath, 0, vHat)     # Compute A0
+        vHat   = randomMatrix(m, l)                     # Take a random VHat
+        A0, A1 = integrate(operator, vHat,              # Compute A_0 and A_1
+                           nNode, radius, origin)
         k, sMax, sMin = rank(A0, rankTol)               # Rank
 
         if(verbose): print       format(sMax, ".2e") + " |",
@@ -84,9 +84,8 @@ def simple(operator, origin, radius,
     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
+    # Compute B
+    B = V0.H * A1 * W0 * S0Inv
 
     # Eigenvalues & eigenvectors (by projecting QHat onto V0)
     if(verbose): print "Solving linear EVP..."
@@ -129,56 +128,37 @@ def rank(A, tol):
     return k, S[0].tolist(), S[nSV - 1].tolist()
 
 
-def path(nodes, origin, radius):
-    """Returns a list with the coordinates of a circular contour
+def integrate(operator, vHat, nNode, radius, origin):
+    """Computes the two first countour integrals of Beyn's method (A_0 and A_1)
+    over a circular contour.
 
     Keyword arguments:
-    nodes -- the number of nodes used to discretise the contour
+    operator -- the solver defining the operator to use
+    vHat -- the RHS matrix defing Beyn's integrals
+    nNode -- the number of nodes used to discretise the circular 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 integrals
+    A0 = np.matlib.zeros(vHat.shape, dtype=complex)
+    A1 = np.matlib.zeros(vHat.shape, dtype=complex)
 
-    # Initialise integration loop
-    F1  = multiSolve(operator, vHat, path[0])
-    F1 *= np.power(path[0], order)
+    # Loop over integation points and integrate
+    for i in range(nNode):
+       t   = 2 * np.pi * i / nNode;
+       phi = origin + radius * np.exp(1j * t)
+       tmp = multiSolve(operator, vHat, phi)
 
-    # 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)
+       A0 += tmp * np.power(phi, 0) * np.exp(1j * t)
+       A1 += tmp * np.power(phi, 1) * np.exp(1j * t)
 
-        I += (F1 + F2) * (path[i + 1] - path[i])
-        F1 = F2
+    # Final scale
+    A0 *= radius/nNode
+    A1 *= radius/nNode
 
     # Done
-    return I
+    return A0, A1
 
 
 def multiSolve(solver, B, w):
@@ -212,10 +192,10 @@ def randomMatrix(n, m):
     return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
 
 
-def display(nodes, maxIt, lStart, lStep, rankTol, origin, radius):
+def display(nNode, maxIt, lStart, lStep, rankTol, origin, radius):
     print "Beyn's contour integral method (simple)"
     print "---------------------------------------"
-    print " # Nodes used for the trapezoidal rule:" + " " + str(nodes)
+    print " # Nodes used for the trapezoidal rule:" + " " + str(nNode)
     print " # Maximum number of iterations:       " + " " + str(maxIt)
     print " # Initial size of col(A0):            " + " " + str(lStart)
     print " # Step size for increasing col(A0):   " + " " + str(lStep)