# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1228375916 -3600
# Node ID 445d27d79f4e65c65faddd896d09584aa95c4fb6
# Parent  8b1a2555c4e23d9b3e89246deb738f84b48c0471
support permutation matrix objects

diff --git a/liboctave/Array2.h b/liboctave/Array2.h
--- a/liboctave/Array2.h
+++ b/liboctave/Array2.h
@@ -116,14 +116,14 @@
       return Array2<T> (tmp, tmp.rows (), tmp.columns ());
     }
 
-  Array2<T> index (idx_vector& i, int resize_ok = 0,
+  Array2<T> index (const idx_vector& i, int resize_ok = 0,
 		   const T& rfv = Array<T>::resize_fill_value ()) const
     {
       Array<T> tmp = Array<T>::index (i, resize_ok, rfv);
       return Array2<T> (tmp, tmp.rows (), tmp.columns ());
     }
 
-  Array2<T> index (idx_vector& i, idx_vector& j, int resize_ok = 0,
+  Array2<T> index (const idx_vector& i, const idx_vector& j, int resize_ok = 0,
 		   const T& rfv = Array<T>::resize_fill_value ()) const
     {
       Array<T> tmp = Array<T>::index (i, j, resize_ok, rfv);
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,25 @@
+2008-12-03  Jaroslav Hajek  <highegg@gmail.com>
+
+	* PermMatrix.h, PermMatrix.cc: New sources.
+	* MDiagArray2.cc (MDiagArray2<T>::is_multiple_of_identity): New method.
+	* MDiagArray2.h (MDiagArray2<T>::is_multiple_of_identity): Declare it.
+	* idx-vector.cc (idx_vector::is_permutation): New method.
+	* idx-vector.h (idx_vector::is_permutation): Declare it.
+	* base-lu.cc (base_lu::getp): New method.
+	(base_lu::P): Call getp.
+	(base_lu::Pvec): Call getp.
+	* base-lu.h (base_lu): Delcare P as PermMatrix. Remove unused template
+	params.
+	* dbleQRP.cc (dbleQRP::dbleQRP): Construct a permutation matrix.
+	(dbleQRP::Pvec): New method.
+	* dbleQRP.h: Declare new method. Declare P as PermMatrix.
+	* CmplxQRP.cc (ComplexQRP): Likewise.
+	* CmplxQRP.h (ComplexQRP): Likewise.
+	* floatQRP.cc (FloatQRP): Likewise.
+	* floatQRP.h (FloatQRP): Likewise.
+	* fCmplxQRP.cc (FloatComplexQRP): Likewise.
+	* fCmplxQRP.h (FloatComplexQRP): Likewise.
+
 2008-12-01  Jaroslav Hajek  <highegg@gmail.com>
 
 	* DiagArray2.h (DiagArray2<T>::DiagArray2<T> (const DiagArray2<U>&)): New template
diff --git a/liboctave/CmplxLU.cc b/liboctave/CmplxLU.cc
--- a/liboctave/CmplxLU.cc
+++ b/liboctave/CmplxLU.cc
@@ -34,7 +34,7 @@
 #include <base-lu.h>
 #include <base-lu.cc>
 
-template class base_lu <ComplexMatrix, Complex, Matrix, double>;
+template class base_lu <ComplexMatrix>;
 
 // Define the constructor for this particular derivation.
 
diff --git a/liboctave/CmplxLU.h b/liboctave/CmplxLU.h
--- a/liboctave/CmplxLU.h
+++ b/liboctave/CmplxLU.h
@@ -30,22 +30,22 @@
 
 class
 OCTAVE_API
-ComplexLU : public base_lu <ComplexMatrix, Complex, Matrix, double>
+ComplexLU : public base_lu <ComplexMatrix>
 {
 public:
 
   ComplexLU (void)
-    : base_lu <ComplexMatrix, Complex, Matrix, double> () { }
+    : base_lu <ComplexMatrix> () { }
 
   ComplexLU (const ComplexMatrix& a);
 
   ComplexLU (const ComplexLU& a)
-    : base_lu <ComplexMatrix, Complex, Matrix, double> (a) { }
+    : base_lu <ComplexMatrix> (a) { }
 
   ComplexLU& operator = (const ComplexLU& a)
     {
       if (this != &a)
-	base_lu <ComplexMatrix, Complex, Matrix, double> :: operator = (a);
+	base_lu <ComplexMatrix> :: operator = (a);
 
       return *this;
     }
diff --git a/liboctave/CmplxQRP.cc b/liboctave/CmplxQRP.cc
--- a/liboctave/CmplxQRP.cc
+++ b/liboctave/CmplxQRP.cc
@@ -86,7 +86,7 @@
   Array<double> rwork (2*n);
   double *prwork = rwork.fortran_vec ();
 
-  Array<octave_idx_type> jpvt (n, 0);
+  MArray<octave_idx_type> jpvt (n, 0);
   octave_idx_type *pjpvt = jpvt.fortran_vec ();
 
   // Code to enforce a certain permutation could go here...
@@ -97,18 +97,8 @@
   // 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;
-    }
+  jpvt -= 1;
+  p = PermMatrix (jpvt, true);
 
   octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
 
@@ -131,6 +121,14 @@
   q.resize (m, n2);
 }
 
+ColumnVector
+ComplexQRP::Pvec (void) const
+{
+  Array<double> pa (p);
+  ColumnVector pv (MArray<double> (pa) + 1.0);
+  return pv;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/CmplxQRP.h b/liboctave/CmplxQRP.h
--- a/liboctave/CmplxQRP.h
+++ b/liboctave/CmplxQRP.h
@@ -27,6 +27,8 @@
 #include <iostream>
 
 #include "CmplxQR.h"
+#include "PermMatrix.h"
+#include "dColVector.h"
 
 class
 OCTAVE_API
@@ -54,13 +56,15 @@
 
   void init (const ComplexMatrix&, QR::type = QR::std);
 
-  Matrix P (void) const { return p; }
+  PermMatrix P (void) const { return p; }
+
+  ColumnVector Pvec (void) const;
 
   friend std::ostream&  operator << (std::ostream&, const ComplexQRP&);
 
 private:
 
-  Matrix p;
+  PermMatrix p;
 };
 
 #endif
diff --git a/liboctave/MDiagArray2.cc b/liboctave/MDiagArray2.cc
--- a/liboctave/MDiagArray2.cc
+++ b/liboctave/MDiagArray2.cc
@@ -31,6 +31,22 @@
 
 #include "MArray-defs.h"
 
+template <class T>
+bool 
+MDiagArray2<T>::is_multiple_of_identity (T val) const
+{
+  bool retval = this->rows () == this->cols ();
+  if (retval)
+    {
+      octave_idx_type len = this->length (), i = 0;
+      for (;i < len; i++) 
+        if (DiagArray2<T>::elem (i, i) != val) break;
+      retval = i == len;
+    }
+
+  return retval;
+}
+
 // Some functions return a reference to this object after a failure.
 template <class T> MDiagArray2<T> MDiagArray2<T>::nil_array;
 
diff --git a/liboctave/MDiagArray2.h b/liboctave/MDiagArray2.h
--- a/liboctave/MDiagArray2.h
+++ b/liboctave/MDiagArray2.h
@@ -104,6 +104,8 @@
   MDiagArray2<T> transpose (void) const { return DiagArray2<T>::transpose (); }
   MDiagArray2<T> hermitian (T (*fcn) (const T&) = 0) const { return DiagArray2<T>::hermitian (fcn); }
 
+  bool is_multiple_of_identity (T val) const;
+
   static MDiagArray2<T> nil_array;
 
   // Currently, the OPS functions don't need to be friends, but that
diff --git a/liboctave/Makefile.in b/liboctave/Makefile.in
--- a/liboctave/Makefile.in
+++ b/liboctave/Makefile.in
@@ -55,7 +55,7 @@
 	Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \
 	sparse-base-chol.h SparseCmplxCHOL.h \
 	SparsedbleCHOL.h SparseCmplxQR.h SparseQR.h Sparse-op-defs.h \
-	MatrixType.h \
+	MatrixType.h PermMatrix.h \
 	int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \
 	int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \
 	intNDArray.h \
@@ -125,7 +125,7 @@
 	dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \
 	MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \
 	SparseCmplxCHOL.cc SparsedbleCHOL.cc \
-	SparseCmplxQR.cc SparseQR.cc MatrixType.cc \
+	SparseCmplxQR.cc SparseQR.cc MatrixType.cc PermMatrix.cc \
 	int8NDArray.cc uint8NDArray.cc int16NDArray.cc uint16NDArray.cc \
 	int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc \
 	fCColVector.cc fCRowVector.cc fCDiagMatrix.cc fCMatrix.cc fCNDArray.cc \
diff --git a/liboctave/PermMatrix.cc b/liboctave/PermMatrix.cc
new file mode 100644
--- /dev/null
+++ b/liboctave/PermMatrix.cc
@@ -0,0 +1,148 @@
+/*
+
+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/>.
+
+*/
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "PermMatrix.h"
+#include "idx-vector.h"
+#include "error.h"
+#include "Array-util.h"
+
+static void
+gripe_invalid_permutation (void)
+{
+  (*current_liboctave_error_handler)
+    ("PermMatrix: invalid permutation vector");
+}
+
+PermMatrix::PermMatrix (const Array<octave_idx_type>& p, bool colp, bool check)
+  : Array<octave_idx_type> (p), _colp(colp)
+{
+  this->dimensions = dim_vector (p.length (), p.length ());
+  if (check)
+    {
+      if (! idx_vector (p).is_permutation (p.length ()))
+        {
+          gripe_invalid_permutation ();
+          Array<octave_idx_type>::operator = (Array<octave_idx_type> ());
+        }
+    }
+}
+
+PermMatrix::PermMatrix (const idx_vector& idx, bool colp, octave_idx_type n)
+  : Array<octave_idx_type> (), _colp(colp)
+{
+  octave_idx_type len = idx.length (n);
+  if (! idx.is_permutation (len))
+    gripe_invalid_permutation ();
+  else
+    {
+      Array<octave_idx_type> idxa (len);
+      for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i);
+      Array<octave_idx_type>::operator = (idxa);
+      this->dimensions = dim_vector (len, len);
+    }
+}
+
+PermMatrix::PermMatrix (octave_idx_type n)
+  : Array<octave_idx_type> (n), _colp (false)
+{
+  this->dimensions = dim_vector (n, n);
+  for (octave_idx_type i = 0; i < n; i++) xelem (i) = i;
+}
+
+octave_idx_type 
+PermMatrix::checkelem (octave_idx_type i, octave_idx_type j) const
+{
+  octave_idx_type len = Array<octave_idx_type>::length ();
+  if (i < 0 || j < 0 || i > len || j > len)
+    {
+      (*current_liboctave_error_handler) ("index out of range");
+      return 0;
+    }
+  else
+    return elem (i, j);
+}
+
+
+PermMatrix 
+PermMatrix::transpose (void) const
+{
+  PermMatrix retval (*this);
+  retval._colp = ! retval._colp;
+  return retval;
+}
+
+PermMatrix 
+PermMatrix::inverse (void) const
+{
+  return transpose ();
+}
+
+octave_idx_type 
+PermMatrix::determinant (void) const
+{
+  Array<octave_idx_type> pa = *this;
+  octave_idx_type len = pa.length (), *p = pa.fortran_vec ();
+  bool neg = false;
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      octave_idx_type j = p[i];
+      if (j != i)
+        {
+          p[i] = p[j];
+          p[j] = j;
+          neg = ! neg;
+        }
+    }
+  
+  return neg ? -1 : 1;
+}
+
+PermMatrix 
+operator *(const PermMatrix& a, const PermMatrix& b)
+{
+  const Array<octave_idx_type>& ia = a, ib = b;
+  PermMatrix r;
+  octave_idx_type n = a.columns ();
+  if (n != b.rows ())
+    gripe_nonconformant ("operator *", n, n, b.rows (), b.rows ());
+  else if (a._colp == b._colp)
+    {
+      r = PermMatrix ((a._colp 
+                       ? ia.index (idx_vector (ib)) 
+                       : ib.index (idx_vector (ia))), a._colp, false);
+    }
+  else
+    {
+      Array<octave_idx_type> ra (n);
+      if (a._colp)
+        ra.assign (idx_vector (ib), ia);
+      else
+        ra.assign (idx_vector (ia), ib);
+      r = PermMatrix (ra, a._colp, false);
+    }
+
+  return r;
+}
diff --git a/liboctave/PermMatrix.h b/liboctave/PermMatrix.h
new file mode 100644
--- /dev/null
+++ b/liboctave/PermMatrix.h
@@ -0,0 +1,91 @@
+/*
+
+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_PermMatrix_h)
+#define octave_PermMatrix_h 1
+
+#include "Array.h"
+#include "mx-defs.h"
+
+class PermMatrix : public Array<octave_idx_type>
+{
+private:
+
+  octave_idx_type get (octave_idx_type i) const { return Array<octave_idx_type>::xelem (i); }
+
+public:
+
+  PermMatrix (void) : Array<octave_idx_type> (), _colp (false) { }
+
+  PermMatrix (octave_idx_type n);
+
+  PermMatrix (const Array<octave_idx_type>& p, bool colp = false, 
+              bool check = true);
+
+  PermMatrix (const PermMatrix& m)
+    : Array<octave_idx_type> (m), _colp(m._colp) 
+    { this->dimensions = m.dims (); }
+  
+  PermMatrix (const idx_vector& idx, bool colp = false, octave_idx_type n = 0); 
+
+  octave_idx_type 
+  elem (octave_idx_type i, octave_idx_type j) const
+    {
+      return (_colp 
+              ? ((get(j) != i) ? 1 : 0)
+              : ((get(i) != j) ? 1 : 0));
+    }
+
+  octave_idx_type 
+  checkelem (octave_idx_type i, octave_idx_type j) const;
+
+  octave_idx_type
+  operator () (octave_idx_type i, octave_idx_type j) const
+    {
+#if defined (BOUNDS_CHECKING)
+      return checkelem (i, j);
+#else
+      return elem (i, j);
+#endif
+    }
+  
+  // These are, in fact, super-fast.
+  PermMatrix transpose (void) const;
+  PermMatrix inverse (void) const;
+
+  // Determinant, i.e. the sign of permutation.
+  octave_idx_type determinant (void) const;
+
+  bool is_col_perm (void) const { return _colp; }
+  bool is_row_perm (void) const { return !_colp; }
+
+  friend PermMatrix operator *(const PermMatrix& a, const PermMatrix& b);
+
+private:
+  bool _colp;
+};
+
+// Multiplying permutations together.
+PermMatrix 
+operator *(const PermMatrix& a, const PermMatrix& b);
+
+#endif
diff --git a/liboctave/base-lu.cc b/liboctave/base-lu.cc
--- a/liboctave/base-lu.cc
+++ b/liboctave/base-lu.cc
@@ -26,9 +26,9 @@
 
 #include "base-lu.h"
 
-template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+template <class lu_type>
 lu_type
-base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: L (void) const
+base_lu <lu_type> :: L (void) const
 {
   octave_idx_type a_nr = a_fact.rows ();
   octave_idx_type a_nc = a_fact.cols ();
@@ -48,9 +48,9 @@
   return l;
 }
 
-template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+template <class lu_type>
 lu_type
-base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: U (void) const
+base_lu <lu_type> :: U (void) const
 {
   octave_idx_type a_nr = a_fact.rows ();
   octave_idx_type a_nc = a_fact.cols ();
@@ -67,9 +67,9 @@
   return u;
 }
 
-template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
-p_type
-base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: P (void) const
+template <class lu_type>
+Array<octave_idx_type>
+base_lu <lu_type> :: getp (void) const
 {
   octave_idx_type a_nr = a_fact.rows ();
 
@@ -90,38 +90,25 @@
 	}
     }
 
-  p_type p (a_nr, a_nr, p_elt_type (0.0));
-
-  for (octave_idx_type i = 0; i < a_nr; i++)
-    p.xelem (i, pvt.xelem (i)) = 1.0;
-
-  return p;
+  return pvt;
 }
 
-template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+template <class lu_type>
+PermMatrix
+base_lu <lu_type> :: P (void) const
+{
+  return PermMatrix (getp (), false);
+}
+
+template <class lu_type>
 ColumnVector
-base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: P_vec (void) const
+base_lu <lu_type> :: P_vec (void) const
 {
   octave_idx_type a_nr = a_fact.rows ();
 
-  Array<octave_idx_type> pvt (a_nr);
-
-  for (octave_idx_type i = 0; i < a_nr; i++)
-    pvt.xelem (i) = i;
-
-  for (octave_idx_type i = 0; i < ipvt.length(); i++)
-    {
-      octave_idx_type k = ipvt.xelem (i);
+  ColumnVector p (a_nr);
 
-      if (k != i)
-	{
-	  octave_idx_type tmp = pvt.xelem (k);
-	  pvt.xelem (k) = pvt.xelem (i);
-	  pvt.xelem (i) = tmp;
-	}
-    }
-
-  ColumnVector p (a_nr);
+  Array<octave_idx_type> pvt = getp ();
 
   for (octave_idx_type i = 0; i < a_nr; i++)
     p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1);
diff --git a/liboctave/base-lu.h b/liboctave/base-lu.h
--- a/liboctave/base-lu.h
+++ b/liboctave/base-lu.h
@@ -25,13 +25,16 @@
 
 #include "MArray.h"
 #include "dColVector.h"
+#include "PermMatrix.h"
 
-template <class lu_type, class lu_elt_type, class p_type, class p_elt_type>
+template <class lu_type>
 class
 base_lu
 {
 public:
 
+  typedef typename lu_type::element_type lu_elt_type;
+
   base_lu (void) { }
 
   base_lu (const base_lu& a) : a_fact (a.a_fact), ipvt (a.ipvt) { }
@@ -54,12 +57,13 @@
 
   lu_type Y (void) const { return a_fact; }
 
-  p_type P (void) const;
+  PermMatrix P (void) const;
 
   ColumnVector P_vec (void) const;
 
 protected:
 
+  Array<octave_idx_type> getp (void) const;
   lu_type a_fact;
   MArray<octave_idx_type> ipvt;
 };
diff --git a/liboctave/dDiagMatrix.cc b/liboctave/dDiagMatrix.cc
--- a/liboctave/dDiagMatrix.cc
+++ b/liboctave/dDiagMatrix.cc
@@ -387,6 +387,7 @@
   return d;
 }
 
+
 std::ostream&
 operator << (std::ostream& os, const DiagMatrix& a)
 {
diff --git a/liboctave/dDiagMatrix.h b/liboctave/dDiagMatrix.h
--- a/liboctave/dDiagMatrix.h
+++ b/liboctave/dDiagMatrix.h
@@ -98,6 +98,8 @@
 
   ColumnVector diag (octave_idx_type k = 0) const;
 
+  bool is_identity (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const DiagMatrix& a);
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -248,6 +248,19 @@
     elem (i, i) = a.elem (i, i);
 }
 
+Matrix::Matrix (const PermMatrix& a)
+  : MArray2<double> (a.rows (), a.cols (), 0.0)
+{
+  const Array<octave_idx_type> ia (a);
+  octave_idx_type len = a.rows ();
+  if (a.is_col_perm ())
+    for (octave_idx_type i = 0; i < len; i++)
+      elem (ia(i), i) = 1.0;
+  else
+    for (octave_idx_type i = 0; i < len; i++)
+      elem (i, ia(i)) = 1.0;
+}
+
 // FIXME -- could we use a templated mixed-type copy function
 // here?
 
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -64,6 +64,8 @@
 
   explicit Matrix (const DiagMatrix& a);
 
+  explicit Matrix (const PermMatrix& a);
+
   explicit Matrix (const boolMatrix& a);
 
   explicit Matrix (const charMatrix& a);
diff --git a/liboctave/dbleLU.cc b/liboctave/dbleLU.cc
--- a/liboctave/dbleLU.cc
+++ b/liboctave/dbleLU.cc
@@ -34,7 +34,7 @@
 #include <base-lu.h>
 #include <base-lu.cc>
 
-template class base_lu <Matrix, double, Matrix, double>;
+template class base_lu <Matrix>;
 
 // Define the constructor for this particular derivation.
 
diff --git a/liboctave/dbleLU.h b/liboctave/dbleLU.h
--- a/liboctave/dbleLU.h
+++ b/liboctave/dbleLU.h
@@ -29,20 +29,20 @@
 
 class
 OCTAVE_API
-LU : public base_lu <Matrix, double, Matrix, double>
+LU : public base_lu <Matrix>
 {
 public:
 
-  LU (void) : base_lu <Matrix, double, Matrix, double> () { }
+  LU (void) : base_lu <Matrix> () { }
 
   LU (const Matrix& a);
 
-  LU (const LU& a) : base_lu <Matrix, double, Matrix, double> (a) { }
+  LU (const LU& a) : base_lu <Matrix> (a) { }
 
   LU& operator = (const LU& a)
     {
       if (this != &a)
-	base_lu <Matrix, double, Matrix, double> :: operator = (a);
+	base_lu <Matrix> :: operator = (a);
 
       return *this;
     }
diff --git a/liboctave/dbleQR.cc b/liboctave/dbleQR.cc
--- a/liboctave/dbleQR.cc
+++ b/liboctave/dbleQR.cc
@@ -2,6 +2,7 @@
 
 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
               John W. Eaton
+Copyright (C) 2008 Jaroslav Hajek
 
 This file is part of Octave.
 
diff --git a/liboctave/dbleQRP.cc b/liboctave/dbleQRP.cc
--- a/liboctave/dbleQRP.cc
+++ b/liboctave/dbleQRP.cc
@@ -81,7 +81,7 @@
 
   double *tmp_data = A_fact.fortran_vec ();
 
-  Array<octave_idx_type> jpvt (n, 0);
+  MArray<octave_idx_type> jpvt (n, 0);
   octave_idx_type *pjpvt = jpvt.fortran_vec ();
 
   // Code to enforce a certain permutation could go here...
@@ -91,18 +91,8 @@
   // 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;
-    }
+  jpvt -= 1;
+  p = PermMatrix (jpvt, true);
 
   octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
 
@@ -125,6 +115,14 @@
   q.resize (m, n2);
 }
 
+ColumnVector
+QRP::Pvec (void) const
+{
+  Array<double> pa (p);
+  ColumnVector pv (MArray<double> (pa) + 1.0);
+  return pv;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/dbleQRP.h b/liboctave/dbleQRP.h
--- a/liboctave/dbleQRP.h
+++ b/liboctave/dbleQRP.h
@@ -27,6 +27,8 @@
 #include <iostream>
 
 #include "dbleQR.h"
+#include "PermMatrix.h"
+#include "dColVector.h"
 
 class
 OCTAVE_API
@@ -55,13 +57,15 @@
 
   void init (const Matrix&, QR::type = QR::std);
 
-  Matrix P (void) const { return p; }
+  PermMatrix P (void) const { return p; }
+
+  ColumnVector Pvec (void) const;
 
   friend std::ostream&  operator << (std::ostream&, const QRP&);
 
 protected:
 
-  Matrix p;
+  PermMatrix p;
 };
 
 #endif
diff --git a/liboctave/fCmplxLU.cc b/liboctave/fCmplxLU.cc
--- a/liboctave/fCmplxLU.cc
+++ b/liboctave/fCmplxLU.cc
@@ -34,7 +34,7 @@
 #include <base-lu.h>
 #include <base-lu.cc>
 
-template class base_lu <FloatComplexMatrix, FloatComplex, Matrix, double>;
+template class base_lu <FloatComplexMatrix>;
 
 // Define the constructor for this particular derivation.
 
diff --git a/liboctave/fCmplxLU.h b/liboctave/fCmplxLU.h
--- a/liboctave/fCmplxLU.h
+++ b/liboctave/fCmplxLU.h
@@ -30,22 +30,22 @@
 
 class
 OCTAVE_API
-FloatComplexLU : public base_lu <FloatComplexMatrix, FloatComplex, Matrix, double>
+FloatComplexLU : public base_lu <FloatComplexMatrix>
 {
 public:
 
   FloatComplexLU (void)
-    : base_lu <FloatComplexMatrix, FloatComplex, Matrix, double> () { }
+    : base_lu <FloatComplexMatrix> () { }
 
   FloatComplexLU (const FloatComplexMatrix& a);
 
   FloatComplexLU (const FloatComplexLU& a)
-    : base_lu <FloatComplexMatrix, FloatComplex, Matrix, double> (a) { }
+    : base_lu <FloatComplexMatrix> (a) { }
 
   FloatComplexLU& operator = (const FloatComplexLU& a)
     {
       if (this != &a)
-	base_lu <FloatComplexMatrix, FloatComplex, Matrix, double> :: operator = (a);
+	base_lu <FloatComplexMatrix> :: operator = (a);
 
       return *this;
     }
diff --git a/liboctave/fCmplxQRP.cc b/liboctave/fCmplxQRP.cc
--- a/liboctave/fCmplxQRP.cc
+++ b/liboctave/fCmplxQRP.cc
@@ -86,7 +86,7 @@
   Array<float> rwork (2*n);
   float *prwork = rwork.fortran_vec ();
 
-  Array<octave_idx_type> jpvt (n, 0);
+  MArray<octave_idx_type> jpvt (n, 0);
   octave_idx_type *pjpvt = jpvt.fortran_vec ();
 
   // Code to enforce a certain permutation could go here...
@@ -97,18 +97,8 @@
   // 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;
-    }
+  jpvt -= 1;
+  p = PermMatrix (jpvt, true);
 
   octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
 
@@ -131,6 +121,14 @@
   q.resize (m, n2);
 }
 
+FloatColumnVector
+FloatComplexQRP::Pvec (void) const
+{
+  Array<float> pa (p);
+  FloatColumnVector pv (MArray<float> (pa) + 1.0f);
+  return pv;
+}
+
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/fCmplxQRP.h b/liboctave/fCmplxQRP.h
--- a/liboctave/fCmplxQRP.h
+++ b/liboctave/fCmplxQRP.h
@@ -27,7 +27,8 @@
 #include <iostream>
 
 #include "fCmplxQR.h"
-#include "fMatrix.h"
+#include "PermMatrix.h"
+#include "fColVector.h"
 
 class
 OCTAVE_API
@@ -55,13 +56,15 @@
 
   void init (const FloatComplexMatrix&, QR::type = QR::std);
 
-  FloatMatrix P (void) const { return p; }
+  PermMatrix P (void) const { return p; }
+
+  FloatColumnVector Pvec (void) const;
 
   friend std::ostream&  operator << (std::ostream&, const FloatComplexQRP&);
 
 private:
 
-  FloatMatrix p;
+  PermMatrix p;
 };
 
 #endif
diff --git a/liboctave/fDiagMatrix.h b/liboctave/fDiagMatrix.h
--- a/liboctave/fDiagMatrix.h
+++ b/liboctave/fDiagMatrix.h
@@ -98,6 +98,8 @@
 
   FloatColumnVector diag (octave_idx_type k = 0) const;
 
+  bool is_identity (void) const;
+
   // i/o
 
   friend OCTAVE_API std::ostream& operator << (std::ostream& os, const FloatDiagMatrix& a);
diff --git a/liboctave/fMatrix.cc b/liboctave/fMatrix.cc
--- a/liboctave/fMatrix.cc
+++ b/liboctave/fMatrix.cc
@@ -247,6 +247,19 @@
     elem (i, i) = a.elem (i, i);
 }
 
+FloatMatrix::FloatMatrix (const PermMatrix& a)
+  : MArray2<float> (a.rows (), a.cols (), 0.0)
+{
+  const Array<octave_idx_type> ia (a);
+  octave_idx_type len = a.rows ();
+  if (a.is_col_perm ())
+    for (octave_idx_type i = 0; i < len; i++)
+      elem (ia(i), i) = 1.0;
+  else
+    for (octave_idx_type i = 0; i < len; i++)
+      elem (i, ia(i)) = 1.0;
+}
+
 // FIXME -- could we use a templated mixed-type copy function
 // here?
 
diff --git a/liboctave/fMatrix.h b/liboctave/fMatrix.h
--- a/liboctave/fMatrix.h
+++ b/liboctave/fMatrix.h
@@ -64,10 +64,13 @@
 
   explicit FloatMatrix (const FloatDiagMatrix& a);
 
+  explicit FloatMatrix (const PermMatrix& a);
+
   explicit FloatMatrix (const boolMatrix& a);
 
   explicit FloatMatrix (const charMatrix& a);
 
+
   FloatMatrix& operator = (const FloatMatrix& a)
     {
       MArray2<float>::operator = (a);
diff --git a/liboctave/floatLU.cc b/liboctave/floatLU.cc
--- a/liboctave/floatLU.cc
+++ b/liboctave/floatLU.cc
@@ -34,7 +34,7 @@
 #include <base-lu.h>
 #include <base-lu.cc>
 
-template class base_lu <FloatMatrix, float, Matrix, double>;
+template class base_lu <FloatMatrix>;
 
 // Define the constructor for this particular derivation.
 
diff --git a/liboctave/floatLU.h b/liboctave/floatLU.h
--- a/liboctave/floatLU.h
+++ b/liboctave/floatLU.h
@@ -30,20 +30,20 @@
 
 class
 OCTAVE_API
-FloatLU : public base_lu <FloatMatrix, float, Matrix, double>
+FloatLU : public base_lu <FloatMatrix>
 {
 public:
 
-  FloatLU (void) : base_lu <FloatMatrix, float, Matrix, double> () { }
+  FloatLU (void) : base_lu <FloatMatrix> () { }
 
   FloatLU (const FloatMatrix& a);
 
-  FloatLU (const FloatLU& a) : base_lu <FloatMatrix, float, Matrix, double> (a) { }
+  FloatLU (const FloatLU& a) : base_lu <FloatMatrix> (a) { }
 
   FloatLU& operator = (const FloatLU& a)
     {
       if (this != &a)
-	base_lu <FloatMatrix, float, Matrix, double> :: operator = (a);
+	base_lu <FloatMatrix> :: operator = (a);
 
       return *this;
     }
diff --git a/liboctave/floatQRP.cc b/liboctave/floatQRP.cc
--- a/liboctave/floatQRP.cc
+++ b/liboctave/floatQRP.cc
@@ -81,7 +81,7 @@
 
   float *tmp_data = A_fact.fortran_vec ();
 
-  Array<octave_idx_type> jpvt (n, 0);
+  MArray<octave_idx_type> jpvt (n, 0);
   octave_idx_type *pjpvt = jpvt.fortran_vec ();
 
   // Code to enforce a certain permutation could go here...
@@ -91,18 +91,8 @@
   // 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;
-    }
+  jpvt -= 1;
+  p = PermMatrix (jpvt, true);
 
   octave_idx_type n2 = (qr_type == QR::economy) ? min_mn : m;
 
@@ -125,6 +115,13 @@
   q.resize (m, n2);
 }
 
+FloatColumnVector
+FloatQRP::Pvec (void) const
+{
+  Array<float> pa (p);
+  FloatColumnVector pv (MArray<float> (pa) + 1.0f);
+  return pv;
+}
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***
diff --git a/liboctave/floatQRP.h b/liboctave/floatQRP.h
--- a/liboctave/floatQRP.h
+++ b/liboctave/floatQRP.h
@@ -27,7 +27,8 @@
 #include <iostream>
 
 #include "floatQR.h"
-#include "fMatrix.h"
+#include "PermMatrix.h"
+#include "fColVector.h"
 
 class
 OCTAVE_API
@@ -56,13 +57,15 @@
 
   void init (const FloatMatrix&, QR::type = QR::std);
 
-  FloatMatrix P (void) const { return p; }
+  PermMatrix P (void) const { return p; }
+
+  FloatColumnVector Pvec (void) const;
 
   friend std::ostream&  operator << (std::ostream&, const FloatQRP&);
 
 protected:
 
-  FloatMatrix p;
+  PermMatrix p;
 };
 
 #endif
diff --git a/liboctave/idx-vector.cc b/liboctave/idx-vector.cc
--- a/liboctave/idx-vector.cc
+++ b/liboctave/idx-vector.cc
@@ -528,6 +528,39 @@
                              dim_vector (1, len), DIRECT);
 }
 
+bool
+idx_vector::is_permutation (octave_idx_type n) const
+{
+  bool retval = false;
+
+  if (is_colon_equiv (n))
+    retval = true;
+  else if (length (n) == n && extent(n) == n)
+    {
+      bool *left = new bool[n];
+
+      std::fill (left, left + n, true);
+
+      retval = true;
+
+      for (octave_idx_type i = 0, len = length (); i < len; i++)
+        { 
+          octave_idx_type k = xelem (i);
+          if (left[k])
+              left[k] = false;
+          else
+            {
+              retval = false;
+              break;
+            }
+        }
+
+      delete left;
+    }
+
+  return retval;
+}
+
 octave_idx_type 
 idx_vector::freeze (octave_idx_type z_len, const char *tag, bool resize_ok)
 {
diff --git a/liboctave/idx-vector.h b/liboctave/idx-vector.h
--- a/liboctave/idx-vector.h
+++ b/liboctave/idx-vector.h
@@ -719,6 +719,8 @@
   idx_vector
   complement (octave_idx_type n) const;
 
+  bool is_permutation (octave_idx_type n) const;
+
   // FIXME -- these are here for compatibility.  They should be removed
   // when no longer in use.
 
diff --git a/liboctave/mx-base.h b/liboctave/mx-base.h
--- a/liboctave/mx-base.h
+++ b/liboctave/mx-base.h
@@ -58,6 +58,9 @@
 #include "fDiagMatrix.h"
 #include "fCDiagMatrix.h"
 
+// Permutation matrix class
+#include "PermMatrix.h"
+
 // Sparse Matrix classes.
 
 #include "boolSparse.h"
diff --git a/liboctave/mx-defs.h b/liboctave/mx-defs.h
--- a/liboctave/mx-defs.h
+++ b/liboctave/mx-defs.h
@@ -55,6 +55,8 @@
 class FloatDiagMatrix;
 class FloatComplexDiagMatrix;
 
+class PermMatrix;
+
 class AEPBALANCE;
 class ComplexAEPBALANCE;
 class FloatAEPBALANCE;
diff --git a/liboctave/mx-op-defs.h b/liboctave/mx-op-defs.h
--- a/liboctave/mx-op-defs.h
+++ b/liboctave/mx-op-defs.h
@@ -1289,6 +1289,60 @@
   extern OCTAVE_API T ## NDArray max (const T ## NDArray& a, \
 				       const T ## NDArray& b);
 
+#define PMM_MULTIPLY_OP(PM, M) \
+M operator * (const PM& p, const M& x) \
+{ \
+  octave_idx_type nr = x.rows (), nc = x.columns (); \
+  M result; \
+  if (p.columns () != nr) \
+    gripe_nonconformant ("operator *", p.rows (), p.columns (), nr, nc); \
+  else \
+    { \
+      if (p.is_col_perm ()) \
+        { \
+          result = M (nr, nc); \
+          result.assign (idx_vector (p), idx_vector::colon, x); \
+        } \
+      else \
+        result = x.index (idx_vector (p), idx_vector::colon); \
+    } \
+  \
+  return result; \
+}
+
+#define MPM_MULTIPLY_OP(M, PM) \
+M operator * (const M& x, const PM& p) \
+{ \
+  octave_idx_type nr = x.rows (), nc = x.columns (); \
+  M result; \
+  if (p.rows () != nc) \
+    gripe_nonconformant ("operator *", nr, nc, p.rows (), p.columns ()); \
+  else \
+    { \
+      if (p.is_col_perm ()) \
+        result = x.index (idx_vector::colon, idx_vector (p)); \
+      else \
+        { \
+          result = M (nr, nc); \
+          result.assign (idx_vector::colon, idx_vector (p), x); \
+        } \
+    } \
+  \
+  return result; \
+}
+
+#define PMM_BIN_OP_DECLS(R, PM, M, API) \
+  BIN_OP_DECL (R, operator *, PM, M, API);
+
+#define PMM_BIN_OPS(R, PM, M) \
+  PMM_MULTIPLY_OP(PM, M);
+
+#define MPM_BIN_OP_DECLS(R, M, PM, API) \
+  BIN_OP_DECL (R, operator *, M, PM, API);
+
+#define MPM_BIN_OPS(R, M, PM) \
+  MPM_MULTIPLY_OP(M, PM);
+
 
 /*
 ;;; Local Variables: ***
diff --git a/liboctave/mx-ops b/liboctave/mx-ops
--- a/liboctave/mx-ops
+++ b/liboctave/mx-ops
@@ -70,6 +70,7 @@
 ui32nda uint32NDArray ND uint32NDArray.h YES octave_uint32::zero uint32_t
 i64nda int64NDArray ND int64NDArray.h YES octave_int64::zero int64_t
 ui64nda uint64NDArray ND uint64NDArray.h YES octave_uint64::zero uint64_t
+pm PermMatrix PM PermMatrix.h YES static_cast<octave_idx_type>(0)
 # ops
 # result_t lhs_t rhs_t op-type lhs_conv rhs_conv headers ...
 #
@@ -141,6 +142,15 @@
 fm fm fdm B
 fm fs fdm B
 #
+m pm m B
+m m pm B
+cm pm cm B
+cm cm pm B
+fm pm fm B
+fm fm pm B
+fcm pm fcm B
+fcm fcm pm B
+#
 i8nda s i8nda BCL NONE NONE boolMatrix.h boolNDArray.h
 i8nda i8nda s BCL NONE NONE boolMatrix.h boolNDArray.h
 ui8nda s ui8nda BCL NONE NONE boolMatrix.h boolNDArray.h
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,29 @@
+2008-12-03  Jaroslav Hajek  <highegg@gmail.com>
+
+	* ov-perm.h: New source.
+	* ov-perm.cc: New source.
+	* ov-flt-perm.h: New source.
+	* ov-flt-perm.cc: New source.
+	* ov-base-diag.cc (octave_base_diag<DMT, MT>::do_index_op):
+	If subscripting an identity matrix by permutation(s), return a
+	permutation matrix object.
+	* ov.cc (octave_value::octave_value (const PermMatrix&)): New
+	constructor.
+	* ov.cc (octave_value::octave_value (const PermMatrix&)): Declare it.
+	* op-cm-pm.cc: New source.
+	* op-fcm-fpm.cc: New source.
+	* op-fm-fpm.cc: New source.
+	* op-fpm-fcm.cc: New source.
+	* op-fpm-fm.cc: New source.
+	* op-fpm-fpm.cc: New source.
+	* op-m-pm.cc: New source.
+	* op-pm-cm.cc: New source.
+	* op-pm-m.cc: New source.
+	* op-pm-pm.cc: New source.
+	* op-pm-template.cc: New source.
+	* Makefile.in: Include new sources.
+	* DLD-FUNCTIONS/qr.cc (Fqr): Reflect interface changes of QR classes.
+
 2008-12-03  John W. Eaton  <jwe@octave.org>
 
 	* debug.cc (bp_table::do_get_breakpoint_list): Style fixes.
diff --git a/src/DLD-FUNCTIONS/lu.cc b/src/DLD-FUNCTIONS/lu.cc
--- a/src/DLD-FUNCTIONS/lu.cc
+++ b/src/DLD-FUNCTIONS/lu.cc
@@ -356,7 +356,7 @@
 
 		    case 2:
 		      {
-			FloatMatrix P = fact.P ();
+			PermMatrix P = fact.P ();
 			FloatMatrix L = P.transpose () * fact.L ();
 			retval(1) = fact.U ();
 			retval(0) = L;
@@ -394,7 +394,7 @@
 
 		    case 2:
 		      {
-			Matrix P = fact.P ();
+			PermMatrix P = fact.P ();
 			Matrix L = P.transpose () * fact.L ();
 			retval(1) = fact.U ();
 			retval(0) = L;
@@ -435,7 +435,7 @@
 
 		    case 2:
 		      {
-			FloatMatrix P = fact.P ();
+			PermMatrix P = fact.P ();
 			FloatComplexMatrix L = P.transpose () * fact.L ();
 			retval(1) = fact.U ();
 			retval(0) = L;
@@ -473,7 +473,7 @@
 
 		    case 2:
 		      {
-			Matrix P = fact.P ();
+			PermMatrix P = fact.P ();
 			ComplexMatrix L = P.transpose () * fact.L ();
 			retval(1) = fact.U ();
 			retval(0) = L;
@@ -515,7 +515,7 @@
 %! [l, u, p] = lu ([1, 2; 3, 4]);
 %! assert(l, [1, 0; 1/3, 1], sqrt (eps));
 %! assert(u, [3, 4; 0, 2/3], sqrt (eps));
-%! assert(p, [0, 1; 1, 0], sqrt (eps));
+%! assert(p(:,:), [0, 1; 1, 0], sqrt (eps));
 
 %!test
 %! [l, u, p] = lu ([1, 2; 3, 4],'vector');
@@ -527,7 +527,7 @@
 %! [l u p] = lu ([1, 2; 3, 4; 5, 6]);
 %! assert(l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
 %! assert(u, [5, 6; 0, 4/5], sqrt (eps));
-%! assert(p, [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
+%! assert(p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
 
 %!assert(lu (single([1, 2; 3, 4])), single([3, 4; 1/3, 2/3]), eps('single'));
 
@@ -540,7 +540,7 @@
 %! [l, u, p] = lu (single([1, 2; 3, 4]));
 %! assert(l, single([1, 0; 1/3, 1]), sqrt (eps('single')));
 %! assert(u, single([3, 4; 0, 2/3]), sqrt (eps('single')));
-%! assert(p, single([0, 1; 1, 0]), sqrt (eps('single')));
+%! assert(p(:,:), single([0, 1; 1, 0]), sqrt (eps('single')));
 
 %!test
 %! [l, u, p] = lu (single([1, 2; 3, 4]),'vector');
@@ -552,7 +552,7 @@
 %! [l u p] = lu (single([1, 2; 3, 4; 5, 6]));
 %! assert(l, single([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps('single')));
 %! assert(u, single([5, 6; 0, 4/5]), sqrt (eps('single')));
-%! assert(p, single([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps('single')));
+%! assert(p(:,:), single([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps('single')));
 
 %!error <Invalid call to lu.*> lu ();
 %!error lu ([1, 2; 3, 4], 2);
diff --git a/src/DLD-FUNCTIONS/qr.cc b/src/DLD-FUNCTIONS/qr.cc
--- a/src/DLD-FUNCTIONS/qr.cc
+++ b/src/DLD-FUNCTIONS/qr.cc
@@ -308,7 +308,10 @@
 		    default:
 		      {
 			FloatQRP fact (m, type);
-			retval(2) = fact.P ();
+                        if (type == QR::economy)
+                          retval(2) = fact.Pvec ();
+                        else
+                          retval(2) = fact.P ();
 			retval(1) = fact.R ();
 			retval(0) = fact.Q ();
 		      }
@@ -343,7 +346,10 @@
 		    default:
 		      {
 			FloatComplexQRP fact (m, type);
-			retval(2) = fact.P ();
+                        if (type == QR::economy)
+                          retval(2) = fact.Pvec ();
+                        else
+                          retval(2) = fact.P ();
 			retval(1) = fact.R ();
 			retval(0) = fact.Q ();
 		      }
@@ -381,7 +387,10 @@
 		    default:
 		      {
 			QRP fact (m, type);
-			retval(2) = fact.P ();
+                        if (type == QR::economy)
+                          retval(2) = fact.Pvec ();
+                        else
+                          retval(2) = fact.P ();
 			retval(1) = fact.R ();
 			retval(0) = fact.Q ();
 		      }
@@ -416,7 +425,10 @@
 		    default:
 		      {
 			ComplexQRP fact (m, type);
-			retval(2) = fact.P ();
+                        if (type == QR::economy)
+                          retval(2) = fact.Pvec ();
+                        else
+                          retval(2) = fact.P ();
 			retval(1) = fact.R ();
 			retval(0) = fact.Q ();
 		      }
diff --git a/src/Makefile.in b/src/Makefile.in
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -102,6 +102,7 @@
 	ov-colon.h ov-base.h ov-base-mat.h ov-base-scalar.h \
 	ov-str-mat.h ov-bool-mat.h ov-null-mat.h ov-bool.h \
 	ov-base-diag.h ov-re-diag.h ov-flt-re-diag.h ov-cx-diag.h ov-flt-cx-diag.h \
+	ov-perm.h ov-flt-perm.h \
 	ov-cell.h ov.h ov-fcn.h ov-builtin.h ov-dld-fcn.h \
 	ov-mex-fcn.h ov-usr-fcn.h ov-fcn-handle.h \
 	ov-fcn-inline.h ov-class.h ov-typeinfo.h ov-type-conv.h \
@@ -171,16 +172,19 @@
 	op-fdm-fcdm.cc op-fdm-fcm.cc op-fdm-fcs.cc op-fdm-fdm.cc \
 	op-fdm-fm.cc op-fdm-fs.cc op-fm-fcdm.cc op-fm-fdm.cc
 
+PERM_OP_XSRC := op-cm-pm.cc op-fcm-fpm.cc op-fm-fpm.cc op-fpm-fcm.cc \
+	op-fpm-fm.cc op-fpm-fpm.cc op-m-pm.cc op-pm-cm.cc op-pm-m.cc op-pm-pm.cc
+
 OP_XSRC :=  op-b-b.cc op-b-bm.cc op-bm-b.cc op-bm-bm.cc op-cell.cc \
 	op-chm.cc op-class.cc op-list.cc op-range.cc op-str-m.cc \
 	op-str-s.cc op-str-str.cc op-struct.cc \
 	$(DOUBLE_OP_XSRC) $(FLOAT_OP_XSRC) $(INTTYPE_OP_XSRC) \
-	$(SPARSE_OP_XSRC) $(DIAG_OP_XSRC) $(FDIAG_OP_XSRC)
+	$(SPARSE_OP_XSRC) $(DIAG_OP_XSRC) $(FDIAG_OP_XSRC) $(PERM_OP_XSRC)
 
 OP_SRC := $(addprefix OPERATORS/, $(OP_XSRC))
 
 OP_INCLUDES := OPERATORS/op-int.h \
-	OPERATORS/op-dm-template.cc OPERATORS/op-dms-template.cc
+	OPERATORS/op-dm-template.cc OPERATORS/op-dms-template.cc OPERATORS/op-pm-template.cc
 
 OV_INTTYPE_SRC := \
 	ov-int8.cc ov-int16.cc ov-int32.cc ov-int64.cc \
@@ -199,6 +203,7 @@
 	ov-class.cc ov-typeinfo.cc \
 	ov-flt-re-mat.cc ov-flt-cx-mat.cc ov-float.cc ov-flt-complex.cc \
 	ov-re-diag.cc ov-flt-re-diag.cc ov-cx-diag.cc ov-flt-cx-diag.cc \
+	ov-perm.cc ov-flt-perm.cc \
 	$(OV_INTTYPE_SRC) \
 	$(OV_SPARSE_SRC)
 
diff --git a/src/OPERATORS/op-cm-pm.cc b/src/OPERATORS/op-cm-pm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-cm-pm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-cx-mat.h"
+
+#define LMATRIX complex_matrix
+#define RMATRIX perm_matrix
+
+#define LSHORT cm
+#define RSHORT pm
+
+#define RIGHT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-fcm-fpm.cc b/src/OPERATORS/op-fcm-fpm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-fcm-fpm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-flt-cx-mat.h"
+
+#define LMATRIX float_complex_matrix
+#define RMATRIX float_perm_matrix
+
+#define LSHORT fcm
+#define RSHORT fpm
+
+#define RIGHT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-fm-fpm.cc b/src/OPERATORS/op-fm-fpm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-fm-fpm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-flt-re-mat.h"
+
+#define LMATRIX float_matrix
+#define RMATRIX float_perm_matrix
+
+#define LSHORT fm
+#define RSHORT fpm
+
+#define RIGHT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-fpm-fcm.cc b/src/OPERATORS/op-fpm-fcm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-fpm-fcm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-flt-cx-mat.h"
+
+#define LMATRIX float_perm_matrix
+#define RMATRIX float_complex_matrix
+
+#define LSHORT fpm
+#define RSHORT fcm
+
+#define LEFT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-fpm-fm.cc b/src/OPERATORS/op-fpm-fm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-fpm-fm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-flt-re-mat.h"
+
+#define LMATRIX float_perm_matrix
+#define RMATRIX float_matrix
+
+#define LSHORT fpm
+#define RSHORT fm
+
+#define LEFT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-fpm-fpm.cc b/src/OPERATORS/op-fpm-fpm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-fpm-fpm.cc
@@ -0,0 +1,91 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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 "gripes.h"
+#include "oct-obj.h"
+#include "ov.h"
+#include "ov-perm.h"
+#include "ov-flt-perm.h"
+#include "ov-flt-re-mat.h"
+#include "ov-typeinfo.h"
+#include "ops.h"
+
+DEFUNOP (transpose, float_perm_matrix)
+{
+  CAST_UNOP_ARG (const octave_float_perm_matrix&);
+  return octave_value (v.perm_matrix_value().transpose (), true);
+}
+
+DEFBINOP (mul, float_perm_matrix, float_perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_float_perm_matrix&, const octave_float_perm_matrix&);
+
+  return octave_value (v1.perm_matrix_value () * v2.perm_matrix_value (), true);
+}
+
+DEFBINOP (div, float_perm_matrix, float_perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_float_perm_matrix&, const octave_float_perm_matrix&);
+  
+  return octave_value (v1.perm_matrix_value () * v2.perm_matrix_value ().inverse (), false);
+}
+
+DEFBINOP (ldiv, float_perm_matrix, float_perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_float_perm_matrix&, const octave_float_perm_matrix&);
+  
+  return octave_value (v1.perm_matrix_value ().inverse () * v2.perm_matrix_value (), false);
+}
+
+CONVDECL (float_perm_matrix_to_float_matrix)
+{
+  CAST_CONV_ARG (const octave_float_perm_matrix&);
+
+  return new octave_float_matrix (v.float_matrix_value ());
+}
+
+CONVDECL (float_perm_matrix_to_perm_matrix)
+{
+  CAST_CONV_ARG (const octave_float_perm_matrix&);
+
+  return new octave_perm_matrix (v.perm_matrix_value ());
+}
+
+void
+install_fpm_fpm_ops (void)
+{
+  INSTALL_UNOP (op_transpose, octave_float_perm_matrix, transpose);
+  INSTALL_UNOP (op_hermitian, octave_float_perm_matrix, transpose);
+
+  INSTALL_BINOP (op_mul, octave_float_perm_matrix, octave_float_perm_matrix, mul);
+  INSTALL_BINOP (op_div, octave_float_perm_matrix, octave_float_perm_matrix, div);
+  INSTALL_BINOP (op_ldiv, octave_float_perm_matrix, octave_float_perm_matrix, ldiv);
+
+  INSTALL_CONVOP (octave_float_perm_matrix, octave_float_matrix, float_perm_matrix_to_float_matrix);
+  INSTALL_CONVOP (octave_float_perm_matrix, octave_perm_matrix, float_perm_matrix_to_perm_matrix);
+  INSTALL_ASSIGNCONV (octave_float_perm_matrix, octave_float_matrix, octave_float_matrix);
+}
+
diff --git a/src/OPERATORS/op-m-pm.cc b/src/OPERATORS/op-m-pm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-m-pm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-re-mat.h"
+
+#define LMATRIX matrix
+#define RMATRIX perm_matrix
+
+#define LSHORT m
+#define RSHORT pm
+
+#define RIGHT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-pm-cm.cc b/src/OPERATORS/op-pm-cm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-pm-cm.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-cx-mat.h"
+
+#define LMATRIX perm_matrix
+#define RMATRIX complex_matrix
+
+#define LSHORT pm
+#define RSHORT cm
+
+#define LEFT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-pm-m.cc b/src/OPERATORS/op-pm-m.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-pm-m.cc
@@ -0,0 +1,33 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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/>.
+
+*/
+
+#define MINCLUDE "ov-re-mat.h"
+
+#define LMATRIX perm_matrix
+#define RMATRIX matrix
+
+#define LSHORT pm
+#define RSHORT m
+
+#define LEFT
+
+#include "op-pm-template.cc"
diff --git a/src/OPERATORS/op-pm-pm.cc b/src/OPERATORS/op-pm-pm.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-pm-pm.cc
@@ -0,0 +1,85 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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 "gripes.h"
+#include "oct-obj.h"
+#include "ov.h"
+#include "ov-perm.h"
+#include "ov-flt-perm.h"
+#include "ov-re-mat.h"
+#include "ov-typeinfo.h"
+#include "ops.h"
+
+DEFUNOP (transpose, perm_matrix)
+{
+  CAST_UNOP_ARG (const octave_perm_matrix&);
+  return octave_value (v.perm_matrix_value().transpose ());
+}
+
+DEFBINOP_OP (mul, perm_matrix, perm_matrix, *)
+
+DEFBINOP (div, perm_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_perm_matrix&);
+  
+  return (v1.perm_matrix_value () * v2.perm_matrix_value ().inverse ());
+}
+
+DEFBINOP (ldiv, perm_matrix, perm_matrix)
+{
+  CAST_BINOP_ARGS (const octave_perm_matrix&, const octave_perm_matrix&);
+  
+  return (v1.perm_matrix_value ().inverse () * v2.perm_matrix_value ());
+}
+
+CONVDECL (perm_matrix_to_matrix)
+{
+  CAST_CONV_ARG (const octave_perm_matrix&);
+
+  return new octave_matrix (v.matrix_value ());
+}
+
+CONVDECL (perm_matrix_to_float_perm_matrix)
+{
+  CAST_CONV_ARG (const octave_perm_matrix&);
+
+  return new octave_float_perm_matrix (v.perm_matrix_value ());
+}
+
+void
+install_pm_pm_ops (void)
+{
+  INSTALL_UNOP (op_transpose, octave_perm_matrix, transpose);
+  INSTALL_UNOP (op_hermitian, octave_perm_matrix, transpose);
+
+  INSTALL_BINOP (op_mul, octave_perm_matrix, octave_perm_matrix, mul);
+  INSTALL_BINOP (op_div, octave_perm_matrix, octave_perm_matrix, div);
+  INSTALL_BINOP (op_ldiv, octave_perm_matrix, octave_perm_matrix, ldiv);
+
+  INSTALL_CONVOP (octave_perm_matrix, octave_matrix, perm_matrix_to_matrix);
+  INSTALL_CONVOP (octave_perm_matrix, octave_float_perm_matrix, perm_matrix_to_float_perm_matrix);
+  INSTALL_ASSIGNCONV (octave_perm_matrix, octave_matrix, octave_matrix);
+}
diff --git a/src/OPERATORS/op-pm-template.cc b/src/OPERATORS/op-pm-template.cc
new file mode 100644
--- /dev/null
+++ b/src/OPERATORS/op-pm-template.cc
@@ -0,0 +1,78 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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 "ov-perm.h"
+#include "ov-flt-perm.h"
+#include MINCLUDE
+#include "ops.h"
+
+#define OCTAVE_LMATRIX CONCAT2(octave_, LMATRIX)
+#define OCTAVE_RMATRIX CONCAT2(octave_, RMATRIX)
+#ifdef LEFT
+#define LMATRIX_VALUE perm_matrix_value
+#define RMATRIX_VALUE CONCAT2(RMATRIX, _value)
+#else
+#define LMATRIX_VALUE CONCAT2(LMATRIX, _value)
+#define RMATRIX_VALUE perm_matrix_value
+#endif
+
+DEFBINOP (mul, LMATRIX, RMATRIX)
+{
+  CAST_BINOP_ARGS (const OCTAVE_LMATRIX&, const OCTAVE_RMATRIX&);
+
+  return v1.LMATRIX_VALUE () * v2.RMATRIX_VALUE ();
+}
+
+#ifdef LEFT
+DEFBINOP (ldiv, LMATRIX, RMATRIX)
+{
+  CAST_BINOP_ARGS (const OCTAVE_LMATRIX&, const OCTAVE_RMATRIX&);
+
+  return v1.perm_matrix_value ().inverse () * v2.RMATRIX_VALUE (); 
+}
+#else
+DEFBINOP (div, LMATRIX, RMATRIX)
+{
+  CAST_BINOP_ARGS (const OCTAVE_LMATRIX&, const OCTAVE_RMATRIX&);
+
+  return v1.LMATRIX_VALUE () * v2.perm_matrix_value ().inverse ();
+}
+#endif
+
+
+#define SHORT_NAME CONCAT3(LSHORT, _, RSHORT)
+#define INST_NAME CONCAT3(install_, SHORT_NAME, _ops)
+
+void
+INST_NAME (void)
+{
+  INSTALL_BINOP (op_mul, OCTAVE_LMATRIX, OCTAVE_RMATRIX, mul);
+#ifdef LEFT
+  INSTALL_BINOP (op_ldiv, OCTAVE_LMATRIX, OCTAVE_RMATRIX, ldiv);
+#else
+  INSTALL_BINOP (op_div, OCTAVE_LMATRIX, OCTAVE_RMATRIX, div);
+#endif
+}
diff --git a/src/ov-base-diag.cc b/src/ov-base-diag.cc
--- a/src/ov-base-diag.cc
+++ b/src/ov-base-diag.cc
@@ -70,18 +70,58 @@
                                         bool resize_ok)
 {
   octave_value retval;
+  typedef typename DMT::element_type el_type;
+  el_type one = 1;
 
-  if (idx.length () == 2 && idx(0).is_scalar_type () 
-      && idx(1).is_scalar_type () && ! resize_ok)
+  octave_idx_type nidx = idx.length ();
+  idx_vector idx0, idx1;
+  if (nidx == 2)
+    {
+      idx0 = idx(0).index_vector ();
+      idx1 = idx(1).index_vector ();
+    }
+
+  // This hack is to allow constructing permutation matrices using
+  // eye(n)(p,:), eye(n)(:,q) && eye(n)(p,q) where p & q are permutation
+  // vectors. 
+  // Note that, for better consistency, eye(n)(:,:) still converts to a full
+  // matrix.
+  // FIXME: This check is probably unnecessary for complex matrices. 
+  if (! error_state && nidx == 2 && matrix.is_multiple_of_identity (one))
     {
-      idx_vector i = idx(0).index_vector (), j = idx(1).index_vector ();
-      // FIXME: the proxy mechanism of DiagArray2 causes problems here.
-      typedef typename DMT::element_type el_type;
-      if (! error_state)
-         retval = el_type (matrix.checkelem (i(0), j(0)));
+      
+      bool left = idx0.is_permutation (matrix.rows ());
+      bool right = idx1.is_permutation (matrix.cols ());
+
+      if (left && right)
+        {
+          if (idx0.is_colon ()) left = false;
+          if (idx1.is_colon ()) right = false;
+          if (left && right)
+              retval = octave_value (PermMatrix (idx0, false) * PermMatrix (idx1, true),
+                                     is_single_type ());
+          else if (left)
+              retval = octave_value (PermMatrix (idx0, false),
+                                     is_single_type ());
+          else if (right)
+              retval = octave_value (PermMatrix (idx1, true),
+                                     is_single_type ());
+        }
     }
-  else
-    retval = to_dense ().do_index_op (idx, resize_ok);
+
+  // if error_state is set, we've already griped.
+  if (! error_state && ! retval.is_defined ())
+    {
+      if (nidx == 2 && ! resize_ok &&
+          idx0.is_scalar () && idx1.is_scalar ())
+        {
+          // FIXME: the proxy mechanism of DiagArray2 causes problems here.
+          retval = el_type (matrix.checkelem (idx0(0), idx1(0)));
+        }
+      else
+        retval = to_dense ().do_index_op (idx, resize_ok);
+    }
+
   return retval;
 }
 
