function [F,G,D] = gsvd(S,C,L)
  ##GSVD	Generalized Singular Value Decomposition
  ##	[F,G,D] = GCVA(S,C,L) finds matrices which fulfill
  ##		F*S*F' = eye(rank(S))
  ##		G*L*G' = eye(rank(L))
  ##		F*C*G' = D
  ##	where D is a diagonal (possibly non-square) matrix
  ##	with real, positive elements between zero and one
  ##	listed in descending order.
  
  [U1,S1,V1] = svd(S,0);
  [U2,S2,V2] = svd(L,0);
  
  ## X1 = sqrt(pinv(S1))
  if min(size(S1)) == 1
    S1 = S1(1);
  else
    S1 = diag(S1);
  endif
  tol = max(size(S)) * S1(1) * eps;
  S1 = sqrt(S1);
  r = sum(S1 > tol);
  [m,n] = size(S);
  if (r == 0)
    X1 = zeros(n,m);
  else
    S1 = diag(ones(r,1)./S1(1:r));
    X1 = eye(n,r)*S1*eye(r,m)';
  endif
  
  ## X2 = sqrt(pinv(S2))
  if min(size(S2)) == 1
    S2 = S2(1);
  else
    S2 = diag(S2);
  endif
  tol = max(size(S)) * S2(1) * eps;
  S2 = sqrt(S2);
  r = sum(S2 > tol);
  [m,n] = size(L);
  if (r == 0)
    X2 = zeros(n,m);
  else
    S2 = diag(ones(r,1)./S2(1:r));
    X2 = eye(n,r)*S2*eye(r,m)';
  endif
  
  X = X1*U1'*C*U2*X2;
  [U3,S3,V3] = svd(X,0);
  F = U3'*X1*U1';
  G = V3'*X2*U2';
  D = S3;
endfunction