| Re: [eigen] Specializing max_coeff_visitor for some number types |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
On Mon, 7 Apr 2014, Christoph Hertzberg wrote:
On 07.04.2014 15:38, Marc Glisse wrote:
On Mon, 7 Apr 2014, Christoph Hertzberg wrote:
I would rather suggest to out-factor a method find_best_pivot which
then can be specialized for special types.
At what point exactly should it split? The easiest would be that LU calls:
DenseBase<Derived>::bestPivotCoeff(IndexType* index) const
which visits with a best_coeff_visitor, which uses
internal::scalar_better_coeff_op, which by default compares the absolute
values. But then we lost the parallelism in computing the absolute values.
I think this is a good point to split. We don't use SIMD-parallelism in abs()
at this point, because the expression is later used by a non-parallel
operation and it is so cheap that evaluating it to a temporary is not worth
it. Theoretically, better_coeff_op could be implemented to also provide a
packet operation (which needs to be finalized with a reduce() operation at
the end), but I guess most of the time finding the pivot is negligible
compared to the heavy matrix operations during decomposition.
The attached patch doesn't change the result of make check (adolc is still
broken, there is a missing -lGL somewhere, and parmetis doesn't chose any
MPI so finds none).
As you can see in LU, there is a subtlety I hadn't thought of:
max(abs) was returning a RealScalar, not a Scalar. That's not a big
problem when we only want to compare it with 0.
I think I didn't break maxPivot() in FullPivLU.
householder, SVD, etc would likely need a similar treatment.
I want the best in my case (well, not necessarily the best, but one good
enough). And the biggest lower-bound is just a simple heuristic, given
the choice between (1,1) and (2,20) I'd pick the first interval as the
better pivot ;-) But it may not be worth doing anything more subtle.
Hm, well the question is for a matrix, e.g.,
[ (1,1) row1 ]
[ (2,20) row2 ]
whether
row2 - (2,20)/(1,1) * row1
or
row1 - (1,1)/(2,20) * row2
is better (simplified example). But I agree that simply using the element
with the bigger lower bound may not always be the best choice.
I believe you need to consider at least a third line. The uncertainty on
the pivot will propagate to the whole matrix (except its own row/col).
Is the patch roughly the right approach? Where do we go from there?
--
Marc Glisse
diff -r 4daecc7950b9 Eigen/src/Core/DenseBase.h
--- a/Eigen/src/Core/DenseBase.h Mon Apr 07 14:14:48 2014 +0100
+++ b/Eigen/src/Core/DenseBase.h Mon Apr 07 20:24:22 2014 +0200
@@ -432,9 +432,13 @@
template<typename IndexType> EIGEN_DEVICE_FUNC
typename internal::traits<Derived>::Scalar maxCoeff(IndexType* row, IndexType* col) const;
template<typename IndexType> EIGEN_DEVICE_FUNC
+ typename internal::traits<Derived>::Scalar bestCoeff(IndexType* row, IndexType* col) const;
+ template<typename IndexType> EIGEN_DEVICE_FUNC
typename internal::traits<Derived>::Scalar minCoeff(IndexType* index) const;
template<typename IndexType> EIGEN_DEVICE_FUNC
typename internal::traits<Derived>::Scalar maxCoeff(IndexType* index) const;
+ template<typename IndexType> EIGEN_DEVICE_FUNC
+ typename internal::traits<Derived>::Scalar bestCoeff(IndexType* index) const;
template<typename BinaryOp>
EIGEN_DEVICE_FUNC
diff -r 4daecc7950b9 Eigen/src/Core/Visitor.h
--- a/Eigen/src/Core/Visitor.h Mon Apr 07 14:14:48 2014 +0100
+++ b/Eigen/src/Core/Visitor.h Mon Apr 07 20:24:22 2014 +0200
@@ -162,6 +162,34 @@
};
};
+/** \internal
+ * \brief Visitor computing the best coefficient (as a pivot) with its value and coordinates
+ *
+ * \sa DenseBase::bestCoeff(Index*, Index*)
+ */
+template <typename Derived>
+struct best_coeff_visitor : coeff_visitor<Derived>
+{
+ typedef typename Derived::Index Index;
+ typedef typename Derived::Scalar Scalar;
+ void operator() (const Scalar& value, Index i, Index j)
+ {
+ if(internal::scalar_better_coeff_op<Scalar>()(value, this->res))
+ {
+ this->res = value;
+ this->row = i;
+ this->col = j;
+ }
+ }
+};
+
+template<typename Scalar>
+struct functor_traits<best_coeff_visitor<Scalar> > {
+ enum {
+ Cost = functor_traits<scalar_better_coeff_op<Scalar> >::Cost
+ };
+};
+
} // end namespace internal
/** \returns the minimum of all coefficients of *this and puts in *row and *col its location.
@@ -232,6 +260,40 @@
return maxVisitor.res;
}
+/** \returns the best (as a pivot) of all coefficients of *this and puts in *row and *col its location.
+ * \warning the result is undefined if \c *this contains NaN.
+ *
+ * \sa DenseBase::minCoeff(IndexType*,IndexType*), DenseBase::maxCoeff(IndexType*,IndexType*), DenseBase::visitor()
+ */
+template<typename Derived>
+template<typename IndexType>
+typename internal::traits<Derived>::Scalar
+DenseBase<Derived>::bestCoeff(IndexType* rowPtr, IndexType* colPtr) const
+{
+ internal::best_coeff_visitor<Derived> bestVisitor;
+ this->visit(bestVisitor);
+ *rowPtr = bestVisitor.row;
+ if (colPtr) *colPtr = bestVisitor.col;
+ return bestVisitor.res;
+}
+
+/** \returns the best (as a pivot) of all coefficients of *this and puts in *index its location.
+ * \warning the result is undefined if \c *this contains NaN.
+ *
+ * \sa DenseBase::bestCoeff(IndexType*,IndexType*), DenseBase::maxCoeff(IndexType*), DenseBase::minCoeff(IndexType*), DenseBase::visitor()
+ */
+template<typename Derived>
+template<typename IndexType>
+typename internal::traits<Derived>::Scalar
+DenseBase<Derived>::bestCoeff(IndexType* index) const
+{
+ EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+ internal::best_coeff_visitor<Derived> bestVisitor;
+ this->visit(bestVisitor);
+ *index = (RowsAtCompileTime==1) ? bestVisitor.col : bestVisitor.row;
+ return bestVisitor.res;
+}
+
} // end namespace Eigen
#endif // EIGEN_VISITOR_H
diff -r 4daecc7950b9 Eigen/src/Core/functors/BinaryFunctors.h
--- a/Eigen/src/Core/functors/BinaryFunctors.h Mon Apr 07 14:14:48 2014 +0100
+++ b/Eigen/src/Core/functors/BinaryFunctors.h Mon Apr 07 20:24:22 2014 +0200
@@ -240,6 +240,27 @@
};
};
+/** \internal
+ * \brief Template functor to determine if a scalar would be a better pivot
+ *
+ * \sa class CwiseBinaryOp, MatrixBase::cwiseMax, class VectorwiseOp, MatrixBase::maxCoeff()
+ */
+template<typename Scalar> struct scalar_better_coeff_op {
+ EIGEN_EMPTY_STRUCT_CTOR(scalar_better_coeff_op)
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool operator() (const Scalar& a, const Scalar& b) const
+ {
+ internal::scalar_abs_op<Scalar> abso;
+ return abso(a) > abso(b);
+ }
+};
+template<typename Scalar>
+struct functor_traits<scalar_better_coeff_op<Scalar> > {
+ enum {
+ Cost = NumTraits<Scalar>::AddCost,
+ PacketAccess = false
+ };
+};
+
/** \internal
diff -r 4daecc7950b9 Eigen/src/LU/FullPivLU.h
--- a/Eigen/src/LU/FullPivLU.h Mon Apr 07 14:14:48 2014 +0100
+++ b/Eigen/src/LU/FullPivLU.h Mon Apr 07 20:24:22 2014 +0200
@@ -442,10 +442,10 @@
// biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
Index row_of_biggest_in_corner, col_of_biggest_in_corner;
+ using std::abs;
RealScalar biggest_in_corner;
- biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
- .cwiseAbs()
- .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+ biggest_in_corner = abs(m_lu.bottomRightCorner(rows-k, cols-k)
+ .bestCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner));
row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
col_of_biggest_in_corner += k; // need to add k to them.
diff -r 4daecc7950b9 Eigen/src/LU/PartialPivLU.h
--- a/Eigen/src/LU/PartialPivLU.h Mon Apr 07 14:14:48 2014 +0100
+++ b/Eigen/src/LU/PartialPivLU.h Mon Apr 07 20:24:22 2014 +0200
@@ -249,13 +249,13 @@
Index rcols = cols-k-1;
Index row_of_biggest_in_col;
- RealScalar biggest_in_corner
- = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
+ Scalar biggest_in_corner
+ = lu.col(k).tail(rows-k).bestCoeff(&row_of_biggest_in_col);
row_of_biggest_in_col += k;
row_transpositions[k] = PivIndex(row_of_biggest_in_col);
- if(biggest_in_corner != RealScalar(0))
+ if(biggest_in_corner != Scalar(0))
{
if(k != row_of_biggest_in_col)
{