Memory Access Error - What to do?

[redpottpotter]

Hey there :)
I am trying to write a programm finding vectorial (3 component) zeros of the Pfaffian determinant of a given matrix. I analytically found an expression and now only need the values of kx, ky, kz where Pfaff is zero. My code gives me a memory access error and I really have no clue what the problem is. I would really appreciate help!

Here my full code:

#include <iostream>
#include <math.h>
#include <numeric>
#include <fstream>
#include <omp.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>

using namespace std;
using namespace boost::numeric::odeint;

ofstream data("data_BFS_quintet.txt");
const double hbar = 6.582119e-01; // Planck constant (meV*ps)
const double Ecut = 250000.0; //cut-off energy (meV)
const double m = 511e3; //electron mass (meV)
const double beta = 2.1e-3; //spin-orbit coupling (meV)
const double alpha = pow(hbar,2)/(2*m)-beta/2;
const double kB = 0.0861733; // Boltzmann constant (meV/K)
const double T = 0.3; // Temperature (K)
const double e = 1.6021766e-07; // elementary charge (A ps)
const double mu = 50.0; //chemical potential (meV)
const double Delta_0 = 2.0; //post-quench gap parameter
const double a = sqrt(2*m*mu/(hbar*hbar))*2;

int counter_e = 377696;

const int Nx = 200;
const int Ny = 200;
const int Nz = 200;
double kx_e[Nx];
double ky_e[Ny];
double kz_e[Nz];
int **sigma_e = new int *[3];


typedef boost::numeric::ublas::matrix< double > state_type;


inline double E_p(double alpha, double beta, double mu, double kx, double ky, double kz){
return (alpha + 9/4*beta)*(kx*kx+ky*ky+kz*kz)-mu;
}

double E_m(double alpha, double beta, double mu, double kx, double ky, double kz){
return (alpha + 1/4*beta)*(kx*kx+ky*ky+kz*kz)-mu;
}

double c_xz(double beta, double kx, double kz){
return sqrt(3)*beta*kx*kz;
}

double c_yz(double beta, double ky, double kz){
return sqrt(3)*beta*ky*kz;
}

double Pfaff(double kx[], double ky[], double kz[], int sigma_x[], int sigma_y[], int sigma_z[], double Delta_0){
double result;
int i;
result += E_p(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)*E_p(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)*E_m(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)*E_m(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)+4*pow(Delta_0,2)*(E_p(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)*E_m(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]], alpha, beta, mu)+c_xz(kx[sigma_x[i]], kz[sigma_z[i]], beta)*c_xz(kx[sigma_x[i]], kz[sigma_z[i]], beta)+c_yz(ky[sigma_y[i]], kz[sigma_z[i]], beta)*c_yz(ky[sigma_y[i]], kz[sigma_z[i]], beta));
return result;
}

void write_data(double x, double y, double z)
{
data << x << '\t' << y << '\t' << z << endl;
cout << x << '\t' << y << '\t' << z << endl;
}


int main(){
double ke_xlimits[2] = {-a, a};
double ke_ylimits[2] = {-a, a};
double ke_zlimits[2] = {-a, a};

kx_e[0] = ke_xlimits[0];
ky_e[0] = ke_ylimits[0];
kz_e[0] = ke_zlimits[0];


for(int i = 1; i < Nx; i++){
kx_e[i] = kx_e[i-1] + (ke_xlimits[1]-ke_xlimits[0])/(Nx-1);
}
for(int i = 1; i < Ny; i++){
ky_e[i] = ky_e[i-1] + (ke_ylimits[1]-ke_ylimits[0])/(Ny-1);
}
for(int i = 1; i < Nz; i++){
kz_e[i] = kz_e[i-1] + (ke_zlimits[1]-ke_zlimits[0])/(Nz-1);
}

int counter_e = 0;

for(int i = 0; i < Nx; i++){
for(int j = 0; j < Ny; j++){
for(int l = 0; l < Nz; l++){
if(abs(E_p(alpha, beta, mu, kx_e[i], ky_e[i], kz_e[i])) < Ecut){
counter_e++;
}
}
}
}
cout << counter_e << endl;

sigma_e[0] = new int[counter_e];
sigma_e[1] = new int[counter_e];
sigma_e[2] = new int[counter_e];


counter_e = 0;

for(int i = 0; i < Nx; i++){
for(int j = 0; j < Ny; j++){
for(int l = 0; l < Nz; l++){
if(abs(E_p(alpha, beta, mu, kx_e[i], ky_e[i], kz_e[i])) < Ecut){
sigma_e[0][counter_e] = i;
sigma_e[1][counter_e] = j;
sigma_e[2][counter_e] = l;

counter_e++;
}
}
}
}
for(int i = 0; i < counter_e; i++){

if (Pfaff(kx_e, ky_e, kz_e, sigma_e[0], sigma_e[1], sigma_e[2], Delta_0) == 0) {
double x = kx_e[sigma_e[0][i]];
double y = ky_e[sigma_e[1][i]];
double z = kz_e[sigma_e[2][i]];
write_data(x, y, z);

}
}
return 0;
}

[lastchance]

PLEASE USE CODE TAGS

double result;
int i;
result += E_p(kx[sigma_x[i]],ky[sigma_y[i]], kz[sigma_z[i]]

What are the values of result and i when that third line is called?

Your code didn't actually fall over in cpp.sh (surprisingly).

[jonnin]

just about every line of that code is hitting an array with an index.
at least one of those, the index has gone wrong and is out of bounds.
you need to find it... either use a debugger or find the line where it crashes (by commenting out code until it runs) and then print the indices in the calls near where it dies to see what went bad.

the c++ header is cmath by the way.