Hi,

I am working on an expression template library for CFD applications. The aim is to support multiple existing vector expression template libraries as backends and write the code in a most generic way. However, Eigen (version 3.2. seems to have some problems when used in combination with the auto return type. The following demo code illustrates the problem. In fact, using auto return type to simply return one of the parameters (but not modify it) works well. However, using auto return type together with arithmetics (see last line of the main routine, where the velocity is computed by dividing the momentum by the density) leads to

Code: Select all

result =
-3.22057e-232
-3.22057e-232
          inf
          inf
          inf
          inf
          inf
          inf
          inf
 1.43815e+308


Any help on combining the Eigen library with nested auto return types and arithmetics is highly appreciated.

Kind regards,
Matthias Möller

Code: Select all

namespace fdbb_mini
{

/// forward declaration
template<int iarg, class ... Ts>
static inline auto getArg(const Ts& ... args)
    -> typename std::tuple_element<iarg, std::tuple<Ts...> >::type;

/// Specialization for iarg==0
template<int iarg, class T0, class ... Ts>
static inline auto getArg(const T0& arg0, const Ts& ... args)
    -> typename std::enable_if<iarg==0, typename std::tuple_element<iarg, std::tuple<T0,Ts...> >::type>::type
{
    return arg0;
}
                   
/// Specialization for iarg!=0                                                             
template<int iarg, class T0, class ... Ts>
static inline auto getArg(const T0& arg0, const Ts& ... args)
    -> typename std::enable_if<iarg!=0, typename std::tuple_element<iarg-1, std::tuple<Ts...> >::type>::type
{
    return getArg<iarg-1>(args...);
}

// Density
template<class ... Ts>
static inline auto density(const Ts& ... vars)
    -> decltype(getArg<0>(vars...))
{
    return getArg<0>(vars...);
}

// Momentum
template<int dim,class ... Ts>
static inline auto momentum(const Ts& ... vars)
    -> typename std::enable_if<dim==0, decltype(getArg<1>(vars...))>::type
{
    return getArg<1>(vars...);
}

// Energy
template<class ... Ts>
static inline auto energy(const Ts& ... vars)
    -> decltype(getArg<2>(vars...))
{
    return getArg<2>(vars...);
}

// Velocity
template<int dim,class ... Ts>
static inline auto velocity(const Ts& ... vars)
    -> typename std::enable_if<dim==0, decltype(momentum<0>(vars...)/density(vars...))>::type
{
    return momentum<0>(vars...)/density(vars...);
}

}

int main(int argc, char** argv){

    const std::size_t N=10;
    Eigen::Array<double,Eigen::Dynamic,1> density(N), momentum(N), energy(N), result(N);
    density.fill(2.0);
    momentum.fill(1.0);
    energy.fill(5.0);
    result.fill(0.0);
   
    // Working examples
    result = fdbb_mini::getArg<0>(density, momentum, energy);
    std::cout << result << std::endl;
    result = fdbb_mini::density(density, momentum, energy);
    std::cout << result << std::endl;

    result = fdbb_mini::getArg<1>(density, momentum, energy);
    std::cout << result << std::endl;
    result = fdbb_mini::momentum<0>(density, momentum, energy);
    std::cout << result << std::endl;
   
    result = fdbb_mini::getArg<2>(density, momentum, energy);
    std::cout << result << std::endl;
    result = fdbb_mini::energy(density, momentum, energy);
    std::cout << result << std::endl;

    // Non-working example
    result = fdbb_mini::velocity<0>(density, momentum, energy, energy);
    std::cout << result << std::endl;
   
    return 0;
}


This is because your momentum() and density() functions returns a copy of the input arrays which are then stored by reference by the quotient expression. Then the copies are deleted before the assignment to result takes place... The solution is to make sure your accessors returns by reference to the original objects.