comparison toolbox/cod.m @ 2:c124219d7bfa draft

Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author Antonio Pino Robles <data.script93@gmail.com>
date Thu, 07 May 2015 18:36:24 +0200
parents 8f23314345f4
children
comparison
equal deleted inserted replaced
1:e471a92d17be 2:c124219d7bfa
1 function [U, R, V] = cod(A, tol)
2 %COD Complete orthogonal decomposition.
3 % [U, R, V] = COD(A, TOL) computes a decomposition A = U*T*V,
4 % where U and V are unitary, T = [R 0; 0 0] has the same dimensions as
5 % A, and R is upper triangular and nonsingular of dimension rank(A).
6 % Rank decisions are made using TOL, which defaults to approximately
7 % MAX(SIZE(A))*NORM(A)*EPS.
8 % By itself, COD(A, TOL) returns R.
9
10 % Reference:
11 % G.H. Golub and C.F. Van Loan, Matrix Computations, Second
12 % Edition, Johns Hopkins University Press, Baltimore, Maryland,
13 % 1989, sec. 5.4.2.
14
15 [m, n] = size(A);
16
17 % QR decomposition.
18 [U, R, P] = qr(A); % AP = UR
19 V = P'; % A = URV;
20 if nargin == 1, tol = max(m,n)*eps*abs(R(1,1)); end % |R(1,1)| approx NORM(A).
21
22 % Determine r = effective rank.
23 r = sum(abs(diag(R)) > tol);
24 r = r(1); % Fix for case where R is vector.
25 R = R(1:r,:); % Throw away negligible rows (incl. all zero rows, m>n).
26
27 if r ~= n
28
29 % Reduce nxr R' = r [L] to lower triangular form: QR' = [Lbar].
30 % n-r [M] [0]
31
32 [Q, R] = trap2tri(R');
33 V = Q*V;
34 R = R';
35
36 end
37
38 if nargout <= 1, U = R; end