собственный автоматический вывод типа в общем продукте

У меня есть следующий фрагмент кода (прошу прощения за немного больший фрагмент кода, это минимальный пример, к которому я смог уменьшить свою проблему):

#include <Eigen/Dense>
#include <complex>
#include <iostream>
#include <typeinfo>

// Dynamic Matrix over Scalar field
template <typename Scalar> 
using DynMat = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;

// Dynamic column vector over Scalar field
template <typename Scalar>
using DynVect = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;

// Returns the D x D Identity matrix over the field Derived::Scalar
// deduced from the expression Eigen::MatrixBase<Derived>& A
template<typename Derived>
DynMat<typename Derived::Scalar> Id(const Eigen::MatrixBase<Derived>& A, std::size_t D)
{   
    DynMat<typename Derived::Scalar> result =
            DynMat<typename Derived::Scalar>::Identity(D, D);

    return result;
}

int main()
{
    //using ScalarField = std::complex<double>; // same issue even if I use complex numbers
    using ScalarField = double; // we use doubles in this example

    // A double dynamic matrix (i.e. MatrixXd)
    DynMat<ScalarField> Foo; // used to deduce the type in Id<>()

    // A double dynamic column vector (i.e. VectorXd)
    DynVect<ScalarField> v(4);
    v << 1., 0. , 0. ,0.; // plug in some values into it

    // Make sure that Id(Foo, 4) correctly deduces the template parameters
    std::cout << "Id(Foo, 4) is indeed the 4 x 4 identiy matrix over the ScalarField of "
              << "typeid().name(): " << typeid(ScalarField).name() << std::endl;
    std::cout << Id(Foo, 4) << std::endl; // Indeed the 4 x 4 complex Identity matrix

    // Use auto type deduction for GenMatProduct, junk is displayed. Why?!
    std::cout << std::endl << "Use auto type deduction for GenMatProduct,\
                 sometimes junk is displayed. Why?!" << std::endl;
    auto autoresult = Id(Foo, 4) * v; // evaluated result must be identically equal to v
    for(int i=0; i<10; i++)
    {
            std::cout << autoresult.transpose(); // thought 1 0 0 0 is the result, but NO, junk
            std::cout << " has norm: " << autoresult.norm() << std::endl; // junk
    }

    // Use implicit cast to Dynamic Matrix, works fine
    std::cout << std::endl << "Use implicit cast to Dynamic Matrix, works fine" << std::endl;
    DynMat<ScalarField> castresult = Id(Foo, 4) * v; // evaluated result must be identically equal to v
    for(int i=0; i<10; i++)
    {
            std::cout << castresult.transpose(); // 1 0 0 0, works ok
            std::cout << " has norm: " << castresult.norm() << std::endl; // ok
    }
}

Основная идея заключается в том, что шаблонная функция Id<>() принимает собственное выражение A в качестве параметра вместе с размером D и создает единичную матрицу над скалярным полем выражения A. Эта функция сама по себе работает нормально. Однако, когда я использую его в продукте Eigen с выведенным типом auto, например, в строке auto autoresult = Id(Foo, 4) * v, я ожидаю умножения вектора v на единичную матрицу, поэтому чистый результат должен быть выражением, которое при оценке должно быть тождественно равно v. Но это не так, см. первый цикл for, всякий раз, когда я отображаю результат и вычисляю его норму, большую часть времени получаю мусор. Если, с другой стороны, я неявно привожу собственное произведение Id(Foo, 4) * v к динамической матрице, все работает нормально, результат оценивается правильно.

Я использую Eigen 3.2.2 в OS X Yosemite и получаю одинаковое странное поведение как с g++ 4.9.1, так и с Apple LLVM версии 6.0 (clang-600.0.54) (на основе LLVM 3.5svn)

ВОПРОС:

• Я не понимаю, что происходит в первом цикле for, почему продукт не оценивается, когда я использую std::cout или даже когда я использую метод norm? Я что-то пропустил? Здесь нет алиасинга, и я действительно озадачен тем, что происходит. Я знаю, что Eigen использует ленивую оценку и вычисляет выражение, когда это необходимо, но здесь, похоже, это не так. Эта проблема чрезвычайно важна для меня, так как у меня есть много функций того же вкуса, что и Id<>(), которые при использовании в auto выводимых выражениях могут дать сбой.

Проблема возникает довольно часто, но не всегда. Однако, если вы запустите программу 3-4 раза, вы обязательно ее увидите.

Команда, которую я использую для компиляции и запуска:

clang++ (g++) -std=c++11 -isystem ./eigen_3.2.2/ testeigen.cpp -otesteigen; ./testeigen

Типичный вывод, который я получил в реальном прогоне:

Id(Foo, 4) is indeed the 4 x 4 identiy matrix over the ScalarField of typeid().name(): d
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

Use GenMatProduct, sometimes junk is displayed. Why?!
1 0 0 0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf
3.10504e+231 3.10504e+231 3.95253e-323            0 has norm: inf

Use implicit cast to Dynamic Matrix, works fine
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1
1 0 0 0 has norm: 1

Даже если я использую eval() в

  std::cout << autoresult.eval().transpose(); // thought 1 0 0 0 is the result, but NO, junk
  std::cout << " has norm: " << autoresult.eval().norm() << std::endl; // junk

Я получаю такое же странное поведение.