Skip to content
Snippets Groups Projects
Commit 99730fba authored by Nicolas Marsic's avatar Nicolas Marsic
Browse files

getdp: pretty print (again)

parent 4c1a2ac5
No related branches found
No related tags found
No related merge requests found
......@@ -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: " +
......
......@@ -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]):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment