# HG changeset patch
# User John W. Eaton <jwe@octave.org>
# Date 1211377006 14400
# Node ID 39c1026191e93605b94bad19be3bae85c80e10ff
# Parent  975e9540be2c29340369f2e4450c6f67d3e329cb
add missing files from single-precision merge

diff --git a/liboctave/CmplxGEPBAL.cc b/liboctave/CmplxGEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/CmplxGEPBAL.cc
@@ -0,0 +1,130 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+#include <vector>
+
+#include "CmplxGEPBAL.h"
+#include "Array-util.h"
+#include "f77-fcn.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (zggbal, ZGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N,
+			     Complex* A, const octave_idx_type& LDA, Complex* B,
+			     const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI,
+			     double* LSCALE, double* RSCALE,
+			     double* WORK, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
+			     F77_CONST_CHAR_ARG_DECL,
+			     const octave_idx_type& N, const octave_idx_type& ILO,
+			     const octave_idx_type& IHI, const double* LSCALE,
+			     const double* RSCALE, octave_idx_type& M, double* V,
+			     const octave_idx_type& LDV, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+
+}
+
+octave_idx_type
+ComplexGEPBALANCE::init (const ComplexMatrix& a, const ComplexMatrix& b, 
+		  const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("ComplexGEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  if (a.dims() != b.dims ())
+    {
+      gripe_nonconformant ("ComplexGEPBALANCE", n, n, b.rows(), b.cols());
+      return -1;
+    } 
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  OCTAVE_LOCAL_BUFFER (double, plscale, n);
+  OCTAVE_LOCAL_BUFFER (double, prscale,  n);
+  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
+
+  balanced_mat = a;
+  Complex *p_balanced_mat = balanced_mat.fortran_vec ();
+  balanced_mat2 = b;
+  Complex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (zggbal, ZGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, p_balanced_mat2,
+			     n, ilo, ihi, plscale, prscale, pwork, info
+			     F77_CHAR_ARG_LEN (1)));
+
+  balancing_mat = Matrix (n, n, 0.0);
+  balancing_mat2 = Matrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+      balancing_mat.elem (i ,i) = 1.0;
+      balancing_mat2.elem (i ,i) = 1.0;
+    }
+
+  double *p_balancing_mat = balancing_mat.fortran_vec ();
+  double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
+
+  // first left
+  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+      
+  // then right
+  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("R", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat2, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/CmplxGEPBAL.h b/liboctave/CmplxGEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/CmplxGEPBAL.h
@@ -0,0 +1,91 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_ComplexGEPBALANCE_h)
+#define octave_ComplexGEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "CMatrix.h"
+#include "dMatrix.h"
+
+class
+OCTAVE_API
+ComplexGEPBALANCE
+{
+public:
+
+  ComplexGEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  ComplexGEPBALANCE (const ComplexMatrix& a, const ComplexMatrix& b, const std::string& balance_job)
+    {
+      init (a, b, balance_job); 
+    }
+
+  ComplexGEPBALANCE (const ComplexGEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2),
+    balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { }
+
+  ComplexGEPBALANCE& operator = (const ComplexGEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balanced_mat2 = a.balanced_mat2;
+	  balancing_mat = a.balancing_mat;
+	  balancing_mat2 = a.balancing_mat2;
+	}
+      return *this;
+    }
+
+  ~ComplexGEPBALANCE (void) { }
+
+  ComplexMatrix balanced_matrix (void) const { return balanced_mat; }
+
+  ComplexMatrix balanced_matrix2 (void) const { return balanced_mat2; }
+
+  Matrix balancing_matrix (void) const { return balancing_mat; }
+
+  Matrix balancing_matrix2 (void) const { return balancing_mat2; }
+
+  friend std::ostream& operator << (std::ostream& os, const ComplexGEPBALANCE& a);
+
+private:
+
+  ComplexMatrix balanced_mat;
+  ComplexMatrix balanced_mat2;
+  Matrix balancing_mat;
+  Matrix balancing_mat2;
+
+  octave_idx_type init (const ComplexMatrix& a, const ComplexMatrix& b, 
+			const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/dbleGEPBAL.cc b/liboctave/dbleGEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/dbleGEPBAL.cc
@@ -0,0 +1,130 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+#include <vector>
+
+#include "dbleGEPBAL.h"
+#include "Array-util.h"
+#include "f77-fcn.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (dggbal, DGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N,
+			     double* A, const octave_idx_type& LDA, double* B,
+			     const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI,
+			     double* LSCALE, double* RSCALE,
+			     double* WORK, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (dggbak, DGGBAK) (F77_CONST_CHAR_ARG_DECL,
+			     F77_CONST_CHAR_ARG_DECL,
+			     const octave_idx_type& N, const octave_idx_type& ILO,
+			     const octave_idx_type& IHI, const double* LSCALE,
+			     const double* RSCALE, octave_idx_type& M, double* V,
+			     const octave_idx_type& LDV, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+
+}
+
+octave_idx_type
+GEPBALANCE::init (const Matrix& a, const Matrix& b, 
+		  const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("GEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  if (a.dims() != b.dims ())
+    {
+      gripe_nonconformant ("GEPBALANCE", n, n, b.rows(), b.cols());
+      return -1;
+    } 
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  OCTAVE_LOCAL_BUFFER (double, plscale, n);
+  OCTAVE_LOCAL_BUFFER (double, prscale, n);
+  OCTAVE_LOCAL_BUFFER (double, pwork, 6 * n);
+
+  balanced_mat = a;
+  double *p_balanced_mat = balanced_mat.fortran_vec ();
+  balanced_mat2 = b;
+  double *p_balanced_mat2 = balanced_mat2.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (dggbal, DGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, p_balanced_mat2,
+			     n, ilo, ihi, plscale, prscale, pwork, info
+			     F77_CHAR_ARG_LEN  (1)));
+
+  balancing_mat = Matrix (n, n, 0.0);
+  balancing_mat2 = Matrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+      balancing_mat.elem (i ,i) = 1.0;
+      balancing_mat2.elem (i ,i) = 1.0;
+    }
+
+  double *p_balancing_mat = balancing_mat.fortran_vec ();
+  double *p_balancing_mat2 = balancing_mat2.fortran_vec ();
+
+  // first left
+  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+      
+  // then right
+  F77_XFCN (dggbak, DGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("R", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat2, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/dbleGEPBAL.h b/liboctave/dbleGEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/dbleGEPBAL.h
@@ -0,0 +1,90 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_GEPBALANCE_h)
+#define octave_GEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "dMatrix.h"
+
+class
+OCTAVE_API
+GEPBALANCE
+{
+public:
+
+  GEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  GEPBALANCE (const Matrix& a, const Matrix& b, const std::string& balance_job)
+    {
+      init (a, b, balance_job); 
+    }
+
+  GEPBALANCE (const GEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2),
+    balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { }
+
+  GEPBALANCE& operator = (const GEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balanced_mat2 = a.balanced_mat2;
+	  balancing_mat = a.balancing_mat;
+	  balancing_mat2 = a.balancing_mat2;
+	}
+      return *this;
+    }
+
+  ~GEPBALANCE (void) { }
+
+  Matrix balanced_matrix (void) const { return balanced_mat; }
+
+  Matrix balanced_matrix2 (void) const { return balanced_mat2; }
+
+  Matrix balancing_matrix (void) const { return balancing_mat; }
+
+  Matrix balancing_matrix2 (void) const { return balancing_mat2; }
+
+  friend std::ostream& operator << (std::ostream& os, const GEPBALANCE& a);
+
+private:
+
+  Matrix balanced_mat;
+  Matrix balanced_mat2;
+  Matrix balancing_mat;
+  Matrix balancing_mat2;
+
+  octave_idx_type init (const Matrix& a, const Matrix& b, 
+			const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxAEPBAL.cc b/liboctave/fCmplxAEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxAEPBAL.cc
@@ -0,0 +1,102 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+
+#include "fCmplxAEPBAL.h"
+#include "fMatrix.h"
+#include "f77-fcn.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 (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
+FloatComplexAEPBALANCE::init (const FloatComplexMatrix& a,
+			 const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  Array<float> scale (n);
+  float *pscale = scale.fortran_vec ();
+
+  balanced_mat = a;
+  FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (cgebal, CGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, ilo, ihi,
+			     pscale, info
+			     F77_CHAR_ARG_LEN (1)));
+
+  balancing_mat = FloatComplexMatrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    balancing_mat.elem (i, i) = 1.0;
+
+  FloatComplex *p_balancing_mat = balancing_mat.fortran_vec ();
+
+  char side = 'R';
+
+  F77_XFCN (cgebak, CGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 (&side, 1),
+			     n, ilo, ihi, pscale, n,
+			     p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxAEPBAL.h b/liboctave/fCmplxAEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxAEPBAL.h
@@ -0,0 +1,80 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatComplexAEPBALANCE_h)
+#define octave_FloatComplexAEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "fCMatrix.h"
+
+class
+OCTAVE_API
+FloatComplexAEPBALANCE
+{
+public:
+
+  FloatComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  FloatComplexAEPBALANCE (const FloatComplexMatrix& a, const std::string& balance_job)
+    {
+      init (a, balance_job); 
+    }
+
+  FloatComplexAEPBALANCE (const FloatComplexAEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { }
+
+  FloatComplexAEPBALANCE& operator = (const FloatComplexAEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balancing_mat = a.balancing_mat;
+	}
+      return *this;
+    }
+
+  ~FloatComplexAEPBALANCE (void) { }
+
+  FloatComplexMatrix balanced_matrix (void) const { return balanced_mat; }
+
+  FloatComplexMatrix balancing_matrix (void) const { return balancing_mat; }
+
+  friend std::ostream& operator << (std::ostream& os, const FloatComplexAEPBALANCE& a);
+
+private:
+
+  FloatComplexMatrix balanced_mat;
+  FloatComplexMatrix balancing_mat;
+
+  octave_idx_type init (const FloatComplexMatrix& a, const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxGEPBAL.cc b/liboctave/fCmplxGEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxGEPBAL.cc
@@ -0,0 +1,130 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+#include <vector>
+
+#include "fCmplxGEPBAL.h"
+#include "Array-util.h"
+#include "f77-fcn.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (cggbal, CGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N,
+			     FloatComplex* A, const octave_idx_type& LDA, FloatComplex* B,
+			     const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI,
+			     float* LSCALE, float* RSCALE,
+			     float* WORK, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sggbak, SGGBAK) (F77_CONST_CHAR_ARG_DECL,
+			     F77_CONST_CHAR_ARG_DECL,
+			     const octave_idx_type& N, const octave_idx_type& ILO,
+			     const octave_idx_type& IHI, const float* LSCALE,
+			     const float* RSCALE, octave_idx_type& M, float* V,
+			     const octave_idx_type& LDV, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+
+}
+
+octave_idx_type
+FloatComplexGEPBALANCE::init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 
+		  const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("FloatComplexGEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  if (a.dims() != b.dims ())
+    {
+      gripe_nonconformant ("FloatComplexGEPBALANCE", n, n, b.rows(), b.cols());
+      return -1;
+    } 
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  OCTAVE_LOCAL_BUFFER (float, plscale, n);
+  OCTAVE_LOCAL_BUFFER (float, prscale, n);
+  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
+
+  balanced_mat = a;
+  FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
+  balanced_mat2 = b;
+  FloatComplex *p_balanced_mat2 = balanced_mat2.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (cggbal, CGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, p_balanced_mat2,
+			     n, ilo, ihi, plscale,prscale, pwork, info
+			     F77_CHAR_ARG_LEN (1)));
+
+  balancing_mat = FloatMatrix (n, n, 0.0);
+  balancing_mat2 = FloatMatrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+      balancing_mat.elem (i ,i) = 1.0;
+      balancing_mat2.elem (i ,i) = 1.0;
+    }
+
+  float *p_balancing_mat = balancing_mat.fortran_vec ();
+  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
+
+  // first left
+  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+      
+  // then right
+  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("R", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat2, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxGEPBAL.h b/liboctave/fCmplxGEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxGEPBAL.h
@@ -0,0 +1,91 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatComplexGEPBALANCE_h)
+#define octave_FloatComplexGEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "fMatrix.h"
+#include "fCMatrix.h"
+
+class
+OCTAVE_API
+FloatComplexGEPBALANCE
+{
+public:
+
+  FloatComplexGEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  FloatComplexGEPBALANCE (const FloatComplexMatrix& a, const FloatComplexMatrix& b, const std::string& balance_job)
+    {
+      init (a, b, balance_job); 
+    }
+
+  FloatComplexGEPBALANCE (const FloatComplexGEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2),
+    balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { }
+
+  FloatComplexGEPBALANCE& operator = (const FloatComplexGEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balanced_mat2 = a.balanced_mat2;
+	  balancing_mat = a.balancing_mat;
+	  balancing_mat2 = a.balancing_mat2;
+	}
+      return *this;
+    }
+
+  ~FloatComplexGEPBALANCE (void) { }
+
+  FloatComplexMatrix balanced_matrix (void) const { return balanced_mat; }
+
+  FloatComplexMatrix balanced_matrix2 (void) const { return balanced_mat2; }
+
+  FloatMatrix balancing_matrix (void) const { return balancing_mat; }
+
+  FloatMatrix balancing_matrix2 (void) const { return balancing_mat2; }
+
+  friend std::ostream& operator << (std::ostream& os, const FloatComplexGEPBALANCE& a);
+
+private:
+
+  FloatComplexMatrix balanced_mat;
+  FloatComplexMatrix balanced_mat2;
+  FloatMatrix balancing_mat;
+  FloatMatrix balancing_mat2;
+
+  octave_idx_type init (const FloatComplexMatrix& a, const FloatComplexMatrix& b, 
+			const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxHESS.cc b/liboctave/fCmplxHESS.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxHESS.cc
@@ -0,0 +1,127 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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 (n);
+  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 (n-1);
+  FloatComplex *ptau = tau.fortran_vec ();
+
+  Array<FloatComplex> work (lwork);
+  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;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxHESS.h b/liboctave/fCmplxHESS.h
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxHESS.h
@@ -0,0 +1,81 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatComplexHESS_h)
+#define octave_FloatComplexHESS_h 1
+
+#include <iostream>
+
+#include "fCMatrix.h"
+
+class
+OCTAVE_API
+FloatComplexHESS
+{
+public:
+
+  FloatComplexHESS (void) : hess_mat (), unitary_hess_mat () { }
+
+  FloatComplexHESS (const FloatComplexMatrix& a) { init (a); }
+
+  FloatComplexHESS (const FloatComplexMatrix& a, octave_idx_type& info) { 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
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxQR.cc b/liboctave/fCmplxQR.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxQR.cc
@@ -0,0 +1,307 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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/>.
+
+*/
+
+// updating/downdating by Jaroslav Hajek 2008
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "fCmplxQR.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+#include "Range.h"
+#include "idx-vector.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (cgeqrf, CGEQRF) (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 (cungqr, CUNGQR) (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&);
+
+  // these come from qrupdate
+
+  F77_RET_T
+  F77_FUNC (cqr1up, CQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             FloatComplex*, FloatComplex*, const FloatComplex*, const FloatComplex*);
+
+  F77_RET_T
+  F77_FUNC (cqrinc, CQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&, const FloatComplex*);
+
+  F77_RET_T
+  F77_FUNC (cqrdec, CQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             FloatComplex*, const FloatComplex*, FloatComplex*, const octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cqrinr, CQRINR) (const octave_idx_type&, const octave_idx_type&, 
+                             const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex*, 
+                             const octave_idx_type&, const FloatComplex*);
+
+  F77_RET_T
+  F77_FUNC (cqrder, CQRDER) (const octave_idx_type&, const octave_idx_type&, 
+                             const FloatComplex*, FloatComplex*, const FloatComplex*, FloatComplex *, 
+                             const octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cqrshc, CQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
+                             FloatComplex*, FloatComplex*, const octave_idx_type&, const octave_idx_type&);
+}
+
+FloatComplexQR::FloatComplexQR (const FloatComplexMatrix& a, QR::type qr_type)
+  : q (), r ()
+{
+  init (a, qr_type);
+}
+
+void
+FloatComplexQR::init (const FloatComplexMatrix& a, QR::type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  if (m == 0 || n == 0)
+    {
+      (*current_liboctave_error_handler)
+	("FloatComplexQR must have non-empty matrix");
+      return;
+    }
+
+  octave_idx_type min_mn = m < n ? m : n;
+
+  Array<FloatComplex> tau (min_mn);
+  FloatComplex *ptau = tau.fortran_vec ();
+
+  octave_idx_type lwork = 32*n;
+  Array<FloatComplex> work (lwork);
+  FloatComplex *pwork = work.fortran_vec ();
+
+  octave_idx_type info = 0;
+
+  FloatComplexMatrix A_fact;
+  if (m > n && qr_type != QR::economy)
+    {
+      A_fact.resize (m, m);
+      A_fact.insert (a, 0, 0);
+    }
+  else
+    A_fact = a;
+
+  FloatComplex *tmp_data = A_fact.fortran_vec ();
+
+  F77_XFCN (cgeqrf, CGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info));
+
+  if (qr_type == QR::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+	{
+	  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+	  for (octave_idx_type i = limit + 1; i < m; i++)
+	    A_fact.elem (i, j) *= tau.elem (j);
+	}
+
+      r = A_fact;
+
+      if (m > n)
+	r.resize (m, n);
+    }
+  else
+    {
+      octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
+
+      if (qr_type == QR::economy && m > n)
+	r.resize (n, n, 0.0);
+      else
+	r.resize (m, n, 0.0);
+
+      for (octave_idx_type j = 0; j < n; j++)
+	{
+	  octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
+	  for (octave_idx_type i = 0; i <= limit; i++)
+	    r.elem (i, j) = A_fact.elem (i, j);
+	}
+
+      lwork = 32 * n2;
+      work.resize (lwork);
+      FloatComplex *pwork2 = work.fortran_vec ();
+
+      F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
+				 pwork2, lwork, info));
+
+      q = A_fact;
+      q.resize (m, n2);
+    }
+}
+
+FloatComplexQR::FloatComplexQR (const FloatComplexMatrix &q, const FloatComplexMatrix& r)
+{
+  if (q.columns () != r.rows ()) 
+    {
+      (*current_liboctave_error_handler) ("QR dimensions mismatch");
+      return;
+    }
+
+  this->q = q;
+  this->r = r;
+}
+
+void
+FloatComplexQR::update (const FloatComplexMatrix& u, const FloatComplexMatrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.length () == m && v.length () == n)
+    F77_XFCN (cqr1up, CQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), 
+			       u.data (), v.data ()));
+  else
+    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+}
+
+void
+FloatComplexQR::insert_col (const FloatComplexMatrix& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("QR insert index out of range");
+  else
+    {
+      FloatComplexMatrix r1 (m,n+1);
+
+      F77_XFCN (cqrinc, CQRINC, (m, n, k, q.fortran_vec (), r.data (),
+				 r1.fortran_vec (), j+1, u.data ()));
+
+      r = r1;
+    }
+}
+
+void
+FloatComplexQR::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (k < m && k < n) 
+    (*current_liboctave_error_handler) ("QR delete dimensions mismatch");
+  else if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("QR delete index out of range");
+  else
+    {
+      FloatComplexMatrix r1 (k, n-1);
+
+      F77_XFCN (cqrdec, CQRDEC, (m, n, k, q.fortran_vec (), r.data (),
+				 r1.fortran_vec (), j+1));
+
+      r = r1;
+    }
+}
+
+void
+FloatComplexQR::insert_row (const FloatComplexMatrix& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("QR insert index out of range");
+  else
+    {
+      FloatComplexMatrix q1 (m+1, m+1);
+      FloatComplexMatrix r1 (m+1, n);
+
+      F77_XFCN (cqrinr, CQRINR, (m, n, q.data (), q1.fortran_vec (), 
+				 r.data (), r1.fortran_vec (), j+1, u.data ()));
+
+      q = q1;
+      r = r1;
+    }
+}
+
+void
+FloatComplexQR::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("QR delete dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("QR delete index out of range");
+  else
+    {
+      FloatComplexMatrix q1 (m-1, m-1);
+      FloatComplexMatrix r1 (m-1, n);
+
+      F77_XFCN (cqrder, CQRDER, (m, n, q.data (), q1.fortran_vec (), 
+				 r.data (), r1.fortran_vec (), j+1 ));
+
+      q = q1;
+      r = r1;
+    }
+}
+
+void
+FloatComplexQR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("QR shift index out of range");
+  else
+    F77_XFCN (cqrshc, CQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1));
+}
+
+void
+FloatComplexQR::economize (void)
+{
+  octave_idx_type r_nc = r.columns ();
+
+  if (r.rows () > r_nc)
+    {
+      q.resize (q.rows (), r_nc);
+      r.resize (r_nc, r_nc);
+    }
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxQR.h b/liboctave/fCmplxQR.h
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxQR.h
@@ -0,0 +1,96 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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/>.
+
+*/
+
+// updating/downdating by Jaroslav Hajek 2008
+
+#if !defined (octave_FloatComplexQR_h)
+#define octave_FloatComplexQR_h 1
+
+#include <iostream>
+
+#include "fCMatrix.h"
+#include "fCColVector.h"
+#include "fCRowVector.h"
+#include "dbleQR.h"
+
+class
+OCTAVE_API
+FloatComplexQR
+{
+public:
+
+  FloatComplexQR (void) : q (), r () { }
+
+  FloatComplexQR (const FloatComplexMatrix&, QR::type = QR::std);
+
+  FloatComplexQR (const FloatComplexMatrix& q, const FloatComplexMatrix& r);
+
+  FloatComplexQR (const FloatComplexQR& a) : q (a.q), r (a.r) { }
+
+  FloatComplexQR& operator = (const FloatComplexQR& a)
+    {
+      if (this != &a)
+	{
+	  q = a.q;
+	  r = a.r;
+	}
+      return *this;
+    }
+
+  ~FloatComplexQR (void) { }
+
+  void init (const FloatComplexMatrix&, QR::type = QR::std);
+
+  FloatComplexMatrix Q (void) const { return q; }
+
+  FloatComplexMatrix R (void) const { return r; }
+
+  void update (const FloatComplexMatrix& u, const FloatComplexMatrix& v);
+
+  void insert_col (const FloatComplexMatrix& u, octave_idx_type j);
+
+  void delete_col (octave_idx_type j);
+
+  void insert_row (const FloatComplexMatrix& u, octave_idx_type j);
+
+  void delete_row (octave_idx_type j);
+
+  void shift_cols (octave_idx_type i, octave_idx_type j);
+
+  void economize();
+
+  friend std::ostream&  operator << (std::ostream&, const FloatComplexQR&);
+
+protected:
+
+  FloatComplexMatrix q;
+  FloatComplexMatrix r;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxQRP.cc b/liboctave/fCmplxQRP.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxQRP.cc
@@ -0,0 +1,138 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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 <cassert>
+
+#include "fCmplxQRP.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (cgeqpf, CGEQPF) (const octave_idx_type&, const octave_idx_type&, FloatComplex*,
+			     const octave_idx_type&, octave_idx_type*, FloatComplex*, FloatComplex*,
+			     float*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (cungqr, CUNGQR) (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&);
+}
+
+// It would be best to share some of this code with FloatComplexQR class...
+
+FloatComplexQRP::FloatComplexQRP (const FloatComplexMatrix& a, QR::type qr_type)
+  : FloatComplexQR (), p ()
+{
+  init (a, qr_type);
+}
+
+void
+FloatComplexQRP::init (const FloatComplexMatrix& a, QR::type qr_type)
+{
+  assert (qr_type != QR::raw);
+
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  if (m == 0 || n == 0)
+    {
+      (*current_liboctave_error_handler)
+	("FloatComplexQR must have non-empty matrix");
+      return;
+    }
+
+  octave_idx_type min_mn = m < n ? m : n;
+  Array<FloatComplex> tau (min_mn);
+  FloatComplex *ptau = tau.fortran_vec ();
+
+  octave_idx_type lwork = 3*n > 32*m ? 3*n : 32*m;
+  Array<FloatComplex> work (lwork);
+  FloatComplex *pwork = work.fortran_vec ();
+
+  octave_idx_type info = 0;
+
+  FloatComplexMatrix A_fact = a;
+  if (m > n && qr_type != QR::economy)
+    A_fact.resize (m, m, 0.0);
+
+  FloatComplex *tmp_data = A_fact.fortran_vec ();
+
+  Array<float> rwork (2*n);
+  float *prwork = rwork.fortran_vec ();
+
+  Array<octave_idx_type> jpvt (n, 0);
+  octave_idx_type *pjpvt = jpvt.fortran_vec ();
+
+  // Code to enforce a certain permutation could go here...
+
+  F77_XFCN (cgeqpf, CGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork,
+			     prwork, info));
+
+  // Form Permutation matrix (if economy is requested, return the
+  // indices only!)
+
+  if (qr_type == QR::economy)
+    {
+      p.resize (1, n, 0.0);
+      for (octave_idx_type j = 0; j < n; j++)
+	p.elem (0, j) = jpvt.elem (j);
+    }
+  else
+    {
+      p.resize (n, n, 0.0);
+      for (octave_idx_type j = 0; j < n; j++)
+	p.elem (jpvt.elem (j) - 1, j) = 1.0;
+    }
+
+  octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
+
+  if (qr_type == QR::economy && m > n)
+    r.resize (n, n, 0.0);
+  else
+    r.resize (m, n, 0.0);
+
+  for (octave_idx_type j = 0; j < n; j++)
+    {
+      octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
+      for (octave_idx_type i = 0; i <= limit; i++)
+	r.elem (i, j) = A_fact.elem (i, j);
+    }
+
+  F77_XFCN (cungqr, CUNGQR, (m, n2, min_mn, tmp_data, m, ptau,
+			     pwork, lwork, info));
+
+  q = A_fact;
+  q.resize (m, n2);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/fCmplxQRP.h b/liboctave/fCmplxQRP.h
new file mode 100644
--- /dev/null
+++ b/liboctave/fCmplxQRP.h
@@ -0,0 +1,73 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatComplexQRP_h)
+#define octave_FloatComplexQRP_h 1
+
+#include <iostream>
+
+#include "fCmplxQR.h"
+#include "fMatrix.h"
+
+class
+OCTAVE_API
+FloatComplexQRP : public FloatComplexQR
+{
+public:
+
+  FloatComplexQRP (void) : FloatComplexQR (), p () { }
+
+  FloatComplexQRP (const FloatComplexMatrix&, QR::type = QR::std);
+
+  FloatComplexQRP (const FloatComplexQRP& a) : FloatComplexQR (a), p (a.p) { }
+
+  FloatComplexQRP& operator = (const FloatComplexQRP& a)
+    {
+      if (this != &a)
+	{
+	  FloatComplexQR::operator = (a);
+	  p = a.p;
+	}
+      return *this;
+    }
+
+  ~FloatComplexQRP (void) { }
+
+  void init (const FloatComplexMatrix&, QR::type = QR::std);
+
+  FloatMatrix P (void) const { return p; }
+
+  friend std::ostream&  operator << (std::ostream&, const FloatComplexQRP&);
+
+private:
+
+  FloatMatrix p;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatAEPBAL.cc b/liboctave/floatAEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/floatAEPBAL.cc
@@ -0,0 +1,99 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+
+#include "floatAEPBAL.h"
+#include "f77-fcn.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 (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
+FloatAEPBALANCE::init (const FloatMatrix& a, const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  Array<float> scale (n);
+  float *pscale = scale.fortran_vec ();
+
+  balanced_mat = a;
+  float *p_balanced_mat = balanced_mat.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (sgebal, SGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, ilo, ihi, pscale, info
+			     F77_CHAR_ARG_LEN (1)));
+
+  balancing_mat = Matrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    balancing_mat.elem (i ,i) = 1.0;
+
+  float *p_balancing_mat = balancing_mat.fortran_vec ();
+
+  char side = 'R';
+
+  F77_XFCN (sgebak, SGEBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 (&side, 1),
+			     n, ilo, ihi, pscale, n,
+			     p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatAEPBAL.h b/liboctave/floatAEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/floatAEPBAL.h
@@ -0,0 +1,80 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatAEPBALANCE_h)
+#define octave_FloatAEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "fMatrix.h"
+
+class
+OCTAVE_API
+FloatAEPBALANCE
+{
+public:
+
+  FloatAEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  FloatAEPBALANCE (const FloatMatrix& a,const std::string& balance_job)
+    {
+      init (a, balance_job); 
+    }
+
+  FloatAEPBALANCE (const FloatAEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { }
+
+  FloatAEPBALANCE& operator = (const FloatAEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balancing_mat = a.balancing_mat;
+	}
+      return *this;
+    }
+
+  ~FloatAEPBALANCE (void) { }
+
+  FloatMatrix balanced_matrix (void) const { return balanced_mat; }
+
+  FloatMatrix balancing_matrix (void) const { return balancing_mat; }
+
+  friend std::ostream& operator << (std::ostream& os, const FloatAEPBALANCE& a);
+
+private:
+
+  FloatMatrix balanced_mat;
+  FloatMatrix balancing_mat;
+
+  octave_idx_type init (const FloatMatrix& a, const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatGEPBAL.cc b/liboctave/floatGEPBAL.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/floatGEPBAL.cc
@@ -0,0 +1,130 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
+              2007 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 <string>
+#include <vector>
+
+#include "floatGEPBAL.h"
+#include "Array-util.h"
+#include "f77-fcn.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (sggbal, SGGBAL) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type& N,
+			     float* A, const octave_idx_type& LDA, float* B,
+			     const octave_idx_type& LDB, octave_idx_type& ILO, octave_idx_type& IHI,
+			     float* LSCALE, float* RSCALE,
+			     float* WORK, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL);
+
+  F77_RET_T
+  F77_FUNC (sggbak, SGGBAK) (F77_CONST_CHAR_ARG_DECL,
+			     F77_CONST_CHAR_ARG_DECL,
+			     const octave_idx_type& N, const octave_idx_type& ILO,
+			     const octave_idx_type& IHI, const float* LSCALE,
+			     const float* RSCALE, octave_idx_type& M, float* V,
+			     const octave_idx_type& LDV, octave_idx_type& INFO
+			     F77_CHAR_ARG_LEN_DECL
+			     F77_CHAR_ARG_LEN_DECL);
+
+}
+
+octave_idx_type
+FloatGEPBALANCE::init (const FloatMatrix& a, const FloatMatrix& b, 
+		  const std::string& balance_job)
+{
+  octave_idx_type n = a.cols ();
+
+  if (a.rows () != n)
+    {
+      (*current_liboctave_error_handler) ("FloatGEPBALANCE requires square matrix");
+      return -1;
+    }
+
+  if (a.dims() != b.dims ())
+    {
+      gripe_nonconformant ("FloatGEPBALANCE", n, n, b.rows(), b.cols());
+      return -1;
+    } 
+
+  octave_idx_type info;
+  octave_idx_type ilo;
+  octave_idx_type ihi;
+
+  OCTAVE_LOCAL_BUFFER (float, plscale, n);
+  OCTAVE_LOCAL_BUFFER (float, prscale, n);
+  OCTAVE_LOCAL_BUFFER (float, pwork, 6 * n);
+
+  balanced_mat = a;
+  float *p_balanced_mat = balanced_mat.fortran_vec ();
+  balanced_mat2 = b;
+  float *p_balanced_mat2 = balanced_mat2.fortran_vec ();
+
+  char job = balance_job[0];
+
+  F77_XFCN (sggbal, SGGBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     n, p_balanced_mat, n, p_balanced_mat2,
+			     n, ilo, ihi, plscale, prscale, pwork, info
+			     F77_CHAR_ARG_LEN  (1)));
+
+  balancing_mat = FloatMatrix (n, n, 0.0);
+  balancing_mat2 = FloatMatrix (n, n, 0.0);
+  for (octave_idx_type i = 0; i < n; i++)
+    {
+      OCTAVE_QUIT;
+      balancing_mat.elem (i ,i) = 1.0;
+      balancing_mat2.elem (i ,i) = 1.0;
+    }
+
+  float *p_balancing_mat = balancing_mat.fortran_vec ();
+  float *p_balancing_mat2 = balancing_mat2.fortran_vec ();
+
+  // first left
+  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("L", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+      
+  // then right
+  F77_XFCN (sggbak, SGGBAK, (F77_CONST_CHAR_ARG2 (&job, 1),
+			     F77_CONST_CHAR_ARG2 ("R", 1),
+			     n, ilo, ihi, plscale, prscale,
+			     n, p_balancing_mat2, n, info
+			     F77_CHAR_ARG_LEN (1)
+			     F77_CHAR_ARG_LEN (1)));
+
+  return info;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatGEPBAL.h b/liboctave/floatGEPBAL.h
new file mode 100644
--- /dev/null
+++ b/liboctave/floatGEPBAL.h
@@ -0,0 +1,90 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatGEPBALANCE_h)
+#define octave_FloatGEPBALANCE_h 1
+
+#include <iostream>
+#include <string>
+
+#include "fMatrix.h"
+
+class
+OCTAVE_API
+FloatGEPBALANCE
+{
+public:
+
+  FloatGEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+
+  FloatGEPBALANCE (const FloatMatrix& a, const FloatMatrix& b, const std::string& balance_job)
+    {
+      init (a, b, balance_job); 
+    }
+
+  FloatGEPBALANCE (const FloatGEPBALANCE& a)
+    : balanced_mat (a.balanced_mat), balanced_mat2 (a.balanced_mat2),
+    balancing_mat (a.balancing_mat), balancing_mat2 (a.balancing_mat2) { }
+
+  FloatGEPBALANCE& operator = (const FloatGEPBALANCE& a)
+    {
+      if (this != &a)
+	{
+	  balanced_mat = a.balanced_mat;
+	  balanced_mat2 = a.balanced_mat2;
+	  balancing_mat = a.balancing_mat;
+	  balancing_mat2 = a.balancing_mat2;
+	}
+      return *this;
+    }
+
+  ~FloatGEPBALANCE (void) { }
+
+  FloatMatrix balanced_matrix (void) const { return balanced_mat; }
+
+  FloatMatrix balanced_matrix2 (void) const { return balanced_mat2; }
+
+  FloatMatrix balancing_matrix (void) const { return balancing_mat; }
+
+  FloatMatrix balancing_matrix2 (void) const { return balancing_mat2; }
+
+  friend std::ostream& operator << (std::ostream& os, const FloatGEPBALANCE& a);
+
+private:
+
+  FloatMatrix balanced_mat;
+  FloatMatrix balanced_mat2;
+  FloatMatrix balancing_mat;
+  FloatMatrix balancing_mat2;
+
+  octave_idx_type init (const FloatMatrix& a, const FloatMatrix& b, 
+			const std::string& balance_job);
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatHESS.cc b/liboctave/floatHESS.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/floatHESS.cc
@@ -0,0 +1,128 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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");
+      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;
+  float *h = hess_mat.fortran_vec ();
+
+  Array<float> scale (n);
+  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 (n-1);
+  float *ptau = tau.fortran_vec ();
+
+  Array<float> work (lwork);
+  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;
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatHESS.h b/liboctave/floatHESS.h
new file mode 100644
--- /dev/null
+++ b/liboctave/floatHESS.h
@@ -0,0 +1,78 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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 <iostream>
+
+#include "fMatrix.h"
+
+class
+OCTAVE_API
+FloatHESS
+{
+public:
+
+  FloatHESS (void) : hess_mat (), unitary_hess_mat () { }
+
+  FloatHESS (const FloatMatrix& a) { init (a); }
+
+  FloatHESS (const FloatMatrix& a, octave_idx_type& info) { 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
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatQR.cc b/liboctave/floatQR.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/floatQR.cc
@@ -0,0 +1,297 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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 "floatQR.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+#include "Range.h"
+#include "idx-vector.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (sgeqrf, SGEQRF) (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 (sorgqr, SORGQR) (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&);
+
+  // these come from qrupdate
+
+  F77_RET_T
+  F77_FUNC (sqr1up, SQR1UP) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             float*, float*, const float*, const float*);
+
+  F77_RET_T
+  F77_FUNC (sqrinc, SQRINC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             float*, const float*, float*, const octave_idx_type&, const float*);
+
+  F77_RET_T
+  F77_FUNC (sqrdec, SQRDEC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             float*, const float*, float*, const octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sqrinr, SQRINR) (const octave_idx_type&, const octave_idx_type&, 
+                             const float*, float*, const float*, float*, 
+                             const octave_idx_type&, const float*);
+
+  F77_RET_T
+  F77_FUNC (sqrder, SQRDER) (const octave_idx_type&, const octave_idx_type&, 
+                             const float*, float*, const float*, float *, 
+                             const octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sqrshc, SQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
+                             float*, float*, const octave_idx_type&, const octave_idx_type&);
+}
+
+FloatQR::FloatQR (const FloatMatrix& a, QR::type qr_type)
+  : q (), r ()
+{
+  init (a, qr_type);
+}
+
+void
+FloatQR::init (const FloatMatrix& a, QR::type qr_type)
+{
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  if (m == 0 || n == 0)
+    {
+      (*current_liboctave_error_handler) ("QR must have non-empty matrix");
+      return;
+    }
+
+  octave_idx_type min_mn = m < n ? m : n;
+  Array<float> tau (min_mn);
+  float *ptau = tau.fortran_vec ();
+
+  octave_idx_type lwork = 32*n;
+  Array<float> work (lwork);
+  float *pwork = work.fortran_vec ();
+
+  octave_idx_type info = 0;
+
+  FloatMatrix A_fact = a;
+  if (m > n && qr_type != QR::economy)
+      A_fact.resize (m, m, 0.0);
+
+  float *tmp_data = A_fact.fortran_vec ();
+
+  F77_XFCN (sgeqrf, SGEQRF, (m, n, tmp_data, m, ptau, pwork, lwork, info));
+
+  if (qr_type == QR::raw)
+    {
+      for (octave_idx_type j = 0; j < min_mn; j++)
+	{
+	  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
+	  for (octave_idx_type i = limit + 1; i < m; i++)
+	    A_fact.elem (i, j) *= tau.elem (j);
+	}
+
+      r = A_fact;
+
+      if (m > n)
+	r.resize (m, n);
+    }
+  else
+    {
+      octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
+
+      if (qr_type == QR::economy && m > n)
+	r.resize (n, n, 0.0);
+      else
+	r.resize (m, n, 0.0);
+
+      for (octave_idx_type j = 0; j < n; j++)
+	{
+	  octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
+	  for (octave_idx_type i = 0; i <= limit; i++)
+	    r.elem (i, j) = tmp_data[m*j+i];
+	}
+
+      lwork = 32 * n2;
+      work.resize (lwork);
+      float *pwork2 = work.fortran_vec ();
+
+      F77_XFCN (sorgqr, SORGQR, (m, n2, min_mn, tmp_data, m, ptau,
+				 pwork2, lwork, info));
+
+      q = A_fact;
+      q.resize (m, n2);
+    }
+}
+
+FloatQR::FloatQR (const FloatMatrix& q, const FloatMatrix& r)
+{
+  if (q.columns () != r.rows ()) 
+    {
+      (*current_liboctave_error_handler) ("QR dimensions mismatch");
+      return;
+    }
+
+  this->q = q;
+  this->r = r;
+}
+
+void
+FloatQR::update (const FloatMatrix& u, const FloatMatrix& v)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.length () == m && v.length () == n)
+    F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), r.fortran_vec (), 
+			       u.data (), v.data ()));
+  else
+    (*current_liboctave_error_handler) ("QR update dimensions mismatch");
+}
+
+void
+FloatQR::insert_col (const FloatMatrix& u, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type n = r.columns ();
+  octave_idx_type k = q.columns ();
+
+  if (u.length () != m)
+    (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
+  else if (j < 0 || j > n) 
+    (*current_liboctave_error_handler) ("QR insert index out of range");
+  else
+    {
+      FloatMatrix r1 (m, n+1);
+
+      F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), r.data (),
+				 r1.fortran_vec (), j+1, u.data ()));
+
+      r = r1;
+    }
+}
+
+void
+FloatQR::delete_col (octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (k < m && k < n) 
+    (*current_liboctave_error_handler) ("QR delete dimensions mismatch");
+  else if (j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("QR delete index out of range");
+  else
+    {
+      FloatMatrix r1 (k, n-1);
+
+      F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), r.data (),
+				 r1.fortran_vec (), j+1));
+
+      r = r1;
+    }
+}
+
+void
+FloatQR::insert_row (const FloatMatrix& u, octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square () || u.length () != n)
+    (*current_liboctave_error_handler) ("QR insert dimensions mismatch");
+  else if (j < 0 || j > m) 
+    (*current_liboctave_error_handler) ("QR insert index out of range");
+  else
+    {
+      FloatMatrix q1 (m+1, m+1);
+      FloatMatrix r1 (m+1, n);
+
+      F77_XFCN (sqrinr, SQRINR, (m, n, q.data (), q1.fortran_vec (), 
+				 r.data (), r1.fortran_vec (), j+1, u.data ()));
+
+      q = q1;
+      r = r1;
+    }
+}
+
+void
+FloatQR::delete_row (octave_idx_type j)
+{
+  octave_idx_type m = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (! q.is_square ())
+    (*current_liboctave_error_handler) ("QR delete dimensions mismatch");
+  else if (j < 0 || j > m-1) 
+    (*current_liboctave_error_handler) ("QR delete index out of range");
+  else
+    {
+      FloatMatrix q1 (m-1, m-1);
+      FloatMatrix r1 (m-1, n);
+
+      F77_XFCN (sqrder, SQRDER, (m, n, q.data (), q1.fortran_vec (), 
+				 r.data (), r1.fortran_vec (), j+1 ));
+
+      q = q1;
+      r = r1;
+    }
+}
+
+void
+FloatQR::shift_cols (octave_idx_type i, octave_idx_type j)
+{
+  octave_idx_type m = q.rows ();
+  octave_idx_type k = r.rows ();
+  octave_idx_type n = r.columns ();
+
+  if (i < 0 || i > n-1 || j < 0 || j > n-1) 
+    (*current_liboctave_error_handler) ("QR shift index out of range");
+  else
+    F77_XFCN (sqrshc, SQRSHC, (m, n, k, q.fortran_vec (), r.fortran_vec (), i+1, j+1));
+}
+
+void
+FloatQR::economize (void)
+{
+  octave_idx_type r_nc = r.columns ();
+
+  if (r.rows () > r_nc)
+    {
+      q.resize (q.rows (), r_nc);
+      r.resize (r_nc, r_nc);
+    }
+}
+
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatQR.h b/liboctave/floatQR.h
new file mode 100644
--- /dev/null
+++ b/liboctave/floatQR.h
@@ -0,0 +1,96 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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/>.
+
+*/
+
+// updating/downdating by Jaroslav Hajek 2008
+
+#if !defined (octave_FloatQR_h)
+#define octave_FloatQR_h 1
+
+#include <iostream>
+
+#include "fMatrix.h"
+#include "fColVector.h"
+#include "fRowVector.h"
+#include "dbleQR.h"
+
+class
+OCTAVE_API
+FloatQR
+{
+public:
+
+  FloatQR (void) : q (), r () { }
+
+  FloatQR (const FloatMatrix&, QR::type = QR::std);
+
+  FloatQR (const FloatMatrix& q, const FloatMatrix& r);
+
+  FloatQR (const FloatQR& a) : q (a.q), r (a.r) { }
+
+  FloatQR& operator = (const FloatQR& a)
+    {
+      if (this != &a)
+	{
+	  q = a.q;
+	  r = a.r;
+	}
+      return *this;
+    }
+
+  ~FloatQR (void) { }
+
+  void init (const FloatMatrix&, QR::type);
+
+  FloatMatrix Q (void) const { return q; }
+
+  FloatMatrix R (void) const { return r; }
+
+  void update (const FloatMatrix& u, const FloatMatrix& v);
+
+  void insert_col (const FloatMatrix& u, octave_idx_type j);
+
+  void delete_col (octave_idx_type j);
+
+  void insert_row (const FloatMatrix& u, octave_idx_type j);
+
+  void delete_row (octave_idx_type j);
+
+  void shift_cols (octave_idx_type i, octave_idx_type j);
+
+  void economize (void);
+
+  friend std::ostream&  operator << (std::ostream&, const FloatQR&);
+
+protected:
+
+  FloatMatrix q;
+  FloatMatrix r;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatQRP.cc b/liboctave/floatQRP.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/floatQRP.cc
@@ -0,0 +1,132 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
+              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 <cassert>
+
+#include "floatQRP.h"
+#include "f77-fcn.h"
+#include "lo-error.h"
+
+extern "C"
+{
+  F77_RET_T
+  F77_FUNC (sgeqpf, SGEQPF) (const octave_idx_type&, const octave_idx_type&, float*,
+			     const octave_idx_type&, octave_idx_type*, float*, float*, octave_idx_type&);
+
+  F77_RET_T
+  F77_FUNC (sorgqr, SORGQR) (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&);
+}
+
+// It would be best to share some of this code with QR class...
+
+FloatQRP::FloatQRP (const FloatMatrix& a, QR::type qr_type)
+  : FloatQR (), p ()
+{
+  init (a, qr_type);
+}
+
+void
+FloatQRP::init (const FloatMatrix& a, QR::type qr_type)
+{
+  assert (qr_type != QR::raw);
+
+  octave_idx_type m = a.rows ();
+  octave_idx_type n = a.cols ();
+
+  if (m == 0 || n == 0)
+    {
+      (*current_liboctave_error_handler) ("QR must have non-empty matrix");
+      return;
+    }
+
+  octave_idx_type min_mn = m < n ? m : n;
+  Array<float> tau (min_mn);
+  float *ptau = tau.fortran_vec ();
+
+  octave_idx_type lwork = 3*n > 32*m ? 3*n : 32*m;
+  Array<float> work (lwork);
+  float *pwork = work.fortran_vec ();
+
+  octave_idx_type info = 0;
+
+  FloatMatrix A_fact = a;
+  if (m > n && qr_type != QR::economy)
+    A_fact.resize (m, m, 0.0);
+
+  float *tmp_data = A_fact.fortran_vec ();
+
+  Array<octave_idx_type> jpvt (n, 0);
+  octave_idx_type *pjpvt = jpvt.fortran_vec ();
+
+  // Code to enforce a certain permutation could go here...
+
+  F77_XFCN (sgeqpf, SGEQPF, (m, n, tmp_data, m, pjpvt, ptau, pwork, info));
+
+  // Form Permutation matrix (if economy is requested, return the
+  // indices only!)
+
+  if (qr_type == QR::economy)
+    {
+      p.resize (1, n, 0.0);
+      for (octave_idx_type j = 0; j < n; j++)
+	p.elem (0, j) = jpvt.elem (j);
+    }
+  else
+    {
+      p.resize (n, n, 0.0);
+      for (octave_idx_type j = 0; j < n; j++)
+	p.elem (jpvt.elem (j) - 1, j) = 1.0;
+    }
+
+  octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
+
+  if (qr_type == QR::economy && m > n)
+    r.resize (n, n, 0.0);
+  else
+    r.resize (m, n, 0.0);
+
+  for (octave_idx_type j = 0; j < n; j++)
+    {
+      octave_idx_type limit = j < min_mn-1 ? j : min_mn-1;
+      for (octave_idx_type i = 0; i <= limit; i++)
+	r.elem (i, j) = A_fact.elem (i, j);
+    }
+
+  F77_XFCN (sorgqr, SORGQR, (m, n2, min_mn, tmp_data, m, ptau,
+			     pwork, lwork, info));
+
+  q = A_fact;
+  q.resize (m, n2);
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/
diff --git a/liboctave/floatQRP.h b/liboctave/floatQRP.h
new file mode 100644
--- /dev/null
+++ b/liboctave/floatQRP.h
@@ -0,0 +1,74 @@
+/*
+
+Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
+              2007 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_FloatQRP_h)
+#define octave_FloatQRP_h 1
+
+#include <iostream>
+
+#include "floatQR.h"
+#include "fMatrix.h"
+
+class
+OCTAVE_API
+FloatQRP : public FloatQR
+{
+public:
+
+  FloatQRP (void) : FloatQR (), p () { }
+
+  FloatQRP (const FloatMatrix&, QR::type = QR::std);
+
+  FloatQRP (const FloatQRP& a) : FloatQR (a), p (a.p) { }
+
+  FloatQRP& operator = (const FloatQRP& a)
+    {
+      if (this != &a)
+	{
+	  FloatQR::operator = (a);
+	  p = a.p;
+	}
+
+      return *this;
+    }
+
+  ~FloatQRP (void) { }
+
+  void init (const FloatMatrix&, QR::type = QR::std);
+
+  FloatMatrix P (void) const { return p; }
+
+  friend std::ostream&  operator << (std::ostream&, const FloatQRP&);
+
+protected:
+
+  FloatMatrix p;
+};
+
+#endif
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/