Can someone please run the following test cases in Matlab?
A = [4 2 6; 3 5 7; 8 7 8; 7 4 1; 8 2 9];
A1 = [4 2 6; 3 5 1; 8 7 9; 7 4 10; 8 2 14];
[U, V, X, C, S] = gsvd(A, diag([1 2 4]))
[U, V, X, C, S] = gsvd(A, diag([1 0 4]))
[U, V, X, C, S] = gsvd(A1, diag([1 2 4]))
[U, V, X, C, S] = gsvd(A1, diag([1 0 4]))
[U, V, X, C, S] = gsvd(A', diag([1 2 4 8 16]))
[U, V, X, C, S] = gsvd(A', diag([1 0 4 8 16]))
[U, V, X, C, S] = gsvd(A1', diag([1 2 4 8 16]))
[U, V, X, C, S] = gsvd(A1', diag([1 0 4 8 16]))
Z = A + i*[8 8 1; 9 9 4; 1 7 5; 8 7 9; 3 1 5];
Z1 = A1 + i*[8 8 8; 9 9 9; 1 7 5; 8 7 9; 3 1 5];
[U, V, X, C, S] = gsvd(Z, diag([1 2 4])+i*diag([4 2 1]))
[U, V, X, C, S] = gsvd(Z, diag([1 0 4])+i*diag([4 0 1]))
[U, V, X, C, S] = gsvd(Z1, diag([1 2 4])+i*diag([4 2 1]))
[U, V, X, C, S] = gsvd(Z1, diag([1 0 4])+i*diag([4 0 1]))
[U, V, X, C, S] = gsvd(Z', diag([1 2 4 8 16])+i*diag([5 4 3 2 1]))
[U, V, X, C, S] = gsvd(Z', diag([1 0 4 8 16])+i*diag([5 0 3 2 1]))
[U, V, X, C, S] = gsvd(Z1', diag([1 2 4 8 16])+i*diag([5 4 3 2 1]))
[U, V, X, C, S] = gsvd(Z1', diag([1 0 4 8 16])+i*diag([5 0 3 2 1]))

It looks like there are differences between Matlab and LAPACK in how GSVD is defined. This is alluded to vaguely in both Matlab's [1] and Octave's [2] documentation.
An interface to LAPACK's GSVD routines is good to have available, regardless of what Matlab is doing.
As far as recovering Matlab's form, we can do that (up to the ordering of singular values and vectors, which seems to be consistently the opposite) for at least their Examples 1 and 2 with
[U, V, C, S, X, R] = gsvd (A, B);
X = (R / X)';
[m, n] = size(A);
if m > n
C = [C; zeros(mn, n)];
elseif m < n
C = [C zeros(m, nm)];
S0 = S; S = eye(n); S(1:m, 1:m) = S0;
endif
[1] cf. https://groups.google.com/forum/#!topic/comp.softsys.matlab/aBgSZNTBcxU
[2] "Other authors, like S. Van Huffel, define this transformation as the simulatenous diagonalisation of the input matrices, this can be achieved by multiplying x by the inverse of [I 0; 0 r]."
