function [U,sm,X,V] = cgsvd(A,L)%CGSVD Compact generalized SVD of a matrix pair in regularization problems.%% sm = cgsvd(A,L)% [U,sm,X,V] = cgsvd(A,L) ,  sm = [sigma,mu]%% Computes the generalized SVD of the matrix pair (A,L). The dimensions% of A and L must be such that [A;L] has more rows than columns.%% If m >= n >= p then the GSVD has the form:%    [ A ] = [ U  0 ]*[ diag(sigma)      0    ]*inv(X)%    [ L ]   [ 0  V ] [      0       eye(n-p) ]%                     [  diag(mu)        0    ]% where%    U  is  m-by-n ,    sigma  is  p-by-1%    V  is  p-by-p ,    mu     is  p-by-1%    X  is  n-by-n .%% Otherwise the GSVD has the form:%    [ A ] = [ U  0 ]*[    0     diag(sigma)    0     ]*inv([X0,X])%    [ L ]   [ 0  V ] [    0          0      eye(n-p) ]%                     [ eye(n-m)      0         0     ]%                     [    0      diag(mu)      0     ] % Reference: C. F. Van Loan, "Computing the CS and the generalized % singular value decomposition", Numer. Math. 46 (1985), 479-491.  % Per Christian Hansen, IMM, April 14, 2003.  % Initialization.[m,n] = size(A); [p,n1] = size(L);if (n1 ~= n)  error('No. columns in A and L must be the same')endif (m+p < n)  error('Dimensions must satisfy m+p >= n')end% Call Matlab's GSVD routine.[U,V,W,C,S] = gsvd(full(A),full(L),0);if (m >= n)  % The overdetermined or square case.  sm = [diag(C(1:p,1:p)),diag(S(1:p,1:p))];   if (nargout < 2)     U = sm;   else     % Full decomposition.     X = inv(W');   endelse  % The underdetermined case.  sm = [diag(C(1:m+p-n,n-m+1:p)),diag(S(n-m+1:p,n-m+1:p))];   if (nargout < 2)     U = sm;   else     % Full decomposition.     X = inv(W');    X = X(:,n-m+1:n);   endend