Confirmed.  I checked the LAPACK routine for DHGEQZ and the Fortran routine requires that aa, bb inputs are in/out.  The documentation definitely refers to returning some results in these bidirectional function parameters.

One solution would be to use a local variable for aa and bb.  However, I note that the call to DHGEQZ is within a loop

while (k < nn)

which runs from 0 to the number of rows in the matrix in steps of 2.  Currently, the local version of 'aa' is initialized just once.  But if we change to initializing the local version every time there will be nr/2 iterations and nr^2 elements to copy per initialization.  This doesn't seem ideal.

Is Octave using the right algorithm?

I notice that for complex aa,bb the code simply calls ZTSGEN.  But for real cases, the code runs this loop which repeatedly calls DHGEQZ and then after all that work is done it calls DTSGEN.  Is there a way to just call DTSGEN once and skip the rest?

there was some changes for octave 9


warning: qz: returns the complex QZ by default on real matrices since version 9
warning: called from
    test33371 at line 26 column 16

ans =

   0.0000 - 2.9278i
  -0.0000 + 2.9278i
   0.5916 - 0.0000i
  -0.5916 +      0i

>> ver

----------------------------------------------------------------------
GNU Octave Version: 9.1.1 (hg id: c59ff2782e04)
GNU Octave License: GNU General Public License
Operating System: Linux 6.8.0-39-generic #39-Ubuntu SMP PREEMPT_DYNAMIC Fri Jul  521:49:14 UTC 2024 x86_64

The function ordqz is used to move eigenvalue blocks in the triangular A and B matrices obtained from qz.

If called repeatedly using the select vector, which is sometimes needed, the result may become unreliable. The code snippet below fails.

tol = 1e-13;

a = [ 0   0   1  -2
      0  -2   0  -1
      1   0  -4   1
     -2  -1   1  -2];

b = [ 0  -1   2   2
      1   0   0  -1
     -2   0   0  -1
     -2   1   1   0];

[aa, bb, q, z] = qz (a, b);
ordeig (aa, bb)

[aa, bb, q, z] = ordqz (aa, bb, q, z, [0 0 1 1]);
assert (norm (aa - q * a * z) < tol);  % succeeds

[aa, bb, q, z] = ordqz (aa, bb, q, z, [0 0 1 1]);
assert (norm (aa - q * a * z) < tol);  % fails

During the check and extraction of eigenvalues in ordqz.cc, 2x2 blocks are passed as arguments to DHGEQZ. This might modify the corresponding 2x2 blocks in aa and bb, but the remaing rows and columns are not affected, neither qq nor zz. This seems to cause the problem.                                           
                                                                               
FIX: copy the 2x2 blocks in aa and bb to local 2x2 matrices before           
calling DHGEQZ, instead of operating on aa and bb directly.

I am using version 7.3.0 but the problem seems to remain in later versions as well.