| [nm,m,n,a,w,matu,u,matv,v,ierr,rv1]=svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1); |
function [nm,m,n,a,w,matu,u,matv,v,ierr,rv1]=svd(nm,m,n,a,w,matu,u,matv,v,ierr,rv1);
%***BEGIN PROLOGUE SVD
%***SUBSIDIARY
%***PURPOSE Perform the singular value decomposition of a rectangular
% matrix.
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (SVD-S)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% This subroutine is a translation of the ALGOL procedure SVD,
% NUM. MATH. 14, 403-420(1970) by Golub and Reinsch.
% HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
%
% This subroutine determines the singular value decomposition
% T
% A=USV of a REAL M by N rectangular matrix. Householder
% bidiagonalization and a variant of the QR algorithm are used.
%
% On Input
%
% NM must be set to the row dimension of the two-dimensional
% array parameters, A, U and V, as declared in the calling
% program dimension statement. NM is an INTEGER variable.
% Note that NM must be at least as large as the maximum
% of M and N.
%
% M is the number of rows of A and U.
%
% N is the number of columns of A and U and the order of V.
%
% A contains the rectangular input matrix to be decomposed. A is
% a two-dimensional REAL array, dimensioned A(NM,N).
%
% MATU should be set to true if the U matrix in the
% decomposition is desired, and to false otherwise.
% MATU is a LOGICAL variable.
%
% MATV should be set to true if the V matrix in the
% decomposition is desired, and to false otherwise.
% MATV is a LOGICAL variable.
%
% On Output
%
% A is unaltered (unless overwritten by U or V).
%
% W contains the N (non-negative) singular values of A (the
% diagonal elements of S). They are unordered. If an
% error exit is made, the singular values should be correct
% for indices IERR+1, IERR+2, ..., N. W is a one-dimensional
% REAL array, dimensioned W(N).
%
% U contains the matrix U (orthogonal column vectors) of the
% decomposition if MATU has been set to true Otherwise,
% U is used as a temporary array. U may coincide with A.
% If an error exit is made, the columns of U corresponding
% to indices of correct singular values should be correct.
% U is a two-dimensional REAL array, dimensioned U(NM,N).
%
% V contains the matrix V (orthogonal) of the decomposition if
% MATV has been set to true Otherwise, V is not referenced.
% V may also coincide with A if U does not. If an error
% exit is made, the columns of V corresponding to indices of
% correct singular values should be correct. V is a two-
% dimensional REAL array, dimensioned V(NM,N).
%
% IERR is an INTEGER flag set to
% Zero for normal return,
% K if the K-th singular value has not been
% determined after 30 iterations.
%
% RV1 is a one-dimensional REAL array used for temporary storage,
% dimensioned RV1(N).
%
% CALLS PYTHAG(A,B) for sqrt(A**2 + B**2).
%
% Questions and comments should be directed to B. S. Garbow,
% APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
% ------------------------------------------------------------------
%
%***SEE ALSO EISDOC
%***ROUTINES CALLED PYTHAG
%***REVISION HISTORY (YYMMDD)
% 811101 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900402 Added TYPE section. (WRB)
%***end PROLOGUE SVD
%
persistent c f g h i i1 igo ii its j k k1 kk l l1 ll mn s s1 scalemlv x y z ;
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(l), l=0; end;
if isempty(ii), ii=0; end;
if isempty(i1), i1=0; end;
if isempty(kk), kk=0; end;
if isempty(k1), k1=0; end;
if isempty(ll), ll=0; end;
if isempty(l1), l1=0; end;
if isempty(mn), mn=0; end;
if isempty(its), its=0; end;
if isempty(igo), igo=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nm])).*prod([nm])-numel(a))],nm,[]);
w_shape=size(w);w=reshape(w,1,[]);
u_shape=size(u);u=reshape([u(:).',zeros(1,ceil(numel(u)./prod([nm])).*prod([nm])-numel(u))],nm,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([nm])).*prod([nm])-numel(v))],nm,[]);
rv1_shape=size(rv1);rv1=reshape(rv1,1,[]);
if isempty(c), c=0; end;
if isempty(f), f=0; end;
if isempty(g), g=0; end;
if isempty(h), h=0; end;
if isempty(s), s=0; end;
if isempty(x), x=0; end;
if isempty(y), y=0; end;
if isempty(z), z=0; end;
if isempty(scalemlv), scalemlv=0; end;
if isempty(s1), s1=0; end;
%
%***FIRST EXECUTABLE STATEMENT SVD
ierr = 0;
%
for i = 1 : m;
%
for j = 1 : n;
u(i,j) = a(i,j);
end; j = fix(n+1);
end; i = fix(m+1);
% .......... HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM ..........
g = 0.0e0;
scalemlv = 0.0e0;
s1 = 0.0e0;
%
for i = 1 : n;
l = fix(i + 1);
rv1(i) = scalemlv.*g;
g = 0.0e0;
s = 0.0e0;
scalemlv = 0.0e0;
if( i<=m )
%
for k = i : m;
scalemlv = scalemlv + abs(u(k,i));
end; k = fix(m+1);
%
if( scalemlv~=0.0e0 )
%
for k = i : m;
u(k,i) = u(k,i)./scalemlv;
s = s + u(k,i).^2;
end; k = fix(m+1);
%
f = u(i,i);
g = -(abs(sqrt(s)).*sign(f));
h = f.*g - s;
u(i,i) = f - g;
if( i~=n )
%
for j = l : n;
s = 0.0e0;
%
for k = i : m;
s = s + u(k,i).*u(k,j);
end; k = fix(m+1);
%
f = s./h;
%
for k = i : m;
u(k,j) = u(k,j) + f.*u(k,i);
end; k = fix(m+1);
end; j = fix(n+1);
end;
%
for k = i : m;
u(k,i) = scalemlv.*u(k,i);
end; k = fix(m+1);
end;
end;
%
w(i) = scalemlv.*g;
g = 0.0e0;
s = 0.0e0;
scalemlv = 0.0e0;
if( i<=m && i~=n )
%
for k = l : n;
scalemlv = scalemlv + abs(u(i,k));
end; k = fix(n+1);
%
if( scalemlv~=0.0e0 )
%
for k = l : n;
u(i,k) = u(i,k)./scalemlv;
s = s + u(i,k).^2;
end; k = fix(n+1);
%
f = u(i,l);
g = -(abs(sqrt(s)).*sign(f));
h = f.*g - s;
u(i,l) = f - g;
%
for k = l : n;
rv1(k) = u(i,k)./h;
end; k = fix(n+1);
%
if( i~=m )
%
for j = l : m;
s = 0.0e0;
%
for k = l : n;
s = s + u(j,k).*u(i,k);
end; k = fix(n+1);
%
for k = l : n;
u(j,k) = u(j,k) + s.*rv1(k);
end; k = fix(n+1);
end; j = fix(m+1);
end;
%
for k = l : n;
u(i,k) = scalemlv.*u(i,k);
end; k = fix(n+1);
end;
end;
%
s1 = max(s1,abs(w(i))+abs(rv1(i)));
end; i = fix(n+1);
% .......... ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS ..........
if( matv )
% .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
for ii = 1 : n;
i = fix(n + 1 - ii);
if( i~=n )
if( g~=0.0e0 )
%
for j = l : n;
% .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
v(j,i) =(u(i,j)./u(i,l))./g;
end; j = fix(n+1);
%
for j = l : n;
s = 0.0e0;
%
for k = l : n;
s = s + u(i,k).*v(k,j);
end; k = fix(n+1);
%
for k = l : n;
v(k,j) = v(k,j) + s.*v(k,i);
end; k = fix(n+1);
end; j = fix(n+1);
end;
%
for j = l : n;
v(i,j) = 0.0e0;
v(j,i) = 0.0e0;
end; j = fix(n+1);
end;
%
v(i,i) = 1.0e0;
g = rv1(i);
l = fix(i);
end; ii = fix(n+1);
end;
% .......... ACCUMULATION OF LEFT-HAND TRANSFORMATIONS ..........
if( matu )
% ..........FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- ..........
mn = fix(n);
if( m<n )
mn = fix(m);
end;
%
for ii = 1 : mn;
i = fix(mn + 1 - ii);
l = fix(i + 1);
g = w(i);
if( i~=n )
%
for j = l : n;
u(i,j) = 0.0e0;
end; j = fix(n+1);
end;
%
if( g==0.0e0 )
%
for j = i : m;
u(j,i) = 0.0e0;
end; j = fix(m+1);
else;
if( i~=mn )
%
for j = l : n;
s = 0.0e0;
%
for k = l : m;
s = s + u(k,i).*u(k,j);
end; k = fix(m+1);
% .......... DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW ..........
f =(s./u(i,i))./g;
%
for k = i : m;
u(k,j) = u(k,j) + f.*u(k,i);
end; k = fix(m+1);
end; j = fix(n+1);
end;
%
for j = i : m;
u(j,i) = u(j,i)./g;
%
end; j = fix(m+1);
end;
%
u(i,i) = u(i,i) + 1.0e0;
end; ii = fix(mn+1);
end;
% .......... DIAGONALIZATION OF THE BIDIAGONAL FORM ..........
% .......... FOR K=N STEP -1 UNTIL 1 DO -- ..........
for kk = 1 : n;
k1 = fix(n - kk);
k = fix(k1 + 1);
its = 0;
% .......... TEST FOR SPLITTING.
% FOR L=K STEP -1 UNTIL 1 DO -- ..........
igo=0;
while( true );
for ll = 1 : k;
l1 = fix(k - ll);
l = fix(l1 + 1);
if( s1+abs(rv1(l))==s1 )
igo=1;
break;
end;
% .......... RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
% THROUGH THE BOTTOM OF THE LOOP ..........
if( s1+abs(w(l1))==s1 )
break;
end;
end;
if(igo==0)
% .......... CANCELLATION OF RV1(L) IF L GREATER THAN 1 ..........
c = 0.0e0;
s = 1.0e0;
%
for i = l : k;
f = s.*rv1(i);
rv1(i) = c.*rv1(i);
if( s1+abs(f)==s1 )
break;
end;
g = w(i);
[h ,f,g]=pythag(f,g);
w(i) = h;
c = g./h;
s = -f./h;
if( matu )
%
%
for j = 1 : m;
y = u(j,l1);
z = u(j,i);
u(j,l1) = y.*c + z.*s;
u(j,i) = -y.*s + z.*c;
end; j = fix(m+1);
end;
end;
end;
% .......... TEST FOR CONVERGENCE ..........
z = w(k);
if( l==k )
break;
end;
%
%
if( its==30 )
% .......... SET ERROR -- NO CONVERGENCE TO A
% SINGULAR VALUE AFTER 30 ITERATIONS ..........
ierr = fix(k);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
return;
else;
its = fix(its + 1);
x = w(l);
y = w(k1);
g = rv1(k1);
h = rv1(k);
f = 0.5e0.*(((g+z)./h).*((g-z)./y)+y./h-h./y);
[g ,f]=pythag(f,1.0e0);
f = x -(z./x).*z +(h./x).*(y./(f+(abs(g).*sign(f)))-h);
% .......... NEXT QR TRANSFORMATION ..........
c = 1.0e0;
s = 1.0e0;
%
for i1 = l : k1;
i = fix(i1 + 1);
g = rv1(i);
y = w(i);
h = s.*g;
g = c.*g;
[z ,f,h]=pythag(f,h);
rv1(i1) = z;
c = f./z;
s = h./z;
f = x.*c + g.*s;
g = -x.*s + g.*c;
h = y.*s;
y = y.*c;
if( matv )
%
for j = 1 : n;
x = v(j,i1);
z = v(j,i);
v(j,i1) = x.*c + z.*s;
v(j,i) = -x.*s + z.*c;
end; j = fix(n+1);
end;
%
[z ,f,h]=pythag(f,h);
w(i1) = z;
% .......... ROTATION CAN BE ARBITRARY IF Z IS ZERO ..........
if( z~=0.0e0 )
c = f./z;
s = h./z;
end;
f = c.*g + s.*y;
x = -s.*g + c.*y;
if( matu )
%
%
for j = 1 : m;
y = u(j,i1);
z = u(j,i);
u(j,i1) = y.*c + z.*s;
u(j,i) = -y.*s + z.*c;
end; j = fix(m+1);
end;
end; i1 = fix(k1+1);
%
rv1(l) = 0.0e0;
rv1(k) = f;
w(k) = x;
end;
end;
% .......... CONVERGENCE ..........
if( z<0.0e0 )
% .......... SHIFT FROM BOTTOM 2 BY 2 MINOR ..........
% .......... W(K) IS MADE NON-NEGATIVE ..........
w(k) = -z;
if( matv )
%
for j = 1 : n;
v(j,k) = -v(j,k);
end; j = fix(n+1);
end;
end;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
rv1_shape=zeros(rv1_shape);rv1_shape(:)=rv1(1:numel(rv1_shape));rv1=rv1_shape;
end %subroutine svd
%DECK SVECS
|
|
