Given a band matrix B, I would like to extract its rows, one at a time.
The monster I've come up with is this:

Code: Select all

B.coeffs().rowwise().reverse().diagonal(k).transpose().reverse();

which for a well-chosen value of k ought to return one of the rows.

My questions are:

1. Is this the best way of doing this?
2. If so, how do I come up with k? It seems like a complicated function of 4 variables (rows, cols, supers, subs).
3. Do I need to start completely from scratch to set one of the rows? I can't seem to get reverse to return a mutable view.

Thank you,

krey

For 2. I found this: http://www.netlib.org/lapack/lug/node124.html

We currently don't have better to offer. Implementing a generic AntiDiagonal expression (similar to Diagonal) would help. Regarding the write issue, the problem is that there is no con-const reverse method in VectorwiseOp, easily fixable:

Code: Select all

diff --git a/Eigen/src/Core/VectorwiseOp.h b/Eigen/src/Core/VectorwiseOp.h
--- a/Eigen/src/Core/VectorwiseOp.h
+++ b/Eigen/src/Core/VectorwiseOp.h
@@ -279,16 +279,17 @@ template<typename ExpressionType, int Di
     typedef typename ReturnType<internal::member_stableNorm,RealScalar>::Type StableNormReturnType;
     typedef typename ReturnType<internal::member_hypotNorm,RealScalar>::Type HypotNormReturnType;
     typedef typename ReturnType<internal::member_sum>::Type SumReturnType;
     typedef typename ReturnType<internal::member_mean>::Type MeanReturnType;
     typedef typename ReturnType<internal::member_all>::Type AllReturnType;
     typedef typename ReturnType<internal::member_any>::Type AnyReturnType;
     typedef PartialReduxExpr<ExpressionType, internal::member_count<Index>, Direction> CountReturnType;
     typedef typename ReturnType<internal::member_prod>::Type ProdReturnType;
+    typedef Reverse<const ExpressionType, Direction> ConstReverseReturnType;
     typedef Reverse<ExpressionType, Direction> ReverseReturnType;

     template<int p> struct LpNormReturnType {
       typedef PartialReduxExpr<ExpressionType, internal::member_lpnorm<p,RealScalar>,Direction> Type;
     };

     /** \returns a row (or column) vector expression of the smallest coefficient
       * of each column (or row) of the referenced expression.
@@ -451,17 +452,25 @@ template<typename ExpressionType, int Di
     /** \returns a matrix expression
       * where each column (or row) are reversed.
       *
       * Example: \include Vectorwise_reverse.cpp
       * Output: \verbinclude Vectorwise_reverse.out
       *
       * \sa DenseBase::reverse() */
     EIGEN_DEVICE_FUNC
-    const ReverseReturnType reverse() const
+    const ConstReverseReturnType reverse() const
+    { return ConstReverseReturnType( _expression() ); }
+
+    /** \returns a writable matrix expression
+      * where each column (or row) are reversed.
+      *
+      * \sa reverse() const */
+    EIGEN_DEVICE_FUNC
+    ReverseReturnType reverse()
     { return ReverseReturnType( _expression() ); }

     typedef Replicate<ExpressionType,(isVertical?Dynamic:1),(isHorizontal?Dynamic:1)> ReplicateReturnType;
     EIGEN_DEVICE_FUNC
     const ReplicateReturnType replicate(Index factor) const;

     /**
       * \return an expression of the replication of each column (or row) of \c *this


Fixed in devel branch: https://bitbucket.org/eigen/eigen/commits/d4db25551a78/

Many thanks for taking care of this, Gael.

For the record, here is the correct code for getting the nth row

Code: Select all

template <typename T>
auto banded_row (BandMatrix<T> & B, const uint n) {
    const uint offset = B.cols()-B.supers()-1;
    return B.coeffs().rowwise().reverse().diagonal(offset-n).transpose().reverse();
}

Please let me know if you can think of something more elegant/faster etc.