--- a/liboctave/numeric/schur.cc	Wed Aug 17 10:37:57 2016 -0400
+++ b/liboctave/numeric/schur.cc	Wed Aug 17 10:55:38 2016 -0400
@@ -34,437 +34,435 @@
 
 namespace octave
 {
-namespace math
-{
-
-// For real types.
+  namespace math
+  {
+    // For real types.
 
-template <typename T>
-static octave_idx_type
-select_ana (const T& a, const T&)
-{
-  return (a < 0.0);
-}
-
-template <typename T>
-static octave_idx_type
-select_dig (const T& a, const T& b)
-{
-  return (hypot (a, b) < 1.0);
-}
-
-// For complex types.
+    template <typename T>
+    static octave_idx_type
+    select_ana (const T& a, const T&)
+    {
+      return (a < 0.0);
+    }
 
-template <typename T>
-static octave_idx_type
-select_ana (const T& a)
-{
-  return a.real () < 0.0;
-}
+    template <typename T>
+    static octave_idx_type
+    select_dig (const T& a, const T& b)
+    {
+      return (hypot (a, b) < 1.0);
+    }
 
-template <typename T>
-static octave_idx_type
-select_dig (const T& a)
-{
-  return (abs (a) < 1.0);
-}
+    // For complex types.
 
-template <>
-octave_idx_type
-schur<Matrix>::init (const Matrix& a, const std::string& ord, bool calc_unitary)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
+    template <typename T>
+    static octave_idx_type
+    select_ana (const T& a)
+    {
+      return a.real () < 0.0;
+    }
 
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("schur: requires square matrix");
-
-  if (a_nr == 0)
+    template <typename T>
+    static octave_idx_type
+    select_dig (const T& a)
     {
-      schur_mat.clear ();
-      unitary_mat.clear ();
-      return 0;
+      return (abs (a) < 1.0);
     }
 
-  // Workspace requirements may need to be fixed if any of the
-  // following change.
+    template <>
+    octave_idx_type
+    schur<Matrix>::init (const Matrix& a, const std::string& ord, bool calc_unitary)
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
+
+      if (a_nr != a_nc)
+        (*current_liboctave_error_handler) ("schur: requires square matrix");
 
-  char jobvs;
-  char sense = 'N';
-  char sort = 'N';
+      if (a_nr == 0)
+        {
+          schur_mat.clear ();
+          unitary_mat.clear ();
+          return 0;
+        }
+
+      // Workspace requirements may need to be fixed if any of the
+      // following change.
 
-  if (calc_unitary)
-    jobvs = 'V';
-  else
-    jobvs = 'N';
+      char jobvs;
+      char sense = 'N';
+      char sort = 'N';
+
+      if (calc_unitary)
+        jobvs = 'V';
+      else
+        jobvs = 'N';
 
-  char ord_char = ord.empty () ? 'U' : ord[0];
+      char ord_char = ord.empty () ? 'U' : ord[0];
+
+      if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
+        sort = 'S';
 
-  if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = 'S';
+      volatile double_selector selector = 0;
+      if (ord_char == 'A' || ord_char == 'a')
+        selector = select_ana<double>;
+      else if (ord_char == 'D' || ord_char == 'd')
+        selector = select_dig<double>;
 
-  volatile double_selector selector = 0;
-  if (ord_char == 'A' || ord_char == 'a')
-    selector = select_ana<double>;
-  else if (ord_char == 'D' || ord_char == 'd')
-    selector = select_dig<double>;
+      octave_idx_type n = a_nc;
+      octave_idx_type lwork = 8 * n;
+      octave_idx_type liwork = 1;
+      octave_idx_type info;
+      octave_idx_type sdim;
+      double rconde;
+      double rcondv;
+
+      schur_mat = a;
+
+      if (calc_unitary)
+        unitary_mat.clear (n, n);
+
+      double *s = schur_mat.fortran_vec ();
+      double *q = unitary_mat.fortran_vec ();
+
+      Array<double> wr (dim_vector (n, 1));
+      double *pwr = wr.fortran_vec ();
+
+      Array<double> wi (dim_vector (n, 1));
+      double *pwi = wi.fortran_vec ();
 
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 8 * n;
-  octave_idx_type liwork = 1;
-  octave_idx_type info;
-  octave_idx_type sdim;
-  double rconde;
-  double rcondv;
+      Array<double> work (dim_vector (lwork, 1));
+      double *pwork = work.fortran_vec ();
+
+      // BWORK is not referenced for the non-ordered Schur routine.
+      octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
+      Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
+      octave_idx_type *pbwork = bwork.fortran_vec ();
+
+      Array<octave_idx_type> iwork (dim_vector (liwork, 1));
+      octave_idx_type *piwork = iwork.fortran_vec ();
 
-  schur_mat = a;
-
-  if (calc_unitary)
-    unitary_mat.clear (n, n);
+      F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
+                                 F77_CONST_CHAR_ARG2 (&sort, 1),
+                                 selector,
+                                 F77_CONST_CHAR_ARG2 (&sense, 1),
+                                 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
+                                 pwork, lwork, piwork, liwork, pbwork, info
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)));
 
