Я использую библиотеку Eigen для решения задачи FEM, в которой я использую несколько подобных разреженных матриц для различных видов производных, которые я вычисляю. Чтобы построить разреженную матрицу для решения системы, я хотел бы использовать инициализатор запятой, но он не поддерживается для разреженных матриц (https://eigen.tuxfamily.org/dox/group__TutorialSparse.html и https://stackoverflow.com/a/43532618/15811117). В этом ответе Анри Менке предлагает неподдерживаемый продукт Кронекера Eigen, хотя я не думаю, что это сработает здесь.
С сеткой NxM я вычисляю несколько производных матриц (которые являются разреженными), Dx2, Dxz, и Dz2 размером N*M x N*M. Моя матрица решателя выглядит следующим образом (K - это просто константы, а 0 - это просто заполнитель - нулевая матрица, N * M x N * M)
K3*Dx2 + K1*Dz2 K2*Dxz 0 0 0
K2*Dxz K3*Dz2 + K1*Dx2 0 0 0
0 0 K1*(Dx2+Dz2) 0 0
0 0 0 K3*Dx2 + K1*Dz2 K2*Dxz
0 0 0 K2*Dxz K3*Dz2 + K1*Dx2
На данный момент я создаю его, преобразовывая разреженный в плотный, используя инициализатор запятой, а затем преобразовывая обратно в разреженный с помощью .sparseView(). Единственная проблема заключается в том, что если сетка намного больше 16x16, я получаю ошибку std::bad_alloc (что вполне ожидаемо, поскольку создается огромная матрица: 5*N*M x 5*N*M), и мне нужна сетка не менее 100x100 (и не всегда будет столько 0 матриц, сколько показано сейчас, хотя она останется разреженной матрицей).
Каков наилучший способ построить эту матрицу? Я подумал о том, чтобы попытаться использовать триплеты для этого, подобно этому (Преобразовать собственную матрицу в триплетную форму C++).
Редактировать: я не думаю, что использовать триплеты будет так просто (хотя я сейчас пытаюсь это сделать), потому что я редактирую несколько значений в каждой матрице D_ _ для граничных условий... если нет способ найти триплет в std::vector с заданным номером строки и столбца?
Обновление: оба метода произведения Триплет и Кронекер возможны, но я столкнулся с несколькими проблемами:
С продуктом Кронекера я получаю следующую ошибку (около сотни раз) при компиляции:
.../eigen/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h:225:97: error: no type named ‘ReturnType’ in ‘struct Eigen::ScalarBinaryOpTraits<int, std::complex<double>, Eigen::internal::scalar_product_op<int, std::complex<double> > >’
С помощью метода Triplet я определил несколько операторов и функцию для «сдвига» матрицы в нужное место (в основном то, что делает произведение Кронекера). Мне просто нужно выяснить проблему выхода индекса за пределы.
