# HG changeset patch # User David Spies # Date 1403141920 21600 # Node ID aa9ca67f09fb95d401ef17b62c317a93d9cdc465 # Parent 3d33fe79816c07aa23a9baedf26c094f1a325a84 make all permutation matrices column permutations (bug #42418) Making all permutation matrices column permutations allows permutation matrices to be iterated over in column-major order just like any other matrix. * PermMatrix.h, PermMatrix.cc (PermMatrix::_colp): Delete member variable and all uses. (PermMatrix::transpose): Create new matrix with internal representation flipped. (PermMatrix::pos_power): New function. (PermMatrix::eye): Call the one-argument constructor which already defaults to the identity matrix. (PermMatrix::is_col_perm): Unconditionally return true. (PermMatrix::is_row_perm): Unconditionally return false. (PermMatrix::data, PermMatrix::fortran_vec, PermMatrix::pvec): Delete. (PermMatrix::col_perm_vec): New function. * lu.cc: New test. * base-lu.cc (base_lu::base_lu): Call transpose in ipvt initialization. * find.cc, kron.cc, pr-output.cc, ov-perm.cc, PermMatrix.cc, PermMatrix.h, Sparse.cc, dMatrix.cc, fMatrix.cc, CmplxQRP.cc, base-lu.cc, dbleQRP.cc, fCmplxQRP.cc, floatQRP.cc, Sparse-perm-op-defs.h, mx-op-defs.h: Adapt to PermMatrix changes. diff -r 3d33fe79816c -r aa9ca67f09fb libinterp/corefcn/find.cc --- a/libinterp/corefcn/find.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/libinterp/corefcn/find.cc Wed Jun 18 19:38:40 2014 -0600 @@ -263,33 +263,15 @@ if (count > 0) { - const octave_idx_type* p = v.data (); - if (v.is_col_perm ()) + const Array& p = v.col_perm_vec (); + for (octave_idx_type k = 0; k < count; k++) { - for (octave_idx_type k = 0; k < count; k++) - { - OCTAVE_QUIT; - const octave_idx_type j = start_nc + k; - const octave_idx_type i = p[j]; - i_idx(k) = static_cast (1+i); - j_idx(k) = static_cast (1+j); - idx(k) = j * nc + i + 1; - } - } - else - { - for (octave_idx_type k = 0; k < count; k++) - { - OCTAVE_QUIT; - const octave_idx_type i = start_nc + k; - const octave_idx_type j = p[i]; - // Scatter into the index arrays according to - // j adjusted by the start point. - const octave_idx_type koff = j - start_nc; - i_idx(koff) = static_cast (1+i); - j_idx(koff) = static_cast (1+j); - idx(koff) = j * nc + i + 1; - } + OCTAVE_QUIT; + const octave_idx_type j = start_nc + k; + const octave_idx_type i = p(j); + i_idx(k) = static_cast (1+i); + j_idx(k) = static_cast (1+j); + idx(k) = j * nc + i + 1; } } else diff -r 3d33fe79816c -r aa9ca67f09fb libinterp/corefcn/kron.cc --- a/libinterp/corefcn/kron.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/libinterp/corefcn/kron.cc Wed Jun 18 19:38:40 2014 -0600 @@ -136,39 +136,18 @@ { octave_idx_type na = a.rows (); octave_idx_type nb = b.rows (); - const octave_idx_type *pa = a.data (); - const octave_idx_type *pb = b.data (); - PermMatrix c(na*nb); // Row permutation. - octave_idx_type *pc = c.fortran_vec (); - - bool cola = a.is_col_perm (); - bool colb = b.is_col_perm (); - if (cola && colb) - { - for (octave_idx_type i = 0; i < na; i++) - for (octave_idx_type j = 0; j < nb; j++) - pc[pa[i]*nb+pb[j]] = i*nb+j; - } - else if (cola) + const Array& pa = a.col_perm_vec (); + const Array& pb = b.col_perm_vec (); + Array res_perm; + octave_idx_type rescol = 0; + for (octave_idx_type i = 0; i < na; i++) { - for (octave_idx_type i = 0; i < na; i++) - for (octave_idx_type j = 0; j < nb; j++) - pc[pa[i]*nb+j] = i*nb+pb[j]; - } - else if (colb) - { - for (octave_idx_type i = 0; i < na; i++) - for (octave_idx_type j = 0; j < nb; j++) - pc[i*nb+pb[j]] = pa[i]*nb+j; - } - else - { - for (octave_idx_type i = 0; i < na; i++) - for (octave_idx_type j = 0; j < nb; j++) - pc[i*nb+j] = pa[i]*nb+pb[j]; + octave_idx_type a_add = pa(i) * nb; + for (octave_idx_type j = 0; j < nb; j++) + res_perm.xelem (rescol++) = a_add + pb(j); } - return c; + return PermMatrix (res_perm, true); } template diff -r 3d33fe79816c -r aa9ca67f09fb libinterp/corefcn/lu.cc --- a/libinterp/corefcn/lu.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/libinterp/corefcn/lu.cc Wed Jun 18 19:38:40 2014 -0600 @@ -842,6 +842,20 @@ %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps); %! %!testif HAVE_QRUPDATE_LUU +%! [L,U,P] = lu (A); +%! [~,ordcols] = max (P,[],1); +%! [~,ordrows] = max (P,[],2); +%! P1 = eye (size(P))(:,ordcols); +%! P2 = eye (size(P))(ordrows,:); +%! assert(P1 == P); +%! assert(P2 == P); +%! [L,U,P] = luupdate (L,U,P,u,v); +%! [L,U,P1] = luupdate (L,U,P1,u,v); +%! [L,U,P2] = luupdate (L,U,P2,u,v); +%! assert(P1 == P); +%! assert(P2 == P); +%! +%!testif HAVE_QRUPDATE_LUU %! [L,U,P] = lu (Ac); %! [L,U,P] = luupdate (L,U,P,uc,vc); %! assert (norm (vec (tril (L)-L), Inf) == 0); diff -r 3d33fe79816c -r aa9ca67f09fb libinterp/corefcn/pr-output.cc --- a/libinterp/corefcn/pr-output.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/libinterp/corefcn/pr-output.cc Wed Jun 18 19:38:40 2014 -0600 @@ -2485,11 +2485,10 @@ if (pr_as_read_syntax) { - Array pvec = m.pvec (); - bool colp = m.is_col_perm (); + Array pvec = m.col_perm_vec (); os << "eye ("; - if (colp) os << ":, "; + os << ":, "; octave_idx_type col = 0; while (col < nc) @@ -2520,7 +2519,6 @@ else os << " ...\n"; } - if (! colp) os << ", :"; os << ")"; } else diff -r 3d33fe79816c -r aa9ca67f09fb libinterp/octave-value/ov-perm.cc --- a/libinterp/octave-value/ov-perm.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/libinterp/octave-value/ov-perm.cc Wed Jun 18 19:38:40 2014 -0600 @@ -261,9 +261,9 @@ octave_perm_matrix::save_ascii (std::ostream& os) { os << "# size: " << matrix.rows () << "\n"; - os << "# orient: " << (matrix.is_col_perm () ? 'c' : 'r') << '\n'; + os << "# orient: c\n"; - Array pvec = matrix.pvec (); + Array pvec = matrix.col_perm_vec (); octave_idx_type n = pvec.length (); ColumnVector tmp (n); for (octave_idx_type i = 0; i < n; i++) tmp(i) = pvec(i) + 1; @@ -314,11 +314,12 @@ { int32_t sz = matrix.rows (); - bool colp = matrix.is_col_perm (); + bool colp = true; os.write (reinterpret_cast (&sz), 4); os.write (reinterpret_cast (&colp), 1); - os.write (reinterpret_cast (matrix.data ()), - matrix.byte_size ()); + const Array& col_perm = matrix.col_perm_vec (); + os.write (reinterpret_cast (col_perm.data ()), + col_perm.byte_size ()); return true; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/array/PermMatrix.cc --- a/liboctave/array/PermMatrix.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/array/PermMatrix.cc Wed Jun 18 19:38:40 2014 -0600 @@ -36,8 +36,8 @@ ("PermMatrix: invalid permutation vector"); } -PermMatrix::PermMatrix (const Array& p, bool colp, bool check) - : Array (p), _colp(colp) +void +PermMatrix::setup (const Array& p, bool colp, bool check) { if (check) { @@ -47,12 +47,28 @@ Array::operator = (Array ()); } } + + if (!colp) + *this = this->transpose (); +} + +PermMatrix::PermMatrix (const Array& p) + : Array (p) +{ + setup (p, false, true); } -PermMatrix::PermMatrix (const idx_vector& idx, bool colp, octave_idx_type n) - : Array (), _colp(colp) +PermMatrix::PermMatrix (const Array& p, bool colp, bool check) + : Array (p) +{ + setup (p, colp, check); +} + +void +PermMatrix::setup (const idx_vector& idx, bool colp, octave_idx_type n) { octave_idx_type len = idx.length (n); + if (! idx.is_permutation (len)) gripe_invalid_permutation (); else @@ -61,12 +77,28 @@ for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i); Array::operator = (idxa); } + + if (!colp) + *this = this->transpose (); +} + +PermMatrix::PermMatrix (const idx_vector& idx) + : Array () +{ + setup (idx, false, 0); +} + +PermMatrix::PermMatrix (const idx_vector& idx, bool colp, octave_idx_type n) + : Array () +{ + setup (idx, colp, n); } PermMatrix::PermMatrix (octave_idx_type n) - : Array (dim_vector (n, 1)), _colp (false) + : Array (dim_vector (n, 1)) { - for (octave_idx_type i = 0; i < n; i++) xelem (i) = i; + for (octave_idx_type i = 0; i < n; i++) + xelem (i) = i; } octave_idx_type @@ -86,8 +118,13 @@ PermMatrix PermMatrix::transpose (void) const { - PermMatrix retval (*this); - retval._colp = ! retval._colp; + octave_idx_type len = Array::length (); + + PermMatrix retval (len); + + for (octave_idx_type i = 0; i < len; ++i) + retval.xelem (xelem (i)) = i; + return retval; } @@ -133,17 +170,20 @@ } PermMatrix -PermMatrix::power (octave_idx_type m) const +PermMatrix::power(octave_idx_type m) const +{ + if (m < 0) + return inverse ().pos_power (-m); + else if (m > 0) + return pos_power (m); + else + return PermMatrix (rows ()); +} + +PermMatrix +PermMatrix::pos_power (octave_idx_type m) const { octave_idx_type n = rows (); - bool res_colp = _colp; - if (m < 0) - { - res_colp = ! res_colp; - m = -m; - } - else if (m == 0) - return PermMatrix (n); const octave_idx_type *p = data (); Array res_pvec (dim_vector (n, 1), -1); @@ -176,43 +216,29 @@ } - return PermMatrix (res_pvec, res_colp, false); + return PermMatrix (res_pvec, true, false); } PermMatrix PermMatrix::eye (octave_idx_type n) { - Array p (dim_vector (n, 1)); - for (octave_idx_type i = 0; i < n; i++) - p(i) = i; - - return PermMatrix (p, false, false); + return PermMatrix (n); } PermMatrix operator *(const PermMatrix& a, const PermMatrix& b) { - const Array ia = a.pvec (); - const Array ib = b.pvec (); PermMatrix r; + + const Array ia = a.col_perm_vec (); + const Array ib = b.col_perm_vec (); + 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 ra (dim_vector (n, 1)); - if (a._colp) - ra.assign (idx_vector (ia), ib); - else - ra.assign (idx_vector (ib), ia); - r = PermMatrix (ra, a._colp, false); - } + r = PermMatrix (ia.index (idx_vector (ib)), true, false); return r; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/array/PermMatrix.h --- a/liboctave/array/PermMatrix.h Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/array/PermMatrix.h Wed Jun 18 19:38:40 2014 -0600 @@ -31,20 +31,21 @@ class OCTAVE_API PermMatrix : protected Array { - public: - PermMatrix (void) : Array (), _colp (false) { } + PermMatrix (void) : Array () { } PermMatrix (octave_idx_type n); - PermMatrix (const Array& p, bool colp = false, - bool check = true); + PermMatrix (const Array& p) GCC_ATTR_DEPRECATED; + + PermMatrix (const Array& p, bool colp, bool check = true); - PermMatrix (const PermMatrix& m) - : Array (m), _colp(m._colp) { } + PermMatrix (const PermMatrix& m) : Array (m) { } - PermMatrix (const idx_vector& idx, bool colp = false, octave_idx_type n = 0); + PermMatrix (const idx_vector& idx) GCC_ATTR_DEPRECATED; + + PermMatrix (const idx_vector& idx, bool colp, octave_idx_type n = 0); octave_idx_type dim1 (void) const { return Array::length (); } @@ -68,15 +69,13 @@ dim_vector dims (void) const { return dim_vector (dim1 (), dim2 ()); } - Array pvec (void) const + const Array& col_perm_vec (void) const { return *this; } octave_idx_type elem (octave_idx_type i, octave_idx_type j) const { - return (_colp - ? ((Array::elem (j) == i) ? 1 : 0) - : ((Array::elem (i) == j) ? 1 : 0)); + return (Array::elem (j) == i) ? 1 : 0; } octave_idx_type @@ -102,20 +101,8 @@ // Efficient integer power of a permutation. PermMatrix power (octave_idx_type n) const; - bool is_col_perm (void) const { return _colp; } - bool is_row_perm (void) const { return !_colp; } - - friend OCTAVE_API PermMatrix operator *(const PermMatrix& a, - const PermMatrix& b); - - const octave_idx_type *data (void) const - { return Array::data (); } - - const octave_idx_type *fortran_vec (void) const - { return Array::fortran_vec (); } - - octave_idx_type *fortran_vec (void) - { return Array::fortran_vec (); } + bool is_col_perm (void) const { return true; } + bool is_row_perm (void) const { return false; } void print_info (std::ostream& os, const std::string& prefix) const { Array::print_info (os, prefix); } @@ -123,12 +110,17 @@ static PermMatrix eye (octave_idx_type n); private: - bool _colp; + + PermMatrix pos_power (octave_idx_type m) const; + + void setup (const Array& p, bool colp, bool check); + + void setup (const idx_vector& idx, bool colp, octave_idx_type n); }; // Multiplying permutations together. PermMatrix OCTAVE_API -operator *(const PermMatrix& a, const PermMatrix& b); +operator * (const PermMatrix& a, const PermMatrix& b); #endif diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/array/Sparse.cc --- a/liboctave/array/Sparse.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/array/Sparse.cc Wed Jun 18 19:38:40 2014 -0600 @@ -61,18 +61,10 @@ for (octave_idx_type i = 0; i <= n; i++) cidx (i) = i; - const Array pv = a.pvec (); - - if (a.is_row_perm ()) - { - for (octave_idx_type i = 0; i < n; i++) - ridx (pv(i)) = i; - } - else - { - for (octave_idx_type i = 0; i < n; i++) - ridx (i) = pv(i); - } + const Array pv = a.col_perm_vec (); + + for (octave_idx_type i = 0; i < n; i++) + ridx (i) = pv(i); for (octave_idx_type i = 0; i < n; i++) data (i) = 1.0; diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/array/dMatrix.cc --- a/liboctave/array/dMatrix.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/array/dMatrix.cc Wed Jun 18 19:38:40 2014 -0600 @@ -259,14 +259,10 @@ Matrix::Matrix (const PermMatrix& a) : MArray (a.dims (), 0.0) { - const Array ia (a.pvec ()); + const Array ia (a.col_perm_vec ()); 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; + for (octave_idx_type i = 0; i < len; i++) + elem (ia(i), i) = 1.0; } // FIXME: could we use a templated mixed-type copy function here? diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/array/fMatrix.cc --- a/liboctave/array/fMatrix.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/array/fMatrix.cc Wed Jun 18 19:38:40 2014 -0600 @@ -259,14 +259,10 @@ FloatMatrix::FloatMatrix (const PermMatrix& a) : MArray (a.dims (), 0.0) { - const Array ia (a.pvec ()); + const Array ia (a.col_perm_vec ()); 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; + for (octave_idx_type i = 0; i < len; i++) + elem (ia(i), i) = 1.0; } // FIXME: could we use a templated mixed-type copy function here? diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/numeric/CmplxQRP.cc --- a/liboctave/numeric/CmplxQRP.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/numeric/CmplxQRP.cc Wed Jun 18 19:38:40 2014 -0600 @@ -103,7 +103,7 @@ RowVector ComplexQRP::Pvec (void) const { - Array pa (p.pvec ()); + Array pa (p.col_perm_vec ()); RowVector pv (MArray (pa) + 1.0); return pv; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/numeric/base-lu.cc --- a/liboctave/numeric/base-lu.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/numeric/base-lu.cc Wed Jun 18 19:38:40 2014 -0600 @@ -30,7 +30,7 @@ template base_lu::base_lu (const lu_type& l, const lu_type& u, const PermMatrix& p) - : a_fact (u), l_fact (l), ipvt (p.pvec ()) + : a_fact (u), l_fact (l), ipvt (p.transpose ().col_perm_vec ()) { if (l.columns () != u.rows ()) (*current_liboctave_error_handler) ("lu: dimension mismatch"); diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/numeric/dbleQRP.cc --- a/liboctave/numeric/dbleQRP.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/numeric/dbleQRP.cc Wed Jun 18 19:38:40 2014 -0600 @@ -100,7 +100,7 @@ RowVector QRP::Pvec (void) const { - Array pa (p.pvec ()); + Array pa (p.col_perm_vec ()); RowVector pv (MArray (pa) + 1.0); return pv; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/numeric/fCmplxQRP.cc --- a/liboctave/numeric/fCmplxQRP.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/numeric/fCmplxQRP.cc Wed Jun 18 19:38:40 2014 -0600 @@ -103,7 +103,7 @@ FloatRowVector FloatComplexQRP::Pvec (void) const { - Array pa (p.pvec ()); + Array pa (p.col_perm_vec ()); FloatRowVector pv (MArray (pa) + 1.0f); return pv; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/numeric/floatQRP.cc --- a/liboctave/numeric/floatQRP.cc Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/numeric/floatQRP.cc Wed Jun 18 19:38:40 2014 -0600 @@ -100,7 +100,7 @@ FloatRowVector FloatQRP::Pvec (void) const { - Array pa (p.pvec ()); + Array pa (p.col_perm_vec ()); FloatRowVector pv (MArray (pa) + 1.0f); return pv; } diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/operators/Sparse-perm-op-defs.h --- a/liboctave/operators/Sparse-perm-op-defs.h Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/operators/Sparse-perm-op-defs.h Wed Jun 18 19:38:40 2014 -0600 @@ -67,17 +67,7 @@ return SM (); } - if (p.is_row_perm ()) - { - // Form the column permutation and then call the colpm_sm routine. - const octave_idx_type *prow = p.pvec ().data (); - OCTAVE_LOCAL_BUFFER (octave_idx_type, pcol, nr); - for (octave_idx_type i = 0; i < nr; ++i) - pcol[prow[i]] = i; - return octinternal_do_mul_colpm_sm (pcol, a); - } - else - return octinternal_do_mul_colpm_sm (p.pvec ().data (), a); + return octinternal_do_mul_colpm_sm (p.col_perm_vec ().data (), a); } template @@ -163,10 +153,7 @@ return SM (); } - if (p.is_row_perm ()) - return octinternal_do_mul_sm_rowpm (a, p.pvec ().data ()); - else - return octinternal_do_mul_sm_colpm (a, p.pvec ().data ()); + return octinternal_do_mul_sm_colpm (a, p.col_perm_vec ().data ()); } #endif // octave_Sparse_perm_op_defs_h diff -r 3d33fe79816c -r aa9ca67f09fb liboctave/operators/mx-op-defs.h --- a/liboctave/operators/mx-op-defs.h Thu Jun 19 10:35:39 2014 -0700 +++ b/liboctave/operators/mx-op-defs.h Wed Jun 18 19:38:40 2014 -0600 @@ -606,13 +606,8 @@ gripe_nonconformant ("operator *", p.rows (), p.columns (), nr, nc); \ else \ { \ - if (p.is_col_perm ()) \ - { \ - result = M (nr, nc); \ - result.assign (p.pvec (), idx_vector::colon, x); \ - } \ - else \ - result = x.index (p.pvec (), idx_vector::colon); \ + result = M (nr, nc); \ + result.assign (p.col_perm_vec (), idx_vector::colon, x); \ } \ \ return result; \ @@ -627,15 +622,7 @@ 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, p.pvec ()); \ - else \ - { \ - result = M (nr, nc); \ - result.assign (idx_vector::colon, p.pvec (), x); \ - } \ - } \ + result = x.index (idx_vector::colon, p.col_perm_vec ()); \ \ return result; \ }