Eigen: selfadjointView не работает, если вектор innerIndexPtr неупорядочен - PullRequest
Новая
       89

Eigen: selfadjointView не работает, если вектор innerIndexPtr неупорядочен

0 голосов
/ 16 июня 2020

Работая в системе реального времени, я использую низкоуровневый API для управления разреженной матрицей Eigen, чтобы повысить производительность на определенном c наборе задач. В частности, я смешиваю Eigen с вариантом алгоритма CSparse cs_symperm , чтобы получить перестановку матрицы A: Ap = PAP '. cs_symperm использовать ту же структуру низкого уровня для разреженной матрицы, но для фиксированного столбца индексы строк могут быть неупорядочены.

Это простой пример, построенный вручную, того, что могло произойти: 

SparseMatrix<double> A( 2, 2 );
A.insert( 0, 0 ) = 1.0;
A.insert( 0, 1 ) = 2.0;
A.insert( 1, 1 ) = 3.0;
A.makeCompressed();
SparseMatrix<double> B = A;
B.innerIndexPtr()[0] = 0;
B.innerIndexPtr()[1] = 1; // <-- Not ordered
B.innerIndexPtr()[2] = 0; // <-- Not ordered
B.valuePtr()[0] = 1.0;
B.valuePtr()[1] = 3.0; // <-- Not ordered
B.valuePtr()[2] = 2.0; // <-- Not ordered

Здесь A и B - одна и та же матрица. Единственное отличие - порядок данных.

Матрично-векторное произведение правильно работает: 

VectorXd x( 2 );
x << 1.0, 2.0;
VectorXd y = A * x;
VectorXd w = B * x;
assert( y( 0 ) == w( 0 ) ); // <-- OK
assert( y( 1 ) == w( 1 ) ); // <-- OK

selfadjointView не работает: 

y = A.selfadjointView<Upper>() * x;
w = B.selfadjointView<Upper>() * x;
assert( y( 0 ) == w( 0 ) ); // <-- Fail!

Пример в документации Eigen (https://eigen.tuxfamily.org/dox/group__TutorialSparse.html) показывает упорядоченные данные, но нет явного указания.

К сожалению, я не могу получить Ap с помощью Eigen из-за динамического c распределения временных библиотек объекты. Есть идеи?

Тест проводился с использованием Eigen v3.3.7.

1 Ответ

2 голосов
/ 16 июня 2020

Чтобы отсортировать элементы матрицы, вы можете дважды транспонировать ее: 

B = B.transpose(); B = B.transpose();

, или вы можете транспонировать ее один раз и использовать selfadjointView<Lower>(), или вы можете назначить ее матрице-строке (это также неявно транспонирует его): 

SparseMatrix<double, RowMajor> C = B;

w = C.selfadjointView<Upper>() * x;
...