Hi all, Just another sqrtm implementation based on Pade iteration... def sqrtm2(A,p): # # Pade Iteration # m,n = shape(A) i=arange(1,p+1) xi = 0.5*(1.+cos(0.5*(2*i-1)*pi/p)) ai2 = 1.0/xi-1. Y0 = A Z0 = identity(n) maxiter = 5 for k in arange(0,maxiter): sum1= zeros((n,n),Complex) sum2= zeros((n,n),Complex) for j in arange(0,p): sum1 = sum1 + linalg.inv(dot(Z0,Y0)+ai2[j]*identity(n))/xi[j] sum2 = sum2 + linalg.inv(dot(Y0,Z0)+ai2[j]*identity(n))/xi[j] Y1 = dot(Y0,sum1)/p Z1 = dot(Z0,sum1)/p Y0 = Y1 Z0 = Z1 return Y1,Z1 # # test matrix # A=rand(5,5) # # Modify the spectrum # A = dot(A,A)/25.0+identity(5) # Y2,Z2 = sqrtm2(A,4) print linalg.norm(dot(Y2,Y2)-A) print linalg.norm(dot(Z2,Z2)-linalg.inv(A)) Any suggestions ? Nils