2
我最近安装了犰狳并尝试了稀疏矩阵的特征值问题。不幸的是,分解失败是参数'N'(下面的代码)太大e.q.我想知道这里发生了什么。矩阵不是很复杂 - 它具有对角结构。对于太大的矩阵,特征值分解失败,armadillo eigs_sym()
UPDATE
数学也有这个矩阵的问题。它告诉我Arnoldi算法不会收敛。也许我需要在arnoldi arpack例程中手动指定一些参数以确保收敛?
这里是我的代码:
#include <armadillo>
int main()
{
double N = 1000.0;
// create matrix
int kmin = 0;
int kmax = static_cast<int>(std::floor(N/2.0));
int dim = (kmax - kmin) + 1;
// locations and values in sparse matrices
arma::umat hc_locations (2, 3*dim-2);
arma::vec hc_values (3*dim-2);
// diagonal part
for (int k=0; k<dim; k++)
{
hc_locations (0,k) = k;
hc_locations (1,k) = k;
hc_values (k) = 2.0/N*static_cast<double>(kmin + k)*(2.0*(N-2.0*static_cast<double>(k + kmin)) - 1.0);
}
// upper and lower diagonal
for (int k=0; k<dim-1; k++)
{
hc_locations (0,k+dim) = k;
hc_locations (1,k+dim) = k+1;
hc_values (k+dim) = 2.0/N*std::sqrt((static_cast<double>(k+1+kmin)) *
(static_cast<double>(k+1+kmin)) *
(N - static_cast<double>(2*(k+1+kmin)) + 1.0) *
(N - static_cast<double>(2*(k+1+kmin)) + 2.0));
hc_locations (0, k+2*dim-1) = k+1;
hc_locations (1, k+2*dim-1) = k;
hc_values (k+2*dim-1) = 2.0/N*std::sqrt ((static_cast<double>(k+1+kmin)) *
(static_cast<double>(k+1+kmin)) *
(N - static_cast<double>(2*(k+kmin))) *
(N - static_cast<double>(2*(k+kmin)) - 1.0));
}
arma::sp_mat Ham(hc_locations, hc_values);
// eigenvalue problem
arma::vec eigval;
arma::mat eigvec;
arma::eigs_sym(eigval, eigvec, Ham, 3, "sa");
}
前段时间我已经解决了这个问题。问题在于ARPACK库。更具体地说,Armadillo库不允许你指定Krylov基础的大小和最大迭代次数,就像原始的ARPACK库一样。该矩阵的最小特征值之间的距离对于阶数5的基础来说太小(默认Armadillo 2 * nev + 1)。我在Fortran中编写了类似的程序,对于矩阵大小N = 10000,我需要90阶的基本大小以便使用Arnoldi迭代进行收敛。 Mathematica 10.0还允许您指定基向量数和最大迭代次数。有用。 – Bociek