Re: [eigen] A not so simple product

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

To: eigen@xxxxxxxxxxxxxxxxxxx
Subject: Re: [eigen] A not so simple product
From: Benjamin Schindler <bschindler@xxxxxxxxxxxxxxx>
Date: Tue, 09 Dec 2008 16:48:06 +0100

Hi

Thanks for the voluminous responses. But I'm running into a segfault itseems if you could give me a hint, that would be great. The code whichis commented out works as expected, but the new one does not

The code is as follows:

vtkFloatArray *_h = vtkFloatArray::SafeDownCast(pointData->GetArray(HName));

vtkFloatArray *_wrho =vtkFloatArray::SafeDownCast(pointData->GetArray(WRhoName));

vtkFloatArray *_v = vtkFloatArray::SafeDownCast(pointData->GetArray(VName));

Map<VectorXf> h(_h->GetPointer(0), _h->GetNumberOfTuples());
Map<VectorXf> wrho(_wrho->GetPointer(0), _wrho->GetNumberOfTuples());
Map<MatrixXf> v(_h->GetPointer(0), _v->GetNumberOfTuples(), 3);

v = wrho.cwise().inverse().asDiagonal() * v;

//v.col(0) = v.col(0).cwise() / wrho;
//v.col(1) = v.col(1).cwise() / wrho;
//v.col(2) = v.col(2).cwise() / wrho;

in gdb, it doesn't just segfault, but it stops. When I force it tocontinue, I get the segfault. Because the output is large, I append itas an attachment. Any ideas what's causing this?



Benoit Jacob wrote:

OK, just had a chat with Gael over IRC.

You should use the diagonal matrix product approach. It will be easy
to add vectorization support for that, I add that to the TODO.

Cheers,
Benoit

2008/12/9 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:

Indeed, the "obvious" way

for(i = 0; i < N; i++) m.row(i) /= rho(i)*w(i);

doesn't vectorize because rows have size 3... so if you want
vectorization you need a complicated trick as Gael shows.

Gael: in order to achieve vectorization of diagonal products, wouldn't
it be simplest to just change DiagonalProduct.h so it takes in account
the fact that coefficients come from a vector?

Something like replacing m_lhs.coeff(row, row) by
m_lhs.m_coeffs.coeff(row), since we know that m_lhs is a
DiagonalMatrix?

Cheers,
Benoit

2008/12/9 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:

Hi,

I guess there is no M and you wanted to write "VectorXf w(N);" instead of
"VectorXf w(M);" ?

Anyway, you can write it as a single expression like that:

m = ( (rho.cwise()*w).cwise().inverse().asDiagonal() ) * m;

the matrix-matrix product will be optimized since Eigen knows that the
left-hand side (lhs) matrix is a diagonal matrix. Unfortunately, Eigen does
not take into account yet the fact that the diagonal of lhs is actually
stored as a vector, and therefore this expression won't be vectorized.

The good think is that this is an optimization we can quite easily add. So I
would recommand to do so and wait a bit that such expressions get vectorized
by Eigen.

Note that in the meantime, if performance matters you can also do:

rho (or w, or a temporary) = rho.cwise()*w).cwise().inverse();
for (int k=0: k<3; ++k)
  m.col(i).cwise() *= rho (or w, or a temporary);

in which case all computations will be vectorized at the cost of additionnal
read/writes operations.



Below are some proposals related to that problem:


Another option would be to add a mechanism to extend a vector to a matrix
with repetitions:

m = m .cwise()* ( (rho.cwise()*w).cwise().inverse().duplicate(3) );

I'm not sure about "duplicate" for the name but I know matlab offers such a
feature. For the vectorization this is a bit tricky because such a
"Duplicate" expression can be vectorized in both direction: you can either
return column packets of row packets (by repeating a scalar in each
component of the packet) to accomodate the storage of m. However, this fine
tunning of the packet generation is something we cannot support without big
changes in Eigen (or ugly workarounds).


Finally, a third solution to implement this kind of thing could be to extend
the rowwise/colwise API:

m = m.colwise() * (rho.cwise()*w).cwise().inverse();

perhaps this last proposal is a bit confusing.


In any case the choice for the traversal of m is a bit tricky. Assuming m is
column major, we can process it column per column that is good wrt to the
cache of m but which requires a temporary to store
"rho.cwise()*w).cwise().inverse()", or we process m per block of packet-size
rows...


cheers,
Gael.


On Tue, Dec 9, 2008 at 2:10 PM, Benjamin Schindler
<bschindler@xxxxxxxxxxxxxxx> wrote:

Hi

Okay, I'm trying to figure the most efficient way to vectorize this
product:

MatrixXf m(N,3);
VectorXf rho(N);
VectorXf w(M);

What I would like to do: The i'th row of m is suppposed to get multiplied
by 1 / (rho(i)*w(i))

I've not yet found a way which actually works as a short
eigen-expression... is there any or do I have to go back to the old rusty
for loop? :)

Thanks
Benjamin

---

---

