function s = SVD3(Y,Z)
  % Y and Z are m-by-3 matrices with full column rank.  Assume m >> 3.
  % s is a column 3-vector where s(1)>=s(2)>=s(3) are the three positive
  % singular values of YZ'
  
  % Note that if QY and QZ are orthogonal, then Y*Z' and (QY'*Y)(QZ*Z)'
  % have the same singular values. 
  
  % If QY'*Y = [ RY;0] and QZ'*Z = [RZ;0] are Qr factorizations with
  % 3-by-3 triangular matrices RY and RZ, then
  
  %                QY'*(Y*Z')QZ = [ RY*RZ'  0 ; 0  0]
  
  % Thus, the three singular values that we want are the three singular
  % values of RY*RZ'.
  
  [QY,RY] = qr(Y,0);
  [QZ,RZ] = qr(Z,0);
  s = svd(RY*RZ');
  
  % Can also do svd instead or QR in the preprocessing step...
  
  [UY,SY,VY] = svd(Y,0);
  [UZ,SZ,VZ] = svd(Z,0);
  s = svd((VY*SY)'*(VZ*SZ));
  
  
 