2017-03-05 53 views
1

我目前正在編寫一個小程序,用C++中的有限差分方法求解微分方程。問題是我正在使用lapack來解決我的三角對角矩陣,但是在例程被調用後我得到了一個Segmentation Fault。該例程仍然退出info = 0,所以我不知道問題是什麼。這是該計劃。Lapack例程中的分割錯誤dstev

#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,所以我不知道什麼是錯誤/正確的關於此代碼,但任何幫助表示讚賞。另外,爲了將一個向量傳入例程,爲什麼我們需要傳入begin()迭代器?這是我學習的方式,但我不認爲我完全理解它。

回答

1

我不知道如果這個代碼的唯一問題,但std::vector< std::vector< double> >不分配連續的內存...因此,LAPACK將無法訪問矩陣ev的所有元素。您應該給出一個大小爲n*n的矢量,其矩陣的條目(i,j)i+j*n