function [U,S,V]=mp_mysvd(A)

[m,n]=size(A);
a_zero=A(1)*0;%a zero of the same data type as input matrix A
a_one=1+a_zero;%a one of the same data type as input matrix A
g=a_zero; 
scale=a_zero;
anorm=a_zero;
for i=1:n 
 L=i+1;
 rv1(i)=scale*g;
 g=a_zero;
 s=a_zero;
 scale=a_zero;
 if(i<=m)
  for k=i:m 
   scale=scale+abs(A(k,i));
  end 
  if (scale~=0)
   for k=i:m
    A(k,i)=A(k,i)/scale;
    s=s+A(k,i)*A(k,i);
   end 
   f=A(i,i);
   g=-signtransf(sqrt(s),f);
   h=f*g-s;
   A(i,i)=f-g;
   if i~=n
    for j=L:n 
     s=a_zero;
     for k=i:m 
      s=s+A(k,i)*A(k,j);
     end 
     f=s/h;
     for k=i:m 
      A(k,j)=A(k,j)+f*A(k,i);
     end 
    end
   end
   for  k=i:m 
    A(k,i)=scale*A(k,i);
   end 
  end
 end
 w(i)=scale *g;
 g=a_zero;
 s=a_zero;
 scale=a_zero;
 if((i<=m) & (i~=n))
  for k=L:n %do 17 k=l,n
   scale=scale+abs(A(i,k));
  end 
  if(scale~=a_zero)
   for k=L:n 
    A(i,k)=A(i,k)/scale;
    s=s+A(i,k)*A(i,k);
   end 
   f=A(i,L);
   g=-signtransf(sqrt(s),f);
   h=f*g-s;
   A(i,L)=f-g;
   for k=L:n 
    rv1(k)=A(i,k)/h;
   end 
   for j=L:m 
    s=a_zero;
    for k=L:n 
     s=s+A(j,k)*A(i,k);
    end 
    for k=L:n 
     A(j,k)=A(j,k)+s*rv1(k);
    end 
   end 
   for k=L:n 
    A(i,k)=scale*A(i,k);
   end 
  end
 end
 anorm=max(anorm,(abs(w(i))+abs(rv1(i))));
end 
v=repmat(a_zero,n,n);
for i=n:-1:1 
 if(i<n)
  if(g~=a_zero)
   for j=L:n 
    v(j,i)=(A(i,j)/A(i,L))/g;
   end 
   for j=L:n 
    s=a_zero;
    for k=L:n 
     s=s+A(i,k)*v(k,j);
    end 
    for k=L:n 
     v(k,j)=v(k,j)+s*v(k,i);
    end 
   end 
  end
  for j=L:n 
   v(i,j)=a_zero;
   v(j,i)=a_zero;
  end 
 end
 v(i,i)=1+a_zero;
 g=rv1(i);
 L=i;
end 
for i=min(m,n):-1:1 
 L=i+1;
 g=w(i);
 for j=L:n 
  A(i,j)=a_zero;
 end 
 if(g~=a_zero)
  g=a_one/g;
  for j=L:n 
   s=a_zero;
   for k=L:m
    s=s+A(k,i)*A(k,j);
   end 
   f=(s/A(i,i))*g;
   for k=i:m 
    A(k,j)=A(k,j)+f*A(k,i);
   end 
  end 
  for j=i:m 
   A(j,i)=A(j,i)*g;
  end 
 else
  for j=i:m 
   A(j,i)=a_zero;
  end 
 end
 A(i,i)=A(i,i)+a_one;
end 
for k=n:-1:1 
 for its=1:30 
  flag=1;
  for L=k:-1:1 
   nm=L-1; 
   if((abs(rv1(L))+anorm)==anorm)                 
    flag=0; 
    break
   elseif((abs(w(nm))+anorm)==anorm) 
    break
   end
  end 
  if flag==1
   
   c=a_zero; 
   s=a_one;
   for i=L:k 
    f=s*rv1(i);
    rv1(i)=c*rv1(i);
    if((abs(f)+anorm)==anorm) 
     break
    end
    g=w(i);
    h=pythag(f,g);
    w(i)=h;
    h=a_one/h;
    c= (g*h);
    s=-(f*h);
    for j=1:m 
     y=A(j,nm);
     z=A(j,i);
     A(j,nm)=(y*c)+(z*s);
     A(j,i)=-(y*s)+(z*c);
    end 
   end 
  end                   
  z=w(k);
  if(L==k)
   if(z<=a_zero)
    w(k)=-z;
    for j=1:n 
     v(j,k)=-v(j,k);
    end 
   end
   break
  end
  if(its==30) 
   error('no convergence in 30 iterations of svdcmp')
  end
  x=w(L);
  nm=k-1;
  y=w(nm);
  g=rv1(nm);
  h=rv1(k);
  f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
  g=pythag(f,a_one);
  f=((x-z)*(x+z)+h*((y/(f+signtransf(g,f)))-h))/x;
  c=a_one; 
  s=a_one;
  for j=L:nm 
   i=j+1;
   g=rv1(i);
   y=w(i);
   h=s*g;
   g=c*g;
   z=pythag(f,h);
   rv1(j)=z;
   c=f/z;
   s=h/z;
   f= (x*c)+(g*s);
   g=-(x*s)+(g*c);
   h=y*s;
   y=y*c;
   for jj=1:n 
    x=v(jj,j);
    z=v(jj,i);
    v(jj,j)= (x*c)+(z*s);
    v(jj,i)=-(x*s)+(z*c);
   end 
   z=pythag(f,h);
   w(j)=z;
   if(z~=a_zero)
    z=a_one/z;
    c=f*z;
    s=h*z;
   end
   f= (c*g)+(s*y);
   x=-(s*g)+(c*y);
   for jj=1:m 
    y=A(jj,j);
    z=A(jj,i);
    A(jj,j)= (y*c)+(z*s);
    A(jj,i)=-(y*s)+(z*c);
   end 
  end 
  rv1(L)=a_zero;
  rv1(k)=f;
  w(k)=x;
 end 
    
end 
[dummy,I]=sort(-w(1:n));%Sort in decreasing order
S=diag(w(I));
if nargin<2
 if n<size(A,1),S(size(A,1),n)=a_zero;end
 if n<size(A,2),S(n,size(A,2))=a_zero;end
 dummy=A(:,I);
 [U,dummy]=qr(dummy);
 q=find(diag(dummy)<0);U(:,q)=-U(:,q);
else
 U=A(:,I);
end
V=v(1:n,I);
if nargout<=1,U=w(I);U=U(:);end
return

function py=pythag(a,b)
absa=abs(a);
absb=abs(b);
if(absa>absb)
 py=absa*sqrt(1.+(absb/absa)^2);
else
 if(absb==0)
  py=absb;
 else
  py=absb*sqrt(1.+(absa/absb)^2);
 end
end

function c=signtransf(A1,A2)
c=abs(A1);
q=find(A2<0);c(q)=-abs(A1(q));
