Комбинация умножения матрицы-вектора с коэффициентом по собственному значению

Question

Я пытаюсь сконденсировать немного кода, используя функции в Eigen, в частности те, которые позволяют вам делать кодовые продукты массивов, но, возможно, используя их неправильно. Основная идея содержится в этом расширен цикл:

for (int dir = 0; dir < NDIM; dir++) 
    { 
    for (int i = 0; i < nlocal; i++) 
    { 
     for (int qp = 0; qp < nVolQuad; qp++) 
     qNewPtr[i] += 
      volQuad.weights(qp)*bigStoredVolMatrices[dir](i,qp)*alpha(dir,qp)*fAtQuad(qp); 
    } 
    } 

, который я пытаюсь конденсируются:

for (int dir = 0; dir < NDIM; dir++) 
    { 
    resultVectorDir[dir].noalias() = bigStoredVolMatrices[dir]* 
     (volQuad.weights.array()*fAtQuad.array()*alpha.row(dir).array()).matrix(); 
    } 

    for (int i = 0; i < nlocal; i++) 
    { 
    for (int dir = 0; dir < NDIM; dir++) 
     qNewPtr[i] += resultVectorDir[dir](i); 
    } 

или без переключения между матрицей и матрицей, используя что-то вроде:

for (int dir = 0; dir < NDIM; dir++) 
    { 
    resultVectorDir[dir].noalias() = bigStoredVolMatrices[dir]* 
     (volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir)))); 
    } 

    for (int i = 0; i < nlocal; i++) 
    { 
    for (int dir = 0; dir < NDIM; dir++) 
     qNewPtr[i] += resultVectorDir[dir](i); 
    } 

Странная вещь для меня иногда работает. Иногда код генерирует желаемый результат, но иногда также генерирует NaN. Я думал, что мне может потребоваться явно обнулить resultVectorDir в сжатой версии, но это не решило проблему. Я подумал, что, возможно, есть что-то тонкое в выполнении этого порядка операций? Любая помощь, которая может быть предоставлена, будет с большой благодарностью.

В качестве дополнения к этому вопросу я напал на эту проблему с помощью добрых старомодных утверждений печати и обнаружил, что я не должен правильно использовать коэффициент полезного действия функции массивов. Например, в ситуации, когда nVolQuad = 9, я побежал в этот сегмент кода:

for (int dir = 0; dir < NDIM; dir++) 
    { 
    tempArray = volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir))); 
    for (int qp = 0; qp < nVolQuad; qp++) 
    { 
     std::cout << std::endl; 
     std::cout << tempArray(qp) << " "; 
     std::cout << volQuad.weights(qp)*alpha(dir,qp)*fAtQuad(qp) << " "; 
     std::cout << std::endl; 
    } 
    std::cout << std::endl; 
    std::cout << std::endl; 
    } 

Кусок вывода выглядит следующим образом:

-2.23064e-05 -2.23064e-05

1.49458e-154 -3.56902e-05

6.94729e-310 -2.23064e-05

-2.68156e + 154 -2.0672e-05

6.94729e-310 -3.30752e-05

6.94729e-310 -2.0672e-05

2.13645e-314 -2.99114e-06

6.94729e-310 -4.78582e-06

6.94729e-310 -2.99114e-06

Другие сегменты вывода похожи. Первая запись верна, но остальные 8 записей tempArray нонсенс. tempArray инициализируется до 0,0 до цикла, поэтому я немного теряю. Я продолжу копаться в документации Эйгена, чтобы убедиться, что я не делаю что-то чрезвычайно глупо в использовании этой функции.

Answer 1

Ключевая ошибка заключается в том, что вы можете немедленно перейти от расширенного цикла к средству с коэффициентом, когда расширенный цикл обращается к данным разных размеров. Здесь:

volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir))) 

С volQuad.weights и fAtQuad инициализируются Эйгеном :: VectorXd, они являются векторы-столбцы, но при использовании alpha.row (DIR), что конкретная структура данных представляет собой вектор-строка. Таким образом, коэффициент полезного продукта не имеет смысла, и вы только получите правильную запись.Это может быть легко исправлено путем изменения синтаксиса:

volQuad.weights.cwiseProduct(fAtQuad.cwiseProduct(alpha.row(dir).transpose()))