[SciPy-user] Matlab to scipy

[Nils Wagner]

Dear experts, 
 
I tried to convert the following Matlab code to scipy. 
(LU-factorization with complete pivoting) 
http://www1.mate.polimi.it/~calnum/programs.html 
 
However it failed. Any help would be appreciated 
 
Nils 
 
function [L,U,P,Q] = LUpivtot(A,n) 
P=eye(n); Q=P; Minv=P; 
for k=1:n-1 
  [Pk,Qk]=pivot(A,k,n); 
  A=Pk*A*Qk; 
  [Mk,Mkinv]=MGauss(A,k,n); 
  A=Mk*A; 
  P=Pk*P; 
  Q=Q*Qk; 
  Minv=Minv*Pk*Mkinv; 
end 
U=triu(A); 
L=P*Minv; 
 
function [Mk,Mkinv]=MGauss(A,k,n) 
Mk=eye(n); 
for i=k+1:n 
   Mk(i,k)=-A(i,k)/A(k,k); 
end 
Mkinv=2*eye(n)-Mk; 
return 
 
function [Pk,Qk]=pivot(A,k,n) 
[y,i]=max(abs(A(k:n,k:n))); 
[piv,jpiv]=max(y); 
ipiv=i(jpiv); 
jpiv=jpiv+k-1; 
ipiv=ipiv+k-1; 
Pk=eye(n); 
Pk(ipiv,ipiv)=0; 
Pk(k,k)=0; 
Pk(k,ipiv)=1; 
Pk(ipiv,k)=1; 
Qk=eye(n); 
Qk(jpiv,jpiv)=0; 
Qk(k,k)=0; 
Qk(k,jpiv)=1; 
Qk(jpiv,k)=1; 
return 
 
------------------------------------------------------------- 
 
from scipy import * 
import MLab 
 
def LUpivot(A,n): 
  P = MLab.eye(n) 
  Q = P 
  Minv = P 
  for k in arange(0,n): 
    Pk,Qk = pivot(A,k,n) 
    A = dot(Pk,dot(A,Qk)) 
    Mk,Mkinv = MGauss(A,k,n) 
    A = dot(Mk,A) 
    P = dot(Pk,P) 
    Q = dot(Q,Qk) 
    Minv = dot(Minv,dot(Pk,Mkinv)) 
  U = MLab.triu(A) 
  L = dot(P,Minv) 
  return L,U,P,Q 
 
def MGauss(A,k,n): 
  Mk = MLab.eye(n) 
  for i in arange(k+1,n): 
     Mk[i,k]=-A[i,k]/A[k,k] 
  Mkinv = 2.0*MLab.eye(n)-Mk 
  return Mk, Mkinv 
 
def pivot(A,k,n): 
  y,i = MLab.max(abs(A[k:n,k:n])) 
  piv,jpiv = MLab.max(y) 
  ipiv=i[jpiv] 
  jpiv=jpiv+k-1 
  ipiv=ipiv+k-1 
  Pk = MLab.eye(n) 
  Pk[ipiv,ipiv]=0 
  Pk[k,k]      =0 
  Pk[k,ipiv]   =1 
  Pk[ipiv,k]   =1 
  Qk=MLab.eye(n) 
  Qk[jpiv,jpiv]=0 
  Qk[k,k]      =0 
  Qk[k,jpiv]   =1 
  Qk[jpiv,k]   =1 
# 
# Example 
# 
A = array(( 
[-4000.0, 2000.0, 2000.0], 
[ 2000.0, 0.78125, 0.0], 
[ 2000.0, 0.0, 0.0])) 
n,m=shape(A) 
L,U,P,Q = LUpivot(A,n)