[PYTHON MATRIX-SIG] Eigenvectors

Doug Heisterkamp drh@oasis.unl.edu
Tue, 13 Aug 1996 08:47:54 -0500 (CDT)


Here's an eigenvector function to add to the updated LinAlg.py module that
Janne Sinkkonen posted a few days ago.  If anyone has other LAPACK access
functions and wishes to share them, email them to me and I'll merge them
in with LinAlg.py.   

Doug Heisterkamp

# Eigenvectors

def eigenvectors(a):
    """eigenvectors(a) returns u,v  where u is the eigenvalues and
v is a matrix of eigenvectors with vector v[i] corresponds to 
eigenvalue u[i].  The eigenvectors are the rows of the matrix v. 
Satisfies the equation dot(a,v[i]) = u[i]*v[i]
    t =_commonType(a)
    real_t = _array_type[0][_array_precision[t]]
    lapack_routine = _findLapackRoutine('geev', t)
    a = _castCopyAndTranspose(t, a)
    n = a.shape[0]
    dummy = Numeric.zeros((1,), t)
    lwork = max(1,8*n) # minimum value: max(1,3*n) real, max(1,2*n) complex
    work = Numeric.zeros((lwork,), t)
    if _array_kind[t] == 1: # Complex routines take different arguments
        w = Numeric.zeros((n,), t)
        rwork = Numeric.zeros((n,),real_t)
        v = Numeric.zeros((n,n), t)
        results = lapack_routine('N', 'V', n, a, n, w,
                                  dummy, 1, v, n, work, lwork, rwork, 0)
        wr = Numeric.zeros((n,), t)
        wi = Numeric.zeros((n,), t)
        vr = Numeric.zeros((n,n), t)
        results = lapack_routine('N', 'V', n, a, n, wr, wi,
                                  dummy, 1, vr, n, work, lwork, 0)
        if Numeric.logicalAnd.reduce(Numeric.equal(wi, 0.)):
            w = wr
            v = vr
            w = wr+1j*wi
            v = Numeric.array(vr,Numeric.Complex())
            ind = Numeric.nonzero(
                              Numeric.equal(wi,0.0) # true for real e-vals
                                       ,0)          # true for complex e-vals
                                 )                  # indices of complex e-vals
            for i in range(len(ind)/2):
                print "processing i=",i
                v[ind[2*i]] = vr[ind[2*i]] + 1j*vr[ind[2*i+1]]
                v[ind[2*i+1]] = vr[ind[2*i]] - 1j*vr[ind[2*i+1]]
                print "New v = ",v
    if results['info'] > 0:
        raise LinAlgError, 'Eigenvalues did not converge'
    return w,v


MATRIX-SIG  - SIG on Matrix Math for Python

send messages to: matrix-sig@python.org
administrivia to: matrix-sig-request@python.org