|
|
Ublas : |
From: Daniel Elliott (danelliottster_at_[hidden])
Date: 2006-12-20 20:00:35
Gunter,
On 12/20/06, Gunter Winkler <guwi17_at_[hidden]> wrote:
> 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);
> }
Let me see if I understand...
One option (especially when dealing with the swap representation of a
permutation matrix) is that we return a matrix_indirect which is a
clever class that looks and feels like a matrix but has the ability to
re-index matrix elements based upon array values?
> permutation_matrix could even be a rectangular selection matrix, that
> doubles, drops or exchanges rows/columns.
> What do you think?
I think that is great. Any idea how the programmer would specify
which type of permutation they want to define?
> 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.
Force them to use map representation for products - I'm not sure I
understand. The swap you use above is like the one in lu.hpp. Isn't
that one best for in-place and the matrix of ones and zeros (the map
representation) is best when assigning to another matrix?
Aside from reducing the amount of code that must be written, what
other benefits are there to restricting the user to one particular
representation? I tend to want to give a user the opportunity to do
something that is not optimal - just for the sake of flexibility.
I need to look at the uBlas code much more closely. How does prod()
currently decide on the type of matrix that will contain the result of
a calculation? Is this different than when you use the overloaded
operators? Can we let the user choose?
> 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).
This idea is very neat and is exactly the kind of optimization that
could really pay off in certain situations.
- dan elliott