Commit 7e026bc0 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Timstat2: OpenMP support.

parent 6d75d1c3
Pipeline #5787 passed with stages
in 14 minutes and 56 seconds
...@@ -29,31 +29,56 @@ ...@@ -29,31 +29,56 @@
// correlation in time // correlation in time
static void static void
correlation_init(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals, correlation_init(bool hasMissValues, size_t gridsize, const Varray<double> &x, const Varray<double> &y, double xmv, double ymv,
double *work0, double *work1, double *work2, double *work3, double *work4) Varray<size_t> &nofvals, Varray<double> &work0, Varray<double> &work1, Varray<double> &work2, Varray<double> &work3, Varray<double> &work4)
{ {
if (hasMissValues)
{
#ifdef _OPENMP
#pragma omp parallel for default(shared)
#endif
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
if ((!DBL_IS_EQUAL(array1[i], missval1)) && (!DBL_IS_EQUAL(array2[i], missval2))) if ((!DBL_IS_EQUAL(x[i], xmv)) && (!DBL_IS_EQUAL(y[i], ymv)))
{ {
work0[i] += array1[i]; work0[i] += x[i];
work1[i] += array2[i]; work1[i] += y[i];
work2[i] += array1[i] * array1[i]; work2[i] += x[i] * x[i];
work3[i] += array2[i] * array2[i]; work3[i] += y[i] * y[i];
work4[i] += array1[i] * array2[i]; work4[i] += x[i] * y[i];
nofvals[i]++;
}
}
}
else
{
#ifdef HAVE_OPENMP4
#pragma omp parallel for simd default(shared)
#elif _OPENMP
#pragma omp parallel for default(shared)
#endif
for (size_t i = 0; i < gridsize; ++i)
{
work0[i] += x[i];
work1[i] += y[i];
work2[i] += x[i] * x[i];
work3[i] += y[i] * y[i];
work4[i] += x[i] * y[i];
nofvals[i]++; nofvals[i]++;
} }
} }
} }
static size_t static size_t
correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals, double *work0, double *work1, double *work2, correlation(size_t gridsize, double xmv, double ymv, const Varray<size_t> &nofvals, Varray<double> &work0, const Varray<double> &work1,
double *work3, double *work4) const Varray<double> &work2, const Varray<double> &work3, const Varray<double> &work4)
{ {
size_t nmiss = 0; size_t nmiss = 0;
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
const auto missval1 = xmv;
const auto missval2 = ymv;
double cor; double cor;
const auto nvals = nofvals[i]; const auto nvals = nofvals[i];
if (nvals > 0) if (nvals > 0)
...@@ -69,12 +94,12 @@ correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals, ...@@ -69,12 +94,12 @@ correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals,
cor = DIVMN(temp1, SQRTMN(temp6)); cor = DIVMN(temp1, SQRTMN(temp6));
cor = std::min(std::max(cor, -1.0), 1.0); cor = std::min(std::max(cor, -1.0), 1.0);
if (DBL_IS_EQUAL(cor, missval1)) nmiss++; if (DBL_IS_EQUAL(cor, xmv)) nmiss++;
} }
else else
{ {
nmiss++; nmiss++;
cor = missval1; cor = xmv;
} }
work0[i] = cor; work0[i] = cor;
...@@ -85,28 +110,51 @@ correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals, ...@@ -85,28 +110,51 @@ correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals,
// covariance in time // covariance in time
static void static void
covariance_init(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals, covariance_init(bool hasMissValues, size_t gridsize, const Varray<double> &x, const Varray<double> &y, double xmv, double ymv,
double *work0, double *work1, double *work2) Varray<size_t> &nofvals, Varray<double> &work0, Varray<double> &work1, Varray<double> &work2)
{ {
if (hasMissValues)
{
#ifdef _OPENMP
#pragma omp parallel for default(shared)
#endif
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
if ((!DBL_IS_EQUAL(array1[i], missval1)) && (!DBL_IS_EQUAL(array2[i], missval2))) if ((!DBL_IS_EQUAL(x[i], xmv)) && (!DBL_IS_EQUAL(y[i], ymv)))
{
work0[i] += x[i];
work1[i] += y[i];
work2[i] += x[i] * y[i];
nofvals[i]++;
}
}
}
else
{ {
work0[i] += array1[i]; #ifdef HAVE_OPENMP4
work1[i] += array2[i]; #pragma omp parallel for simd default(shared)
work2[i] += array1[i] * array2[i]; #elif _OPENMP
#pragma omp parallel for default(shared)
#endif
for (size_t i = 0; i < gridsize; ++i)
{
work0[i] += x[i];
work1[i] += y[i];
work2[i] += x[i] * y[i];
nofvals[i]++; nofvals[i]++;
} }
} }
} }
static size_t static size_t
covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, double *work0, double *work1, double *work2) covariance(size_t gridsize, double xmv, double ymv, const Varray<size_t> &nofvals, Varray<double> &work0, const Varray<double> &work1, const Varray<double> &work2)
{ {
size_t nmiss = 0; size_t nmiss = 0;
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
const auto missval1 = xmv;
const auto missval2 = ymv;
double covar; double covar;
const auto nvals = nofvals[i]; const auto nvals = nofvals[i];
if (nvals > 0) if (nvals > 0)
...@@ -114,12 +162,12 @@ covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, d ...@@ -114,12 +162,12 @@ covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, d
double dnvals = nvals; double dnvals = nvals;
const auto temp = DIVMN(MULMN(work0[i], work1[i]), dnvals * dnvals); const auto temp = DIVMN(MULMN(work0[i], work1[i]), dnvals * dnvals);
covar = SUBMN(DIVMN(work2[i], dnvals), temp); covar = SUBMN(DIVMN(work2[i], dnvals), temp);
if (DBL_IS_EQUAL(covar, missval1)) nmiss++; if (DBL_IS_EQUAL(covar, xmv)) nmiss++;
} }
else else
{ {
nmiss++; nmiss++;
covar = missval1; covar = xmv;
} }
work0[i] = covar; work0[i] = covar;
...@@ -131,12 +179,12 @@ covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, d ...@@ -131,12 +179,12 @@ covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, d
// rms in time // rms in time
static void static void
rmsd_init(size_t gridsize, const double *x, const double *y, double xmissval, double ymissval, size_t *nofvals, rmsd_init(size_t gridsize, const Varray<double> &x, const Varray<double> &y, double xmv, double ymv, Varray<size_t> &nofvals,
double *rmsd) Varray<double> &rmsd)
{ {
for (size_t i = 0; i < gridsize; ++i) for (size_t i = 0; i < gridsize; ++i)
{ {
if ((!DBL_IS_EQUAL(x[i], xmissval)) && (!DBL_IS_EQUAL(y[i], ymissval))) if ((!DBL_IS_EQUAL(x[i], xmv)) && (!DBL_IS_EQUAL(y[i], ymv)))
{ {
rmsd[i] += ((x[i] - y[i]) * (x[i] - y[i])); rmsd[i] += ((x[i] - y[i]) * (x[i] - y[i]));
nofvals[i]++; nofvals[i]++;
...@@ -145,7 +193,7 @@ rmsd_init(size_t gridsize, const double *x, const double *y, double xmissval, do ...@@ -145,7 +193,7 @@ rmsd_init(size_t gridsize, const double *x, const double *y, double xmissval, do
} }
static size_t static size_t
rmsd_compute(size_t gridsize, double missval, size_t *nofvals, double *rmsd) rmsd_compute(size_t gridsize, double missval, const Varray<size_t> &nofvals, Varray<double> &rmsd)
{ {
size_t nmiss = 0; size_t nmiss = 0;
...@@ -170,8 +218,6 @@ Timstat2(void *process) ...@@ -170,8 +218,6 @@ Timstat2(void *process)
{ {
int64_t vdate = 0; int64_t vdate = 0;
int vtime = 0; int vtime = 0;
int varID, levelID;
size_t nmiss = 0;
cdoInitialize(process); cdoInitialize(process);
...@@ -210,7 +256,7 @@ Timstat2(void *process) ...@@ -210,7 +256,7 @@ Timstat2(void *process)
const auto taxisID3 = taxisDuplicate(taxisID1); const auto taxisID3 = taxisDuplicate(taxisID1);
if (timeIsConst) if (timeIsConst)
for (varID = 0; varID < nvars; ++varID) for (int varID = 0; varID < nvars; ++varID)
vlistDefVarTimetype(vlistID3, varID, TIME_CONSTANT); vlistDefVarTimetype(vlistID3, varID, TIME_CONSTANT);
vlistDefTaxis(vlistID3, taxisID3); vlistDefTaxis(vlistID3, taxisID3);
...@@ -223,7 +269,7 @@ Timstat2(void *process) ...@@ -223,7 +269,7 @@ Timstat2(void *process)
Varray4D<double> work(nvars); Varray4D<double> work(nvars);
Varray3D<size_t> nofvals(nvars); Varray3D<size_t> nofvals(nvars);
for (varID = 0; varID < nvars; varID++) for (int varID = 0; varID < nvars; varID++)
{ {
const auto gridsize = varList1[varID].gridsize; const auto gridsize = varList1[varID].gridsize;
const auto nlevs = varList1[varID].nlevels; const auto nlevs = varList1[varID].nlevels;
...@@ -231,7 +277,7 @@ Timstat2(void *process) ...@@ -231,7 +277,7 @@ Timstat2(void *process)
work[varID].resize(nlevs); work[varID].resize(nlevs);
nofvals[varID].resize(nlevs); nofvals[varID].resize(nlevs);
for (levelID = 0; levelID < nlevs; levelID++) for (int levelID = 0; levelID < nlevs; levelID++)
{ {
nofvals[varID][levelID].resize(gridsize, 0); nofvals[varID][levelID].resize(gridsize, 0);
work[varID][levelID].resize(nwork); work[varID][levelID].resize(nwork);
...@@ -253,6 +299,7 @@ Timstat2(void *process) ...@@ -253,6 +299,7 @@ Timstat2(void *process)
for (int recID = 0; recID < nrecs; recID++) for (int recID = 0; recID < nrecs; recID++)
{ {
int varID, levelID;
cdoInqRecord(streamID1, &varID, &levelID); cdoInqRecord(streamID1, &varID, &levelID);
cdoInqRecord(streamID2, &varID, &levelID); cdoInqRecord(streamID2, &varID, &levelID);
...@@ -266,25 +313,27 @@ Timstat2(void *process) ...@@ -266,25 +313,27 @@ Timstat2(void *process)
const auto missval1 = varList1[varID].missval; const auto missval1 = varList1[varID].missval;
const auto missval2 = varList2[varID].missval; const auto missval2 = varList2[varID].missval;
cdoReadRecord(streamID1, &array1[0], &nmiss); size_t nmiss1 = 0, nmiss2 = 0;
cdoReadRecord(streamID2, &array2[0], &nmiss); cdoReadRecord(streamID1, &array1[0], &nmiss1);
cdoReadRecord(streamID2, &array2[0], &nmiss1);
bool hasMissValues = (nmiss1 > 0 || nmiss2 > 0);
auto &rwork = work[varID][levelID]; auto &rwork = work[varID][levelID];
auto &rnofvals = nofvals[varID][levelID]; auto &rnofvals = nofvals[varID][levelID];
if (operfunc == func_cor) if (operfunc == func_cor)
{ {
correlation_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(), correlation_init(hasMissValues, gridsize, array1, array2, missval1, missval2, rnofvals,
rwork[0].data(), rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data()); rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
} }
else if (operfunc == func_covar) else if (operfunc == func_covar)
{ {
covariance_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(), covariance_init(hasMissValues, gridsize, array1, array2, missval1, missval2, rnofvals,
rwork[0].data(), rwork[1].data(), rwork[2].data()); rwork[0], rwork[1], rwork[2]);
} }
else if (operfunc == func_rmsd) else if (operfunc == func_rmsd)
{ {
rmsd_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(), rwork[0].data()); rmsd_init(gridsize, array1, array2, missval1, missval2, rnofvals, rwork[0]);
} }
} }
...@@ -298,8 +347,8 @@ Timstat2(void *process) ...@@ -298,8 +347,8 @@ Timstat2(void *process)
for (int recID = 0; recID < nrecs1; recID++) for (int recID = 0; recID < nrecs1; recID++)
{ {
varID = recVarID[recID]; const auto varID = recVarID[recID];
levelID = recLevelID[recID]; const auto levelID = recLevelID[recID];
const auto gridsize = varList1[varID].gridsize; const auto gridsize = varList1[varID].gridsize;
const auto missval1 = varList1[varID].missval; const auto missval1 = varList1[varID].missval;
...@@ -307,20 +356,19 @@ Timstat2(void *process) ...@@ -307,20 +356,19 @@ Timstat2(void *process)
auto &rwork = work[varID][levelID]; auto &rwork = work[varID][levelID];
auto &rnofvals = nofvals[varID][levelID]; auto &rnofvals = nofvals[varID][levelID];
size_t nmiss = 0;
if (operfunc == func_cor) if (operfunc == func_cor)
{ {
nmiss = correlation(gridsize, missval1, missval2, rnofvals.data(), rwork[0].data(), nmiss = correlation(gridsize, missval1, missval2, rnofvals, rwork[0], rwork[1], rwork[2], rwork[3], rwork[4]);
rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data());
} }
else if (operfunc == func_covar) else if (operfunc == func_covar)
{ {
nmiss = covariance(gridsize, missval1, missval2, rnofvals.data(), rwork[0].data(), nmiss = covariance(gridsize, missval1, missval2, rnofvals, rwork[0], rwork[1], rwork[2]);
rwork[1].data(), rwork[2].data());
} }
else if (operfunc == func_rmsd) else if (operfunc == func_rmsd)
{ {
nmiss = rmsd_compute(gridsize, missval1, rnofvals.data(), rwork[0].data()); nmiss = rmsd_compute(gridsize, missval1, rnofvals, rwork[0]);
} }
cdoDefRecord(streamID3, varID, levelID); cdoDefRecord(streamID3, varID, levelID);
......
...@@ -62,10 +62,10 @@ CheckVector ...@@ -62,10 +62,10 @@ CheckVector
T& operator[](size_t pos) { (void) pos; return ptr[0]; } T& operator[](size_t pos) { (void) pos; return ptr[0]; }
const T & operator[](size_t pos) const { (void) pos; return ptr[0]; } const T & operator[](size_t pos) const { (void) pos; return ptr[0]; }
CheckVector& operator=(const CheckVector& other) { (void)other; ptr = dummy; m_count = 0; return *this; } CheckVector& operator=(const CheckVector& other) { *this = other; ptr = dummy; m_count = 0; return *this; }
CheckVector& operator=(CheckVector&& other) { (void)other; ptr = dummy; m_count = 0; return *this; } CheckVector& operator=(CheckVector&& other) { *this = other; ptr = dummy; m_count = 0; return *this; }
// Copy constructor // Copy constructor
CheckVector(const CheckVector &v2) { m_count = 0; ptr[0] = v2.ptr[0]; } CheckVector(const CheckVector &v2) { m_count = v2.count; ptr[0] = v2.ptr[0]; }
bool operator==(const CheckVector& other) { (void)other; return true; } bool operator==(const CheckVector& other) { (void)other; return true; }
......
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