--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/cod.m	Wed May 06 14:56:53 2015 +0200
@@ -0,0 +1,38 @@
+function [U, R, V] = cod(A, tol)
+%COD    Complete orthogonal decomposition.
+%       [U, R, V] = COD(A, TOL) computes a decomposition A = U*T*V,
+%       where U and V are unitary, T = [R 0; 0 0] has the same dimensions as
+%       A, and R is upper triangular and nonsingular of dimension rank(A).
+%       Rank decisions are made using TOL, which defaults to approximately
+%       MAX(SIZE(A))*NORM(A)*EPS.
+%       By itself, COD(A, TOL) returns R.
+
+%       Reference:
+%       G.H. Golub and C.F. Van Loan, Matrix Computations, Second
+%       Edition, Johns Hopkins University Press, Baltimore, Maryland,
+%       1989, sec. 5.4.2.
+
+[m, n] = size(A);
+
+% QR decomposition.
+[U, R, P] = qr(A);    % AP = UR
+V = P';               % A = URV;
+if nargin == 1, tol = max(m,n)*eps*abs(R(1,1)); end  % |R(1,1)| approx NORM(A).
+
+% Determine r = effective rank.
+r = sum(abs(diag(R)) > tol);
+r = r(1);             % Fix for case where R is vector.
+R = R(1:r,:);         % Throw away negligible rows (incl. all zero rows, m>n).
+
+if r ~= n
+
+   % Reduce nxr R' =  r  [L]  to lower triangular form: QR' = [Lbar].
+   %                 n-r [M]                                  [0]
+
+   [Q, R] = trap2tri(R');
+   V = Q*V;
+   R = R';
+
+end
+
+if nargout <= 1, U = R; end