# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1228813699 -3600
# Node ID a5e080076778ad789b1305be4384325f0d609aee
# Parent  6e9660cd3bf211360d937f1fe5bca5a9c2c58fe8
make balance more Matlab-compatible

diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/ChangeLog	Tue Dec 09 10:08:19 2008 +0100
@@ -1,3 +1,11 @@
+2008-12-09  Jaroslav Hajek  <highegg@gmail.com>
+
+	* base-aepbal.h: New source.
+	* dbleAEPBAL.h, dbleAEPBAL.cc: Rebase AEPBAL on base_aepbal.
+	* floatAEPBAL.h, floatAEPBAL.cc: Rebase FloatAEPBAL on base_aepbal.
+	* CmplxAEPBAL.h, CmplxAEPBAL.cc: Rebase ComplexAEPBAL on base_aepbal.
+	* fCmplxAEPBAL.h, fCmplxAEPBAL.cc: Rebase FloatComplexAEPBAL on base_aepbal.
+
 2008-12-08  Jaroslav Hajek  <highegg@gmail.com>
 
 	* idx-vector.cc (idx_vector::idx_vector_rep::idx_vector_rep (const 
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/CmplxAEPBAL.cc
--- a/liboctave/CmplxAEPBAL.cc	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/CmplxAEPBAL.cc	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -42,46 +43,52 @@
   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&, const 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
-ComplexAEPBALANCE::init (const ComplexMatrix& a,
-			 const std::string& balance_job)
+ComplexAEPBALANCE::ComplexAEPBALANCE (const ComplexMatrix& a, 
+                                      bool noperm, bool noscal)
+  : base_aepbal<ComplexMatrix, ColumnVector> ()
 {
   octave_idx_type n = a.cols ();
 
   if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
-      return -1;
+      return;
     }
 
   octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
 
-  Array<double> scale (n);
+  scale = ColumnVector (n);
   double *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
   Complex *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  char job = balance_job[0];
+  job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
 
   F77_XFCN (zgebal, ZGEBAL, (F77_CONST_CHAR_ARG2 (&job, 1),
 			     n, p_balanced_mat, n, ilo, ihi,
 			     pscale, info
 			     F77_CHAR_ARG_LEN (1)));
+}
 
-  balancing_mat = ComplexMatrix (n, n, 0.0);
+ComplexMatrix
+ComplexAEPBALANCE::balancing_matrix (void) const
+{
+  octave_idx_type n = balanced_mat.rows ();
+  ComplexMatrix balancing_mat (n, n, 0.0);
   for (octave_idx_type i = 0; i < n; i++)
     balancing_mat.elem (i, i) = 1.0;
 
   Complex *p_balancing_mat = balancing_mat.fortran_vec ();
+  const double *pscale = scale.fortran_vec ();
+
+  octave_idx_type info;
 
   char side = 'R';
 
@@ -92,7 +99,7 @@
 			     F77_CHAR_ARG_LEN (1)
 			     F77_CHAR_ARG_LEN (1)));
 