diff --git a/src/ov-flt-perm.cc b/src/ov-flt-perm.cc
new file mode 100644
--- /dev/null
+++ b/src/ov-flt-perm.cc
@@ -0,0 +1,77 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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 "ov-flt-perm.h"
+#include "ov-flt-re-mat.h"
+#include "ov-float.h"
+#include "ops.h"
+
+DEFINE_OCTAVE_ALLOCATOR (octave_float_perm_matrix);
+
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_float_perm_matrix, 
+                                     "float permutation matrix", "single");
+
+static octave_base_value *
+default_float_numeric_conversion_function (const octave_base_value& a)
+{
+  CAST_CONV_ARG (const octave_float_perm_matrix&);
+
+  return new octave_float_matrix (v.float_matrix_value ());
+}
+
+octave_base_value::type_conv_info
+octave_float_perm_matrix::numeric_conversion_function (void) const
+{
+  return octave_base_value::type_conv_info (default_float_numeric_conversion_function,
+                                            octave_float_matrix::static_type_id ());
+}
+
+octave_base_value::type_conv_info
+octave_float_perm_matrix::numeric_demotion_function (void) const
+{
+  return octave_base_value::type_conv_info (0);
+}
+
+octave_base_value *
+octave_float_perm_matrix::try_narrowing_conversion (void)
+{
+  octave_base_value *retval = 0;
+
+  if (matrix.nelem () == 1)
+    retval = new octave_float_scalar (matrix (0, 0));
+
+  return retval;
+}
+
+octave_value
+octave_float_perm_matrix::to_dense (void) const
+{
+  if (! dense_cache.is_defined ())
+      dense_cache = FloatMatrix (matrix);
+
+  return dense_cache;
+}
+
diff --git a/src/ov-flt-perm.h b/src/ov-flt-perm.h
new file mode 100644
--- /dev/null
+++ b/src/ov-flt-perm.h
@@ -0,0 +1,62 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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_float_perm_matrix_h)
+#define octave_float_perm_matrix_h 1
+
+#include "ov-perm.h"
+#include "ov-perm.h"
+
+class 
+OCTINTERP_API
+octave_float_perm_matrix : public octave_perm_matrix
+{
+public:
+  octave_float_perm_matrix (void) : octave_perm_matrix () { }
+
+  octave_float_perm_matrix (const PermMatrix& p) : octave_perm_matrix (p) { }
+
+  octave_base_value *clone (void) const { return new octave_float_perm_matrix (*this); }
+  octave_base_value *empty_clone (void) const { return new octave_float_perm_matrix (); }
+
+  bool is_double_type (void) const { return false; }
+
+  bool is_single_type (void) const { return true; }
+
+  type_conv_info numeric_conversion_function (void) const;
+
+  type_conv_info numeric_demotion_function (void) const;
+
+  octave_base_value *try_narrowing_conversion (void);
+
+protected:
+
+  virtual octave_value to_dense (void) const;
+
+private:
+
+  DECLARE_OCTAVE_ALLOCATOR
+
+  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
+};
+
+#endif
diff --git a/src/ov-perm.cc b/src/ov-perm.cc
new file mode 100644
--- /dev/null
+++ b/src/ov-perm.cc
@@ -0,0 +1,406 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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 "ov-perm.h"
+#include "ov-flt-perm.h"
+#include "ov-re-mat.h"
+#include "ov-scalar.h"
+#include "error.h"
+#include "gripes.h"
+#include "ops.h"
+
+octave_value
+octave_perm_matrix::subsref (const std::string& type,
+                             const std::list<octave_value_list>& idx)
+{
+  octave_value retval;
+
+  switch (type[0])
+    {
+    case '(':
+      retval = do_index_op (idx.front ());
+      break;
+
+    case '{':
+    case '.':
+      {
+	std::string nm = type_name ();
+	error ("%s cannot be indexed with %c", nm.c_str (), type[0]);
+      }
+      break;
+
+    default:
+      panic_impossible ();
+    }
+
+  return retval.next_subsref (type, idx);
+}
+
+octave_value
+octave_perm_matrix::do_index_op (const octave_value_list& idx,
+                                 bool resize_ok)
+{
+  octave_value retval;
+  octave_idx_type nidx = idx.length ();
+  idx_vector idx0, idx1;
+  if (nidx == 2)
+    {
+      idx0 = idx(0).index_vector ();
+      idx1 = idx(1).index_vector ();
+    }
+
+  // This hack is to allow constructing permutation matrices using
+  // eye(n)(p,:), eye(n)(:,q) && eye(n)(p,q) where p & q are permutation
+  // vectors. 
+  // Note that, for better consistency, eye(n)(:,:) still converts to a full
+  // matrix.
+  if (! error_state && nidx == 2)
+    {
+      bool left = idx0.is_permutation (matrix.rows ());
+      bool right = idx1.is_permutation (matrix.cols ());
+
+      if (left && right)
+        {
+          if (idx0.is_colon ()) left = false;
+          if (idx1.is_colon ()) right = false;
+          if (left || right)
+            {
+              PermMatrix p = matrix;
+              if (left)
+                p = PermMatrix (idx0, false) * p;
+              if (right)
+                p = p * PermMatrix (idx1, true);
+              retval = octave_value (p, is_single_type ());
+            }
+        }
+    }
+
+  // if error_state is set, we've already griped.
+  if (! error_state && ! retval.is_defined ())
+    {
+      if (nidx == 2 && ! resize_ok &&
+          idx0.is_scalar () && idx1.is_scalar ())
+        {
+          retval = matrix.checkelem (idx0(0), idx1(0));
+        }
+      else
+        retval = to_dense ().do_index_op (idx, resize_ok);
+    }
+
+  return retval;
+}
+
+bool
+octave_perm_matrix::is_true (void) const
+{
+  return to_dense ().is_true ();
+}
+
+bool
+octave_perm_matrix::valid_as_scalar_index (void) const
+{
+  return false;
+}
+
+double
+octave_perm_matrix::double_value (bool) const
+{
+  double retval = lo_ieee_nan_value ();
+
+  if (numel () > 0)
+    {
+      gripe_implicit_conversion ("Octave:array-as-scalar",
+				 type_name (), "real scalar");
+
+      retval = matrix (0, 0);
+    }
+  else
+    gripe_invalid_conversion (type_name (), "real scalar");
+
+  return retval;
+}
+
+float
+octave_perm_matrix::float_value (bool) const
+{
+  float retval = lo_ieee_float_nan_value ();
+
+  if (numel () > 0)
+    {
+      gripe_implicit_conversion ("Octave:array-as-scalar",
+				 type_name (), "real scalar");
+
+      retval = matrix (0, 0);
+    }
+  else
+    gripe_invalid_conversion (type_name (), "real scalar");
+
+  return retval;
+}
+
+Complex
+octave_perm_matrix::complex_value (bool) const
+{
+  double tmp = lo_ieee_nan_value ();
+
+  Complex retval (tmp, tmp);
+
+  if (rows () > 0 && columns () > 0)
+    {
+      gripe_implicit_conversion ("Octave:array-as-scalar",
+				 type_name (), "complex scalar");
+
+      retval = matrix (0, 0);
+    }
+  else
+    gripe_invalid_conversion (type_name (), "complex scalar");
+
+  return retval;
+}
+
+FloatComplex
+octave_perm_matrix::float_complex_value (bool) const
+{
+  float tmp = lo_ieee_float_nan_value ();
+
+  FloatComplex retval (tmp, tmp);
+
+  if (rows () > 0 && columns () > 0)
+    {
+      gripe_implicit_conversion ("Octave:array-as-scalar",
+				 type_name (), "complex scalar");
+
+      retval = matrix (0, 0);
+    }
+  else
+    gripe_invalid_conversion (type_name (), "complex scalar");
+
+  return retval;
+}
+
+#define FORWARD_MATRIX_VALUE(TYPE, PREFIX) \
+TYPE \
+octave_perm_matrix::PREFIX ## _value (bool frc_str_conv) const \
+{ \
+  return to_dense ().PREFIX ## _value (frc_str_conv); \
+}
+
+FORWARD_MATRIX_VALUE (Matrix, matrix)
+FORWARD_MATRIX_VALUE (FloatMatrix, float_matrix)
+FORWARD_MATRIX_VALUE (ComplexMatrix, complex_matrix)
+FORWARD_MATRIX_VALUE (FloatComplexMatrix, float_complex_matrix)
+
+FORWARD_MATRIX_VALUE (NDArray, array)
+FORWARD_MATRIX_VALUE (FloatNDArray, float_array)
+FORWARD_MATRIX_VALUE (ComplexNDArray, complex_array)
+FORWARD_MATRIX_VALUE (FloatComplexNDArray, float_complex_array)
+
+FORWARD_MATRIX_VALUE (boolNDArray, bool_array)
+FORWARD_MATRIX_VALUE (charNDArray, char_array)
+
+FORWARD_MATRIX_VALUE (SparseMatrix, sparse_matrix)
+FORWARD_MATRIX_VALUE (SparseComplexMatrix, sparse_complex_matrix)
+
+idx_vector
+octave_perm_matrix::index_vector (void) const
+{
+  return to_dense ().index_vector ();
+}
+
+octave_value
+octave_perm_matrix::convert_to_str_internal (bool pad, bool force, char type) const
+{
+  return to_dense ().convert_to_str_internal (pad, force, type);
+}
+
+bool 
+octave_perm_matrix::save_ascii (std::ostream& os)
+{
+  // FIXME: this should probably save the matrix as permutation.
+  return to_dense ().save_ascii (os);
+}
+
+bool 
+octave_perm_matrix::save_binary (std::ostream& os, bool& save_as_floats)
+{
+  return to_dense ().save_binary (os, save_as_floats);
+}
+
+#if defined (HAVE_HDF5)
+
+bool
+octave_perm_matrix::save_hdf5 (hid_t loc_id, const char *name, bool save_as_floats)
+{
+  return to_dense ().save_hdf5 (loc_id, name, save_as_floats);
+}
+
+#endif
+
+void
+octave_perm_matrix::print_raw (std::ostream& os,
+			  bool pr_as_read_syntax) const
+{
+  return to_dense ().print_raw (os, pr_as_read_syntax);
+}
+
+mxArray *
+octave_perm_matrix::as_mxArray (void) const
+{
+  return to_dense ().as_mxArray ();
+}
+
+bool
+octave_perm_matrix::print_as_scalar (void) const
+{
+  dim_vector dv = dims ();
+
+  return (dv.all_ones () || dv.any_zero ());
+}
+
+void
+octave_perm_matrix::print (std::ostream& os, bool pr_as_read_syntax) const
+{
+  to_dense ().print (os, pr_as_read_syntax);
+}
+
+int
+octave_perm_matrix::write (octave_stream& os, int block_size,
+                                oct_data_conv::data_type output_type, int skip,
+                                oct_mach_info::float_format flt_fmt) const
+{ 
+  return to_dense ().write (os, block_size, output_type, skip, flt_fmt); 
+}
+
+void
+octave_perm_matrix::print_info (std::ostream& os,
+				    const std::string& prefix) const
+{
+  matrix.print_info (os, prefix);
+}
+
+
+octave_value
+octave_perm_matrix::to_dense (void) const
+{
+  if (! dense_cache.is_defined ())
+      dense_cache = Matrix (matrix);
+
+  return dense_cache;
+}
+
+#define FORWARD_MAPPER(MAP) \
+  octave_value \
+  octave_perm_matrix::MAP (void) const \
+  { \
+    return to_dense ().MAP (); \
+  }
+
+FORWARD_MAPPER (erf)
+FORWARD_MAPPER (erfc)
+FORWARD_MAPPER (gamma)
+FORWARD_MAPPER (lgamma)
+FORWARD_MAPPER (abs)
+FORWARD_MAPPER (acos)
+FORWARD_MAPPER (acosh)
+FORWARD_MAPPER (angle)
+FORWARD_MAPPER (arg)
+FORWARD_MAPPER (asin)
+FORWARD_MAPPER (asinh)
+FORWARD_MAPPER (atan)
+FORWARD_MAPPER (atanh)
+FORWARD_MAPPER (ceil)
+FORWARD_MAPPER (conj)
+FORWARD_MAPPER (cos)
+FORWARD_MAPPER (cosh)
+FORWARD_MAPPER (exp)
+FORWARD_MAPPER (expm1)
+FORWARD_MAPPER (fix)
+FORWARD_MAPPER (floor)
+FORWARD_MAPPER (imag)
+FORWARD_MAPPER (log)
+FORWARD_MAPPER (log2)
+FORWARD_MAPPER (log10)
+FORWARD_MAPPER (log1p)
+FORWARD_MAPPER (real)
+FORWARD_MAPPER (round)
+FORWARD_MAPPER (roundb)
+FORWARD_MAPPER (signum)
+FORWARD_MAPPER (sin)
+FORWARD_MAPPER (sinh)
+FORWARD_MAPPER (sqrt)
+FORWARD_MAPPER (tan)
+FORWARD_MAPPER (tanh)
+FORWARD_MAPPER (finite)
+FORWARD_MAPPER (isinf)
+FORWARD_MAPPER (isna)
+FORWARD_MAPPER (isnan)
+
+DEFINE_OCTAVE_ALLOCATOR (octave_perm_matrix);
+
+DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_perm_matrix, 
+                                     "permutation matrix", "double");
+
+static octave_base_value *
+default_numeric_conversion_function (const octave_base_value& a)
+{
+  CAST_CONV_ARG (const octave_perm_matrix&);
+
+  return new octave_matrix (v.matrix_value ());
+}
+
+octave_base_value::type_conv_info
+octave_perm_matrix::numeric_conversion_function (void) const
+{
+  return octave_base_value::type_conv_info (default_numeric_conversion_function,
+                                            octave_matrix::static_type_id ());
+}
+
+static octave_base_value *
+default_numeric_demotion_function (const octave_base_value& a)
+{
+  CAST_CONV_ARG (const octave_perm_matrix&);
+
+  return new octave_float_perm_matrix (v.perm_matrix_value ());
+}
+
+octave_base_value::type_conv_info
+octave_perm_matrix::numeric_demotion_function (void) const
+{
+  return octave_base_value::type_conv_info (default_numeric_demotion_function,
+                                            octave_float_perm_matrix::static_type_id ());
+}
+
+octave_base_value *
+octave_perm_matrix::try_narrowing_conversion (void)
+{
+  octave_base_value *retval = 0;
+
+  if (matrix.nelem () == 1)
+    retval = new octave_scalar (matrix (0, 0));
+
+  return retval;
+}
+
diff --git a/src/ov-perm.h b/src/ov-perm.h
new file mode 100644
--- /dev/null
+++ b/src/ov-perm.h
@@ -0,0 +1,257 @@
+/*
+
+Copyright (C) 2008 Jaroslav Hajek <highegg@gmail.com>
+
+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_perm_matrix_h)
+#define octave_perm_matrix_h 1
+
+#include "mx-base.h"
+#include "str-vec.h"
+
+#include "ov-base.h"
+#include "ov-typeinfo.h"
+#include "oct-obj.h"
+
+class 
+OCTINTERP_API
+octave_perm_matrix : public octave_base_value
+{
+public:
+  octave_perm_matrix (void) : matrix () { }
+
+  octave_perm_matrix (const PermMatrix& p) : matrix (p) { }
+
+  octave_base_value *clone (void) const { return new octave_perm_matrix (*this); }
+  octave_base_value *empty_clone (void) const { return new octave_perm_matrix (); }
+
+  type_conv_info numeric_conversion_function (void) const;
+
+  type_conv_info numeric_demotion_function (void) const;
+
+  octave_base_value *try_narrowing_conversion (void);
+
+  size_t byte_size (void) const { return matrix.byte_size (); }
+
+  octave_value squeeze (void) const { return matrix; }
+
+  octave_value subsref (const std::string& type,
+			const std::list<octave_value_list>& idx);
+
+  octave_value_list subsref (const std::string& type,
+			     const std::list<octave_value_list>& idx, int)
+    { return subsref (type, idx); }
+
+  octave_value do_index_op (const octave_value_list& idx,
+			    bool resize_ok = false);
+
+  dim_vector dims (void) const { return matrix.dims (); }
+
+  octave_idx_type nnz (void) const { return matrix.rows (); }
+
+  octave_value reshape (const dim_vector& new_dims) const
+    { return to_dense ().reshape (new_dims); }
+
+  octave_value permute (const Array<int>& vec, bool inv = false) const
+    { return to_dense ().permute (vec, inv); }
+
+  octave_value resize (const dim_vector& dv, bool fill = false) const
+    { return to_dense ().resize (dv, fill); }
+
+  octave_value all (int dim = 0) const { return to_dense ().all (dim); }
+  octave_value any (int dim = 0) const { return to_dense ().any (dim); }
+
+  MatrixType matrix_type (void) const { return MatrixType::Permuted_Diagonal; }
+  MatrixType matrix_type (const MatrixType&) const
+    { return matrix_type (); }
+
+  octave_value diag (octave_idx_type k = 0) const
+    { return to_dense () .diag (k); }
+
+  octave_value sort (octave_idx_type dim = 0, sortmode mode = ASCENDING) const
+    { return to_dense ().sort (dim, mode); }
+  octave_value sort (Array<octave_idx_type> &sidx, octave_idx_type dim = 0,
+		     sortmode mode = ASCENDING) const
+    { return to_dense ().sort (sidx, dim, mode); }
+
+  bool is_matrix_type (void) const { return true; }
+
+  bool is_numeric_type (void) const { return true; }
+
+  bool is_defined (void) const { return true; }
+
+  bool is_constant (void) const { return true; }
+
+  bool is_real_matrix (void) const { return true; }
+
+  bool is_real_type (void) const { return true; }
+
+  bool is_double_type (void) const { return true; }
+
+  bool is_float_type (void) const { return true; }
+
+  bool is_true (void) const;
+
+  bool valid_as_scalar_index (void) const;
+
+  double double_value (bool = false) const;
+
+  float float_value (bool = false) const;
+
+  double scalar_value (bool frc_str_conv = false) const
+    { return double_value (frc_str_conv); }
+
+  idx_vector index_vector (void) const;
+
+  PermMatrix perm_matrix_value (void) const
+    { return matrix; }
+
+  Matrix matrix_value (bool = false) const;
+
+  FloatMatrix float_matrix_value (bool = false) const;
+
+  Complex complex_value (bool = false) const;
+
+  FloatComplex float_complex_value (bool = false) const;
+
+  ComplexMatrix complex_matrix_value (bool = false) const;
+
+  FloatComplexMatrix float_complex_matrix_value (bool = false) const;
+
+  ComplexNDArray complex_array_value (bool = false) const;
+   
+  FloatComplexNDArray float_complex_array_value (bool = false) const;
+   
+  boolNDArray bool_array_value (bool warn = false) const;
+
+  charNDArray char_array_value (bool = false) const;
+  
+  NDArray array_value (bool = false) const; 
+
+  FloatNDArray float_array_value (bool = false) const;
+
+  SparseMatrix sparse_matrix_value (bool = false) const;
+
+  SparseComplexMatrix sparse_complex_matrix_value (bool = false) const;
+
+  int8NDArray
+  int8_array_value (void) const { return to_dense ().int8_array_value (); }
+
+  int16NDArray
+  int16_array_value (void) const { return to_dense ().int16_array_value (); }
+
+  int32NDArray
+  int32_array_value (void) const { return to_dense ().int32_array_value (); }
+
+  int64NDArray
+  int64_array_value (void) const { return to_dense ().int64_array_value (); }
+
+  uint8NDArray
+  uint8_array_value (void) const { return to_dense ().uint8_array_value (); }
+
+  uint16NDArray
+  uint16_array_value (void) const { return to_dense ().uint16_array_value (); }
+
+  uint32NDArray
+  uint32_array_value (void) const { return to_dense ().uint32_array_value (); }
+
+  uint64NDArray
+  uint64_array_value (void) const { return to_dense ().uint64_array_value (); }
+
+  octave_value convert_to_str_internal (bool pad, bool force, char type) const;
+
+  void print_raw (std::ostream& os, bool pr_as_read_syntax = false) const;
+
+  bool save_ascii (std::ostream& os);
+
+  bool save_binary (std::ostream& os, bool& save_as_floats);
+
+#if defined (HAVE_HDF5)
+  bool save_hdf5 (hid_t loc_id, const char *name, bool save_as_floats);
+#endif
+
+  int write (octave_stream& os, int block_size,
+	     oct_data_conv::data_type output_type, int skip,
+	     oct_mach_info::float_format flt_fmt) const;
+
+  mxArray *as_mxArray (void) const;
+
+  bool print_as_scalar (void) const;
+
+  void print (std::ostream& os, bool pr_as_read_syntax = false) const;
+
+  void print_info (std::ostream& os, const std::string& prefix) const;
+
+  octave_value erf (void) const;
+  octave_value erfc (void) const;
+  octave_value gamma (void) const;
+  octave_value lgamma (void) const;
+  octave_value abs (void) const;
+  octave_value acos (void) const;
+  octave_value acosh (void) const;
+  octave_value angle (void) const;
+  octave_value arg (void) const;
+  octave_value asin (void) const;
+  octave_value asinh (void) const;
+  octave_value atan (void) const;
+  octave_value atanh (void) const;
+  octave_value ceil (void) const;
+  octave_value conj (void) const;
+  octave_value cos (void) const;
+  octave_value cosh (void) const;
+  octave_value exp (void) const;
+  octave_value expm1 (void) const;
+  octave_value fix (void) const;
+  octave_value floor (void) const;
+  octave_value imag (void) const;
+  octave_value log (void) const;
+  octave_value log2 (void) const;
+  octave_value log10 (void) const;
+  octave_value log1p (void) const;
+  octave_value real (void) const;
+  octave_value round (void) const;
+  octave_value roundb (void) const;
+  octave_value signum (void) const;
+  octave_value sin (void) const;
+  octave_value sinh (void) const;
+  octave_value sqrt (void) const;
+  octave_value tan (void) const;
+  octave_value tanh (void) const;
+  octave_value finite (void) const;
+  octave_value isinf (void) const;
+  octave_value isna (void) const;
+  octave_value isnan (void) const;
+
+protected:
+
+  PermMatrix matrix;  
+
+  virtual octave_value to_dense (void) const;
+
+  mutable octave_value dense_cache;
+
+private:
+
+  DECLARE_OCTAVE_ALLOCATOR
+
+  DECLARE_OV_TYPEID_FUNCTIONS_AND_DATA
+};
+
+#endif
diff --git a/src/ov.cc b/src/ov.cc
--- a/src/ov.cc
+++ b/src/ov.cc
@@ -42,6 +42,8 @@
 #include "ov-flt-re-mat.h"
 #include "ov-re-diag.h"
 #include "ov-flt-re-diag.h"
+#include "ov-perm.h"
+#include "ov-flt-perm.h"
 #include "ov-bool-sparse.h"
 #include "ov-cx-sparse.h"
 #include "ov-re-sparse.h"
@@ -691,6 +693,15 @@
   maybe_mutate ();
 }
 
+octave_value::octave_value (const PermMatrix& p, bool single)
+  : rep ()
+{
+  if (single)
+    rep = new octave_float_perm_matrix (p);
+  else
+    rep = new octave_perm_matrix (p);
+}
+
 octave_value::octave_value (bool b)
   : rep (new octave_bool (b))
 {
@@ -2367,6 +2378,8 @@
   octave_float_diag_matrix::register_type ();
   octave_float_complex_matrix::register_type ();
   octave_float_complex_diag_matrix::register_type ();
+  octave_perm_matrix::register_type ();
+  octave_float_perm_matrix::register_type ();
   octave_null_matrix::register_type ();
   octave_null_str::register_type ();
   octave_null_sq_str::register_type ();
diff --git a/src/ov.h b/src/ov.h
--- a/src/ov.h
+++ b/src/ov.h
@@ -207,6 +207,7 @@
   octave_value (const FloatComplexRowVector& v);
   octave_value (const ComplexColumnVector& v);
   octave_value (const FloatComplexColumnVector& v);
+  octave_value (const PermMatrix& p, bool single = false);
   octave_value (bool b);
   octave_value (const boolMatrix& bm, const MatrixType& t = MatrixType());
   octave_value (const boolNDArray& bnda);