import numpy as np
import laff
import flame

Lambda = np.matrix( ' 4., 0., 0., 0;\
                      0., 3., 0., 0;\
                      0., 0., 2., 0;\
                      0., 0., 0., 1' )

lambda0 = Lambda[ 0,0 ]

V = np.matrix( np.random.rand( 4,4 ) )

# normalize the columns of V to equal one

for j in range( 0, 4 ):
    V[ :, j ] = V[ :, j ] / np.sqrt( np.transpose( V[:,j] ) * V[:, j ] )

A = V * Lambda * np.linalg.inv( V )

print( 'Lambda = ' )
print( Lambda)

print( 'V = ' )
print( V )

print( 'A = ' )
print( A )


# Pick a random starting vector

x = np.matrix( np.random.rand( 4,1 ) )


mu = 0.    # Let's start by not shifting, so hopefully we hone in on the smallest eigenvalue

for i in range(0,10):
    # We should really compute a factorization of A, but let's be lazy, and compute the inverse
    # explicitly
    Ainv = np.linalg.inv( A - mu * np.eye( 4, 4 ) )
    
    x = Ainv * x 
    
    # normalize x to length one
    x = x / np.sqrt( np.transpose( x ) * x )
    
    # Notice we compute the Rayleigh quotient with matrix A, not Ainv.  This is because
    # the eigenvector of A is an eigenvector of Ainv
    
    mu = np.transpose( x ) * A * x
    
    # The above returns a 1 x 1 matrix.  Let's set mu to the scalar
    
    mu = mu[ 0, 0 ]
    
    print( 'Rayleigh quotient with vector x:', np.transpose( x ) * A * x / ( np.transpose( x ) * x ))
    print( 'inner product of x with v3     :', np.transpose( x ) * V[ :, 3 ] )
    print( ' ' )