Следуя тексту метода конечных элементов Ценкевича, я решаю проблему линейной эластичности со связанными подобластями. Короче говоря, мне нужно вычислить следующий продукт тройной разреженной матрицы:
Q * V^-1 * Q^T
, где Q - это разреженная матрица размером n на m, а V - это симметрия размером m на m * 1034. * положительно определенная разреженная матрица. Как правило, n >> m. В настоящее время следующий код дает правильные результаты, но медленен относительно моих ожиданий:
using namespace Eigen;
typedef double Scalar_t;
typedef SparseMatrix<Scalar_t,ColMajor> Sparse_t;
int n, m;
Sparse_t result(n,n);
Sparse_t V(m,m);
Sparse_t Q(n,m);
//begin
SimplicialLDLT<Sparse_t> solver(V);
assert(solver.info() == Success);
//this is required because all my sparse matrices are column-major
Sparse_t QT(Q.transpose());
result = QT.transpose() * solver.solve(QT);
Учитывая разреженность, симметрию и общий идеальный символ Q и V, я ожидаю, что лучший способ сделать это. Одна идея, которую я имею, состоит в том, чтобы вместо этого использовать декомпозицию LLT.
V = L * L^T
V^-1 = L^-T * L^-1
Q * V^-1 * Q^T = (Q * L^-T) * (L^-1 * Q^T)
Поскольку (Q * L^-T)^T = L^-1 * Q^T Я должен иметь возможность использовать обновление ранга симметрии c для конечного матричного продукта. Это работает, но даже медленнее, чем описанный выше метод!
SimplicialLLT<Sparse_t> V_solver(V);
assert(V_solver.info() == Success);
Sparse_t L(V_solver.matrixL());
SparseLU<Sparse_t> L_solver(L);
assert(L_solver.info() == Success);
Sparse_t QT(Q.transpose());
Sparse_t L_inv_QT = L_solver.solve(QT);
result.selfadjointView<Lower>().rankUpdate(L_inv_QT.transpose());
С точки зрения линейной алгебры, метод LLT должен включать меньше работы и использовать все приятные свойства моих матриц. Я хотел использовать решатель triangularview, то есть: V_solver.matrixL().solve(QT), но он принимает только аргументы с плотной правой частью. Версия решателя представления tri angular не работает, вероятно, потому, что мой аргумент rhs является прямоугольным angular (может включать изменение размеров матрицы).
Есть идеи? Это узкое место моего приложения. Спасибо!
РЕДАКТИРОВАТЬ:
По сообщению Ави Гинсберга Понимание операции solveInPlace в Eigen
Я не использовал solveInPlace правильно. Это позволяет мне напрямую вычислять L^-1 * Q^T с помощью решателя tri angular и выполнять обновление ранга симметрии c. Следующий код дает на 67% меньше времени выполнения. Это то, чего я ожидал!
SimplicialLLT<Sparse_t> solver(V);
assert(solver.info() == Success);
Sparse_t tmp(Q.transpose());
solver.matrixL().solveInPlace(tmp);
result.selfadjointView<Lower>().rankUpdate(tmp.transpose());
Ави, не стесняйся копировать этот фрагмент, и я отмечу твой пост как ответ.