function [Q, T, S, r] = cod(A, tol, rnk)
% [Q, T, S, r] = cod(A, tol, rnk)
%
% Compute the complete orthogonal decomposition (COD) of A:
%   A = Q  [ T 0 ] S^T
%          [ 0 0 ]
% If rnk is specified, it is the rank of A. If rnk is not specified,
% then tol is used to determine the numerical rank of A.
%
% Input:
%   A : m-by-n matrix
%   tol : real << 1. Optional, default value is ||A||_F * machine epsilon.
%   rnk : integer >= 1. Optional, rank of A. No default value.
% Output:
%   Q : m-by-m orthogonal matrix
%   T : r-by-r upper triangular, nonsingular matrix
%   S : n-by-n orthogonal matrix
%   r : integer, rank of A

% Compensate for missing arguments.
if nargin < 2
   tol = norm(A, 'fro') * eps;
end

% Compute QR factorization with column pivoting of A^T = S R1 E^T.
[S, R1, E] = qr(A');

% Find rank of A^T.
if nargin < 3
  r = sum( abs(diag(R1)) > tol);
else
  r = rnk;
end

% Compute QR factorization of R1(1:r, :)^T.
[Q2, R2] = qr( (R1(1:r, :))' );

% Then Q = Q2 E and T = R2(1:r, :).
Q = E * Q2;
T = R2(1:r, :);