-  return info;
+  return balancing_mat;
 }
 
 /*
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/CmplxAEPBAL.h
--- a/liboctave/CmplxAEPBAL.h	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/CmplxAEPBAL.h	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -27,48 +28,25 @@
 #include <iostream>
 #include <string>
 
+#include "base-aepbal.h"
 #include "CMatrix.h"
+#include "dColVector.h"
 
 class
 OCTAVE_API
-ComplexAEPBALANCE
+ComplexAEPBALANCE : public base_aepbal<ComplexMatrix, ColumnVector>
 {
 public:
 
-  ComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+  ComplexAEPBALANCE (void) : base_aepbal<ComplexMatrix, ColumnVector> () { }
 
-  ComplexAEPBALANCE (const ComplexMatrix& a, const std::string& balance_job)
-    {
-      init (a, balance_job); 
-    }
-
-  ComplexAEPBALANCE (const ComplexAEPBALANCE& a)
-    : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { }
+  ComplexAEPBALANCE (const ComplexMatrix& a, bool noperm = false,
+                     bool noscal = false);
 
-  ComplexAEPBALANCE& operator = (const ComplexAEPBALANCE& a)
-    {
-      if (this != &a)
-	{
-	  balanced_mat = a.balanced_mat;
-	  balancing_mat = a.balancing_mat;
-	}
-      return *this;
-    }
-
-  ~ComplexAEPBALANCE (void) { }
+  ComplexAEPBALANCE (const ComplexAEPBALANCE& a) 
+    : base_aepbal<ComplexMatrix, ColumnVector> (a) { }
 
-  ComplexMatrix balanced_matrix (void) const { return balanced_mat; }
-
-  ComplexMatrix balancing_matrix (void) const { return balancing_mat; }
-
-  friend std::ostream& operator << (std::ostream& os, const ComplexAEPBALANCE& a);
-
-private:
-
-  ComplexMatrix balanced_mat;
-  ComplexMatrix balancing_mat;
-
-  octave_idx_type init (const ComplexMatrix& a, const std::string& balance_job);
+  ComplexMatrix balancing_matrix (void) const;
 };
 
 #endif
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/Makefile.in
--- a/liboctave/Makefile.in	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/Makefile.in	Tue Dec 09 10:08:19 2008 +0100
@@ -43,7 +43,7 @@
 MATRIX_INC := Array.h Array2.h Array3.h ArrayN.h DiagArray2.h \
 	Array-util.h ArrayN-idx.h MArray-defs.h \
 	MArray.h MArray2.h MDiagArray2.h Matrix.h MArrayN.h \
-	base-lu.h dim-vector.h mx-base.h mx-op-defs.h \
+	base-lu.h base-aepbal.h dim-vector.h mx-base.h mx-op-defs.h \
 	mx-defs.h mx-ext.h CColVector.h CDiagMatrix.h CMatrix.h \
 	CNDArray.h CRowVector.h CmplxAEPBAL.h CmplxCHOL.h \
 	CmplxGEPBAL.h CmplxHESS.h CmplxLU.h CmplxQR.h CmplxQRP.h \
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/base-aepbal.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/liboctave/base-aepbal.h	Tue Dec 09 10:08:19 2008 +0100
@@ -0,0 +1,96 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek
+
+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_base_aepbal_h)
+#define octave_base_aepbal_h 1
+
+#include "oct-types.h"
+
+template <class MatrixT, class VectorT>
+class base_aepbal
+{
+protected:
+  MatrixT balanced_mat;
+  VectorT scale;
+  octave_idx_type ilo, ihi;
+  char job;
+  base_aepbal (void) : balanced_mat (), scale (), ilo (), ihi (), job () { }
+
+public:
+
+  base_aepbal (const base_aepbal& a) 
+    : balanced_mat (a.balanced_mat), scale (a.scale), 
+      ilo(a.ilo), ihi(a.ihi), job(a.job)
+  { 
+  }
+
+  base_aepbal& operator = (const base_aepbal& a)
+    {
+      balanced_mat = a.balanced_mat;
+      scale = a.scale;
+      ilo = a.ilo;
+      ihi = a.ihi;
+      job = a.job;
+    }
+
+  MatrixT balanced_matrix (void) const { return balanced_mat; }
+
+  VectorT permuting_vector (void) const
+    {
+      octave_idx_type n = balanced_mat.rows ();
+      VectorT pv (n);
+      for (octave_idx_type i = 0; i < n; i++)
+        pv(i) = i+1;
+      for (octave_idx_type i = n-1; i >= ihi; i--)
+        {
+          octave_idx_type j = scale(i) - 1;
+          octave_idx_type k = pv(j);
+          pv(j) = pv(i);
+          pv(i) = k;
+        }
+      for (octave_idx_type i = ilo-2; i >= 0; i--)
+        {
+          octave_idx_type j = scale(i) - 1;
+          octave_idx_type k = pv(j);
+          pv(j) = pv(i);
+          pv(i) = k;
+        }
+      
+      return pv;
+    }
+
+  VectorT scaling_vector (void) const
+    {
+      octave_idx_type n = balanced_mat.rows ();
+      VectorT scv (n);
+      for (octave_idx_type i = 0; i < ilo-1; i++)
+        scv(i) = 1;
+      for (octave_idx_type i = ilo-1; i < ihi; i++)
+        scv(i) = scale(i);
+      for (octave_idx_type i = ihi; i < n; i++)
+        scv(i) = 1;
+
+      return scv;
+    }
+};
+
+#endif
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/dbleAEPBAL.cc
--- a/liboctave/dbleAEPBAL.cc	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/dbleAEPBAL.cc	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -41,44 +42,51 @@
   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&
+			     const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             const 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
-AEPBALANCE::init (const Matrix& a, const std::string& balance_job)
+AEPBALANCE::AEPBALANCE (const Matrix& a, bool noperm, bool noscal)
+  : base_aepbal<Matrix, ColumnVector> ()
 {
   octave_idx_type n = a.cols ();
 
   if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
-      return -1;
+      return;
     }
 
   octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
 
-  Array<double> scale (n);
+  scale = ColumnVector (n);
   double *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
   double *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  char job = balance_job[0];
+  job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
 
   F77_XFCN (dgebal, DGEBAL, (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);
+Matrix
+AEPBALANCE::balancing_matrix (void) const
+{
+  octave_idx_type n = balanced_mat.rows ();
+  Matrix balancing_mat (n, n, 0.0);
   for (octave_idx_type i = 0; i < n; i++)
     balancing_mat.elem (i ,i) = 1.0;
 
   double *p_balancing_mat = balancing_mat.fortran_vec ();
+  const double *pscale = scale.fortran_vec ();
+
+  octave_idx_type info;
 
   char side = 'R';
 
@@ -89,7 +97,7 @@
 			     F77_CHAR_ARG_LEN (1)
 			     F77_CHAR_ARG_LEN (1)));
 
-  return info;
+  return balancing_mat;
 }
 
 /*
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/dbleAEPBAL.h
--- a/liboctave/dbleAEPBAL.h	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/dbleAEPBAL.h	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -27,48 +28,25 @@
 #include <iostream>
 #include <string>
 
+#include "base-aepbal.h"
 #include "dMatrix.h"
+#include "dColVector.h"
 
 class
 OCTAVE_API
-AEPBALANCE
+AEPBALANCE : public base_aepbal<Matrix, ColumnVector>
 {
 public:
 
-  AEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+  AEPBALANCE (void) : base_aepbal<Matrix, ColumnVector> () { }
 
-  AEPBALANCE (const Matrix& a,const std::string& balance_job)
-    {
-      init (a, balance_job); 
-    }
-
-  AEPBALANCE (const AEPBALANCE& a)
-    : balanced_mat (a.balanced_mat), balancing_mat (a.balancing_mat) { }
+  AEPBALANCE (const Matrix& a, bool noperm = false,
+              bool noscal = false);
 
-  AEPBALANCE& operator = (const AEPBALANCE& a)
-    {
-      if (this != &a)
-	{
-	  balanced_mat = a.balanced_mat;
-	  balancing_mat = a.balancing_mat;
-	}
-      return *this;
-    }
-
-  ~AEPBALANCE (void) { }
+  AEPBALANCE (const AEPBALANCE& a) 
+    : base_aepbal<Matrix, ColumnVector> (a) { }
 
-  Matrix balanced_matrix (void) const { return balanced_mat; }
-
-  Matrix balancing_matrix (void) const { return balancing_mat; }
-
-  friend std::ostream& operator << (std::ostream& os, const AEPBALANCE& a);
-
-private:
-
-  Matrix balanced_mat;
-  Matrix balancing_mat;
-
-  octave_idx_type init (const Matrix& a, const std::string& balance_job);
+  Matrix balancing_matrix (void) const;
 };
 
 #endif
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/fCmplxAEPBAL.cc
--- a/liboctave/fCmplxAEPBAL.cc	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/fCmplxAEPBAL.cc	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -42,46 +43,52 @@
   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&, const 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)
+FloatComplexAEPBALANCE::FloatComplexAEPBALANCE (const FloatComplexMatrix& a, 
+                                                bool noperm, bool noscal)
+ : base_aepbal<FloatComplexMatrix, FloatColumnVector> ()
 {
   octave_idx_type n = a.cols ();
 
   if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
-      return -1;
+      return;
     }
 
   octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
 
-  Array<float> scale (n);
+  scale = FloatColumnVector (n);
   float *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
   FloatComplex *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  char job = balance_job[0];
+  job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
 
   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);
+FloatComplexMatrix
+FloatComplexAEPBALANCE::balancing_matrix (void) const
+{
+  octave_idx_type n = balanced_mat.rows ();
+  FloatComplexMatrix balancing_mat (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 ();
+  const float *pscale = scale.fortran_vec ();
+
+  octave_idx_type info;
 
   char side = 'R';
 
@@ -92,7 +99,7 @@
 			     F77_CHAR_ARG_LEN (1)
 			     F77_CHAR_ARG_LEN (1)));
 
-  return info;
+  return balancing_mat;
 }
 
 /*
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/fCmplxAEPBAL.h
--- a/liboctave/fCmplxAEPBAL.h	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/fCmplxAEPBAL.h	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -27,48 +28,25 @@
 #include <iostream>
 #include <string>
 
+#include "base-aepbal.h"
 #include "fCMatrix.h"
+#include "fColVector.h"
 
 class
 OCTAVE_API
-FloatComplexAEPBALANCE
+FloatComplexAEPBALANCE : public base_aepbal<FloatComplexMatrix, FloatColumnVector>
 {
 public:
 
-  FloatComplexAEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+  FloatComplexAEPBALANCE (void) : base_aepbal<FloatComplexMatrix, FloatColumnVector> () { }
 
-  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 (const FloatComplexMatrix& a, bool noperm = false,
+                          bool noscal = false);
 
-  FloatComplexAEPBALANCE& operator = (const FloatComplexAEPBALANCE& a)
-    {
-      if (this != &a)
-	{
-	  balanced_mat = a.balanced_mat;
-	  balancing_mat = a.balancing_mat;
-	}
-      return *this;
-    }
-
-  ~FloatComplexAEPBALANCE (void) { }
+  FloatComplexAEPBALANCE (const FloatComplexAEPBALANCE& a) 
+    : base_aepbal<FloatComplexMatrix, FloatColumnVector> (a) { }
 
-  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);
+  FloatComplexMatrix balancing_matrix (void) const;
 };
 
 #endif
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/floatAEPBAL.cc
--- a/liboctave/floatAEPBAL.cc	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/floatAEPBAL.cc	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, 2005,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -41,44 +42,52 @@
   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&
+			     const octave_idx_type&, const octave_idx_type&, const octave_idx_type&, 
+                             const 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)
+FloatAEPBALANCE::FloatAEPBALANCE (const FloatMatrix& a, 
+                                  bool noperm, bool noscal)
+  : base_aepbal<FloatMatrix, FloatColumnVector> ()
 {
   octave_idx_type n = a.cols ();
 
   if (a.rows () != n)
     {
       (*current_liboctave_error_handler) ("AEPBALANCE requires square matrix");
-      return -1;
+      return;
     }
 
   octave_idx_type info;
-  octave_idx_type ilo;
-  octave_idx_type ihi;
 
-  Array<float> scale (n);
+  scale = FloatColumnVector (n);
   float *pscale = scale.fortran_vec ();
 
   balanced_mat = a;
   float *p_balanced_mat = balanced_mat.fortran_vec ();
 
-  char job = balance_job[0];
+  job = noperm ? (noscal ? 'N' : 'S') : (noscal ? 'P' : 'B');
 
   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 = FloatMatrix (n, n, 0.0);
+FloatMatrix
+FloatAEPBALANCE::balancing_matrix (void) const
+{
+  octave_idx_type n = balanced_mat.rows ();
+  FloatMatrix balancing_mat (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 ();
+  const float *pscale = scale.fortran_vec ();
+
+  octave_idx_type info;
 
   char side = 'R';
 
@@ -89,7 +98,7 @@
 			     F77_CHAR_ARG_LEN (1)
 			     F77_CHAR_ARG_LEN (1)));
 
-  return info;
+  return balancing_mat;
 }
 
 /*
diff -r 6e9660cd3bf2 -r a5e080076778 liboctave/floatAEPBAL.h
--- a/liboctave/floatAEPBAL.h	Tue Dec 09 03:25:31 2008 +0100
+++ b/liboctave/floatAEPBAL.h	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2000, 2002, 2004, 2005, 2006,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -27,48 +28,25 @@
 #include <iostream>
 #include <string>
 
+#include "base-aepbal.h"
 #include "fMatrix.h"
+#include "fColVector.h"
 
 class
 OCTAVE_API
-FloatAEPBALANCE
+FloatAEPBALANCE : public base_aepbal<FloatMatrix, FloatColumnVector>
 {
 public:
 
-  FloatAEPBALANCE (void) : balanced_mat (), balancing_mat () { }
+  FloatAEPBALANCE (void) : base_aepbal<FloatMatrix, FloatColumnVector> () { }
 
-  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 (const FloatMatrix& a, bool noperm = false,
+                   bool noscal = false);
 
-  FloatAEPBALANCE& operator = (const FloatAEPBALANCE& a)
-    {
-      if (this != &a)
-	{
-	  balanced_mat = a.balanced_mat;
-	  balancing_mat = a.balancing_mat;
-	}
-      return *this;
-    }
-
-  ~FloatAEPBALANCE (void) { }
+  FloatAEPBALANCE (const FloatAEPBALANCE& a) 
+    : base_aepbal<FloatMatrix, FloatColumnVector> (a) { }
 
-  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);
+  FloatMatrix balancing_matrix (void) const;
 };
 
 #endif
diff -r 6e9660cd3bf2 -r a5e080076778 src/ChangeLog
--- a/src/ChangeLog	Tue Dec 09 03:25:31 2008 +0100
+++ b/src/ChangeLog	Tue Dec 09 10:08:19 2008 +0100
@@ -1,3 +1,7 @@
+2008-12-09  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/balance.cc (Fbalance): Exploit the new AEPBAL functionality.
+
 2008-12-08  Jaroslav Hajek  <highegg@gmail.com>
 	
 	* xpow.cc ( xpow (const DiagMatrix& a, double b), 
diff -r 6e9660cd3bf2 -r a5e080076778 src/DLD-FUNCTIONS/balance.cc
--- a/src/DLD-FUNCTIONS/balance.cc	Tue Dec 09 03:25:31 2008 +0100
+++ b/src/DLD-FUNCTIONS/balance.cc	Tue Dec 09 10:08:19 2008 +0100
@@ -2,6 +2,7 @@
 
 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006,
               2007 John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
@@ -50,6 +51,7 @@
   "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {@var{aa} =} balance (@var{a}, @var{opt})\n\
 @deftypefnx {Loadable Function} {[@var{dd}, @var{aa}] =} balance (@var{a}, @var{opt})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{d}, @var{aa}] =} balance (@var{a}, @var{opt})\n\
 @deftypefnx {Loadable Function} {[@var{cc}, @var{dd}, @var{aa}, @var{bb}] =} balance (@var{a}, @var{b}, @var{opt})\n\
 \n\
 Compute @code{aa = dd \\ a * dd} in which @code{aa} is a matrix whose\n\
@@ -59,6 +61,11 @@
 the equilibration to be computed without roundoff.  Results of\n\
 eigenvalue calculation are typically improved by balancing first.\n\
 \n\
+If two output values are requested, @code{balance} returns \n\
+the permutation @code{p} and the diagonal @code{d} separately as vectors. \n\
+In this case, @code{dd = eye(n)(p,:) * diag (d)}, where @code{n} is the matrix \n\
+size. \n\
+\n\
 If four output values are requested, compute @code{aa = cc*a*dd} and\n\
 @code{bb = cc*b*dd)}, in which @code{aa} and @code{bb} have non-zero\n\
 elements of approximately the same magnitude and @code{cc} and @code{dd}\n\
@@ -68,19 +75,11 @@
 The eigenvalue balancing option @code{opt} may be one of:\n\
 \n\
 @table @asis\n\
-@item @code{\"N\"}, @code{\"n\"}\n\
-No balancing; arguments copied, transformation(s) set to identity.\n\
-\n\
-@item @code{\"P\"}, @code{\"p\"}\n\
-Permute argument(s) to isolate eigenvalues where possible.\n\
+@item @code{\"noperm\"}, @code{\"S\"}\n\
+Scale only; do not permute.\n\
 \n\
-@item @code{\"S\"}, @code{\"s\"}\n\
-Scale to improve accuracy of computed eigenvalues.\n\
-\n\
-@item @code{\"B\"}, @code{\"b\"}\n\
-Permute and scale, in that order. Rows/columns of a (and b)\n\
-that are isolated by permutation are not scaled.  This is the default\n\
-behavior.\n\
+@item @code{\"noscal\"}, @code{\"P\"}\n\
+Permute only; do not scale.\n\
 @end table\n\
 \n\
 Algebraic eigenvalue balancing uses standard @sc{Lapack} routines.\n\
@@ -100,20 +99,11 @@
     }
 
   // determine if it's AEP or GEP
-  int AEPcase = nargin == 1 ? 1 : args(1).is_string ();
-  std::string bal_job;
+  bool AEPcase = nargin == 1 || args(1).is_string ();
 
   // problem dimension
   octave_idx_type nn = args(0).rows ();
 
-  octave_idx_type arg_is_empty = empty_arg ("balance", nn, args(0).columns());
-
-  if (arg_is_empty < 0)
-    return retval;
-
-  if (arg_is_empty > 0)
-    return octave_value_list (2, Matrix ());
-
   if (nn != args(0).columns())
     {
       gripe_square_matrix_required ("balance");
@@ -154,75 +144,98 @@
   if (AEPcase)
     {  
       // Algebraic eigenvalue problem.
-
-      if (nargin == 1)
-	bal_job = "B";
-      else if (args(1).is_string ())
-	bal_job = args(1).string_value ();
-      else
-	{
-	  error ("balance: AEP argument 2 must be a string");
-	  return retval;
-	}
+      bool noperm = false, noscal = false;
+      if (nargin > 1)
+        {
+          std::string a1s = args(1).string_value ();
+          noperm = a1s == "noperm" || a1s == "S";
+          noscal = a1s == "noscal" || a1s == "P";
+        }
 
       // balance the AEP
       if (isfloat)
 	{
 	  if (complex_case)
 	    {
-	      FloatComplexAEPBALANCE result (fcaa, bal_job);
+	      FloatComplexAEPBALANCE result (fcaa, noperm, noscal);
 
 	      if (nargout == 0 || nargout == 1)
 		retval(0) = result.balanced_matrix ();
-	      else
+	      else if (nargout == 2)
 		{
 		  retval(1) = result.balanced_matrix ();
 		  retval(0) = result.balancing_matrix ();
 		}
+              else
+                {
+                  retval(2) = result.balanced_matrix ();
+                  retval(1) = result.scaling_vector ();
+                  retval(0) = result.permuting_vector ();
+                }
+
 	    }
 	  else
 	    {
-	      FloatAEPBALANCE result (faa, bal_job);
+	      FloatAEPBALANCE result (faa, noperm, noscal);
 
 	      if (nargout == 0 || nargout == 1)
 		retval(0) = result.balanced_matrix ();
-	      else
+	      else if (nargout == 2)
 		{
 		  retval(1) = result.balanced_matrix ();
 		  retval(0) = result.balancing_matrix ();
 		}
+              else
+                {
+                  retval(2) = result.balanced_matrix ();
+                  retval(1) = result.scaling_vector ();
+                  retval(0) = result.permuting_vector ();
+                }
 	    }
 	}
       else
 	{
 	  if (complex_case)
 	    {
-	      ComplexAEPBALANCE result (caa, bal_job);
+	      ComplexAEPBALANCE result (caa, noperm, noscal);
 
 	      if (nargout == 0 || nargout == 1)
 		retval(0) = result.balanced_matrix ();
-	      else
+	      else if (nargout == 2)
 		{
 		  retval(1) = result.balanced_matrix ();
 		  retval(0) = result.balancing_matrix ();
 		}
+              else
+                {
+                  retval(2) = result.balanced_matrix ();
+                  retval(1) = result.scaling_vector ();
+                  retval(0) = result.permuting_vector ();
+                }
 	    }
 	  else
 	    {
-	      AEPBALANCE result (aa, bal_job);
+	      AEPBALANCE result (aa, noperm, noscal);
 
 	      if (nargout == 0 || nargout == 1)
 		retval(0) = result.balanced_matrix ();
-	      else
+	      else if (nargout == 2)
 		{
 		  retval(1) = result.balanced_matrix ();
 		  retval(0) = result.balancing_matrix ();
 		}
+              else
+                {
+                  retval(2) = result.balanced_matrix ();
+                  retval(1) = result.scaling_vector ();
+                  retval(0) = result.permuting_vector ();
+                }
 	    }
 	}
     }
   else
     {
+      std::string bal_job;
       if (nargout == 1)
 	warning ("balance: used GEP, should have two output arguments");