bschindl@foppa /local/viewer-build $ gdb console/parallel
GNU gdb Red Hat Linux (6.5-37.el5_2.2rh)
Copyright (C) 2006 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB.  Type "show warranty" for details.
This GDB was configured as "i386-redhat-linux-gnu"...Using host libthread_db library "/lib/libthread_db.so.1".

(gdb) run output.vtk
Starting program: /local/viewer-build/console/parallel output.vtk
[Thread debugging using libthread_db enabled]
0.00501 0 9.922e-07 / 0.0009312
0.00501 0 9.835e-07 / 0.0009339
[New Thread -1208998192 (LWP 11631)]
Cannot find new threads: debugger service failed
(gdb) where
#0  0x00fcf859 in Eigen::MatrixBase<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000> >::copyCoeff<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, row=306062, col=538, 
    other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Coeffs.h:329
#1  0x00fcf8c2 in Eigen::ei_assign_impl<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, 3, 0>::run (dst=@0xbf916168, src=@0xbf916178)
    at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:223
#2  0x00fcfbdd in Eigen::MatrixBase<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000> >::lazyAssign<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, other=@0xbf916178)
    at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:405
#3  0x00fcfc1c in Eigen::ei_assign_selector<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, false, false>::run (dst=@0xbf916168, 
    other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:420
#4  0x00fcffb5 in Matrix<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Matrix.h:417
#5  0x00fcfff9 in Product<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, Eigen::Map<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, 1> > (this=0xbf916168, lhs=@0xbf916178, 
    rhs=@0xbf916124) at /home/bschindl/env/include/eigen2/Eigen/src/Core/DiagonalProduct.h:79
#6  0x00fd0095 in Eigen::MatrixBase<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > >::operator*<Eigen::Map<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, 1> > (this=0xbf916178, 
    other=@0xbf916124) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Product.h:289
#7  0x00fce89d in vtkDamBreakSource::RequestData (this=0x813e5f8, request=0x81414f8, inputVector=0x0, outputVector=0x813f298)
    at /home/bschindl/masterarbeit/viewer/libsph/vtkDamBreakSource.cpp:59
---- snip -----
Quit
(gdb) c
Continuing.

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread -1208998192 (LWP 11631)]
0x00fcf859 in Eigen::MatrixBase<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000> >::copyCoeff<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, row=306062, col=538, 
    other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Coeffs.h:329
329       derived().coeffRef(row, col) = other.derived().coeff(row, col);
(gdb) where
#0  0x00fcf859 in Eigen::MatrixBase<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000> >::copyCoeff<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, row=306062, col=538, 
    other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Coeffs.h:329
#1  0x00fcf8c2 in Eigen::ei_assign_impl<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, 3, 0>::run (dst=@0xbf916168, src=@0xbf916178)
    at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:223
#2  0x00fcfbdd in Eigen::MatrixBase<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000> >::lazyAssign<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, other=@0xbf916178)
    at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:405
#3  0x00fcfc1c in Eigen::ei_assign_selector<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, false, false>::run (dst=@0xbf916168, 
    other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Assign.h:420
#4  0x00fcffb5 in Matrix<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > > (this=0xbf916168, other=@0xbf916178) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Matrix.h:417
#5  0x00fcfff9 in Product<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > >, Eigen::Map<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, 1> > (this=0xbf916168, lhs=@0xbf916178, 
    rhs=@0xbf916124) at /home/bschindl/env/include/eigen2/Eigen/src/Core/DiagonalProduct.h:79
#6  0x00fd0095 in Eigen::MatrixBase<Eigen::DiagonalMatrix<Eigen::CwiseUnaryOp<Eigen::ei_scalar_inverse_op<float>, Eigen::Map<Eigen::Matrix<float, 10000, 1, 0, 10000, 1>, 1> > > >::operator*<Eigen::Map<Eigen::Matrix<float, 10000, 10000, 0, 10000, 10000>, 1> > (this=0xbf916178, 
    other=@0xbf916124) at /home/bschindl/env/include/eigen2/Eigen/src/Core/Product.h:289
#7  0x00fce89d in vtkDamBreakSource::RequestData (this=0x813e5f8, request=0x81414f8, inputVector=0x0, outputVector=0x813f298)
    at /home/bschindl/masterarbeit/viewer/libsph/vtkDamBreakSource.cpp:59
----- snip -------

Follow-Ups:
- Re: [eigen] A not so simple product
  - From: Benoit Jacob

References:
- [eigen] A not so simple product
  - From: Benjamin Schindler
- Re: [eigen] A not so simple product
  - From: Gael Guennebaud
- Re: [eigen] A not so simple product
  - From: Benoit Jacob
- Re: [eigen] A not so simple product
  - From: Benoit Jacob

Messages sorted by: [ date | thread ]
Prev by Date: Re: [eigen] A not so simple product
Next by Date: Re: [eigen] A not so simple product
Previous by thread: Re: [eigen] A not so simple product
Next by thread: Re: [eigen] A not so simple product

Mail converted by MHonArc 2.6.19+

http://listengine.tuxfamily.org/