# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1455593625 18000
# Node ID f5b8c3aca5f866e0f14492405dc370c2f70a8c5c
# Parent  e69eaee287374979d7d9378cdde4f5e47bca5067
make better use of templates for Hessenberg decomposition

* liboctave/numeric/hess.h, liboctave/numeric/hess.cc:
New files generated from HESS.h, HESS.cc, CmplxHESS.h,
CmplxHESS.cc, dbleHESS.h, dbleHESS.cc, fCmplxHESS.h,
fCmplxHESS.cc, floatHESS.h, and floatHESS.cc and making them
templates.
* liboctave/numeric/module.mk: Update.
* libinterp/corefcn/hess.cc, mx-defs.h, mx-ext.h: Use new template
classes and header file.

diff -r e69eaee28737 -r f5b8c3aca5f8 libinterp/corefcn/hess.cc
--- a/libinterp/corefcn/hess.cc	Mon Feb 15 20:06:12 2016 -0500
+++ b/libinterp/corefcn/hess.cc	Mon Feb 15 22:33:45 2016 -0500
@@ -24,10 +24,7 @@
 #  include <config.h>
 #endif
 
-#include "CmplxHESS.h"
-#include "dbleHESS.h"
-#include "fCmplxHESS.h"
-#include "floatHESS.h"
+#include "hess.h"
 
 #include "defun.h"
 #include "error.h"
