亲爱的先生/女士,这是MCVE
#include
#include
#include
using namespace Eigen;
Eigen::SparseMatrix spdiags(const MatrixXd& B,
const Eigen::Matrix& d, size_t m, size_t n)
{
Eigen::SparseMatrix A(m, n);
typedef Eigen::Triplet T;
std::vector triplets;
triplets.reserve(std::min(m, n)*d.size());
for (int k = 0; k < d.size(); k++)
{
int i_min = std::max(0, -d(k));
int i_max = std::min(m - 1, n - d(k) - 1);
int B_idx_start = m >= n ? d(k) : 0;
for (int i = i_min; i <= i_max; i++)
triplets.push_back(T(i, i + k, B(B_idx_start + i, k)));
}
A.setFromTriplets(triplets.begin(), triplets.end());
std::cout << "Row\tCol\tVal" << std::endl;
for (int k = 0; k < A.outerSize(); ++k)
{
for (SparseMatrix::InnerIterator it(A, k); it; ++it)
{
std::cout << it.row() << "\t"; // row index
std::cout << it.col() << "\t";
std::cout << it.value() << std::endl;
}
}
return A;
}
int main()
{
Matrix d1; d1(0) = 0;
MatrixXd d0; d0.resize(1, 5);
// Note that you *have* to use (x,y) indices on a MatrixXd
// Otherwise, you get a different assertion failure
d0(0,0) = 10; d0(0,1) = 20;
d0(0,2) = 30; d0(0,3) = 30;
d0(0,4) = 40;
// d0(0,5) = 50; // OUT OF BOUNDS!!!
Eigen::SparseMatrix Diag_laplacian = spdiags(d0, d1, 5, 5);
}
预期的结果是(如您所述):
Row Col Val
0 0 10
1 1 20
2 2 30
3 3 30
4 4 40
为了重现结果,我可以使用VS(在我的情况下为2013)或g ++(即,不是Linux vs. Mac)。 当您使用g ++时,我也会这样做。
为了重现您在Linux构建中描述的行为,我使用
g++ -O3 -I"C:\usr\include" Source.cpp -o a.exe
运行a.exe给了我(如您所说)
断言失败:行> = 0 &&行 = 0 && col
调试它告诉我它失败了
triplets.push_back(T(i, i + k, B(B_idx_start + i, k)));
当i == 1 。 为什么? 就像@marc和我所说的一样。 B在使用时未定型/定型。 用B(k, B_idx_start + i)更改B(B_idx_start + i, k) B(k, B_idx_start + i)解决了此问题。
现在,为什么它在Mac上可以使用? 答案与错误本身有关。 这是一个断言错误。 定义NDEBUG时不检查断言 。 所以你可能使用类似
g++ -DNDEBUG -O3 -I"C:\usr\include" Source.cpp -o a.exe
在Mac上运行正常,因为断言被忽略了:
#ifdef NDEBUG
#define assert(_Expression) ((void)0)
#else
因此,如果存在断言失败,为什么在定义NDEBUG时它起作用? 答案是数据指针指向五个分配的doubles的第一个。 使用正确的索引,我们应该得到index = k*1 + (B_idx_start + i) ,并且由于在这种情况下k==0和B_idx_start==0 ,我们得到index=i 。 这在范围之内,因此我们不会出现超出范围的异常。 使用不正确的索引,我们得到index = (B_idx_start + i)*1 + k ,这又导致index=i 。 如果矩阵的大小是(例如)2x5,那么我们将获得一个超出范围的异常。