Re: [ublas] permutation matrix class

From: Gunter Winkler (guwi17_at_[hidden])
Date: 2006-12-20 03:37:45

Daniel Elliott schrieb:
> Yes, I understand your confusion.  I am not being nearly clear enough.
> I should emphasize that this all began with a desire to use
> permutation matrices like I use other matrix types in uBlas (with
> operators and functions like prod() ).
In this case the transformation (or selection) matrix behaves exactly 
like matrix_indirect. Thus one should store the mapping i->p(i). Then a 
specialized prod() could look like
matrix_indirect< MATRIX >
prod(permutation_matrix& P, matrix& A) {
    return P.make_row_indirect_matrix(A);
}
(Therefore P stores the data in an indirect_array, which already 
provides some useful functions)
matrix_indirect< MATRIX >
prod(matrix& A, permutation_matrix& P) {
    return P.make_column_indirect_matrix(A);
}
matrix_indirect< MATRIX >
prod3(permutation_matrix& P, matrix& A, permutation_matrix& Q) {
    return P.make_indirect_matrix(A,Q);
}

permutation_matrix could even be a rectangular selection matrix, that 
doubles, drops or exchanges rows/columns.
What do you think?

The second application is the actual permutation of a vector:
x = prod(P, y) // similar to matrix permutation
x = prod(P, x) // needs temporary or swap representation of P

    // this should work, but its from a quite old file ...
    template <class VEC>
    void permutation_to_swap(VEC & x)
    {
      typedef typename VEC::size_type  size_type;

      for (size_type i=0; i<x.size(); ++i) {
        size_type j = x(i);
        if (j<i) {
          while (j<i) {
            j = x(j);
          }
          x(i) = j;
        }
      }
    }

  template <class VEC>
  VEC inverse_permutation(VEC & p)
  {
    typedef typename VEC::size_type  size_type;
    VEC     res(p.size());

    for (size_type i=0; i<p.size(); ++i) {
      res(p(i)) = i;
    }

    return res;
  }

> With respect to permutation matrices, I would prefer it if the user
> could choose which representation would fit them the best.
The swap representation is bad for prod(P, A), because it needs a 
temporary copy of A. So IMO, one should force the user to only use the 
map representation for products.

btw. for the inverse permutation one could use some "reverse indirect 
matrix" proxy. That internally iterates from the first row/column of A 
to the last, but it.index() always returns p(i) as index. So one saves 
an inversion of the permutation, but one gets a nonconformant matrix 
proxy (the elements in rows and columns are not sorted and thus cannot 
(efficiently) be accessed by index).

mfg
Gunter