@@ -90,7 +87,7 @@
         {
           FloatMatrix tmp = arg.float_matrix_value ();
 
-          FloatHESS result (tmp);
+          hess<FloatMatrix> result (tmp);
 
           if (nargout <= 1)
             retval = ovl (result.hess_matrix ());
@@ -102,7 +99,7 @@
         {
           FloatComplexMatrix ctmp = arg.float_complex_matrix_value ();
 
-          FloatComplexHESS result (ctmp);
+          hess<FloatComplexMatrix> result (ctmp);
 
           if (nargout <= 1)
             retval = ovl (result.hess_matrix ());
@@ -117,7 +114,7 @@
         {
           Matrix tmp = arg.matrix_value ();
 
-          HESS result (tmp);
+          hess<Matrix> result (tmp);
 
           if (nargout <= 1)
             retval = ovl (result.hess_matrix ());
@@ -129,7 +126,7 @@
         {
           ComplexMatrix ctmp = arg.complex_matrix_value ();
 
-          ComplexHESS result (ctmp);
+          hess<ComplexMatrix> result (ctmp);
 
           if (nargout <= 1)
             retval = ovl (result.hess_matrix ());
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/CmplxHESS.cc
--- a/liboctave/numeric/CmplxHESS.cc	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,121 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "CmplxHESS.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, octave_idx_type&,
-                             octave_idx_type&, double*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (zgehrd, ZGEHRD) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*, Complex*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (zunghr, ZUNGHR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, Complex*, Complex*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, Complex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-octave_idx_type
-ComplexHESS::init (const ComplexMatrix& a)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("ComplexHESS requires square matrix");
-
-  char job = 'N';
-  char side = 'R';
-
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 32 * n;
-  octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
-
-  hess_mat = a;
-  Complex *h = hess_mat.fortran_vec ();
-
-  Array<double> scale (dim_vector (n, 1));
-  double *pscale = scale.fortran_vec ();
-
-  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             n, h, n, ilo, ihi, pscale, info
-                             F77_CHAR_ARG_LEN (1)));
-
-  Array<Complex> tau (dim_vector (n-1, 1));
-  Complex *ptau = tau.fortran_vec ();
-
-  Array<Complex> work (dim_vector (lwork, 1));
-  Complex *pwork = work.fortran_vec ();
-
-  F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
-
-  unitary_hess_mat = hess_mat;
-  Complex *z = unitary_hess_mat.fortran_vec ();
-
-  F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
-                             lwork, info));
-
-  F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             F77_CONST_CHAR_ARG2 (&side, 1),
-                             n, ilo, ihi, pscale, n, z, n, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of
-  // doing this (or faster for that matter :-)),
-  // please let me know!
-
-  if (n > 2)
-    for (octave_idx_type j = 0; j < a_nc; j++)
-      for (octave_idx_type i = j+2; i < a_nr; i++)
-        hess_mat.elem (i, j) = 0;
-
-  return info;
-}
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/CmplxHESS.h
--- a/liboctave/numeric/CmplxHESS.h	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,84 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_CmplxHESS_h)
-#define octave_CmplxHESS_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "CMatrix.h"
-
-class
-OCTAVE_API
-ComplexHESS
-{
-public:
-
-  ComplexHESS (void) : hess_mat (), unitary_hess_mat () { }
-
-  ComplexHESS (const ComplexMatrix& a)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    init (a);
-  }
-
-  ComplexHESS (const ComplexMatrix& a, octave_idx_type& info)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    info = init (a);
-  }
-
-  ComplexHESS (const ComplexHESS& a)
-    : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { }
-
-  ComplexHESS& operator = (const ComplexHESS& a)
-  {
-    if (this != &a)
-      {
-        hess_mat = a.hess_mat;
-        unitary_hess_mat = a.unitary_hess_mat;
-      }
-    return *this;
-  }
-
-  ~ComplexHESS (void) { }
-
-  ComplexMatrix hess_matrix (void) const { return hess_mat; }
-
-  ComplexMatrix unitary_hess_matrix (void) const
-  {
-    return unitary_hess_mat;
-  }
-
-  friend std::ostream& operator << (std::ostream& os, const ComplexHESS& a);
-
-private:
-
-  ComplexMatrix hess_mat;
-  ComplexMatrix unitary_hess_mat;
-
-  octave_idx_type init (const ComplexMatrix& a);
-};
-
-#endif
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/dbleHESS.cc
--- a/liboctave/numeric/dbleHESS.cc	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,123 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "dbleHESS.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type&,
-                             octave_idx_type&, double*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (dgehrd, DGEHRD) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dorghr, DORGHR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*, double*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, double*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-octave_idx_type
-HESS::init (const Matrix& a)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("HESS requires square matrix");
-
-  char job = 'N';
-  char side = 'R';
-
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 32 * n;
-  octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
-
-  hess_mat = a;
-  double *h = hess_mat.fortran_vec ();
-
-  Array<double> scale (dim_vector (n, 1));
-  double *pscale = scale.fortran_vec ();
-
-  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             n, h, n, ilo, ihi, pscale, info
-                             F77_CHAR_ARG_LEN (1)));
-
-  Array<double> tau (dim_vector (n-1, 1));
-  double *ptau = tau.fortran_vec ();
-
-  Array<double> work (dim_vector (lwork, 1));
-  double *pwork = work.fortran_vec ();
-
-  F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
-                             lwork, info));
-
-  unitary_hess_mat = hess_mat;
-  double *z = unitary_hess_mat.fortran_vec ();
-
-  F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork,
-                             lwork, info));
-
-  F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             F77_CONST_CHAR_ARG2 (&side, 1),
-                             n, ilo, ihi, pscale, n, z,
-                             n, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of doing
-  // this (or faster for that matter :-)), please let
-  // me know!
-
-  if (n > 2)
-    for (octave_idx_type j = 0; j < a_nc; j++)
-      for (octave_idx_type i = j+2; i < a_nr; i++)
-        hess_mat.elem (i, j) = 0;
-
-  return info;
-}
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/dbleHESS.h
--- a/liboctave/numeric/dbleHESS.h	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,77 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_dbleHESS_h)
-#define octave_dbleHESS_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "dMatrix.h"
-
-class
-OCTAVE_API
-HESS
-{
-public:
-
-  HESS (void) : hess_mat (), unitary_hess_mat () { }
-
-  HESS (const Matrix& a) : hess_mat (), unitary_hess_mat () { init (a); }
-
-  HESS (const Matrix& a, octave_idx_type& info)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    info = init (a);
-  }
-
-  HESS (const HESS& a)
-    : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { }
-
-  HESS& operator = (const HESS& a)
-  {
-    if (this != &a)
-      {
-        hess_mat = a.hess_mat;
-        unitary_hess_mat = a.unitary_hess_mat;
-      }
-    return *this;
-  }
-
-  ~HESS (void) { }
-
-  Matrix hess_matrix (void) const { return hess_mat; }
-
-  Matrix unitary_hess_matrix (void) const { return unitary_hess_mat; }
-
-  friend std::ostream& operator << (std::ostream& os, const HESS& a);
-
-private:
-
-  Matrix hess_mat;
-  Matrix unitary_hess_mat;
-
-  octave_idx_type init (const Matrix& a);
-};
-
-#endif
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/fCmplxHESS.cc
--- a/liboctave/numeric/fCmplxHESS.cc	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,127 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "fCmplxHESS.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, octave_idx_type&,
-                             octave_idx_type&, float*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (cgehrd, CGEHRD) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             FloatComplex*, const octave_idx_type&,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cunghr, CUNGHR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, FloatComplex*,
-                             FloatComplex*, const octave_idx_type&,
-                             octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, FloatComplex*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-octave_idx_type
-FloatComplexHESS::init (const FloatComplexMatrix& a)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    {
-      (*current_liboctave_error_handler)
-        ("FloatComplexHESS requires square matrix");
-      return -1;
-    }
-
-  char job = 'N';
-  char side = 'R';
-
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 32 * n;
-  octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
-
-  hess_mat = a;
-  FloatComplex *h = hess_mat.fortran_vec ();
-
-  Array<float> scale (dim_vector (n, 1));
-  float *pscale = scale.fortran_vec ();
-
-  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             n, h, n, ilo, ihi, pscale, info
-                             F77_CHAR_ARG_LEN (1)));
-
-  Array<FloatComplex> tau (dim_vector (n-1, 1));
-  FloatComplex *ptau = tau.fortran_vec ();
-
-  Array<FloatComplex> work (dim_vector (lwork, 1));
-  FloatComplex *pwork = work.fortran_vec ();
-
-  F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
-
-  unitary_hess_mat = hess_mat;
-  FloatComplex *z = unitary_hess_mat.fortran_vec ();
-
-  F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
-                             lwork, info));
-
-  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             F77_CONST_CHAR_ARG2 (&side, 1),
-                             n, ilo, ihi, pscale, n, z, n, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of
-  // doing this (or faster for that matter :-)),
-  // please let me know!
-
-  if (n > 2)
-    for (octave_idx_type j = 0; j < a_nc; j++)
-      for (octave_idx_type i = j+2; i < a_nr; i++)
-        hess_mat.elem (i, j) = 0;
-
-  return info;
-}
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/fCmplxHESS.h
--- a/liboctave/numeric/fCmplxHESS.h	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,85 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_fCmplxHESS_h)
-#define octave_fCmplxHESS_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fCMatrix.h"
-
-class
-OCTAVE_API
-FloatComplexHESS
-{
-public:
-
-  FloatComplexHESS (void) : hess_mat (), unitary_hess_mat () { }
-
-  FloatComplexHESS (const FloatComplexMatrix& a)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    init (a);
-  }
-
-  FloatComplexHESS (const FloatComplexMatrix& a, octave_idx_type& info)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    info = init (a);
-  }
-
-  FloatComplexHESS (const FloatComplexHESS& a)
-    : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { }
-
-  FloatComplexHESS& operator = (const FloatComplexHESS& a)
-  {
-    if (this != &a)
-      {
-        hess_mat = a.hess_mat;
-        unitary_hess_mat = a.unitary_hess_mat;
-      }
-    return *this;
-  }
-
-  ~FloatComplexHESS (void) { }
-
-  FloatComplexMatrix hess_matrix (void) const { return hess_mat; }
-
-  FloatComplexMatrix unitary_hess_matrix (void) const
-  {
-    return unitary_hess_mat;
-  }
-
-  friend std::ostream& operator << (std::ostream& os,
-                                    const FloatComplexHESS& a);
-
-private:
-
-  FloatComplexMatrix hess_mat;
-  FloatComplexMatrix unitary_hess_mat;
-
-  octave_idx_type init (const FloatComplexMatrix& a);
-};
-
-#endif
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/floatHESS.cc
--- a/liboctave/numeric/floatHESS.cc	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,123 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#  include <config.h>
-#endif
-
-#include "floatHESS.h"
-#include "f77-fcn.h"
-#include "lo-error.h"
-
-extern "C"
-{
-  F77_RET_T
-  F77_FUNC (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type&,
-                             octave_idx_type&, float*, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL);
-
-  F77_RET_T
-  F77_FUNC (sgehrd, SGEHRD) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (sorghr, SORGHR) (const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*, float*,
-                             const octave_idx_type&, octave_idx_type&);
-
-  F77_RET_T
-  F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL,
-                             F77_CONST_CHAR_ARG_DECL,
-                             const octave_idx_type&, const octave_idx_type&,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, float*,
-                             const octave_idx_type&, octave_idx_type&
-                             F77_CHAR_ARG_LEN_DECL
-                             F77_CHAR_ARG_LEN_DECL);
-}
-
-octave_idx_type
-FloatHESS::init (const FloatMatrix& a)
-{
-  octave_idx_type a_nr = a.rows ();
-  octave_idx_type a_nc = a.cols ();
-
-  if (a_nr != a_nc)
-    (*current_liboctave_error_handler) ("FloatHESS requires square matrix");
-
-  char job = 'N';
-  char side = 'R';
-
-  octave_idx_type n = a_nc;
-  octave_idx_type lwork = 32 * n;
-  octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
-
-  hess_mat = a;
-  float *h = hess_mat.fortran_vec ();
-
-  Array<float> scale (dim_vector (n, 1));
-  float *pscale = scale.fortran_vec ();
-
-  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             n, h, n, ilo, ihi, pscale, info
-                             F77_CHAR_ARG_LEN (1)));
-
-  Array<float> tau (dim_vector (n-1, 1));
-  float *ptau = tau.fortran_vec ();
-
-  Array<float> work (dim_vector (lwork, 1));
-  float *pwork = work.fortran_vec ();
-
-  F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
-                             lwork, info));
-
-  unitary_hess_mat = hess_mat;
-  float *z = unitary_hess_mat.fortran_vec ();
-
-  F77_XFCN (sorghr, SORGHR, (n, ilo, ihi, z, n, ptau, pwork,
-                             lwork, info));
-
-  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
-                             F77_CONST_CHAR_ARG2 (&side, 1),
-                             n, ilo, ihi, pscale, n, z,
-                             n, info
-                             F77_CHAR_ARG_LEN (1)
-                             F77_CHAR_ARG_LEN (1)));
-
-  // If someone thinks of a more graceful way of doing
-  // this (or faster for that matter :-)), please let
-  // me know!
-
-  if (n > 2)
-    for (octave_idx_type j = 0; j < a_nc; j++)
-      for (octave_idx_type i = j+2; i < a_nr; i++)
-        hess_mat.elem (i, j) = 0;
-
-  return info;
-}
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/floatHESS.h
--- a/liboctave/numeric/floatHESS.h	Mon Feb 15 20:06:12 2016 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,81 +0,0 @@
-/*
-
-Copyright (C) 1994-2015 John W. Eaton
-
-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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#if ! defined (octave_floatHESS_h)
-#define octave_floatHESS_h 1
-
-#include "octave-config.h"
-
-#include <iosfwd>
-
-#include "fMatrix.h"
-
-class
-OCTAVE_API
-FloatHESS
-{
-public:
-
-  FloatHESS (void) : hess_mat (), unitary_hess_mat () { }
-
-  FloatHESS (const FloatMatrix& a)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    init (a);
-  }
-
-  FloatHESS (const FloatMatrix& a, octave_idx_type& info)
-    : hess_mat (), unitary_hess_mat ()
-  {
-    info = init (a);
-  }
-
-  FloatHESS (const FloatHESS& a)
-    : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat) { }
-
-  FloatHESS& operator = (const FloatHESS& a)
-  {
-    if (this != &a)
-      {
-        hess_mat = a.hess_mat;
-        unitary_hess_mat = a.unitary_hess_mat;
-      }
-    return *this;
-  }
-
-  ~FloatHESS (void) { }
-
-  FloatMatrix hess_matrix (void) const { return hess_mat; }
-
-  FloatMatrix unitary_hess_matrix (void) const { return unitary_hess_mat; }
-
-  friend std::ostream& operator << (std::ostream& os, const FloatHESS& a);
-
-private:
-
-  FloatMatrix hess_mat;
-  FloatMatrix unitary_hess_mat;
-
-  octave_idx_type init (const FloatMatrix& a);
-};
-
-#endif
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/hess.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/hess.cc	Mon Feb 15 22:33:45 2016 -0500
@@ -0,0 +1,406 @@
+/*
+
+Copyright (C) 1994-2015 John W. Eaton
+
+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
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#  include <config.h>
+#endif
+
+
+#include "CMatrix.h"
+#include "dMatrix.h"
+#include "f77-fcn.h"
+#include "fCMatrix.h"
+#include "fMatrix.h"
+#include "hess.h"
+#include "lo-error.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dgebal, DGEBAL) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, double*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dgehrd, DGEHRD) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dorghr, DORGHR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*, double*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (dgebak, DGEBAK) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sgebal, SGEBAL) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, float*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sgehrd, SGEHRD) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sorghr, SORGHR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*, float*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sgebak, SGEBAK) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zgebal, ZGEBAL) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, double*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (zgehrd, ZGEHRD) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*, Complex*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zunghr, ZUNGHR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, Complex*, Complex*,
+                             const octave_idx_type&, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (zgebak, ZGEBAK) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, double*,
+                             const octave_idx_type&, Complex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cgebal, CGEBAL) (F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&,
+                             octave_idx_type&, float*, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (cgehrd, CGEHRD) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             FloatComplex*, const octave_idx_type&,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cunghr, CUNGHR) (const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, FloatComplex*,
+                             FloatComplex*, const octave_idx_type&,
+                             octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cgebak, CGEBAK) (F77_CONST_CHAR_ARG_DECL,
+                             F77_CONST_CHAR_ARG_DECL,
+                             const octave_idx_type&, const octave_idx_type&,
+                             const octave_idx_type&, float*,
+                             const octave_idx_type&, FloatComplex*,
+                             const octave_idx_type&, octave_idx_type&
+                             F77_CHAR_ARG_LEN_DECL
+                             F77_CHAR_ARG_LEN_DECL);
+}
+
+template <>
+octave_idx_type
+hess<Matrix>::init (const Matrix& a)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("hess: requires square matrix");
+
+  char job = 'N';
+  char side = 'R';
+
+  octave_idx_type n = a_nc;
+  octave_idx_type lwork = 32 * n;
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  hess_mat = a;
+  double *h = hess_mat.fortran_vec ();
+
+  Array<double> scale (dim_vector (n, 1));
+  double *pscale = scale.fortran_vec ();
+
+  F77_XFCN (dgebal, DGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             n, h, n, ilo, ihi, pscale, info
+                             F77_CHAR_ARG_LEN (1)));
+
+  Array<double> tau (dim_vector (n-1, 1));
+  double *ptau = tau.fortran_vec ();
+
+  Array<double> work (dim_vector (lwork, 1));
+  double *pwork = work.fortran_vec ();
+
+  F77_XFCN (dgehrd, DGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
+                             lwork, info));
+
+  unitary_hess_mat = hess_mat;
+  double *z = unitary_hess_mat.fortran_vec ();
+
+  F77_XFCN (dorghr, DORGHR, (n, ilo, ihi, z, n, ptau, pwork,
+                             lwork, info));
+
+  F77_XFCN (dgebak, DGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             F77_CONST_CHAR_ARG2 (&side, 1),
+                             n, ilo, ihi, pscale, n, z,
+                             n, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of doing
+  // this (or faster for that matter :-)), please let
+  // me know!
+
+  if (n > 2)
+    for (octave_idx_type j = 0; j < a_nc; j++)
+      for (octave_idx_type i = j+2; i < a_nr; i++)
+        hess_mat.elem (i, j) = 0;
+
+  return info;
+}
+
+template <>
+octave_idx_type
+hess<FloatMatrix>::init (const FloatMatrix& a)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("hess: requires square matrix");
+
+  char job = 'N';
+  char side = 'R';
+
+  octave_idx_type n = a_nc;
+  octave_idx_type lwork = 32 * n;
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  hess_mat = a;
+  float *h = hess_mat.fortran_vec ();
+
+  Array<float> scale (dim_vector (n, 1));
+  float *pscale = scale.fortran_vec ();
+
+  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             n, h, n, ilo, ihi, pscale, info
+                             F77_CHAR_ARG_LEN (1)));
+
+  Array<float> tau (dim_vector (n-1, 1));
+  float *ptau = tau.fortran_vec ();
+
+  Array<float> work (dim_vector (lwork, 1));
+  float *pwork = work.fortran_vec ();
+
+  F77_XFCN (sgehrd, SGEHRD, (n, ilo, ihi, h, n, ptau, pwork,
+                             lwork, info));
+
+  unitary_hess_mat = hess_mat;
+  float *z = unitary_hess_mat.fortran_vec ();
+
+  F77_XFCN (sorghr, SORGHR, (n, ilo, ihi, z, n, ptau, pwork,
+                             lwork, info));
+
+  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             F77_CONST_CHAR_ARG2 (&side, 1),
+                             n, ilo, ihi, pscale, n, z,
+                             n, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of doing
+  // this (or faster for that matter :-)), please let
+  // me know!
+
+  if (n > 2)
+    for (octave_idx_type j = 0; j < a_nc; j++)
+      for (octave_idx_type i = j+2; i < a_nr; i++)
+        hess_mat.elem (i, j) = 0;
+
+  return info;
+}
+
+template <>
+octave_idx_type
+hess<ComplexMatrix>::init (const ComplexMatrix& a)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    (*current_liboctave_error_handler) ("hess: requires square matrix");
+
+  char job = 'N';
+  char side = 'R';
+
+  octave_idx_type n = a_nc;
+  octave_idx_type lwork = 32 * n;
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  hess_mat = a;
+  Complex *h = hess_mat.fortran_vec ();
+
+  Array<double> scale (dim_vector (n, 1));
+  double *pscale = scale.fortran_vec ();
+
+  F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             n, h, n, ilo, ihi, pscale, info
+                             F77_CHAR_ARG_LEN (1)));
+
+  Array<Complex> tau (dim_vector (n-1, 1));
+  Complex *ptau = tau.fortran_vec ();
+
+  Array<Complex> work (dim_vector (lwork, 1));
+  Complex *pwork = work.fortran_vec ();
+
+  F77_XFCN (zgehrd, ZGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
+
+  unitary_hess_mat = hess_mat;
+  Complex *z = unitary_hess_mat.fortran_vec ();
+
+  F77_XFCN (zunghr, ZUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
+                             lwork, info));
+
+  F77_XFCN (zgebak, ZGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             F77_CONST_CHAR_ARG2 (&side, 1),
+                             n, ilo, ihi, pscale, n, z, n, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of
+  // doing this (or faster for that matter :-)),
+  // please let me know!
+
+  if (n > 2)
+    for (octave_idx_type j = 0; j < a_nc; j++)
+      for (octave_idx_type i = j+2; i < a_nr; i++)
+        hess_mat.elem (i, j) = 0;
+
+  return info;
+}
+
+template <>
+octave_idx_type
+hess<FloatComplexMatrix>::init (const FloatComplexMatrix& a)
+{
+  octave_idx_type a_nr = a.rows ();
+  octave_idx_type a_nc = a.cols ();
+
+  if (a_nr != a_nc)
+    {
+      (*current_liboctave_error_handler) ("hess: requires square matrix");
+      return -1;
+    }
+
+  char job = 'N';
+  char side = 'R';
+
+  octave_idx_type n = a_nc;
+  octave_idx_type lwork = 32 * n;
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  hess_mat = a;
+  FloatComplex *h = hess_mat.fortran_vec ();
+
+  Array<float> scale (dim_vector (n, 1));
+  float *pscale = scale.fortran_vec ();
+
+  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             n, h, n, ilo, ihi, pscale, info
+                             F77_CHAR_ARG_LEN (1)));
+
+  Array<FloatComplex> tau (dim_vector (n-1, 1));
+  FloatComplex *ptau = tau.fortran_vec ();
+
+  Array<FloatComplex> work (dim_vector (lwork, 1));
+  FloatComplex *pwork = work.fortran_vec ();
+
+  F77_XFCN (cgehrd, CGEHRD, (n, ilo, ihi, h, n, ptau, pwork, lwork, info));
+
+  unitary_hess_mat = hess_mat;
+  FloatComplex *z = unitary_hess_mat.fortran_vec ();
+
+  F77_XFCN (cunghr, CUNGHR, (n, ilo, ihi, z, n, ptau, pwork,
+                             lwork, info));
+
+  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+                             F77_CONST_CHAR_ARG2 (&side, 1),
+                             n, ilo, ihi, pscale, n, z, n, info
+                             F77_CHAR_ARG_LEN (1)
+                             F77_CHAR_ARG_LEN (1)));
+
+  // If someone thinks of a more graceful way of
+  // doing this (or faster for that matter :-)),
+  // please let me know!
+
+  if (n > 2)
+    for (octave_idx_type j = 0; j < a_nc; j++)
+      for (octave_idx_type i = j+2; i < a_nr; i++)
+        hess_mat.elem (i, j) = 0;
+
+  return info;
+}
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/hess.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/numeric/hess.h	Mon Feb 15 22:33:45 2016 -0500
@@ -0,0 +1,85 @@
+/*
+
+Copyright (C) 1994-2015 John W. Eaton
+
+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
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#if ! defined (octave_hess_h)
+#define octave_hess_h 1
+
+#include "octave-config.h"
+
+#include <iosfwd>
+
+template <typename T>
+class
+hess
+{
+public:
+
+  hess (void)
+    : hess_mat (), unitary_hess_mat ()
+  { }
+
+  hess (const T& a)
+    : hess_mat (), unitary_hess_mat ()
+  {
+    init (a);
+  }
+
+  hess (const T& a, octave_idx_type& info)
+    : hess_mat (), unitary_hess_mat ()
+  {
+    info = init (a);
+  }
+
+  hess (const hess& a)
+    : hess_mat (a.hess_mat), unitary_hess_mat (a.unitary_hess_mat)
+  { }
+
+  hess& operator = (const hess& a)
+  {
+    if (this != &a)
+      {
+        hess_mat = a.hess_mat;
+        unitary_hess_mat = a.unitary_hess_mat;
+      }
+
+    return *this;
+  }
+
+  ~hess (void) { }
+
+  T hess_matrix (void) const { return hess_mat; }
+
+  T unitary_hess_matrix (void) const { return unitary_hess_mat; }
+
+private:
+
+  T hess_mat;
+  T unitary_hess_mat;
+
+  octave_idx_type init (const T& a);
+};
+
+template <typename T>
+extern std::ostream&
+operator << (std::ostream& os, const hess<T>& a);
+
+#endif
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/numeric/module.mk
--- a/liboctave/numeric/module.mk	Mon Feb 15 20:06:12 2016 -0500
+++ b/liboctave/numeric/module.mk	Mon Feb 15 22:33:45 2016 -0500
@@ -11,7 +11,6 @@
   liboctave/numeric/CmplxAEPBAL.h \
   liboctave/numeric/CmplxCHOL.h \
   liboctave/numeric/CmplxGEPBAL.h \
-  liboctave/numeric/CmplxHESS.h \
   liboctave/numeric/CmplxLU.h \
   liboctave/numeric/CmplxQR.h \
   liboctave/numeric/CmplxQRP.h \
@@ -43,7 +42,6 @@
   liboctave/numeric/dbleAEPBAL.h \
   liboctave/numeric/dbleCHOL.h \
   liboctave/numeric/dbleGEPBAL.h \
-  liboctave/numeric/dbleHESS.h \
   liboctave/numeric/dbleLU.h \
   liboctave/numeric/dbleQR.h \
   liboctave/numeric/dbleQRP.h \
@@ -52,7 +50,6 @@
   liboctave/numeric/fCmplxAEPBAL.h \
   liboctave/numeric/fCmplxCHOL.h \
   liboctave/numeric/fCmplxGEPBAL.h \
-  liboctave/numeric/fCmplxHESS.h \
   liboctave/numeric/fCmplxLU.h \
   liboctave/numeric/fCmplxQR.h \
   liboctave/numeric/fCmplxQRP.h \
@@ -61,11 +58,11 @@
   liboctave/numeric/floatAEPBAL.h \
   liboctave/numeric/floatCHOL.h \
   liboctave/numeric/floatGEPBAL.h \
-  liboctave/numeric/floatHESS.h \
   liboctave/numeric/floatLU.h \
   liboctave/numeric/floatQR.h \
   liboctave/numeric/floatQRP.h \
   liboctave/numeric/floatSVD.h \
+  liboctave/numeric/hess.h \
   liboctave/numeric/lo-mappers.h \
   liboctave/numeric/lo-specfun.h \
   liboctave/numeric/oct-convn.h \
@@ -91,7 +88,6 @@
   liboctave/numeric/CmplxAEPBAL.cc \
   liboctave/numeric/CmplxCHOL.cc \
   liboctave/numeric/CmplxGEPBAL.cc \
-  liboctave/numeric/CmplxHESS.cc \
   liboctave/numeric/CmplxLU.cc \
   liboctave/numeric/CmplxQR.cc \
   liboctave/numeric/CmplxQRP.cc \
@@ -107,7 +103,6 @@
   liboctave/numeric/dbleAEPBAL.cc \
   liboctave/numeric/dbleCHOL.cc \
   liboctave/numeric/dbleGEPBAL.cc \
-  liboctave/numeric/dbleHESS.cc \
   liboctave/numeric/dbleLU.cc \
   liboctave/numeric/dbleQR.cc \
   liboctave/numeric/dbleQRP.cc \
@@ -116,7 +111,6 @@
   liboctave/numeric/fCmplxAEPBAL.cc \
   liboctave/numeric/fCmplxCHOL.cc \
   liboctave/numeric/fCmplxGEPBAL.cc \
-  liboctave/numeric/fCmplxHESS.cc \
   liboctave/numeric/fCmplxLU.cc \
   liboctave/numeric/fCmplxQR.cc \
   liboctave/numeric/fCmplxQRP.cc \
@@ -125,11 +119,11 @@
   liboctave/numeric/floatAEPBAL.cc \
   liboctave/numeric/floatCHOL.cc \
   liboctave/numeric/floatGEPBAL.cc \
-  liboctave/numeric/floatHESS.cc \
   liboctave/numeric/floatLU.cc \
   liboctave/numeric/floatQR.cc \
   liboctave/numeric/floatQRP.cc \
   liboctave/numeric/floatSVD.cc \
+  liboctave/numeric/hess.cc \
   liboctave/numeric/lo-mappers.cc \
   liboctave/numeric/lo-specfun.cc \
   liboctave/numeric/oct-convn.cc \
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/operators/mx-defs.h
--- a/liboctave/operators/mx-defs.h	Mon Feb 15 20:06:12 2016 -0500
+++ b/liboctave/operators/mx-defs.h	Mon Feb 15 22:33:45 2016 -0500
@@ -75,10 +75,7 @@
 
 class EIG;
 
-class HESS;
-class ComplexHESS;
-class FloatHESS;
-class FloatComplexHESS;
+template <typename T> class hess;
 
 template <typename T> class schur;
 
diff -r e69eaee28737 -r f5b8c3aca5f8 liboctave/operators/mx-ext.h
--- a/liboctave/operators/mx-ext.h	Mon Feb 15 20:06:12 2016 -0500
+++ b/liboctave/operators/mx-ext.h	Mon Feb 15 22:33:45 2016 -0500
@@ -43,8 +43,7 @@
 
 // Result of a Hessenberg Decomposition
 
-#include "dbleHESS.h"
-#include "CmplxHESS.h"
+#include "hess.h"
 
 // Result of a Schur Decomposition