c++ - EIGEN C++ 中的 Matlab Spdiags 等价物

Question

我正在寻找 C++ 中 A=Spdiags(B,d,N,N) 的等价物。此函数提取矩阵 B 的对角线元素，方法是获取 B 的列并将它们放置在 vector d 指定的对角线上。 N N 是输出矩阵A的大小。

在Eigen里找过，好像没有。

有什么想法吗？

Answer 1

据我所知，没有内置方法，但通过索引构建新矩阵并不难。请注意 k第 th 条对角线从索引 ( 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>如果你愿意，和以前一样