Re: [eigen] Iterators with dense matrices: (was: Areas of eigen3 variation and using eigen2/3 in generic code)

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

2010/1/4 Paul C. Leopardi <paul.leopardi@xxxxxxxxxxxx>:
> On Monday 04 January 2010 12:59:27 Benoit Jacob wrote:
>> Let me explain the Eigen approach to this problem. We try to let the user
>>  avoid as much as possible writing his own loops iterating over the
>>  coefficients of a dense matrix. The reason is because it is extremely
>> complicated for the user to write such loops in an efficient way.
>
>> As I said, the user can provide his own expressions and functors, so the
>> only limitation is about the built-in loop types. We have essentially
>> two loop types:
>>  - assignment
>>  - reduction
>>
>> For example, the above sum() is reduction with a "sum" functor. Again, the
>>  user can provide his own functors here.
>
> Benoit,
> Thanks for the explanation. I have two more examples where I don't understand
> what the Eigen approach would be:
>
> 1. Left inverse of Kronecker product. nork(A, kron(A, B)) is close to B.
>
>  /// Left inverse of Kronecker product
>  template< typename LHS_T, typename RHS_T >
>  const
>  RHS_T
>  nork(const LHS_T& lhs, const RHS_T& rhs, const bool mono)
>  {
>    // nork(A, kron(A, B)) is close to B
>    typedef RHS_T matrix_t;
>    typedef typename LHS_T::size_type lhs_index_t;
>    typedef typename matrix_t::size_type matrix_index_t;
>    const lhs_index_t lhs_s1 = lhs.size1();
>    const lhs_index_t lhs_s2 = lhs.size2();
>    const matrix_index_t rhs_s1 = rhs.size1();
>    const matrix_index_t rhs_s2 = rhs.size2();
>    const matrix_index_t res_s1 = rhs_s1 / lhs_s1;
>    const matrix_index_t res_s2 = rhs_s2 / lhs_s2;
>    typedef typename matrix_t::value_type Scalar_T;
>    const Scalar_T& nnz_lhs = Scalar_T( double(mono ? lhs_s1 : nnz(lhs)) );
>    if (!mono)
>    {
>      typedef error<matrix_t> error_t;
>      if (rhs_s1 == 0)
>        throw error_t("matrix", "nork: number of rows must not be 0");
>      if (rhs_s2 == 0)
>        throw error_t("matrix", "nork: number of cols must not be 0");
>      if (res_s1 * lhs_s1 != rhs_s1)
>        throw error_t("matrix", "nork: incompatible numbers of rows");
>      if (res_s2 * lhs_s2 != rhs_s2)
>        throw error_t("matrix", "nork: incompatible numbers of cols");
>      if (nnz_lhs == Scalar_T(0))
>        throw error_t("matrix", "nork: LHS must not be 0");
>    }
>    matrix_t result(res_s1, res_s2);
>    result.clear();
>    for (typename LHS_T::const_iterator1
>        lhs_it1 = lhs.begin1();
>        lhs_it1 != lhs.end1();
>        ++lhs_it1)
>      for (typename LHS_T::const_iterator2
>          lhs_it2 = lhs_it1.begin();
>          lhs_it2 != lhs_it1.end();
>          ++lhs_it2)
>        if (*lhs_it2 != Scalar_T(0))
>        {
>          typedef ublas::matrix_range<const matrix_t> matrix_range_t;
>          const Scalar_T& lhs_scale = *lhs_it2 * nnz_lhs;
>          const matrix_index_t start1 = res_s1 * lhs_it2.index1();
>          const matrix_index_t start2 = res_s2 * lhs_it2.index2();
>          using ublas::range;
>          const range& range1 = range(start1, start1 + res_s1);
>          const range& range2 = range(start2, start2 + res_s2);
>          const matrix_range_t& rhs_range =
>                   matrix_range_t(rhs, range1, range2);
>          for (typename matrix_range_t::const_iterator1
>              rhs_it1 = rhs_range.begin1();
>              rhs_it1 != rhs_range.end1();
>              ++rhs_it1)
>            for (typename matrix_range_t::const_iterator2
>                rhs_it2 = rhs_it1.begin();
>                rhs_it2 != rhs_it1.end();
>                ++rhs_it2)
>                result(rhs_it2.index1(), rhs_it2.index2()) +=
>                        *rhs_it2 / lhs_scale;
>        }
>    return result;
>  }

So, I am already convinced by the std::binary_search example of the
usefulness of STL algorithms, but just for reference, here's how you
write this with Eigen for dense matrices:

matrix_t result = matrix_t::Zero(res_rows, res_cols);
for(int row = 0; row < lhs.rows(); ++row)
{
  for(int col = 0; col < lhs.cols(); ++col)
  {
    Scalar_T x = lhs(row,col);
    if(x != Scalar_T(0))
    {
      const Scalar_T lhs_scale = x * nnz_lhs;
      const int start_row = res_rows * row;
      const int start_col = res_cols * col;
      result += rhs.block(start_row, start_col, res_rows, res_cols) / lhs_scale;
    }
  }
}

As you can see: the two innermost for loops have been replaced by a
single Eigen expression assignment. Eigen will then easily vectorize
and/or unroll that.

For sparse matrices, of course Eigen already works with iterators.

I agree that the two outermost loops would benefit from using
iterators over nonzero elements, so this code applies readily to both
dense and sparse matrices.

>
> 2. Product of monomial matrices. Called with RHS_T a sparse matrix type.
>
>  /// Product of monomial matrices
>  template< typename LHS_T, typename RHS_T >
>  const typename RHS_T::expression_type
>  mono_prod(const ublas::matrix_expression<LHS_T>& lhs,
>            const ublas::matrix_expression<RHS_T>& rhs)
>  {
>    typedef const LHS_T lhs_expression_t;
>    typedef const RHS_T rhs_expression_t;
>    typedef typename RHS_T::expression_type matrix_t;
>    typedef typename matrix_t::size_type  matrix_index_t;
>    typedef typename matrix_t::value_type Scalar_T;
>    typedef typename lhs_expression_t::const_iterator1   lhs_const_iterator1;
>    typedef typename lhs_expression_t::const_iterator2   lhs_const_iterator2;
>    typedef typename ublas::matrix_row<rhs_expression_t> matrix_row_t;
>    typedef typename matrix_row_t::const_iterator        row_const_iterator;
>
>    const matrix_index_t dim = lhs().size1();
>    // The following assumes that matrix_t is a sparse matrix type.
>    matrix_t result = matrix_t(dim, dim, dim);
>    for (lhs_const_iterator1
>        lhs_row = lhs().begin1();
>        lhs_row != lhs().end1();
>        ++lhs_row)
>    {
>      const lhs_const_iterator2& lhs_it = lhs_row.begin();
>      if (lhs_it != lhs_row.end())
>      {
>        const matrix_row_t& rhs_row = matrix_row_t(rhs(), lhs_it..index2());
>        const row_const_iterator& rhs_it = rhs_row.begin();
>        if (rhs_it != rhs_row.end())
>          result(lhs_it.index1(), rhs_it.index()) = (*lhs_it) * (*rhs_it);
>      }
>    }
>    return result;
>  }

Eigen doesn't do 3D matrices at the moment, although the soon-coming
Array class can be used as a Scalar type so you could use a Matrix of
Arrays.

>
> To use these routines with Eigen matrices, should I try to wrap the Eigen
> classes into classes which have the required iterators, or is there a way to
> use the built-in loop types with these routines?

I understand that we need to add iterators to cover your needs, so
let's do that.

Benoit