-  double *s = schur_mat.fortran_vec ();
-  double *q = unitary_mat.fortran_vec ();
-
-  Array<double> wr (dim_vector (n, 1));
-  double *pwr = wr.fortran_vec ();
+      return info;
+    }
 
-  Array<double> wi (dim_vector (n, 1));
-  double *pwi = wi.fortran_vec ();
+    template <>
+    octave_idx_type
+    schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord,
+                              bool calc_unitary)
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
 
-  Array<double> work (dim_vector (lwork, 1));
-  double *pwork = work.fortran_vec ();
+      if (a_nr != a_nc)
+        (*current_liboctave_error_handler) ("SCHUR requires square matrix");
+
+      if (a_nr == 0)
+        {
+          schur_mat.clear ();
+          unitary_mat.clear ();
+          return 0;
+        }
 
-  // BWORK is not referenced for the non-ordered Schur routine.
-  octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
-  Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
-  octave_idx_type *pbwork = bwork.fortran_vec ();
+      // Workspace requirements may need to be fixed if any of the
+      // following change.
+
+      char jobvs;
+      char sense = 'N';
+      char sort = 'N';
 
-  Array<octave_idx_type> iwork (dim_vector (liwork, 1));
-  octave_idx_type *piwork = iwork.fortran_vec ();
+      if (calc_unitary)
+        jobvs = 'V';
+      else
+        jobvs = 'N';
+
+      char ord_char = ord.empty () ? 'U' : ord[0];
+
+      if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
+        sort = 'S';
+
+      volatile float_selector selector = 0;
+      if (ord_char == 'A' || ord_char == 'a')
+        selector = select_ana<float>;
+      else if (ord_char == 'D' || ord_char == 'd')
+        selector = select_dig<float>;
 
-  F77_XFCN (dgeesx, DGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
-                             F77_CONST_CHAR_ARG2 (&sort, 1),
-                             selector,
-                             F77_CONST_CHAR_ARG2 (&sense, 1),
-                             n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
-                             pwork, lwork, piwork, liwork, pbwork, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
+      octave_idx_type n = a_nc;
+      octave_idx_type lwork = 8 * n;
+      octave_idx_type liwork = 1;
+      octave_idx_type info;
+      octave_idx_type sdim;
+      float rconde;
+      float rcondv;
+
+      schur_mat = a;
 
-  return info;
-}
+      if (calc_unitary)
+        unitary_mat.clear (n, n);
+
+      float *s = schur_mat.fortran_vec ();
+      float *q = unitary_mat.fortran_vec ();
+
+      Array<float> wr (dim_vector (n, 1));
+      float *pwr = wr.fortran_vec ();
+
+      Array<float> wi (dim_vector (n, 1));
+      float *pwi = wi.fortran_vec ();
 
-template <>
-octave_idx_type
-schur<FloatMatrix>::init (const FloatMatrix& a, const std::string& ord,
-                          bool calc_unitary)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
+      Array<float> work (dim_vector (lwork, 1));
+      float *pwork = work.fortran_vec ();
+
+      // BWORK is not referenced for the non-ordered Schur routine.
+      octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
+      Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
+      octave_idx_type *pbwork = bwork.fortran_vec ();
+
+      Array<octave_idx_type> iwork (dim_vector (liwork, 1));
+      octave_idx_type *piwork = iwork.fortran_vec ();
 
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("SCHUR requires square matrix");
+      F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
+                                 F77_CONST_CHAR_ARG2 (&sort, 1),
+                                 selector,
+                                 F77_CONST_CHAR_ARG2 (&sense, 1),
+                                 n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
+                                 pwork, lwork, piwork, liwork, pbwork, info
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)));
 
