Mercurial > hg > octave-max
changeset 8367:445d27d79f4e
support permutation matrix objects
line wrap: on
line diff
--- 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);
--- 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
--- 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.
--- 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; }
--- 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++ ***
--- 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
--- 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;
--- 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
--- 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 \
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; +}
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
--- 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);
--- 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; };
--- a/liboctave/dDiagMatrix.cc +++ b/liboctave/dDiagMatrix.cc @@ -387,6 +387,7 @@ return d; } + std::ostream& operator << (std::ostream& os, const DiagMatrix& a) {
--- 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);
--- 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?
--- 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);
--- 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.
--- 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; }
--- 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.
--- 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++ ***
--- 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
--- 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.
--- 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; }
--- 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++ ***
--- 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
--- 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);
--- 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?
--- 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);
--- 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.
--- 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; }
--- 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++ ***
--- 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
--- 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) {
--- 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.
--- 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"
--- 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;
--- 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: ***
--- 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
--- 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.
--- 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);
--- 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 (); }
--- 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)
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"
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"
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"
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"
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"
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); +} +
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"
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"
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"
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); +}
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 +}
--- 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; }
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; +} +
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
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; +} +
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
--- 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 ();
--- 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);