Segmentation Fault in Lapack routine dstev

Question

I am currently just writing a small program to solve a differential equation using a finite difference method in C++. The problem is I am using lapack to solve my tridiagonal matrix but I am getting a Segmentation Fault after the routine has been called. The routine still exits with info=0 so I don't know what the issue is. Here's the program.

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

This is my first time using lapack so I don't know what is wrong/right about this code but any help is appreciated. Also, for passing a vector into the routine, why do we need to pass in the begin() iterator? This is how I learned it but I don't think I fully understand it.

Answer 1

I don't know if its the only issue with this code, but std::vector< std::vector< double> > does not allocate contiguous memory... Therefore, Lapack won't be able to access all element of the matrix ev. You should give a vector of size n*n for which the entry (i,j) of a matrix is i+j*n.