import numpy  as np      # fundamental package needed for scientific computing.
import scipy.linalg as splin

def PuissItereInv(MatA,MatM,mu,eps):

#   Algo de la puissance iteree inverse

    [n1,n2] = np.shape(MatA)
    
    w1 = np.zeros(n1) + 1.
    w1 = w1/np.sqrt( np.dot(w1,w1) )
    
    [LU,P] = splin.lu_factor(MatA - mu*MatM)
    
    Mw1 = np.dot(MatM,w1)
    
    w2 = splin.lu_solve( (LU,P), Mw1 )
    
    scal1 = np.dot(w2,w2)
    scal2 = np.dot(w1,w2)
    
    w1 = w2 / np.sqrt(scal1)
    
    lambda1= scal2/scal1  + mu
    lambda2 = 2*lambda1
    
    tmp0 = np.dot( np.dot(MatA,w1) - lambda1*np.dot(MatM,w1),w1)
    tmp = tmp0

    k = 0
    while ( np.abs( tmp/tmp0 ) > eps ) :
    
          lambda2 = lambda1
          
          Mw1 = np.dot(MatM,w1)
          
          w2 = splin.lu_solve( (LU,P), Mw1 )
          
          scal1 = np.dot(w2,w2)
          scal2 = np.dot(w1,w2)
          
          w1 = w2 / np.sqrt(scal1)
          
          lambda1= scal2/scal1  + mu
          
          tmp = np.dot( np.dot(MatA,w1) - lambda1*np.dot(MatM,w1),w1)
          
	  k = k + 1
	  
    print "   (Aw - lambda Mw,w)/Init = %8.1e / %8.1e = %8.1e ==> VP min = %12.5e" % ( tmp, tmp0, tmp/tmp0, lambda1 )

    return lambda1