从几个较小的稀疏矩阵构建大型稀疏矩阵



我正在使用特征库来解决一个FEM问题,在该问题中,我使用几个类似的稀疏矩阵来计算不同类型的导数。为了构建稀疏矩阵来解决系统,我想使用逗号初始化器,但稀疏矩阵不支持这种方式(https://eigen.tuxfamily.org/dox/group__TutorialSparse.html和https://stackoverflow.com/a/43532618/15811117)。在这个答案中,Henri Menke确实提出了Eigen不支持的Kronecker产品,尽管我认为这在这里不起作用。

使用NxM网格,我计算了几个导数矩阵(它们是稀疏的(,尺寸为N*MxN*M的Dx2, Dxz,Dz2。我的求解器矩阵如下,(K只是常数,0只是填充物——一个零矩阵,N*MxN*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*Mx5*N*M(,并且我需要一个至少100x100的网格(并且不会总是有如当前所示的那么多0矩阵,尽管它将保持稀疏矩阵(。

构建此矩阵的最佳方法是什么我曾想过尝试使用Triplet来实现这一点,类似于此(将特征矩阵转换为Triplet形式C++(。

编辑:我认为使用Triplets不会那么容易(尽管我目前正在尝试(,因为我为边界条件编辑了每个D_矩阵中的几个值。。。除非有一种方法可以在给定行和列编号的std::向量中找到Triplet?

更新:Triplet和Kronecker乘积方法都是可行的,但我遇到了几个问题:

使用Kronecker产品,我在编译时会出现以下错误(大约一百次(:

.../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方法,我定义了一些运算符和一个函数来将矩阵"移动"到正确的位置(基本上就是Kronecker乘积所做的(。我只需要找出一个索引越界的问题。

问题来自于我在Kronecker乘积中使用的矩阵类型为SparseMatrix<int>SparseMatrix<complex<double>>intcomplex<double>没有运算符重载,因此通过使这两个矩阵的类型兼容,代码可以按预期编译和运行。

最新更新