Skip to content

GitLab

Commit 10eef8a7 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Diff: Added memory support for 32-bit float data.

parent e30d5b0a
Hide whitespace changes
Inline Side-by-side
ChangeLog
View file @ 10eef8a7
2020-03-21 Uwe Schulzweida
* Diff: Added memory support for 32-bit float data.
2020-03-17 Uwe Schulzweida 2020-03-17 Uwe Schulzweida
* Vertintap: Added memory support for 32-bit float data. * Vertintap: Added memory support for 32-bit float data.
......
src/Diff.cc
View file @ 10eef8a7
...@@ -36,6 +36,79 @@ ...@@ -36,6 +36,79 @@
#include "pmlist.h" #include "pmlist.h"
#include "cdo_zaxis.h" #include "cdo_zaxis.h"
static void
varListSynchronizeMemtype(VarList &varList1, VarList &varList2, const std::map<int, int> &mapOfVarIDs)
{
int nvars = varList1.size();
for (int varID1 = 0; varID1 < nvars; varID1++)
{
auto it = mapOfVarIDs.find(varID1);
if (it == mapOfVarIDs.end()) continue;
auto varID2 = it->second;
if (varList1[varID1].memType == MemType::Float && varList2[varID2].memType == MemType::Double)
varList1[varID1].memType = MemType::Double;
else if (varList1[varID1].memType == MemType::Double && varList2[varID2].memType == MemType::Float)
varList2[varID2].memType = MemType::Double;
}
}
template <typename T>
static void
diff_kernel(bool lmiss, const T v1, const T v2, const T missval1, const T missval2,
size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
{
const auto v1isnan = std::isnan(v1);
const auto v2isnan = std::isnan(v2);
const auto v1ismissval = lmiss ? DBL_IS_EQUAL(v1, missval1) : false;
const auto v2ismissval = lmiss ? DBL_IS_EQUAL(v2, missval2) : false;
if ((v1isnan && !v2isnan) || (!v1isnan && v2isnan))
{
ndiff++;
relm = 1.0;
}
else if (lmiss == false || (!v1ismissval && !v2ismissval))
{
const double absdiff = std::fabs(v1 - v2);
if (absdiff > 0.) ndiff++;
absm = std::max(absm, absdiff);
const auto vv = v1 * v2;
if (vv < 0.)
dsgn = true;
else if (IS_EQUAL(vv, 0.))
zero = true;
else
relm = std::max(relm, absdiff / std::max(std::fabs(v1), std::fabs(v2)));
}
else if ((v1ismissval && !v2ismissval) || (!v1ismissval && v2ismissval))
{
ndiff++;
relm = 1.0;
}
}
static void
diff_kernel(size_t i, bool lmiss, const Field &field1, const Field &field2, size_t &ndiff, bool &dsgn, bool &zero, double &absm, double &relm)
{
if (field1.memType == MemType::Float)
diff_kernel(lmiss, field1.vec_f[i], field2.vec_f[i], (float)field1.missval, (float)field2.missval, ndiff, dsgn, zero, absm, relm);
else
diff_kernel(lmiss, field1.vec_d[i], field2.vec_d[i], field1.missval, field2.missval, ndiff, dsgn, zero, absm, relm);
}
static void
use_real_part(size_t gridsize, Field &field)
{
if (field.memType == MemType::Float)
for (size_t i = 0; i < gridsize; ++i) field.vec_f[i] = field.vec_f[i * 2];
else
for (size_t i = 0; i < gridsize; ++i) field.vec_d[i] = field.vec_d[i * 2];
}
static void static void
diffGetParameter(double &abslim, double &abslim2, double &rellim, int &mapflag, int &maxcount) diffGetParameter(double &abslim, double &abslim2, double &rellim, int &mapflag, int &maxcount)
{ {
...@@ -81,9 +154,7 @@ Diff(void *process) ...@@ -81,9 +154,7 @@ Diff(void *process)
int nrecs, nrecs2; int nrecs, nrecs2;
int varID1, varID2 = -1; int varID1, varID2 = -1;
int levelID; int levelID;
size_t nmiss1, nmiss2;
int ndrec = 0, nd2rec = 0, ngrec = 0; int ndrec = 0, nd2rec = 0, ngrec = 0;
char varname[CDI_MAX_NAME];
char paramstr[32]; char paramstr[32];
cdoInitialize(process); cdoInitialize(process);
...@@ -124,10 +195,12 @@ Diff(void *process) ...@@ -124,10 +195,12 @@ Diff(void *process)
vlistMap(vlistID1, vlistID2, CMP_ALL, mapflag, mapOfVarIDs); vlistMap(vlistID1, vlistID2, CMP_ALL, mapflag, mapOfVarIDs);
} }
const auto gridsizemax = vlistGridsizeMax(vlistID1); VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
varListSynchronizeMemtype(varList1, varList2, mapOfVarIDs);
Varray<double> array1(gridsizemax), array2(gridsizemax); Field field1, field2;
Varray<double> work(vlistNumber(vlistID1) != CDI_REAL ? 2 * gridsizemax : 0);
const auto taxisID = vlistInqTaxis(vlistID1); const auto taxisID = vlistInqTaxis(vlistID1);
...@@ -138,10 +211,8 @@ Diff(void *process) ...@@ -138,10 +211,8 @@ Diff(void *process)
bool lstop = false; bool lstop = false;
nrecs = cdoStreamInqTimestep(streamID1, tsID); nrecs = cdoStreamInqTimestep(streamID1, tsID);
const auto vdate = taxisInqVdate(taxisID); const auto vdateString = dateToString(taxisInqVdate(taxisID));
const auto vtime = taxisInqVtime(taxisID); const auto vtimeString = timeToString(taxisInqVtime(taxisID));
const auto vdateString = dateToString(vdate);
const auto vtimeString = timeToString(vtime);
nrecs2 = cdoStreamInqTimestep(streamID2, tsID); nrecs2 = cdoStreamInqTimestep(streamID2, tsID);
...@@ -175,67 +246,29 @@ Diff(void *process) ...@@ -175,67 +246,29 @@ Diff(void *process)
indg += 1; indg += 1;
const auto param = vlistInqVarParam(vlistID1, varID1); const auto gridsize = varList1[varID1].gridsize;
const auto code = vlistInqVarCode(vlistID1, varID1);
const auto gridID = vlistInqVarGrid(vlistID1, varID1);
const auto zaxisID = vlistInqVarZaxis(vlistID1, varID1);
const auto gridsize = gridInqSize(gridID);
const auto missval1 = vlistInqVarMissval(vlistID1, varID1);
const auto missval2 = vlistInqVarMissval(vlistID2, varID2);
// checkrel = gridInqType(gridID) != GRID_SPECTRAL; // checkrel = gridInqType(gridID) != GRID_SPECTRAL;
const bool checkrel = true; const bool checkrel = true;
cdiParamToString(param, paramstr, sizeof(paramstr)); cdiParamToString(varList1[varID1].param, paramstr, sizeof(paramstr));
if (vlistInqVarNumber(vlistID1, varID1) == CDI_COMP) field1.init(varList1[varID1]);
{ cdoReadRecord(streamID1, field1);
cdoReadRecord(streamID1, work.data(), &nmiss1); if (varList1[varID1].nwpv == CDI_COMP) use_real_part(gridsize, field1);
for (size_t i = 0; i < gridsize; ++i) array1[i] = work[i * 2];
}
else
cdoReadRecord(streamID1, array1.data(), &nmiss1);
if (vlistInqVarNumber(vlistID2, varID2) == CDI_COMP) field2.init(varList2[varID2]);
{ cdoReadRecord(streamID2, field2);
cdoReadRecord(streamID2, work.data(), &nmiss2); if (varList2[varID2].nwpv == CDI_COMP) use_real_part(gridsize, field2);
for (size_t i = 0; i < gridsize; ++i) array2[i] = work[i * 2];
}
else
cdoReadRecord(streamID2, array2.data(), &nmiss2);
int ndiff = 0; bool lmiss = field1.nmiss || field2.nmiss;
size_t ndiff = 0;
bool dsgn = false, zero = false; bool dsgn = false, zero = false;
double absm = 0.0, relm = 0.0; double absm = 0.0, relm = 0.0;
double absdiff;
for (size_t i = 0; i < gridsize; i++) for (size_t i = 0; i < gridsize; i++)
{ {
if ((std::isnan(array1[i]) && !std::isnan(array2[i])) || (!std::isnan(array1[i]) && std::isnan(array2[i]))) diff_kernel(i, lmiss, field1, field2, ndiff, dsgn, zero, absm, relm);
{
ndiff++;
relm = 1.0;
}
else if (!DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2))
{
absdiff = std::fabs(array1[i] - array2[i]);
if (absdiff > 0.) ndiff++;
absm = std::max(absm, absdiff);
if (array1[i] * array2[i] < 0.)
dsgn = true;
else if (IS_EQUAL(array1[i] * array2[i], 0.))
zero = true;
else
relm = std::max(relm, absdiff / std::max(std::fabs(array1[i]), std::fabs(array2[i])));
}
else if ((DBL_IS_EQUAL(array1[i], missval1) && !DBL_IS_EQUAL(array2[i], missval2))
|| (!DBL_IS_EQUAL(array1[i], missval1) && DBL_IS_EQUAL(array2[i], missval2)))
{
ndiff++;
relm = 1.0;
}
} }
if (!Options::silentMode || Options::cdoVerbose) if (!Options::silentMode || Options::cdoVerbose)
...@@ -260,8 +293,6 @@ Diff(void *process) ...@@ -260,8 +293,6 @@ Diff(void *process)
fprintf(stdout, "\n"); fprintf(stdout, "\n");
} }
if (operatorID == DIFFN) vlistInqVarName(vlistID1, varID1, varname);
fprintf(stdout, "%6d ", indg); fprintf(stdout, "%6d ", indg);
fprintf(stdout, ":"); fprintf(stdout, ":");
...@@ -269,9 +300,9 @@ Diff(void *process) ...@@ -269,9 +300,9 @@ Diff(void *process)
fprintf(stdout, "%s %s ", vdateString.c_str(), vtimeString.c_str()); fprintf(stdout, "%s %s ", vdateString.c_str(), vtimeString.c_str());
reset_text_color(stdout); reset_text_color(stdout);
set_text_color(stdout, GREEN); set_text_color(stdout, GREEN);
fprintf(stdout, "%7g ", cdoZaxisInqLevel(zaxisID, levelID)); fprintf(stdout, "%7g ", cdoZaxisInqLevel(varList1[varID1].zaxisID, levelID));
fprintf(stdout, "%8zu %7zu ", gridsize, std::max(nmiss1, nmiss2)); fprintf(stdout, "%8zu %7zu ", gridsize, std::max(field1.nmiss, field2.nmiss));
fprintf(stdout, "%7d ", ndiff); fprintf(stdout, "%7zu ", ndiff);
reset_text_color(stdout); reset_text_color(stdout);
fprintf(stdout, ":"); fprintf(stdout, ":");
...@@ -283,11 +314,11 @@ Diff(void *process) ...@@ -283,11 +314,11 @@ Diff(void *process)
set_text_color(stdout, BRIGHT, GREEN); set_text_color(stdout, BRIGHT, GREEN);
if (operatorID == DIFFN) if (operatorID == DIFFN)
fprintf(stdout, "%-11s", varname); fprintf(stdout, "%-11s", varList1[varID1].name);
else if (operatorID == DIFF || operatorID == DIFFP) else if (operatorID == DIFF || operatorID == DIFFP)
fprintf(stdout, "%-11s", paramstr); fprintf(stdout, "%-11s", paramstr);
else if (operatorID == DIFFC) else if (operatorID == DIFFC)
fprintf(stdout, "%4d", code); fprintf(stdout, "%4d", varList1[varID1].code);
reset_text_color(stdout); reset_text_color(stdout);
fprintf(stdout, "\n"); fprintf(stdout, "\n");
......
Markdown is supported
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

Imprint | Privacy Policy