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

["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;
  }

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;
  }

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?
Best, Paul