WSL/SLF GitLab Repository

Commit 44388749 authored by Mathias Bavay's avatar Mathias Bavay
Browse files

The array memory alignement has been fixed for Array2D (this was the only that...

The array memory alignement has been fixed for Array2D (this was the only that was wrongly aligned). This is a fix to issue 43
parent 11b3fa7d
......@@ -166,10 +166,8 @@ template<class T> T& Array2D<T>::operator()(const unsigned int& x, const unsigne
throw IndexOutOfBoundsException("", AT);
}
#endif
//the 2D array is stored by columns in a 1D vector. Each column follows the previous one.
//This matches the data usage in Alpine3D, data access being done by columns.
return vecData[x*ny + y];
//ROW-MAJOR alignment of the vector: fully C-compatible memory layout
return vecData[x + y*nx];
}
template<class T> const T Array2D<T>::operator()(const unsigned int& x, const unsigned int& y) const {
......@@ -178,7 +176,7 @@ template<class T> const T Array2D<T>::operator()(const unsigned int& x, const un
throw IndexOutOfBoundsException("", AT);
}
#endif
return vecData[x*ny + y];
return vecData[x + y*nx];
}
template<class T> Array2DProxy<T> Array2D<T>::operator[](const unsigned int& i) {
......@@ -207,10 +205,9 @@ template<class T> void Array2D<T>::subset(const Array2D<T>& _array2D, const unsi
resize(_ncols, _nrows); //create new Array2D object
//Copy by value subspace
for (unsigned int ii=0; ii<ny; ii++) {
for (unsigned int jj=0; jj<nx; jj++) {
//Running through the vector in order of memory alignment HACK
operator()(jj,ii) = _array2D(_nx+jj, _ny+ii);
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(jj,ii) = _array2D(_nx+ii, _ny+jj);
}
}
}
......@@ -240,8 +237,8 @@ template<class T> void Array2D<T>::resize(const unsigned int& anx, const unsigne
template<class T> void Array2D<T>::resize(const unsigned int& anx, const unsigned int& any, const T& init) {
resize(anx, any);
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) = init;
}
}
......@@ -259,8 +256,8 @@ template<class T> void Array2D<T>::clear() {
template<class T> std::ostream& operator<<(std::ostream& os, const Array2D<T>& array) {
os << "<array2d>\n";
for(unsigned int ii=0; ii<array.nx; ii++) {
for (unsigned int jj=0; jj<array.ny; jj++) {
for(unsigned int jj=0; jj<array.ny; jj++) {
for (unsigned int ii=0; ii<array.nx; ii++) {
os << array(ii,jj) << " ";
}
os << "\n";
......@@ -274,16 +271,16 @@ template<class T> T Array2D<T>::getMin(const IOUtils::nodata_handling flag_nodat
T min = std::numeric_limits<T>::max();
if(flag_nodata==IOUtils::RAW_NODATA) {
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
if(val<min) min=val;
}
}
return min;
} else if(flag_nodata==IOUtils::PARSE_NODATA) {
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
if(val!=IOUtils::nodata && val<min) min=val;
}
......@@ -300,16 +297,16 @@ template<class T> T Array2D<T>::getMax(const IOUtils::nodata_handling flag_nodat
T max = -std::numeric_limits<T>::max();
if(flag_nodata==IOUtils::RAW_NODATA) {
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
if(val>max) max=val;
}
}
return max;
} else if(flag_nodata==IOUtils::PARSE_NODATA) {
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
if(val!=IOUtils::nodata && val>max) max=val;
}
......@@ -326,8 +323,8 @@ template<class T> T Array2D<T>::getMean(const IOUtils::nodata_handling flag_noda
T mean = 0;
if(flag_nodata==IOUtils::RAW_NODATA) {
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
mean += val;
}
......@@ -337,8 +334,8 @@ template<class T> T Array2D<T>::getMean(const IOUtils::nodata_handling flag_noda
else return (T)0;
} else if(flag_nodata==IOUtils::PARSE_NODATA) {
unsigned int count = 0;
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
const T val = operator()(ii,jj);
if(val!=IOUtils::nodata) {
mean += val;
......@@ -370,8 +367,8 @@ template<class T> Array2D<T>& Array2D<T>::operator+=(const Array2D<T>& rhs)
throw IOException("Trying to add two Array2D objects with different dimensions", AT);
//Add to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) += rhs(ii,jj);
}
}
......@@ -390,8 +387,8 @@ template<class T> const Array2D<T> Array2D<T>::operator+(const Array2D<T>& rhs)
template<class T> Array2D<T>& Array2D<T>::operator+=(const T& rhs)
{
//Add to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) += rhs;
}
}
......@@ -414,8 +411,8 @@ template<class T> Array2D<T>& Array2D<T>::operator-=(const Array2D<T>& rhs)
throw IOException("Trying to substract two Array2D objects with different dimensions", AT);
//Substract to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) -= rhs(ii,jj);
}
}
......@@ -434,8 +431,8 @@ template<class T> const Array2D<T> Array2D<T>::operator-(const Array2D<T>& rhs)
template<class T> Array2D<T>& Array2D<T>::operator-=(const T& rhs)
{
//Substract to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) -= rhs;
}
}
......@@ -458,8 +455,8 @@ template<class T> Array2D<T>& Array2D<T>::operator*=(const Array2D<T>& rhs)
throw IOException("Trying to multiply two Array2D objects with different dimensions", AT);
//Add to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) *= rhs(ii,jj);
}
}
......@@ -478,8 +475,8 @@ template<class T> const Array2D<T> Array2D<T>::operator*(const Array2D<T>& rhs)
template<class T> Array2D<T>& Array2D<T>::operator*=(const T& rhs)
{
//Add to every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) *= rhs;
}
}
......@@ -502,8 +499,8 @@ template<class T> Array2D<T>& Array2D<T>::operator/=(const Array2D<T>& rhs)
throw IOException("Trying to divide two Array2D objects with different dimensions", AT);
//Divide every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) /= rhs(ii,jj);
}
}
......@@ -522,8 +519,8 @@ template<class T> const Array2D<T> Array2D<T>::operator/(const Array2D<T>& rhs)
template<class T> Array2D<T>& Array2D<T>::operator/=(const T& rhs)
{
//Divide every single member of the Array2D<T>
for (unsigned int ii=0; ii<nx; ii++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int jj=0; jj<ny; jj++) {
for (unsigned int ii=0; ii<nx; ii++) {
operator()(ii,jj) /= rhs;
}
}
......
......@@ -295,8 +295,8 @@ void DEMObject::printFailures() {
bool header=true;
if(update_flag&SLOPE) {
for ( unsigned int i = 0; i < ncols; i++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int i = 0; i < ncols; i++ ) {
if((slope(i,j)==IOUtils::nodata) && (grid2D(i,j)!=IOUtils::nodata)) {
if(header==true) {
std::cout << "[i] DEM slope could not be computed at the following points \n";
......@@ -310,8 +310,8 @@ void DEMObject::printFailures() {
}
if(update_flag&CURVATURE) {
for ( unsigned int i = 0; i < ncols; i++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int i = 0; i < ncols; i++ ) {
if((curvature(i,j)==IOUtils::nodata) && (grid2D(i,j)!=IOUtils::nodata)) {
if(header==true) {
std::cout << "[i] DEM curvature could not be computed at the following points \n";
......@@ -339,8 +339,8 @@ void DEMObject::printFailures() {
*/
void DEMObject::sanitize() {
if(slope_failures>0 || curvature_failures>0) {
for ( unsigned int i = 0; i < ncols; i++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int i = 0; i < ncols; i++ ) {
if(update_flag&SLOPE) {
if((slope(i,j)==IOUtils::nodata) && (grid2D(i,j)!=IOUtils::nodata)) {
grid2D(i,j) = IOUtils::nodata;
......@@ -523,8 +523,8 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
}
//Now, calculate the parameters using the previously defined function pointer
for ( unsigned int i = 0; i < ncols; i++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int i = 0; i < ncols; i++ ) {
if( grid2D(i,j) == IOUtils::nodata ) {
if(update_flag&SLOPE) {
_slope = _azi = IOUtils::nodata;
......@@ -557,8 +557,8 @@ void DEMObject::CalculateAziSlopeCurve(slope_type algorithm) {
}
if((update_flag&SLOPE) && (algorithm==D8)) { //extra processing required: discretization
for ( unsigned int i = 0; i < ncols; i++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int j = 0; j < nrows; j++ ) {
for ( unsigned int i = 0; i < ncols; i++ ) {
//TODO: process flats by an extra algorithm
if(azi(i,j)!=IOUtils::nodata)
azi(i,j) = fmod(floor( (azi(i,j)+22.5)/45. )*45., 360.);
......
......@@ -371,8 +371,8 @@ void RHAlgorithm::calculate(const MeteoData::Parameters& param, Grid2DObject& gr
Interpol2D::LapseIDW(vecTd, vecMeta, dem, vecCoefficients, &Interpol2D::LinProject, grid);
//Recompute Rh from the interpolated td
for (unsigned int ii=0; ii<grid.ncols; ii++) {
for (unsigned int jj=0; jj<grid.nrows; jj++) {
for (unsigned int jj=0; jj<grid.nrows; jj++) {
for (unsigned int ii=0; ii<grid.ncols; ii++) {
grid.grid2D(ii,jj) = Interpol2D::DewPointtoRh(grid.grid2D(ii,jj), ta.grid2D(ii,jj), 1);
}
}
......
......@@ -228,8 +228,8 @@ void Interpol2D::stdPressureGrid2DFill(const DEMObject& dem, Grid2DObject& grid)
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
//provide each point with an altitude dependant pressure... it is worth what it is...
for (unsigned int i=0; i<grid.ncols; i++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
if (dem.grid2D(i,j)!=IOUtils::nodata) {
grid.grid2D(i,j) = lw_AirPressure(dem.grid2D(i,j));
} else {
......@@ -251,8 +251,8 @@ void Interpol2D::constantGrid2DFill(const double& value, const DEMObject& dem, G
grid.set(dem.ncols, dem.nrows, dem.cellsize, dem.llcorner);
//fills a data table with constant values
for (unsigned int i=0; i<grid.ncols; i++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
if (dem.grid2D(i,j)!=IOUtils::nodata) {
grid.grid2D(i,j) = value;
} else {
......@@ -281,8 +281,8 @@ void Interpol2D::constantLapseGrid2DFill(const double& value, const double& alti
//fills a data table with constant values and then reprojects it to the DEM's elevation from a given altitude
//the laspe rate parameters must have been set before
for (unsigned int i=0; i<grid.ncols; i++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
if (dem.grid2D(i,j)!=IOUtils::nodata) {
grid.grid2D(i,j) = funcptr(value, altitude, dem.grid2D(i,j), vecCoefficients);
} else {
......@@ -338,8 +338,8 @@ void Interpol2D::LapseIDW(const std::vector<double>& vecData_in, const std::vect
const double xllcorner = dem.llcorner.getEasting();
const double yllcorner = dem.llcorner.getNorthing();
const double cellsize = dem.cellsize;
for (unsigned int i=0; i<grid.ncols; i++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
if (dem.grid2D(i,j)!=IOUtils::nodata) {
grid.grid2D(i,j) = IDWCore((xllcorner+i*cellsize),
(yllcorner+j*cellsize), vecTref, vecStations_in);
......@@ -371,8 +371,8 @@ void Interpol2D::IDW(const std::vector<double>& vecData_in, const std::vector<St
const double xllcorner = dem.llcorner.getEasting();
const double yllcorner = dem.llcorner.getNorthing();
const double cellsize = dem.cellsize;
for (unsigned int i=0; i<grid.ncols; i++) {
for(unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int j=0; j<grid.nrows; j++) {
for (unsigned int i=0; i<grid.ncols; i++) {
if (dem.grid2D(i,j)!=IOUtils::nodata) {
grid.grid2D(i,j) = IDWCore((xllcorner+i*cellsize), (yllcorner+j*cellsize),
vecData_in, vecStations_in);
......@@ -409,8 +409,8 @@ void Interpol2D::SimpleDEMWindInterpolate(const DEMObject& dem, Grid2DObject& VW
double Ww; // Wind weighting
double Od; // Diverting factor
for (unsigned int i=0;i<VW.ncols-1;i++) {
for (unsigned int j=0;j<VW.nrows-1;j++){
for (unsigned int j=0;j<VW.nrows-1;j++) {
for (unsigned int i=0;i<VW.ncols-1;i++){
// Get input data
speed = VW.grid2D(i,j);
dir = DW.grid2D(i,j);
......
......@@ -46,6 +46,8 @@ namespace mio {
* @endcode
*/
//TODO: keep a pointer to last read position, so we can fseek for the next read?
const std::string SMETIO::smet_version = "0.99";
const unsigned int SMETIO::buffer_reserve = 23*24*2; //kind of average size of a buffer for optimizing vectors
map<string, MeteoData::Parameters> SMETIO::mapParameterByName;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment