I am writing a C++ code for some Monte Carlo simulations. In the simulations, I am storing matrices and updating them based on some rules. So to do this, I created two classes : a 5D matrix class, and a 2D matrix class. The 5D matrix is a big 5D table which can store data types of the 2D matrix class, which is just supposed to mimic a 2D array. Both these classes are essentially just hidden vector classes, which I have defined in a way to give me the structure of the matrices required. This was the background.

Now, when I go through my Monte Carlo code updating the 5D matrix with different 2D matrices based on my requirements, for some reason, I always get a std::out_of_range: vector exception in the 5D matrix class. Moreover, this happens when I get to an element which is supposed to updated at index (size/4) of the underlying 5Dmatrix vector. When I print the size of the 5DMatrix, it gives me the same big size required though. So I am not sure where I am making a mistake? Could someone help me? My MonteCarlo Code is attached below, and I have tried to make it readable (the loops are over various dimensions of the 5D matrix)

```
void MonteCarlo(){
for(int i = 0; i<4; i++){
for(int j = 0; j<8; j++){
for(int k = 0; k<8; k++){
for(int l = 0; l<8; l++){
for(int m = 0; m<8; m++){
//Generate a random matrix for the MonetCarlo procedure
complex<double>** trialU = GetRandomMatrix();
TwoDMatrix< complex<double> > trial(3,3);
TwoDMatrix< complex<double> > current = U.get(i,j,k,l,m);
complex<double>** currentU = new complex<double>* (3);//Convert from MyMatrixType to C++ Matrix
for(int a = 0; a<3; a++){
currentU(a) = new complex<double>(3);
for(int b = 0; b<3; b++){
currentU(a)(b) = current.get(a,b);
}
}
//Calculate action (weights)
double act_new = action(currentU, i,j,k,l,m);
double act_old = action(trialU,i,j,k,l,m);
if(act_new < act_old){
U.replace(trial,i,j,k,l,m);
}
else {
double prob = (rand()%100)/100.0;
if(prob < exp(beta*(act_old - act_new))){
U.replace(trial,i,j,k,l,m);
}
}
for(int c1 =0; c1<3; c1++){
delete() trialU(c1);
delete() currentU(c1);
}
delete() trialU;
delete() currentU;
}
}
}
}
}
```

}

In the following, U.replace(…) is the method which gives me troubles and the error. The 5D class is defined below :

```
class FiveDMatrix{
int dim1;
int dim2;
int dim3;
int dim4;
int dim5;
int size;
vector<T> FDmat;
public:
FiveDMatrix(int _dim1, int _dim2, int _dim3, int _dim4,int _dim5){
this->dim1 = _dim1;
this->dim2 = _dim2;
this->dim3 = _dim3;
this->dim4 = _dim4;
this->dim5 = _dim5;
size = 0;
FDmat.reserve(dim1*dim2*dim3*dim4*dim5);
}
void insert(T value){
size++;
FDmat.push_back(value);
}
void replace(T value, int pos1, int pos2, int pos3, int pos4, int pos5){
//Check this weird result : If I include this line, I get an error pretty quickly FDmat(size/4) = value;
int index = pos5 + dim5*pos4 + dim5*dim4*pos3 + dim5*dim4*dim3*pos2 + dim5*dim4*dim3*dim2*pos1;
FDmat(index) = value;
}
T get(int pos1, int pos2, int pos3, int pos4, int pos5){
int index = pos5 + dim5*pos4 + dim5*dim4*pos3 + dim5*dim4*dim3*pos2 + dim5*dim4*dim3*dim2*pos1;
return FDmat.at(index);
}
int num(){
return FDmat.size();
}
```

};

In fact, something which is even more weird is that if I include a line in the replace() method like I noted in the comment, I get the out of bounds error even more quickly. I do not understand why this happens. Is dynamic memory allocation playing a role? Because if I do not run this in a MonteCarlo loop but just separately run the values which seem to give error, I do not get an error. Somehow looping this turns it into an error.

Sorry for the long post, and hopefully my question is understandable. Thanks!