Function SOAR

python/eigenvalues.py:201–239

(A, B, maxiter=200)

# author: A. Schoefl

def SOAR (A, B, maxiter=200):
# first version without memory saving and breakdown 

    q = A.CreateVector()
    p = A.CreateVector()
    s = A.CreateVector()
    r = A.CreateVector()

    Q = MultiVector(q,0)
    P = MultiVector(p,0)
    p[:] = 0
    q.SetRandom()
    q.FV().imag = 0      
    q /= Norm(q)
    T = Matrix(maxiter, complex=A.is_complex)
    T[:,:] = 0
    
    for j in range(maxiter):
        Q.Append(q)
        P.Append(p)
        r.data = A*q + B*p
        s.data = q

        for i in range(j+1):

            T[i,j] = InnerProduct(r, Q[i])
            r -= T[i,j]*Q[i]
            s -= T[i,j]*P[i]
        
        if j+1 < maxiter:
            T[j+1,j] = Norm(r)
            if T[j+1,j] == 0:
                print("SOAR stopped at iteration j = ", j)
                break

            q.data = 1/T[j+1,j]*r
            p.data = 1/T[j+1,j]*s

    return Q