Correct vector to matrix mapping when reading configuration data from HDF5

src/libnptm/Configuration.cpp

@@ -475,9 +475,10 @@ ScattererConfiguration* ScattererConfiguration::from_dedfb(string dedfb_file_nam
 
   complex<double> ***dc0m = new complex<double>**[configurations];
   complex<double> *dc0 = new complex<double>[configurations]();
+  int dim3 = (_idfc == 0) ? nxi : 1;
   for (int di = 0; di < configurations; di++) {
     dc0m[di] = new complex<double>*[nsph];
-    for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[nxi]();
+    for (int dj = 0; dj < nsph; dj++) dc0m[di][dj] = new complex<double>[dim3]();
   } // di loop
   for (int jxi468 = 1; jxi468 <= nxi; jxi468++) {
     if (_idfc != 0 && jxi468 > 1) continue; // jxi468 loop
@@ -574,7 +575,7 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
       str_type = "FLOAT64_(" + to_string(nsh) + ")";
       status = hdf_file->read(str_name, str_type, (rcf_vector[i115 - 1]));
     }
-    xi_vec = new double[nxi];
+    xi_vec = new double[nxi]();
     str_name = "XIVEC";
     str_type = "FLOAT64_(" + to_string(nxi) + ")";
     status = hdf_file->read(str_name, str_type, xi_vec);
@@ -586,13 +587,15 @@ ScattererConfiguration* ScattererConfiguration::from_hdf5(string file_name) {
     str_type = "FLOAT64_(" + to_string(element_size) + ")";
     status = hdf_file->read(str_name, str_type, elements);
     dc0m = new complex<double>**[configuration_count];
+    int dc_index = 0;
     for (int di = 0; di < configuration_count; di++) {
       dc0m[di] = new complex<double>*[nsph];
       for (int dj = 0; dj < nsph; dj++) {
 	dc0m[di][dj] = new complex<double>[dim3]();
 	for (int dk = 0; dk < dim3; dk++) {
-	  double rval = elements[(di * nsph) + (dj * dim3) + 2 * dk];
-	  double ival = elements[(di * nsph) + (dj * dim3) + 2 * dk + 1];
+	  double rval = elements[2 * dc_index];
+	  double ival = elements[2 * dc_index + 1];
+	  dc_index++;
 	  dc0m[di][dj][dk] = complex<double>(rval, ival);
 	}
       } // dj loop
@@ -665,9 +668,10 @@ ScattererConfiguration* ScattererConfiguration::from_legacy(string file_name) {
       input.read(reinterpret_cast<char *>(&(_rcf_vec[i115 - 1][nsi])), sizeof(double));
   }
   _dc0m = new complex<double>**[configurations];
+  int dim3 = (_idfc == 0) ? _nxi : 1;
   for (int di = 0; di < configurations; di++) {
     _dc0m[di] = new complex<double>*[_nsph];
-    for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[_nxi]();
+    for (int dj = 0; dj < _nsph; dj++) _dc0m[di][dj] = new complex<double>[dim3]();
   } // di loop
   for (int jxi468 = 1; jxi468 <= _nxi; jxi468++) {
     if (_idfc != 0 && jxi468 > 1) continue;
@@ -822,7 +826,7 @@ void ScattererConfiguration::write_hdf5(string file_name) {
 
   int dim3 = (idfc == 0) ? number_of_scales : 1;
   int dc0m_size = 2 * dim3 * number_of_spheres * configurations;
-  double *dc0m = new double[dc0m_size];
+  double *dc0m = new double[dc0m_size]();
   int dc0_index = 0;
   for (int jxi468 = 1; jxi468 <= number_of_scales; jxi468++) {
     if (idfc != 0 && jxi468 > 1) continue;
@@ -835,8 +839,9 @@ void ScattererConfiguration::write_hdf5(string file_name) {
 	double dc0_real, dc0_imag;
 	dc0_real = dc0_matrix[i157][i162 - 1][jxi468 - 1].real();
 	dc0_imag = dc0_matrix[i157][i162 - 1][jxi468 - 1].imag();
-	dc0m[dc0_index++] = dc0_real;
-	dc0m[dc0_index++] = dc0_imag;
+	dc0m[2 * dc0_index] = dc0_real;
+	dc0m[2 * dc0_index + 1] = dc0_imag;
+	dc0_index++;
       }
     }
   }
@@ -1109,6 +1114,7 @@ bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
       return false;
     }
   int dj_index = 0;
+  int dim3 = (idfc == 0) ? number_of_scales : 1;
   for (int dj = 0; dj < number_of_spheres; dj++) {
     bool check_matrixes = false;
     if (iog_vec[dj] >= dj + 1) {
@@ -1130,7 +1136,7 @@ bool ScattererConfiguration::operator ==(ScattererConfiguration &other) {
 	if (rcf[dj_index][di] != other.rcf[dj_index][di]) {
 	  return false;
 	}
-	for (int dk = 0; dk < number_of_scales; dk++) {
+	for (int dk = 0; dk < dim3; dk++) {
 	  if (dc0_matrix[di][dj_index][dk] != other.dc0_matrix[di][dj_index][dk]) {
 	    return false;
 	  }