Lapack 例程 dstev 中的分段错误
Segmentation Fault in Lapack routine dstev
我目前正在编写一个小程序,使用 C++ 中的有限差分法求解微分方程。问题是我正在使用 lapack 来求解我的三对角矩阵,但是在调用例程后我遇到了分段错误。该例程仍然以 info=0 退出,所以我不知道问题出在哪里。这是程序。
#include <iostream>
#include <cmath>
#include <vector>
#define pi 3.14159265358979323846
using namespace std;
extern "C" void dstev_(char* job, int* N, double* D, double* OFFD, double* EV, int* VDIM, double* WORK, int* INFO);
int main(){
int i = 0;
// systems parameters
double a = 2e10-6;
double k = pi/1.55e-6;
double n1_sq = pow(3.0, 2), n2_sq = pow(3.1, 2);
double step = 1e-6;
double factor = 1/(k*k*step*step);
// parameters for lapack
char job = 'V';
int N = 80, vdim = 100, info;
// create and initialize vector arrays for lapack routine
vector<double> d, offd, work;
vector< vector<double> > ev(vdim, vector<double>(N));
// set up diagonal elements
d.push_back(n1_sq - factor);
for(i=1; i<N-1; i++){
if(i<N/4 || i>=3*N/4){
d.push_back(n1_sq - 2*factor);
}
else{
d.push_back(n2_sq - 2*factor);
}
}
d.push_back(n1_sq - factor);
// set up off diagonal elements
for(i=0; i<N-1; i++){
offd.push_back(factor);
}
// initialize other arrays
for(i=0; i<vdim; i++){
work.push_back(0.0);
}
cout << "Before routine" << endl;
dstev_(&job, &N, &*d.begin(), &*offd.begin(), &*(ev.begin()->begin()), &vdim, &*work.begin(), &info);
cout << "After routine"<< endl;
return 0;
}
这是我第一次使用 lapack,所以我不知道这段代码 wrong/right 有什么用,但如有任何帮助,我们将不胜感激。另外,为了将向量传递到例程中,为什么我们需要传递 begin() 迭代器?我是这样学的,但我不认为我完全理解它。
我不知道这是不是这段代码的唯一问题,但是 std::vector< std::vector< double> >
没有分配连续的内存...因此,Lapack 将无法访问矩阵的所有元素 ev
。您应该给出一个大小为 n*n
的向量,其矩阵的条目 (i,j)
为 i+j*n
.
我目前正在编写一个小程序,使用 C++ 中的有限差分法求解微分方程。问题是我正在使用 lapack 来求解我的三对角矩阵,但是在调用例程后我遇到了分段错误。该例程仍然以 info=0 退出,所以我不知道问题出在哪里。这是程序。
#include <iostream>
#include <cmath>
#include <vector>
#define pi 3.14159265358979323846
using namespace std;
extern "C" void dstev_(char* job, int* N, double* D, double* OFFD, double* EV, int* VDIM, double* WORK, int* INFO);
int main(){
int i = 0;
// systems parameters
double a = 2e10-6;
double k = pi/1.55e-6;
double n1_sq = pow(3.0, 2), n2_sq = pow(3.1, 2);
double step = 1e-6;
double factor = 1/(k*k*step*step);
// parameters for lapack
char job = 'V';
int N = 80, vdim = 100, info;
// create and initialize vector arrays for lapack routine
vector<double> d, offd, work;
vector< vector<double> > ev(vdim, vector<double>(N));
// set up diagonal elements
d.push_back(n1_sq - factor);
for(i=1; i<N-1; i++){
if(i<N/4 || i>=3*N/4){
d.push_back(n1_sq - 2*factor);
}
else{
d.push_back(n2_sq - 2*factor);
}
}
d.push_back(n1_sq - factor);
// set up off diagonal elements
for(i=0; i<N-1; i++){
offd.push_back(factor);
}
// initialize other arrays
for(i=0; i<vdim; i++){
work.push_back(0.0);
}
cout << "Before routine" << endl;
dstev_(&job, &N, &*d.begin(), &*offd.begin(), &*(ev.begin()->begin()), &vdim, &*work.begin(), &info);
cout << "After routine"<< endl;
return 0;
}
这是我第一次使用 lapack,所以我不知道这段代码 wrong/right 有什么用,但如有任何帮助,我们将不胜感激。另外,为了将向量传递到例程中,为什么我们需要传递 begin() 迭代器?我是这样学的,但我不认为我完全理解它。
我不知道这是不是这段代码的唯一问题,但是 std::vector< std::vector< double> >
没有分配连续的内存...因此,Lapack 将无法访问矩阵的所有元素 ev
。您应该给出一个大小为 n*n
的向量,其矩阵的条目 (i,j)
为 i+j*n
.