# HG changeset patch
# User Martin Koehler <koehlerm@mpi-magdeburg.mpg.de>
# Date 1587723938 -7200
#      Fri Apr 24 12:25:38 2020 +0200
# Node ID 50b2b31abb4b5c6bdbfbe86d9da307f6c2ae6dee
# Parent  a6b7085c8ac258a23d7dce0631b2504869e390bb
Add function ordqz (patch #9897).

* libinterp/corefcn/ordqz.cc: New file.
* libinterp/corefcn/module.mk: Add new file to build system.
* scripts/help/__unimplemented__.m: Remove file from list.
* liboctave/numeric/lo-lapack-proto.h (ztgsen): Add new function prototype.
* NEWS: Announce new function.
* doc/interpreter/linalg.txi: Add docstring to manual.
* libinterp/corefcn/ordschur.cc (Fordschur), scripts/linear-algebra/ordeig.m:
  Add seealso link to new function.

diff -r a6b7085c8ac2 -r 50b2b31abb4b doc/interpreter/linalg.txi
--- a/doc/interpreter/linalg.txi	Thu Oct 22 18:59:34 2020 +0200
+++ b/doc/interpreter/linalg.txi	Fri Apr 24 12:25:38 2020 +0200
@@ -183,6 +183,8 @@
 
 @DOCSTRING(ordschur)
 
+@DOCSTRING(ordqz)
+
 @DOCSTRING(ordeig)
 
 @DOCSTRING(subspace)
diff -r a6b7085c8ac2 -r 50b2b31abb4b libinterp/corefcn/module.mk
--- a/libinterp/corefcn/module.mk	Thu Oct 22 18:59:34 2020 +0200
+++ b/libinterp/corefcn/module.mk	Fri Apr 24 12:25:38 2020 +0200
@@ -221,6 +221,7 @@
   %reldir%/oct-tex-lexer.ll \
   %reldir%/oct-tex-parser.h \
   %reldir%/oct-tex-parser.yy \
+  %reldir%/ordqz.cc \
   %reldir%/ordschur.cc \
   %reldir%/pager.cc \
   %reldir%/pinv.cc \
diff -r a6b7085c8ac2 -r 50b2b31abb4b libinterp/corefcn/ordqz.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libinterp/corefcn/ordqz.cc	Fri Apr 24 12:25:38 2020 +0200
@@ -0,0 +1,661 @@
+////////////////////////////////////////////////////////////////////////
+//
+// Copyright (C) 2020 The Octave Project Developers
+//
+// See the file COPYRIGHT.md in the top-level directory of this
+// distribution or <https://octave.org/copyright/>.
+//
+// This file is part of Octave.
+//
+// Octave is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// Octave is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Octave; see the file COPYING.  If not, see
+// <https://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////////////
+
+// Generalized eigenvalue reordering via LAPACK
+
+// Originally written by M. Koehler <koehlerm(AT)mpi-magdeburg.mpg.de>
+
+#undef DEBUG
+
+#if defined (HAVE_CONFIG_H)
+#  include "config.h"
+#endif
+
+#include <cctype>
+#include <cmath>
+
+#include "f77-fcn.h"
+#include "lo-lapack-proto.h"
+#include "qr.h"
+#include "quit.h"
+
+#include "defun.h"
+#include "error.h"
+#include "errwarn.h"
+#include "ovl.h"
+
+
+#if defined (DEBUG)
+#  include "pager.h"
+#  include "pr-output.h"
+#endif
+
+
+DEFUN (ordqz, args, nargout,
+       doc: /* -*- texinfo -*-
+@deftypefn   {} {[@var{AR}, @var{BR}, @var{QR}, @var{ZR}] =} ordqz (@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{keyword})
+@deftypefnx  {} {[@var{AR}, @var{BR}, @var{QR}, @var{ZR}] =} ordqz (@var{AA}, @var{BB}, @var{Q}, @var{Z}, @var{select})
+Reorder the QZ@tie{}decomposition of a generalized eigenvalue problem.
+
+The generalized eigenvalue problem is defined as
+
+@tex
+$$A x = \lambda B x$$
+@end tex
+@ifnottex
+
+@math{A x = @var{lambda} B x}
+
+@end ifnottex
+
+Its generalized Schur decomposition is computed using the @code{qz} algorithm:
+
+@code{[@var{AA}, @var{BB}, @var{Q}, @var{Z}] = qz ( @var{A}, @var{B})}
+
+where @var{AA}, @var{BB}, @var{QQ}, and @var{ZZ} fulfill
+@tex
+$$ AA = QAZ, BB = QBZ $$
+@end tex
+@ifnottex
+
+@example
+@group
+
+@var{AA} = @var{Q} * @var{A} * @var{Z}, @var{BB} = @var{Q} * @var{B} * @var{Z}
+
+@end group
+@end example
+
+@end ifnottex
+
+The @code{ordqz} function computes a unitary transformation @var{QR} and
+@var{ZR} such that the order of the eigenvalue on the diagonal of @var{AA} and
+@var{BB} is changed. The resulting reordered matrices @var{AR} and @var{BR}
+fulfill:
+
+@tex
+$$ A_R = Q_R \cdot A \cdot Z_R, B_R = Q_R \cdot B \cdot Z_R $$
+@end tex
+@ifnottex
+
+@example
+@group
+
+@var{AR} = @var{QR} * @var{A} * @var{ZR}, @var{BR} = @var{QR} * @var{B} * @var{ZR}
+
+@end group
+@end example
+
+@end ifnottex
+
+The function can either be called with the @code{keyword} argument which
+selects the eigenvalues in the top left block of @var{AR} and @var{BR} in the
+following way:
+
+@table @asis
+@item @qcode{"S"}, @qcode{"udi"}
+small: leading block has all
+@tex
+$|\lambda| < 1$
+@end tex
+@ifnottex
+|@var{lambda}| < 1
+@end ifnottex
+
+@item @qcode{"B"}, @qcode{"udo"}
+big: leading block has all
+@tex
+$|\lambda| \geq 1$
+@end tex
+@ifnottex
+|@var{lambda}| @geq{} 1
+@end ifnottex
+
+@item @qcode{"-"}, @qcode{"lhp"}
+negative real part: leading block has all eigenvalues in the open left
+half-plane
+
+@item @qcode{"+"}, @qcode{"rhp"}
+non-negative real part: leading block has all eigenvalues in the closed right
+half-plane
+@end table
+
+If a logical vector @var{select} is given instead of a keyword the @code{ordqz}
+function reorders all eigenvalues @code{k} to the to left block for which
+@code{select(k)} is true.
+
+Note: The keywords are compatible to the ones of @code{qr} as well as
+to the @sc{matlab} implementation of @code{ordqz}.
+
+@seealso{eig, ordeig, qz, schur, ordschur}
+@end deftypefn */)
+{
+  int nargin = args.length ();
+  enum { LHP, RHP, UDI, UDO, VEC, NONE } select_mode = NONE;
+
+  if (nargin != 5)
+    print_usage ();
+
+  // Check select argument
+  if (args(4).is_string())
+    {
+      std::string opts = args(4).string_value ();
+      std::for_each (opts.begin (), opts.end (),
+                     [] (char & c) { c = std::tolower (c); });
+      if (opts == "lhp" || opts == "-")
+        select_mode = LHP;
+      else if (opts == "rhp" || opts == "+")
+        select_mode = RHP;
+      else if (opts == "udi" || opts == "s")
+        select_mode = UDI;
+      else if (opts == "udo" || opts == "b")
+        select_mode = UDO;
+      else
+        error_with_id ("Octave:ordqz:unknown-keyword",
+                       "ordqz: unknown keyword, possible values: "
+                       "lhp, rhp, udi, udo");
+    }
+  else if (args(4).isreal () || args(4).isinteger () || args(4).islogical ())
+    {
+      if (args(4).rows () > 1 && args(4).columns () > 1)
+        error_with_id ("Octave:ordqz:select-not-vector",
+                       "ordqz: select argument must be a vector.");
+      select_mode = VEC;
+    }
+  else
+    error_with_id ("Octave:ordqz:unknown-arg",
+                   "ordqz: OPT must be string or a logicial vector");
+
+  if (nargout > 4)
+    error_with_id ("Octave:ordqz:nargout",
+                   "ordqz: at most four output arguments required.");
+
+
+
+  // Matrix A: check dimensions.
+  F77_INT nn = octave::to_f77_int (args(0).rows ());
+  F77_INT nc = octave::to_f77_int (args(0).columns ());
+
+  if (args(0).isempty ())
+    {
+      warn_empty_arg ("qz: A");
+      return octave_value_list (2, Matrix ());
+    }
+  else if (nc != nn)
+    err_square_matrix_required ("qz", "A");
+
+  // Matrix A: get value.
+  Matrix aa;
+  ComplexMatrix caa;
+
+  if (args(0).iscomplex ())
+    caa = args(0).complex_matrix_value ();
+  else
+    aa = args(0).matrix_value ();
+
+
+  // Extract argument 2 (bb, or cbb if complex).
+  F77_INT b_nr = octave::to_f77_int (args(1).rows ());
+  F77_INT b_nc = octave::to_f77_int (args(1).columns ());
+
+  if (nn != b_nc || nn != b_nr)
+    err_nonconformant ();
+
+  Matrix bb;
+  ComplexMatrix cbb;
+
+  if (args(1).iscomplex ())
+    cbb = args(1).complex_matrix_value ();
+  else
+    bb = args(1).matrix_value ();
+
+  // Extract argument 3 (qq, or cqq if complex).
+  F77_INT q_nr = octave::to_f77_int (args(2).rows ());
+  F77_INT q_nc = octave::to_f77_int (args(2).columns ());
+
+  if (nn != q_nc || nn != q_nr)
+    err_nonconformant ();
+
+  Matrix qq;
+  ComplexMatrix cqq;
+
+  if (args(2).iscomplex ())
+    cqq = args(2).complex_matrix_value ().hermitian ();
+  else
+    qq = args(2).matrix_value ().transpose ();
+
+  // Extract argument 4 (zz, or czz if complex).
+  F77_INT z_nr = octave::to_f77_int (args(3).rows ());
+  F77_INT z_nc = octave::to_f77_int (args(3).columns ());
+
+  if (nn != z_nc || nn != z_nr)
+    err_nonconformant ();
+
+  Matrix zz;
+  ComplexMatrix czz;
+
+  if (args(3).iscomplex ())
+    czz = args(3).complex_matrix_value ();
+  else
+    zz = args(3).matrix_value ();
+
+  bool complex_case = (args(0).iscomplex () || args(1).iscomplex ())
+                      || args(2).iscomplex () || args(3).iscomplex ();
+
+  if (select_mode == VEC && args(4).rows () != nn &&  args(4).columns () != nn)
+    error_with_id ("Octave:ordqz:numel_select",
+                   "ordqz: selection vector has the wrong number of elements.");
+
+  Array<double> select_array (dim_vector (nn,1));
+  if (select_mode == VEC)
+    select_array = args(4).vector_value ();
+
+  Array<F77_LOGICAL> select (dim_vector (nn, 1));
+
+  if (complex_case)
+    {
+      // Complex
+      if (args(0).isreal ())
+        caa = ComplexMatrix (aa);
+      if (args(1).isreal ())
+        cbb = ComplexMatrix (bb);
+      if (args(2).isreal ())
+        cqq = ComplexMatrix (qq);
+      if (args(3).isreal ())
+        czz = ComplexMatrix (zz);
+
+      ComplexRowVector alpha (dim_vector (nn, 1));
+      ComplexRowVector beta  (dim_vector (nn, 1));
+      octave_idx_type k;
+
+      for (k = 0; k < nn-1; k++)
+        {
+          if (caa(k+1,k) != 0.0)
+            error_with_id ("Octave:ordqz:unsupported_AA",
+                           "ordqz: quasi upper triangular matrices are not "
+                           "allowed with complex data.");
+        }
+
+
+      for (k = 0; k < nn; k++)
+        {
+          alpha(k) = caa(k,k);
+          beta(k)  = cbb(k,k);
+        }
+
+      for (k = 0; k < nn; k++)
+        {
+          switch( select_mode)
+            {
+              case LHP:
+                select(k) = real (alpha(k) * beta(k)) < 0;
+                break;
+              case RHP:
+                select(k) = real (alpha(k) * beta(k)) > 0;
+                break;
+              case UDI:
+                if (beta(k) != 0.0)
+                  select(k) = abs (alpha(k)/beta(k)) < 1.0;
+                else
+                  select(k) = false;
+                break;
+              case UDO:
+                if (beta(k) != 0.0)
+                  select(k) = abs (alpha(k)/beta(k)) > 1.0;
+                else
+                  select(k) = true;
+                break;
+              case VEC:
+                if (select_array(k) != 0.0)
+                  select(k) = true;
+                else
+                  select(k) = false;
+                break;
+              default:
+                // Already checked
+                panic_impossible();
+            }
+        }
+      F77_LOGICAL wantq = 1, wantz = 1;
+      F77_INT ijob = 0, mm, lrwork3 = 1, liwork = 1, info;
+      F77_DBLE pl, pr;
+      ComplexRowVector work3 (lrwork3);
+      Array<F77_INT> iwork (dim_vector (liwork, 1));
+
+      F77_XFCN (ztgsen, ZTGSEN,
+                (ijob, wantq, wantz,
+                 select.fortran_vec (), nn,
+                 F77_DBLE_CMPLX_ARG (caa.fortran_vec ()), nn,
+                 F77_DBLE_CMPLX_ARG (cbb.fortran_vec ()), nn,
+                 F77_DBLE_CMPLX_ARG (alpha.fortran_vec ()),
+                 F77_DBLE_CMPLX_ARG (beta.fortran_vec ()),
+                 F77_DBLE_CMPLX_ARG (cqq.fortran_vec ()), nn,
+                 F77_DBLE_CMPLX_ARG (czz.fortran_vec ()), nn,
+                 mm,
+                 pl, pr,
+                 nullptr,
+                 F77_DBLE_CMPLX_ARG (work3.fortran_vec ()), lrwork3,
+                 iwork.fortran_vec (), liwork,
+                 info));
+      if (info != 0)
+        error_with_id ("Octave:ordqz:ztgsen_failed",
+                       "ordqz: failed to reorder eigenvalues");
+    }
+  else
+    {
+      // Extract eigenvalues
+      RowVector alphar (dim_vector (nn, 1));
+      RowVector alphai (dim_vector (nn, 1));
+      RowVector beta (dim_vector (nn, 1));
+
+      octave_idx_type k;
+
+      k = 0;
+      while (k < nn)
+        {
+#ifdef DEBUG
+          octave_stdout << "ordqz: k = " << k  << " nn = " << nn << " \n";
+#endif
+          if ((k < nn-1 && aa(k+1,k) == 0.0) || k == nn-1)
+            {
+              alphar(k) = aa(k,k);
+              alphai(k) = 0.0;
+              beta(k)   = bb(k,k);
+              k++;
+            }
+          else
+            {
+              double ar[2], ai[2], b[2], work[4];
+              char qz_job = 'E';
+              char comp_q = 'N';
+              char comp_z = 'N';
+              F77_INT nl = 2;
+              F77_INT ilo = 1;
+              F77_INT ihi = 2;
+              F77_INT lwork = 4;
+              F77_INT info = 0;
+              double * aa_vec = aa.fortran_vec ();
+              double * bb_vec = bb.fortran_vec ();
+
+              F77_XFCN (dhgeqz, DHGEQZ,
+                        (F77_CONST_CHAR_ARG2 (&qz_job, 1),
+                         F77_CONST_CHAR_ARG2 (&comp_q, 1),
+                         F77_CONST_CHAR_ARG2 (&comp_z, 1),
+                         nl, ilo, ihi,
+                         &aa_vec[k+k*nn] , nn,
+                         &bb_vec[k+k*nn], nn,
+                         ar, ai, b,
+                         nullptr, nn,
+                         nullptr, nn, work, lwork, info
+                         F77_CHAR_ARG_LEN (1)
+                         F77_CHAR_ARG_LEN (1)
+                         F77_CHAR_ARG_LEN (1)));
+              if (info != 0)
+                error("ordqz: failed to extract eigenvalues");
+
+              alphar(k)   = ar[0];
+              alphar(k+1) = ar[1];
+              alphai(k)   = ai[0];
+              alphai(k+1) = ai[1];
+              beta(k) = b[0];
+              beta(k+1) = b[1];
+
+              k += 2;
+            }
+
+        }
+
+      for (k = 0; k < nn; k++)
+        {
+          switch (select_mode)
+            {
+              case LHP:
+                select(k) = alphar(k) * beta(k) < 0;
+                break;
+              case RHP:
+                select(k) = alphar(k) * beta(k) > 0;
+                break;
+              case UDI:
+                select(k) = alphar(k)*alphar(k)
+                            + alphai(k)*alphai(k) < beta(k)*beta(k);
+                break;
+              case UDO:
+                select(k) = alphar(k)*alphar(k)
+                            + alphai(k)*alphai(k) > beta(k)*beta(k);
+                break;
+              case VEC:
+                if (select_array(k) != 0.0)
+                  select(k) = true;
+                else
+                  select(k) = false;
+                break;
+              default:
+                // Already checked
+                panic_impossible();
+            }
+        }
+
+      F77_LOGICAL wantq = 1, wantz = 1;
+      F77_INT ijob = 0, mm, lrwork3 = 4*nn+16, liwork = nn, info;
+      F77_DBLE pl, pr;
+      RowVector rwork3(lrwork3);
+      Array<F77_INT> iwork (dim_vector (liwork, 1));
+
+      F77_XFCN (dtgsen, DTGSEN,
+                (ijob, wantq, wantz,
+                 select.fortran_vec (), nn,
+                 aa.fortran_vec (), nn,
+                 bb.fortran_vec (), nn,
+                 alphar.fortran_vec (),
+                 alphai.fortran_vec (),
+                 beta.fortran_vec (),
+                 qq.fortran_vec (), nn,
+                 zz.fortran_vec (), nn,
+                 mm,
+                 pl, pr,
+                 nullptr,
+                 rwork3.fortran_vec (), lrwork3,
+                 iwork.fortran_vec (), liwork,
+                 info));
+      if ( info != 0 )
+        error("ordqz: failed to reorder eigenvalues");
+    }
+
+  octave_value_list retval(nargout);
+  switch (nargout)
+    {
+      case 4:
+        if (complex_case)
+          retval(3) = czz;
+        else
+          retval(3) = zz;
+        OCTAVE_FALLTHROUGH;
+      case 3:
+        if (complex_case)
+          retval(2) = cqq.hermitian();
+        else
+          retval(2) = qq.transpose();
+        OCTAVE_FALLTHROUGH;
+      case 2:
+        if (complex_case)
+          retval(1) = cbb;
+        else
+          retval(1) = bb;
+        OCTAVE_FALLTHROUGH;
+      case 1:
+        if (complex_case)
+          retval(0) = caa;
+        else
+          retval(0) = aa;
+        break;
+      case 0:
+        if (complex_case)
+          retval(0) = caa;
+        else
+          retval(0) = aa;
+        break;
+      default:
+        // Already catched
+        error ("ordqz: too many return arguments");
+    }
+  return retval;
+}
+
+
+/*
+%!shared A, B, AA, BB, QQ, ZZ, AC, BC, AAC, BBC, QQC, ZZC, select, selectc
+%! A = [ -1.03428  0.24929  0.43205 -0.12860;
+%!        1.16228  0.27870  2.12954  0.69250;
+%!       -0.51524 -0.34939 -0.77820  2.13721;
+%!       -1.32941  2.11870  0.72005  1.00835 ];
+%! B = [  1.407302 -0.632956 -0.360628  0.068534;
+%!        0.149898  0.298248  0.991777  0.023652;
+%!        0.169281 -0.405205 -1.775834  1.511730;
+%!        0.717770  1.291390 -1.766607 -0.531352 ];
+%! AC = [ 0.4577 + 0.7199i   0.1476 + 0.6946i   0.6202 + 0.2092i   0.7559 + 0.2759i;
+%!        0.5868 + 0.7275i   0.9174 + 0.8781i   0.6741 + 0.1985i   0.4320 + 0.7023i;
+%!        0.2408 + 0.6359i   0.2959 + 0.8501i   0.3904 + 0.5613i   0.5000 + 0.1428i;
+%!        0.8177 + 0.8581i   0.2583 + 0.8970i   0.7706 + 0.5451i   0.1068 + 0.1650i];
+%! BC = [ 0.089898 + 0.209257i   0.157769 + 0.311387i   0.018926 + 0.622517i   0.058825 + 0.374647i;
+%!        0.009367 + 0.098211i   0.736087 + 0.095797i   0.973192 + 0.583765i   0.434018 + 0.461909i;
+%!        0.880784 + 0.868215i   0.032839 + 0.569461i   0.873437 + 0.266081i   0.739426 + 0.362017i;
+%!        0.121649 + 0.115111i   0.426695 + 0.492222i   0.247670 + 0.034414i   0.771629 + 0.078153i];
+%! [AA, BB, QQ, ZZ] = qz (A, B);
+%! [AAC, BBC, QQC, ZZC] = qz (AC, BC);
+%! select = [0 0 1 1];
+%! selectc = [0 0 0 1];
+
+%!error id=Octave:invalid-fun-call ordqz ()
+%!error id=Octave:invalid-fun-call ordqz (eye (2))
+%!error id=Octave:invalid-fun-call ordqz (eye (2), eye (2))
+%!error id=Octave:invalid-fun-call ordqz (eye (2), eye (2), eye (2))
+%!error id=Octave:invalid-fun-call ordqz (eye (2), eye (2), eye (2), eye (2))
+%!error id=Octave:ordqz:unknown-keyword ordqz (eye (2), eye (2), eye (2), eye (2), "test")
+%!error id=Octave:ordqz:select-not-vector ordqz (eye (2), eye (2), eye (2), eye (2), eye (2))
+%!error id=Octave:ordqz:unknown-arg ordqz (eye (2), eye (2), eye (2), eye (2), eye (2) + eye (2)*1i)
+%!error <nonconformant matrices> ordqz (eye (2), eye (3), eye (2), eye (2), "udi")
+%!error <nonconformant matrices> ordqz (eye (2), eye (2), eye (3), eye (2), "udi")
+%!error <nonconformant matrices> ordqz (eye (2), eye (2), eye (2), eye (3), "udi")
+%!error <nonconformant matrices> ordqz (eye (3), eye (2), eye (2), eye (2), "udi")
+%!error id=Octave:ordqz:nargout [a,b,c,d,e] =  ordqz (eye (2), eye (2), eye (2), eye (2), "udi")
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "rhp");
+%! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
+%! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
+%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
+%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "+");
+%! assert (all (real (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 0));
+%! assert (all (real (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 0));
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "lhp");
+%! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
+%! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
+%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
+%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "-");
+%! assert (all (real (eig (AAX(2:4,2:4), BBX(2:4,2:4))) >= 0));
+%! assert (all (real (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 0));
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udi");
+%! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
+%! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
+%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
+%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "S");
+%! assert (all (abs (eig (AAX(1:1,1:1), BBX(1:1,1:1))) < 1));
+%! assert (all (abs (eig (AAX(2:4,2:4), BBX(2:4,2:4))) > 1));
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "udo");
+%! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
+%! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
+%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
+%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, "B");
+%! assert (all (abs (eig (AAX(1:3,1:3), BBX(1:3,1:3))) >= 1));
+%! assert (all (abs (eig (AAX(4:4,4:4), BBX(4:4,4:4))) < 1));
+
+%!test
+%! [AAX, BBX, QQX, ZZX] = ordqz (AA, BB, QQ, ZZ, select);
+%! assert (all (iscomplex (eig (AAX(1:2,1:2), BBX(1:2,1:2)))));
+%! assert (norm (QQX'*AAX*ZZX' - A, "fro"), 0, 1e-12);
+%! assert (norm (QQX'*BBX*ZZX' - B, "fro"), 0, 1e-12);
+
+%!test
+%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "rhp");
+%! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 0));
+%! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 0));
+%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
+%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
+
+%!test
+%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "lhp");
+%! assert (all (real (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 0));
+%! assert (all (real (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 0));
+%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
+%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
+
+%!test
+%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udi");
+%! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) < 1));
+%! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) >= 1));
+%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
+%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
+
+%!test
+%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, "udo");
+%! assert (all (abs (eig (AACX(1:2,1:2), BBCX(1:2,1:2))) >= 1));
+%! assert (all (abs (eig (AACX(3:4,3:4), BBCX(3:4,3:4))) < 1));
+%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
+%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0,  1e-12);
+
+%!test
+%! [AACX, BBCX, QQCX, ZZCX] = ordqz (AAC, BBC, QQC, ZZC, selectc);
+%! ev = abs (eig (AACX(1:1,1:1), BBCX(1:1,1:1)));
+%! assert(ev > 0.6 && ev < 0.7);
+%! assert (norm (QQCX'*AACX*ZZCX' - AC, "fro"), 0, 1e-12);
+%! assert (norm (QQCX'*BBCX*ZZCX' - BC, "fro"), 0, 1e-12);
+
+%!test
+%! A = magic (5);
+%! B = hilb (5);
+%! [AA, BB, Q, Z] = qz (A, B);
+%! [AAS, BBS, QS, ZS] = ordqz (AA, BB, Q, Z, "rhp");
+%! E2 = ordeig (AAS, BBS);
+%! ECOMP = [2.887137073068260e+06; 4.873038494509828e+02; 4.363814079693997e+01; -2.573723338161374e+04; -8.557816814087249e+02];
+%! assert (norm (ECOMP - E2, "Inf"), 0, 1e-8);
+*/
diff -r a6b7085c8ac2 -r 50b2b31abb4b libinterp/corefcn/ordschur.cc
--- a/libinterp/corefcn/ordschur.cc	Thu Oct 22 18:59:34 2020 +0200
+++ b/libinterp/corefcn/ordschur.cc	Fri Apr 24 12:25:38 2020 +0200
@@ -74,7 +74,7 @@
 [@var{U}, @var{S}] = ordschur (@var{U}, @var{S}, [0,1])
 @end example
 
-@seealso{schur, ordeig}
+@seealso{schur, ordeig, ordqz}
 @end deftypefn */)
 {
   if (args.length () != 3)
diff -r a6b7085c8ac2 -r 50b2b31abb4b liboctave/numeric/lo-lapack-proto.h
--- a/liboctave/numeric/lo-lapack-proto.h	Thu Oct 22 18:59:34 2020 +0200
+++ b/liboctave/numeric/lo-lapack-proto.h	Fri Apr 24 12:25:38 2020 +0200
@@ -1640,6 +1640,32 @@
                              const F77_INT& LIWORK,
                              F77_INT& INFO);
 