-  if (a_nr == 0)
-    {
-      schur_mat.clear ();
-      unitary_mat.clear ();
-      return 0;
+      return info;
     }
 
-  // Workspace requirements may need to be fixed if any of the
-  // following change.
+    template <>
+    octave_idx_type
+    schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord,
+                                bool calc_unitary)
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
 
-  char jobvs;
-  char sense = 'N';
-  char sort = 'N';
+      if (a_nr != a_nc)
+        (*current_liboctave_error_handler) ("SCHUR requires square matrix");
 
-  if (calc_unitary)
-    jobvs = 'V';
-  else
-    jobvs = 'N';
+      if (a_nr == 0)
+        {
+          schur_mat.clear ();
+          unitary_mat.clear ();
+          return 0;
+        }
 
-  char ord_char = ord.empty () ? 'U' : ord[0];
+      // Workspace requirements may need to be fixed if any of the
+      // following change.
 
-  if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = 'S';
+      char jobvs;
+      char sense = 'N';
+      char sort = 'N';
 
-  volatile float_selector selector = 0;
-  if (ord_char == 'A' || ord_char == 'a')
-    selector = select_ana<float>;
-  else if (ord_char == 'D' || ord_char == 'd')
-    selector = select_dig<float>;
+      if (calc_unitary)
+        jobvs = 'V';
+      else
+        jobvs = 'N';
+
+      char ord_char = ord.empty () ? 'U' : ord[0];
+
+      if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
+        sort = 'S';
 
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 8 * n;
-  octave_idx_type liwork = 1;
-  octave_idx_type info;
-  octave_idx_type sdim;
-  float rconde;
-  float rcondv;
+      volatile complex_selector selector = 0;
+      if (ord_char == 'A' || ord_char == 'a')
+        selector = select_ana<Complex>;
+      else if (ord_char == 'D' || ord_char == 'd')
+        selector = select_dig<Complex>;
 
-  schur_mat = a;
-
-  if (calc_unitary)
-    unitary_mat.clear (n, n);
+      octave_idx_type n = a_nc;
+      octave_idx_type lwork = 8 * n;
+      octave_idx_type info;
+      octave_idx_type sdim;
+      double rconde;
+      double rcondv;
 
-  float *s = schur_mat.fortran_vec ();
-  float *q = unitary_mat.fortran_vec ();
-
-  Array<float> wr (dim_vector (n, 1));
-  float *pwr = wr.fortran_vec ();
+      schur_mat = a;
+      if (calc_unitary)
+        unitary_mat.clear (n, n);
 
-  Array<float> wi (dim_vector (n, 1));
-  float *pwi = wi.fortran_vec ();
+      Complex *s = schur_mat.fortran_vec ();
+      Complex *q = unitary_mat.fortran_vec ();
+
+      Array<double> rwork (dim_vector (n, 1));
+      double *prwork = rwork.fortran_vec ();
 
-  Array<float> work (dim_vector (lwork, 1));
-  float *pwork = work.fortran_vec ();
+      Array<Complex> w (dim_vector (n, 1));
+      Complex *pw = w.fortran_vec ();
+
+      Array<Complex> work (dim_vector (lwork, 1));
+      Complex *pwork = work.fortran_vec ();
 
-  // BWORK is not referenced for the non-ordered Schur routine.
-  octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
-  Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
-  octave_idx_type *pbwork = bwork.fortran_vec ();
+      // BWORK is not referenced for non-ordered Schur.
+      octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
+      Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
+      octave_idx_type *pbwork = bwork.fortran_vec ();
 
-  Array<octave_idx_type> iwork (dim_vector (liwork, 1));
-  octave_idx_type *piwork = iwork.fortran_vec ();
+      F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
+                                 F77_CONST_CHAR_ARG2 (&sort, 1),
+                                 selector,
+                                 F77_CONST_CHAR_ARG2 (&sense, 1),
+                                 n, F77_DBLE_CMPLX_ARG (s), n, sdim, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (q), n, rconde, rcondv,
+                                 F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)));
 
