function [S U V it t t1 F]=svd_GolRein_1(A)
% Se calcula al descomposición en valores singulares, [U S V]=U*S*V=A,
%    de una matriz cualquiera m x n.
% Sigue casi exactamente el "paper" de Golub y Reinsch, "Singular Value 
%  Decomposition and Least Squares Solutions", Numeri. Math., 14, 1970.
%
 [m,n]=size(A); tol=norm(A,'fro')*sqrt(eps); itr=0; tic;   
 if n>m, A=A'; [m,~]=size(A); itr=1; end      % Es necesario que m >= n

 % Bidiagonaliza A por Househoulder: U*B_n*V' = A
 [U,B,n,V] = bidiag_Hansen(A); t1=toc; tic;
 q=B_n(:,1); e=[0; B_n(1:n-1,2)]; U1=eye(m); V1=eye(n);
 k=n; it=0;
 while k>0                         % Bucle principal proceso iterativo
    caso=1;
    for l=k:-1:1
       if abs(e(l))  <=tol, e(l)=0;   caso=0; break, end
       if abs(q(l-1))<=tol, q(l-1)=0;         break, end
    end
    if caso
       c=0; s=1; l1=l-1;  % Cancelación columna l-1 con Givens
       for i=l:k
          f=s*e(i); e(i)=c*e(i);
          if abs(f)<=tol, break, end
          g=q(i); h=sqrt(f*f+g*g); c=g/h; s=-f/h; q(i)=h;
          U1(:,[l1 i])=U1(:,[l1 i])*[c -s;s c];
       end
    end
    z=q(k);                        % Comprobar convergencia                       
    if l==k
       if z<0, q(k)=-z; V1(:,k)=-V1(:,k); end
       k=k-1; continue
    end
    x=q(l); y=q(k-1); g=e(k-1); h=e(k);  % Clculo desplazami. con 2x2
    f=((y-z)*(y+z)+(g-h)*(g+h))/(2*h*y); g=sqrt(f*f+1);
    f=((x-z)*(x+z)+h*(y/(f+(sign(f)+(f==0))*g)-h))/x;
    c=1; s=1;                            %  Comienzo iteración QR
    for i=l+1:k
       g=e(i); y=q(i); h=s*g; g=c*g;
       e(i-1)=sqrt(f*f+h*h); z=e(i-1); c=f/z; s=h/z;
       f=x*c+g*s; g=-x*s+g*c; h=y*s; y=y*c;
       V1(:,[i-1 i])=V1(:,[i-1 i])*[c -s;s c];
       q(i-1)=sqrt(f*f+h*h); z=q(i-1); c=f/z; s=h/z;
       f=c*g+s*y; x=-s*g+c*y;
       U1(:,[i-1 i])=U1(:,[i-1 i])*[c -s;s c];
    end
    e(l)=0; e(k)=f; q(k)=x;
    it=it+1;
 end
 U=U*U1; V=V1'*V';
 S=zeros(m,n);        % Ordenar valores singulares en S y salida resul.
 for i=1:n                        
    ssn=q(i); S(i,i)=abs(ssn); if ssn<0, V(i,:)=-V(i,:); end, 
 end
 [q,ix]=sort(diag(S),'descend');         
 for i=1:n, S(i,i)=q(i); end              % Valores sing. en S de + a -
 ix1=1:m; P1=eye(m); ix1(1:n)=ix; P1=P1(:,ix1); U=U*P1; V=V(ix,:);
 F=max(max(U*S*V-A)); t=toc;              % Chequeo precisión resultados
 if itr==1, S=S'; Te=U; U=V'; V=Te'; end  % Si m<n, trasponer resultados
 if nargout==1, S=diag(S); end
end