diff --git a/beyn.py b/beyn.py index 1e97763659aea0cad803afa192cd59c4b4277032..45606f35a5ed79a2057655adf98841cd3944c1aa 100644 --- a/beyn.py +++ b/beyn.py @@ -2,22 +2,26 @@ 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): +def simple(operator, origin, radius, + nodes=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 - radius -- the radius of the circular contour used to search the eigenvalues 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) 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 + verbose -- should I be verbose? (optional) Returns the computed eigenvalues """ + # Display the parameter used + if(verbose): display(nodes, maxIt, lStart, lStep, rankTol, origin, radius) + # Initialise A0 search myPath = path(nodes, origin, radius) hasK = False @@ -27,8 +31,9 @@ def simple(operator, radius, origin, k = -1 # Search A0 + if(verbose): print "Searching A0..." while(not hasK and it != maxIt): - print "Iteration: " + str(it+1) + if(verbose): print " # Iteration: " + str(it+1) vHat = randomMatrix(m ,l) # Take a random VHat A0 = integrate(operator, myPath, 0, vHat) # Compute A0 @@ -49,7 +54,7 @@ def simple(operator, radius, origin, if(it == maxIt): raise RuntimeError(maxItErrorStr()) else: - print "Constructing linear EVP..." + if(verbose): print "Constructing linear EVP..." # Compute V, S and Wh # NB: For SVD(A) = V*S*Wh, numpy computes {v, s, w}, such that: @@ -66,10 +71,11 @@ def simple(operator, radius, origin, B = V0.H * A1 * W0 * S0Inv # Eigenvalues of B - print "Solving linear EVP..." + if(verbose): print "Solving linear EVP..." myLambda, QHat = numpy.linalg.eig(B) # Done + if(verbose): print "Done!" return myLambda @@ -161,6 +167,25 @@ 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): + print "Beyn's contour integral method (simple)" + print "---------------------------------------" + print " # Nodes used for the trapezoidal rule:" + " " + str(nodes) + print " # Maximum number of iterations: " + " " + str(maxIt) + print " # Initial size of col(A0): " + " " + str(lStart) + print " # Step size for col(A0): " + " " + str(lStep) + print " # Rank test tolerance: ", + print format(rankTol, '.2e') + print "---------------------------------------" + print " # Cirular path origin: ", + print "(" + format(np.real(origin).tolist(), '+.2e') + ")", + print "+", + print "(" + format(np.imag(origin).tolist(), '+.2e') + ")j" + print " # Cirular path radius: ", + print format(radius, '+.2e') + print "---------------------------------------" + + def zeroRankErrorStr(): """Returns a string explaining the probable reason of a zero rank""" return ("Found a rank of zero: " + diff --git a/main.py b/main.py index 96c9833866d53ece3b5f38383c522ebbc46b5d9b..cbb1e2551edcd4ada4fa645db5505ddf7cebfe62 100644 --- a/main.py +++ b/main.py @@ -11,8 +11,9 @@ resolution = "Maxwell" operator = Solver.GetDPWave(pro, mesh, resolution) # Compute -v = Beyn.simple(operator, 1e8, complex(9e8, 0)) +v = Beyn.simple(operator, complex(9e8, 0), 1e8) +print print "Eigenvalues:" print "----------- " for i in range(v.shape[0]):