-  F77_XFCN (sgeesx, SGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
-                             F77_CONST_CHAR_ARG2 (&sort, 1),
-                             selector,
-                             F77_CONST_CHAR_ARG2 (&sense, 1),
-                             n, s, n, sdim, pwr, pwi, q, n, rconde, rcondv,
-                             pwork, lwork, piwork, liwork, pbwork, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
+      return info;
+    }
 
-  return info;
-}
+    template <>
+    schur<ComplexMatrix>
+    rsf2csf<ComplexMatrix, Matrix> (const Matrix& s_arg, const Matrix& u_arg)
+    {
+      ComplexMatrix s (s_arg);
+      ComplexMatrix u (u_arg);
+
+      octave_idx_type n = s.rows ();
 
-template <>
-octave_idx_type
-schur<ComplexMatrix>::init (const ComplexMatrix& a, const std::string& ord,
-                            bool calc_unitary)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
+      if (s.columns () != n || u.rows () != n || u.columns () != n)
+        (*current_liboctave_error_handler)
+          ("rsf2csf: inconsistent matrix dimensions");
 
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("SCHUR requires square matrix");
+      if (n > 0)
+        {
+          OCTAVE_LOCAL_BUFFER (double, c, n-1);
+          OCTAVE_LOCAL_BUFFER (double, sx, n-1);
 
-  if (a_nr == 0)
-    {
-      schur_mat.clear ();
-      unitary_mat.clear ();
-      return 0;
+          F77_XFCN (zrsf2csf, ZRSF2CSF, (n, F77_DBLE_CMPLX_ARG (s.fortran_vec ()),
+                                         F77_DBLE_CMPLX_ARG (u.fortran_vec ()), c, sx));
+        }
+
+      return schur<ComplexMatrix> (s, u);
     }
 
-  // Workspace requirements may need to be fixed if any of the
-  // following change.
-
-  char jobvs;
-  char sense = 'N';
-  char sort = 'N';
+    template <>
+    octave_idx_type
+    schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
+                                     const std::string& ord, bool calc_unitary)
+    {
+      octave_idx_type a_nr = a.rows ();
+      octave_idx_type a_nc = a.cols ();
 
-  if (calc_unitary)
-    jobvs = 'V';
-  else
-    jobvs = 'N';
-
-  char ord_char = ord.empty () ? 'U' : ord[0];
-
-  if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = 'S';
+      if (a_nr != a_nc)
+        (*current_liboctave_error_handler) ("SCHUR requires square matrix");
 
-  volatile complex_selector selector = 0;
-  if (ord_char == 'A' || ord_char == 'a')
-    selector = select_ana<Complex>;
-  else if (ord_char == 'D' || ord_char == 'd')
-    selector = select_dig<Complex>;
+      if (a_nr == 0)
+        {
+          schur_mat.clear ();
+          unitary_mat.clear ();
+          return 0;
+        }
+
+      // Workspace requirements may need to be fixed if any of the
+      // following change.
 
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 8 * n;
-  octave_idx_type info;
-  octave_idx_type sdim;
-  double rconde;
-  double rcondv;
+      char jobvs;
+      char sense = 'N';
+      char sort = 'N';
 
-  schur_mat = a;
-  if (calc_unitary)
-    unitary_mat.clear (n, n);
+      if (calc_unitary)
+        jobvs = 'V';
+      else
+        jobvs = 'N';
 
-  Complex *s = schur_mat.fortran_vec ();
-  Complex *q = unitary_mat.fortran_vec ();
+      char ord_char = ord.empty () ? 'U' : ord[0];
 
-  Array<double> rwork (dim_vector (n, 1));
-  double *prwork = rwork.fortran_vec ();
+      if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
+        sort = 'S';
 
-  Array<Complex> w (dim_vector (n, 1));
-  Complex *pw = w.fortran_vec ();
-
-  Array<Complex> work (dim_vector (lwork, 1));
-  Complex *pwork = work.fortran_vec ();
+      volatile float_complex_selector selector = 0;
+      if (ord_char == 'A' || ord_char == 'a')
+        selector = select_ana<FloatComplex>;
+      else if (ord_char == 'D' || ord_char == 'd')
+        selector = select_dig<FloatComplex>;
 
-  // BWORK is not referenced for non-ordered Schur.
-  octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
-  Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
-  octave_idx_type *pbwork = bwork.fortran_vec ();
+      octave_idx_type n = a_nc;
+      octave_idx_type lwork = 8 * n;
+      octave_idx_type info;
+      octave_idx_type sdim;
+      float rconde;
+      float rcondv;
 
-  F77_XFCN (zgeesx, ZGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
-                             F77_CONST_CHAR_ARG2 (&sort, 1),
-                             selector,
-                             F77_CONST_CHAR_ARG2 (&sense, 1),
-                             n, F77_DBLE_CMPLX_ARG (s), n, sdim, F77_DBLE_CMPLX_ARG (pw), F77_DBLE_CMPLX_ARG (q), n, rconde, rcondv,
-                             F77_DBLE_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
+      schur_mat = a;
+      if (calc_unitary)
+        unitary_mat.clear (n, n);
+
+      FloatComplex *s = schur_mat.fortran_vec ();
+      FloatComplex *q = unitary_mat.fortran_vec ();
+
+      Array<float> rwork (dim_vector (n, 1));
+      float *prwork = rwork.fortran_vec ();
 
-  return info;
-}
+      Array<FloatComplex> w (dim_vector (n, 1));
+      FloatComplex *pw = w.fortran_vec ();
 
-template <>
-schur<ComplexMatrix>
-rsf2csf<ComplexMatrix, Matrix> (const Matrix& s_arg, const Matrix& u_arg)
-{
-  ComplexMatrix s (s_arg);
-  ComplexMatrix u (u_arg);
+      Array<FloatComplex> work (dim_vector (lwork, 1));
+      FloatComplex *pwork = work.fortran_vec ();
 
-  octave_idx_type n = s.rows ();
+      // BWORK is not referenced for non-ordered Schur.
+      octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
+      Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
+      octave_idx_type *pbwork = bwork.fortran_vec ();
 
-  if (s.columns () != n || u.rows () != n || u.columns () != n)
-    (*current_liboctave_error_handler)
-      ("rsf2csf: inconsistent matrix dimensions");
+      F77_XFCN (cgeesx, CGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
+                                 F77_CONST_CHAR_ARG2 (&sort, 1),
+                                 selector,
+                                 F77_CONST_CHAR_ARG2 (&sense, 1),
+                                 n, F77_CMPLX_ARG (s), n, sdim, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (q), n, rconde, rcondv,
+                                 F77_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)
+                                 F77_CHAR_ARG_LEN (1)));
 
