title: eigen: expression or function to set/return full column
body: <p>I have several instances in my code, where I have a condition based on coefficients of 1xN arrays, and need to set whole columns of MxN arrays depending on these conditions. In my case, N is <code>Eigen::Dynamic</code> and M ranges from 2 to 4, but is a compile-time constant in each instance.</p>
  <p>Here's a simple function illustrating what I mean, with <code>a</code> and <code>b</code> being the 1xN arrays which form the condition, <code>c</code> being a 2xN array with additional data, and <code>res</code> being an out-parameter, whose columns are always set as a whole:</p>
  <pre><code>#include &lt;iostream&gt;
  #include &lt;Eigen/Dense&gt;
  
  using namespace Eigen;
  
  template&lt;Index nRows&gt;
  using ArrayNXd = Array&lt;double, nRows, Dynamic&gt;;
  
  using Array1Xd = ArrayNXd&lt;1&gt;;
  using Array2Xd = ArrayNXd&lt;2&gt;;
  using Array3Xd = ArrayNXd&lt;3&gt;;
  
  void asFunction(
      Array3Xd&amp; res,
      const Array1Xd&amp; a, const Array1Xd&amp; b, const Array2Xd&amp; c
  ){
      for (Index col{0}; col&lt;a.cols(); ++col){
          if ( a[col] &gt; b[col] )
              res.col(col) = Array3d{
                   a[col] + b[col],
                  (a[col] + b[col]) * c(0, col),
                  (a[col] - b[col]) * c(1, col)
              };
          else
              res.col(col) = Array3d{
                   a[col] - b[col],
                   a[col] + b[col],
                  (a[col] + b[col]) * (a[col] - b[col])
              };
      }
  }
  
  
  int main(){
      Array1Xd a (3), b(3);
      Array2Xd c (2, 3);
      
      a &lt;&lt; 1, 2, 3;
      b &lt;&lt; 0, 1, 2;
      c &lt;&lt;
          0, 1, 2,
          1, 2, 3;
  
      Array3Xd res (3,3);
      
      asFunction(res, a, b, c);
  
      std::cout &lt;&lt; &quot;as function:\n&quot; &lt;&lt; res &lt;&lt; &quot;\n&quot;;
  
      return 0;
  }
  </code></pre>
  <p>Functions similar to this are used in a performance cricital section* of my code, and I feel like I'm leaving performance on the table, because using loops with <code>Eigen</code> types is typically not the optimal solution.</p>
  <p>*yes, I profiled it.</p>
  <p>I wrote the same function as a <code>NullaryExpr</code>, but that was a bit slower. I guess that makes sense, given the additional evaluations of the condition(s) and the branching for each row:</p>
  <pre><code>#include &lt;iostream&gt;
  #include &lt;Eigen/Dense&gt;
  
  using namespace Eigen;
  
  template&lt;Index nRows&gt;
  using ArrayNXd = Array&lt;double, nRows, Dynamic&gt;;
  
  using Array1Xd = ArrayNXd&lt;1&gt;;
  using Array2Xd = ArrayNXd&lt;2&gt;;
  using Array3Xd = ArrayNXd&lt;3&gt;;
  
  class MyFunctor
  {
  public:
      using Scalar = double;
  
      static constexpr Index
          RowsAtCompileTime { 3 },
          MaxRowsAtCompileTime { 3 },
          ColsAtCompileTime { Dynamic },
          MaxColsAtCompileTime { Dynamic };
  
      using DenseType = Array&lt;
          Scalar  ,    RowsAtCompileTime,    ColsAtCompileTime,
          ColMajor, MaxRowsAtCompileTime, MaxColsAtCompileTime
      &gt;;
  
  private:
      typename Array1Xd::Nested m_a;
      typename Array1Xd::Nested m_b;
      typename Array2Xd::Nested m_c;
  
  public:
      MyFunctor(
          const Array1Xd&amp; a,
          const Array1Xd&amp; b,
          const Array2Xd&amp; c
      ) : m_a {a}, m_b {b}, m_c{c}
      {}
  
      bool cond(Index col) const {
          return m_a[col] &gt; m_b[col];
      }
  
      Scalar func1(Index col) const {
          return m_a[col] + m_b[col];
      }
  
      Scalar func2(Index col) const {
          return m_a[col] - m_b[col];
      }
  
      Scalar func3(Index row, Index col) const {
          switch(row){
              case 0: return func1(col);
              case 1: return func1(col) * m_c(0, col);
              case 2: return func2(col) * m_c(1, col);
              default: __builtin_unreachable();
          }
      }
  
      Scalar func4(Index row, Index col) const {
          switch (row){
              case 0: return func2(col);
              case 1: return func1(col);
              case 2: return func1(col) / func2(col);
              default: __builtin_unreachable();
          }
      }
  
      Scalar operator() (Index row, Index col) const {
          if ( cond(col) )
              return func3(row, col);
          else
              return func4(row, col);
      }
  };
  
  using MyReturnType = Eigen::CwiseNullaryOp&lt;
      MyFunctor, typename MyFunctor::DenseType
  &gt;;
  
  
  MyReturnType asFunctor(
      const Array1Xd&amp; a,
      const Array1Xd&amp; b,
      const Array2Xd&amp; c
  ){
      using DenseType = typename MyFunctor::DenseType;
      return DenseType::NullaryExpr(
          3, a.cols(),
          MyFunctor(a, b, c)
      );
  }
  
  
  int main(){
      Array1Xd a (3), b(3);
      Array2Xd c (2, 3);
      
      a &lt;&lt; 1, 2, 3;
      b &lt;&lt; 0, 1, 2;
      c &lt;&lt;
          0, 1, 2,
          1, 2, 3;
  
      std::cout &lt;&lt; &quot;as functor:\n&quot; &lt;&lt; asFunctor(a,b,c) &lt;&lt; &quot;\n&quot;;
  
      return 0;
  }
  </code></pre>
  <p>My question is: <strong>Is there a more efficient way of implementing logic similar to the above</strong> (evaluate scalar condition for each column of a matrix, return values for the whole column based on the condition) <strong>using the <code>eigen</code> library?</strong></p>
  <p>Note: using an expression would be slightly preferred, because I don't need to worry about memory allocation, out-parameters, etc., and the code can be written with scalars in mind, which makes it much more easily understandable.</p>
  <p>Edit: Note2: I tried using <code>&lt;Condition&gt;.template replicate&lt;nRows,1&gt;().select(..., ...)</code> as well, but it was slower and harder to read.</p>