+  F77_RET_T
+  F77_FUNC (ztgsen, ZTGSEN) (const F77_INT& IJOB,
+                             const F77_LOGICAL& WANTQ,
+                             const F77_LOGICAL& WANTZ,
+                             const F77_LOGICAL *SELECT,
+                             const F77_INT& N,
+                             F77_DBLE_CMPLX *A,
+                             const F77_INT& LDA,
+                             F77_DBLE_CMPLX *B,
+                             const F77_INT& LDB,
+                             F77_DBLE_CMPLX *ALPHA,
+                             F77_DBLE_CMPLX *BETA,
+                             F77_DBLE_CMPLX *Q,
+                             const F77_INT& LDQ,
+                             F77_DBLE_CMPLX *Z,
+                             const F77_INT& LDZ,
+                             F77_INT& M,
+                             F77_DBLE& PL,
+                             F77_DBLE& PR,
+                             F77_DBLE *DIF,
+                             F77_DBLE_CMPLX *WORK,
+                             const F77_INT& LWORK,
+                             F77_INT *IWORK,
+                             const F77_INT& LIWORK,
+                             F77_INT& INFO);
+
   // TRCON
 
   F77_RET_T
diff -r a6b7085c8ac2 -r 50b2b31abb4b scripts/help/__unimplemented__.m
--- a/scripts/help/__unimplemented__.m	Thu Oct 22 18:59:34 2020 +0200
+++ b/scripts/help/__unimplemented__.m	Fri Apr 24 12:25:38 2020 +0200
@@ -1044,7 +1044,6 @@
   "odextend",
   "openFile",
   "opengl",
-  "ordqz",
   "outdegree",
   "outerjoin",
   "pad",
diff -r a6b7085c8ac2 -r 50b2b31abb4b scripts/linear-algebra/ordeig.m
--- a/scripts/linear-algebra/ordeig.m	Thu Oct 22 18:59:34 2020 +0200
+++ b/scripts/linear-algebra/ordeig.m	Fri Apr 24 12:25:38 2020 +0200
@@ -35,7 +35,7 @@
 ## appearance of the matrix @code{@var{A}-@var{lambda}*@var{B}}.  The pair
 ## @var{A}, @var{B} is usually the result of a QZ decomposition.
 ##
-## @seealso{ordschur, eig, schur, qz}
+## @seealso{ordschur, ordqz, eig, schur, qz}
 ## @end deftypefn
 
 function lambda = ordeig (A, B)
