我在C++中搜索A = Spdiags(B,d,N,N)的等效项。此函数通过取B的列并将其沿 vector d指定的对角线放置来提取矩阵B的对角线元素。 N N是输出矩阵A的大小。
我在Eigen中进行了搜索,但似乎不存在。
有任何想法吗?
最佳答案
据我所知,还没有内置方法,但是通过索引构建新矩阵并不难。请注意,对角线k
从索引(max(1, 1-k)
,max(1, 1-k)+k
)到(min(m, n-k)
,min(m, n-k)+k
)
template <typename Scalar>
Eigen::SparseMatrix<Scalar> spdiags(const Matrix<Scalar, -1, -1>& B, const Eigen::Matrix<int, -1, 1>& d, size_t m, size_t n) {
Eigen::SparseMatrix<Scalar> A(m,n);
typedef Eigen::Triplet<Scalar> T;
std::vector<T> 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());
return A;
}
注意,我还没有测试过,但是您知道了。 B的第一个索引有些奇怪,但我认为是正确的。
其他版本,spdiags(A):
Eigen::MatrixXd spdiags(const Eigen::SparseMatrix<double>& A) {
// find nonzero diagonals by iterating over all nonzero elements
// d(i) = 1 if the ith diagonal of A contains a nonzero, 0 else
Eigen::VectorXi d = Eigen::VectorXi::Zero(A.rows() + A.cols() - 1);
for (int k=0; k < A.outerSize(); ++k) {
for (SparseMatrix<double>::InnerIterator it(A,k); it; ++it) {
d(it.col() - it.row() + A.rows() - 1) = 1;
}
}
int num_diags = d.sum();
Eigen::MatrixXd B(std::min(A.cols(), A.rows()), num_diags);
// fill B with diagonals
int B_col_idx = 0;
int B_row_sign = A.rows() >= A.cols() ? 1 : -1;
for (int i = 1 - A.rows(); i <= A.cols() - 1; i++) {
if (d(i + A.rows() - 1)) {
const auto& diag = A.diagonal(i);
int B_row_start = std::max(0, B_row_sign * i);
B.block(B_row_start, B_col_idx, diag.size(), 1) = diag;
B_col_idx++;
}
}
return B;
}
同样的免责声明:未经测试,但应该可以使用。如果需要,请像以前一样用
double
替换template <typename Scalar>
关于c++ - EIGEN C++中的Matlab Spdiags等效项,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/30103824/