Erreur de segmentation dans la routine Lapack dstev

Question

Je suis actuellement en train d'écrire un petit programme pour résoudre une équation différentielle en utilisant une méthode de différences finies en C++. Le problème est que j'utilise lapack pour résoudre ma matrice tridiagonale mais que je reçois une erreur de segmentation après que la routine a été appelée. La routine se termine toujours avec info = 0 donc je ne sais pas quel est le problème. Voici le programme.

#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; 
} 

C'est la première fois que j'utilise le lapack, donc je ne sais pas ce qui est faux/correct à propos de ce code mais toute aide est appréciée. De même, pour passer un vecteur dans la routine, pourquoi devons-nous passer dans l'itérateur begin()? C'est ainsi que je l'ai appris mais je ne pense pas le comprendre complètement.

Answer 1

Je ne sais pas si c'est le seul problème avec ce code, mais std::vector< std::vector< double> > n'alloue pas de mémoire contiguë ... Par conséquent, Lapack ne sera pas en mesure d'accéder à tous les éléments de la matrice ev. Vous devez donner un vecteur de taille n*n pour lequel l'entrée (i,j) d'une matrice est i+j*n.