c ++ Собственная копия?

Question

c ++ Собственная копия?

Пусть A — симметричная матрица, а v — вектор. Я извлекаю из A блок из n столбцов, начиная с j, и умножаю его на v, используя

VectorXd a;
a = A.middleCols(j,n).selfadjointView<Lower>() * v // does not compile

Как это не компилируется, тогда как это

a = MatrixXd(A.middleCols(j,n).selfadjointView<Lower>()) * v

интересно, делает ли вторая версия копию

A.middleCols(j,n).selfadjointView<Lower>()

или выполняет вычисления напрямую?

Спасибо за любую подсказку.

РЕДАКТИРОВАТЬ: Я подозреваю, что проблема как-то связана с типами аргументов, так как я получаю сообщение об ошибке:

invalid argument type 'typename ConstSelfAdjointViewReturnType.... to unary expression'

Действительно, A является аргументом функции, переданной по константной ссылке с использованием любого из

const MatrixXd& A
const Ref<const MatrixXd>& A

Вот пример:

// this version doesn't compile
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = A.middleCols(3, 4).selfadjointView<Lower>() * v;
MatrixXd M(code_fill(){...});
// code_fill is the function filling the lower triangular part of a symmetric matrix
M.block(1, 2, 3, 4).triangularView<Lower>() += B.selfadjointView<Lower>();

return M;
}

// this version compiles
MatrixXd returnSymmetricMatrix(const MatrixXd& A, const VectorXd& v, const MatrixXd& B){
// B is a symmetric matrix

VectorXd a;
a = MatrixXd(A.middleCols(3, 4).selfadjointView<Lower>()) * v;
MatrixXd M(code_fill(){...});
// code_fill is the function filling the lower triangular part of a symmetric matrix
Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

return M;
}

EDIT2 Что касается моего первоначального вопроса и примера, который я добавил в разделе «Правка», я немного запутался в отношении копирования. Как я понимаю разницу между рабочей и нерабочей версиями, строка

Matrix(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

работает, потому что его lhs сообщает Eigen, что M.block (1, 2, 3, 4) .triangularView () на самом деле является матрицей, а не ссылкой на матрицу. В противном случае оператор + = из-за ошибки не будет перегружен для .block (). Таким образом, мой первоначальный вопрос заключается в том, говорит ли Matrix (…), что это Matrix, чтобы включить вычисления, или, скорее, скопировать … в Matrix? Спасибо!

0

c++eigen symmetric

Решение

Другие решения

Других решений пока нет …

Источник

Accepted Answer

Следующее выражение:

A.middleCols(j,n).selfadjointView<Lower>()

не создает никакой копии.

С другой стороны, чтобы избежать временного для результата продукта, можно добавить .noalias():

a.noalias() = M.middleCols(j,n).selfadjointView<Lower>() * v;

Это необходимо только для немедленного назначения плотного продукта, чтобы разрешить такой код:

a = M * a;

вести себя как ожидалось.

РЕДАКТИРОВАТЬ:

Что касается вашей проблемы компиляции, следующие компилируется нормально:

#include <Eigen/Dense>
using namespace Eigen;
int main()
{
int n = 10;
MatrixXd M = MatrixXd::Random(n,2*n);
VectorXd v = VectorXd::Random(n);
VectorXd a;
a.noalias() = M.middleCols(2,n).selfadjointView<Upper>() * v;

const Ref<const MatrixXd>& A = M;
a.noalias() = A.middleCols(2,n).selfadjointView<Upper>() * v;
}

EDIT2

Следующая строка:

MatrixXd(M.block(1, 2, 3, 4).triangularView<Lower>()) += B.selfadjointView<Lower>();

не имеет смысла, поскольку вы создаете временную копию, которую вы назначаете. Напомним, что здесь MatrixXd(whatever) вызывает Eigen::Matrix конструктор. Также не имеет смысла присваивать вид самосопряженного треугольному виду. Я даже не могу думать о том, что может быть разумным поведением для этого. Если вы хотите обновить только нижнюю треугольную часть, выполните:

M.block(1, 2, 3, 4).triangularView<Lower>() += B;

В этом случае, triangularView ведет себя как пишущая маска.

1