hurvdemo.m

UTV工具包提供46个Matlab函数

function hurvdemo
%  hurvdemo --> Demonstrates high-rank URV algorithms.

%  <Revision>
%    Ricardo D. Fierro, California State University San Marcos
%    Per Christian Hansen, IMM, Technical University of Denmark
%    Peter S. K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: June 22, 1999
%-----------------------------------------------------------------------

fprintf(1,'.                                                           \n');
fprintf(1,'.  The purpose of HURVDEMO is to demonstrate                \n');
fprintf(1,'.  the high-rank URV algorithms in HURV and                 \n');
fprintf(1,'.  HURV_A.                                                  \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

format short e;

% Test matrix generation.
m = 50;
n = 20;
randn('seed',100);
A = randn(m,n);
numrank = 13;

[u,s,v] = svd(A);

s1 = 2*logspace(1,-3,numrank);
s2 = 5*logspace(-4,-6,n-numrank);
s(1:n,:) = diag([s1 s2]);

A = u*s*v';

% Define input parameters.
tol_rank = 0.001;
tol_ref  = 1e-04;
max_ref  = 1;

fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to HURV                                 \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    no. rows m of A                     = %3.0f \n',m);
fprintf(1,'.    no. cols n of A                     = %3.0f \n',n);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'.    refinement tolerance (tol_ref)      = %6.4e \n',tol_ref);
fprintf(1,'.    max refinement steps (max_ref)      = %3.0f \n',max_ref);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Compute the rank-revealing URV decomposition using       \n');
fprintf(1,'.  Stewart''s algorithm:                                    \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.    [p,R,V,U,vec] = hurv(A,tol_rank,tol_ref,max_ref)       \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Compute URV.
flops(0);
t0 = cputime;

[p,R,V,U,vec] = hurv(A,tol_rank,tol_ref,max_ref);

time = cputime - t0;
noflops = flops;

rnorm = norm(R(1:p,p+1:n));
theta = norm(V(:,1:p)'*v(:,p+1:n));
phi   = norm(U(:,1:p)'*u(:,p+1:m));

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from HURV                                 \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',p);
fprintf(1,'.    upper bound of ||R(1:p,p+1:n)||     = %6.4e \n',vec(1));
fprintf(1,'.    estimate of pth singular value      = %6.4e \n',vec(2));
fprintf(1,'.    estimate of (p+1)th singular value  = %6.4e \n',vec(3));
fprintf(1,'.    estimate of num. nullspace angle    = %6.4e \n',vec(4));
fprintf(1,'.    estimate of num. range angle        = %6.4e \n',vec(5));
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  SVD based results                                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',numrank);
fprintf(1,'.    ||R(1:p,p+1:n)||                    = %6.4e \n',rnorm);
fprintf(1,'.    pth singular value                  = %6.4e \n',s1(numrank));
fprintf(1,'.    (p+1)th singular value              = %6.4e \n',s2(1));
fprintf(1,'.    numerical nullspace angle           = %6.4e \n',theta);
fprintf(1,'.    numerical range angle               = %6.4e \n',phi);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Computational complexity of HURV                         \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    floating point operations           = %6.4e \n',noflops);
fprintf(1,'.    CPU time                            = %6.4e \n',time);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

% Define input parameters.
tol_rank = 0.001;
max_iter = 5;
tol_ref  = 1e-04;
max_ref  = 1;

fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to HURV_A                               \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    no. rows m of A                     = %3.0f \n',m);
fprintf(1,'.    no. cols n of A                     = %3.0f \n',n);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'.    no. inverse iterations (max_iter)   = %3.0f \n',max_iter);
fprintf(1,'.    refinement tolerance (tol_ref)      = %6.4e \n',tol_ref);
fprintf(1,'.    max refinement steps (max_ref)      = %3.0f \n',max_ref);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Compute the rank-revealing URV decomposition using the   \n');
fprintf(1,'.  alternative algorithm:                                   \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.    [p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter, ...        \n');
fprintf(1,'.                                         tol_ref,max_ref)  \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Compute URV.
flops(0);
t0 = cputime;

[p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter,tol_ref,max_ref);

time = cputime - t0;
noflops = flops;

rnorm = norm(R(1:p,p+1:n));
theta = norm(V(:,1:p)'*v(:,p+1:n));
phi   = norm(U(:,1:p)'*u(:,p+1:m));

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from HURV_A                               \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',p);
fprintf(1,'.    upper bound of ||R(1:p,p+1:n)||     = %6.4e \n',vec(1));
fprintf(1,'.    estimate of pth singular value      = %6.4e \n',vec(2));
fprintf(1,'.    estimate of (p+1)th singular value  = %6.4e \n',vec(3));
fprintf(1,'.    estimate of num. nullspace angle    = %6.4e \n',vec(4));
fprintf(1,'.    estimate of num. range angle        = %6.4e \n',vec(5));
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  SVD based results                                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',numrank);
fprintf(1,'.    ||R(1:p,p+1:n)||                    = %6.4e \n',rnorm);
fprintf(1,'.    pth singular value                  = %6.4e \n',s1(numrank));
fprintf(1,'.    (p+1)th singular value              = %6.4e \n',s2(1));
fprintf(1,'.    numerical nullspace angle           = %6.4e \n',theta);
fprintf(1,'.    numerical range angle               = %6.4e \n',phi);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Computational complexity of HURV_A                       \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    floating point operations           = %6.4e \n',noflops);
fprintf(1,'.    CPU time                            = %6.4e \n',time);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');

%-----------------------------------------------------------------------
% End of function hurvdemo
%-----------------------------------------------------------------------