-  if (n > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (double, c, n-1);
-      OCTAVE_LOCAL_BUFFER (double, sx, n-1);
-
-      F77_XFCN (zrsf2csf, ZRSF2CSF, (n, F77_DBLE_CMPLX_ARG (s.fortran_vec ()),
-                                     F77_DBLE_CMPLX_ARG (u.fortran_vec ()), c, sx));
+      return info;
     }
 
-  return schur<ComplexMatrix> (s, u);
-}
+    template <>
+    schur<FloatComplexMatrix>
+    rsf2csf<FloatComplexMatrix, FloatMatrix> (const FloatMatrix& s_arg, const FloatMatrix& u_arg)
+    {
+      FloatComplexMatrix s (s_arg);
+      FloatComplexMatrix u (u_arg);
+
+      octave_idx_type n = s.rows ();
 
-template <>
-octave_idx_type
-schur<FloatComplexMatrix>::init (const FloatComplexMatrix& a,
-                                 const std::string& ord, bool calc_unitary)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
+      if (s.columns () != n || u.rows () != n || u.columns () != n)
+        (*current_liboctave_error_handler)
+          ("rsf2csf: inconsistent matrix dimensions");
 
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("SCHUR requires square matrix");
+      if (n > 0)
+        {
+          OCTAVE_LOCAL_BUFFER (float, c, n-1);
+          OCTAVE_LOCAL_BUFFER (float, sx, n-1);
 
-  if (a_nr == 0)
-    {
-      schur_mat.clear ();
-      unitary_mat.clear ();
-      return 0;
+          F77_XFCN (crsf2csf, CRSF2CSF, (n, F77_CMPLX_ARG (s.fortran_vec ()),
+                                         F77_CMPLX_ARG (u.fortran_vec ()), c, sx));
+        }
+
+      return schur<FloatComplexMatrix> (s, u);
     }
 
-  // Workspace requirements may need to be fixed if any of the
-  // following change.
-
-  char jobvs;
-  char sense = 'N';
-  char sort = 'N';
-
-  if (calc_unitary)
-    jobvs = 'V';
-  else
-    jobvs = 'N';
-
-  char ord_char = ord.empty () ? 'U' : ord[0];
-
-  if (ord_char == 'A' || ord_char == 'D' || ord_char == 'a' || ord_char == 'd')
-    sort = 'S';
-
-  volatile float_complex_selector selector = 0;
-  if (ord_char == 'A' || ord_char == 'a')
-    selector = select_ana<FloatComplex>;
-  else if (ord_char == 'D' || ord_char == 'd')
-    selector = select_dig<FloatComplex>;
+    // Instantiations we need.
 
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 8 * n;
-  octave_idx_type info;
-  octave_idx_type sdim;
-  float rconde;
-  float rcondv;
-
-  schur_mat = a;
-  if (calc_unitary)
-    unitary_mat.clear (n, n);
-
-  FloatComplex *s = schur_mat.fortran_vec ();
-  FloatComplex *q = unitary_mat.fortran_vec ();
+    template class schur<ComplexMatrix>;
 
-  Array<float> rwork (dim_vector (n, 1));
-  float *prwork = rwork.fortran_vec ();
-
-  Array<FloatComplex> w (dim_vector (n, 1));
-  FloatComplex *pw = w.fortran_vec ();
-
-  Array<FloatComplex> work (dim_vector (lwork, 1));
-  FloatComplex *pwork = work.fortran_vec ();
-
-  // BWORK is not referenced for non-ordered Schur.
-  octave_idx_type ntmp = (ord_char == 'N' || ord_char == 'n') ? 0 : n;
-  Array<octave_idx_type> bwork (dim_vector (ntmp, 1));
-  octave_idx_type *pbwork = bwork.fortran_vec ();
+    template class schur<FloatComplexMatrix>;
 
-  F77_XFCN (cgeesx, CGEESX, (F77_CONST_CHAR_ARG2 (&jobvs, 1),
-                             F77_CONST_CHAR_ARG2 (&sort, 1),
-                             selector,
-                             F77_CONST_CHAR_ARG2 (&sense, 1),
-                             n, F77_CMPLX_ARG (s), n, sdim, F77_CMPLX_ARG (pw), F77_CMPLX_ARG (q), n, rconde, rcondv,
-                             F77_CMPLX_ARG (pwork), lwork, prwork, pbwork, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
-
-  return info;
-}
-
-template <>
-schur<FloatComplexMatrix>
-rsf2csf<FloatComplexMatrix, FloatMatrix> (const FloatMatrix& s_arg, const FloatMatrix& u_arg)
-{
-  FloatComplexMatrix s (s_arg);
-  FloatComplexMatrix u (u_arg);
-
-  octave_idx_type n = s.rows ();
+    template class schur<FloatMatrix>;
 
-  if (s.columns () != n || u.rows () != n || u.columns () != n)
-    (*current_liboctave_error_handler)
-      ("rsf2csf: inconsistent matrix dimensions");
-
-  if (n > 0)
-    {
-      OCTAVE_LOCAL_BUFFER (float, c, n-1);
-      OCTAVE_LOCAL_BUFFER (float, sx, n-1);
-
-      F77_XFCN (crsf2csf, CRSF2CSF, (n, F77_CMPLX_ARG (s.fortran_vec ()),
-                                     F77_CMPLX_ARG (u.fortran_vec ()), c, sx));
-    }
-
-  return schur<FloatComplexMatrix> (s, u);
+    template class schur<Matrix>;
+  }
 }
-
-// Instantiations we need.
-
-template class schur<ComplexMatrix>;
-
-template class schur<FloatComplexMatrix>;
-
-template class schur<FloatMatrix>;
-
-template class schur<Matrix>;
-
-}
-}