Re: [eigen] A not so simple product

["Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>]

Did you leave asserts? If not, can you retry with asserts enabled?

If this isn't catched by an assert, then there are basically two possibilities:
1) you discovered a really bad bug in Eigen; OR
2) you did something wrong with those vtk pointers.

Cheers,
Benoit

2008/12/9 Benjamin Schindler <bschindler@xxxxxxxxxxxxxxx>:
> Hi
>
> Thanks for the voluminous responses. But I'm running into a segfault it
> seems if you could give me a hint, that would be great. The code which is
> 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 to continue,
> I get the segfault. Because the output is large, I append it as 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 -------
>
>

---