The following illustrates how the routine gls.m fails when the second output is desired:

s = [0.5 1.0 1.5];
B = [0.5 0 0.5; 0 0.5 0.5; 0 0 0.5];
X = randn(1000,3) * [1 0 0; 0 2 0; 0 0 3];
O = kron(diag(s)*diag(s), eye(size(X,1)));
E = randn(size(X)) * diag(s);
Y = X*B + E;
[BETA, SIGMA, R] = gls(Y,X,O);
plot(X,Y,'.');
BETA
SIGMA
sqrt(SIGMA)

resulting in

octave:14> [BETA, SIGMA, R] = gls(Y,X,O);
error: gls: operator /: nonconformant arguments (op1 is 1x1, op2 is 1000x3)
error: called from
    gls at line 120 column 9

However, the routine completes if there is only the first output computed, i.e.,

octave:15> [BETA,~,~] = gls(Y,X,O);

The issue is that the gls.m routine is mistakenly using the variable 'r' for two things, the rank and the residuals.  Here is a diff hunk to fix it, and I will post a changeset after thinking about this a bit because I'm wondering what to do in the case the observation length minus rank is zero ((rx*cy - rnk) = 0 ... probably isn't very likely).

diff --git a/scripts/statistics/base/gls.m b/scripts/statistics/base/gls.m
--- a/scripts/statistics/base/gls.m
+++ b/scripts/statistics/base/gls.m
@@ -104,9 +104,9 @@ function [beta, v, r] = gls (y, x, o)
   z = o * z;
   y1 = o * reshape (y, ry*cy, 1);
   u = z' * z;
-  r = rank (u);
+  rnk = rank (u);

-  if (r == cx*cy)
+  if (rnk == cx*cy)
     b = inv (u) * z' * y1;
   else
     b = pinv (z) * y1;
@@ -117,7 +117,7 @@ function [beta, v, r] = gls (y, x, o)
   if (isargout (2) || isargout (3))
     r = y - x * beta;
     if (isargout (2))
-      v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy - r);
+      v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy - rnk);
     endif
   endif