Skip to content
Snippets Groups Projects
Commit 044c386b authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

eof: cleanup

parent 44c7e28e
No related branches found
No related tags found
1 merge request!282M214003/develop
Hide whitespace changes
Inline Side-by-side
src/EOFs.cc
+ 67
− 66
View file @ 044c386b
......@@ -53,7 +53,6 @@ scale_eigvec_time(Varray<double> &out, int tsID, int nts, size_t npack, const st
{
double sum = 0.0;
for (int j = 0; j < nts; ++j) sum += data[j][i] * covar[tsID][j];
out[pack[i]] = sum;
}
/*
......@@ -118,7 +117,7 @@ public:
{
bool init = false;
bool first_call = true;
Varray<double> eig_val;
Varray<double> eigenValues;
Varray2D<double> covar;
Varray2D<double> data;
};
......@@ -126,7 +125,7 @@ public:
size_t numMissVals = 0;
int varID, levelID;
int nts = 1;
int grid_space = 0, time_space = 0;
int gridSpace = 0, timeSpace = 0;
double sum_w = 1.0;
CdoStreamID streamID1;
......@@ -143,20 +142,20 @@ public:
int calendar = CALENDAR_STANDARD;
T_WEIGHT_MODE weight_mode;
T_EIGEN_MODE eigen_mode;
T_WEIGHT_MODE weightMode;
T_EIGEN_MODE eigenMode;
size_t num_eigen_functions = 0;
size_t numEigenFunctions = 0;
int numVars;
int n_eig;
int numEigen;
int ngrids;
size_t npack = SIZE_MAX;
VarList varList1;
std::vector<size_t> pack;
Varray<double> in;
std::vector<std::vector<eofdata_t>> eofdata;
std::vector<std::vector<eofdata_t>> eofData2D;
Varray<double> weights;
......@@ -168,10 +167,10 @@ public:
auto operfunc = cdo_operator_f1(operatorID);
operator_input_arg("Number of eigen functions to write out");
n_eig = parameter_to_int(cdo_operator_argv(0));
numEigen = parameter_to_int(cdo_operator_argv(0));
eigen_mode = get_eigenmode();
weight_mode = get_weightmode();
eigenMode = get_eigenmode();
weightMode = get_weightmode();
streamID1 = cdo_open_read(0);
vlistID1 = cdo_stream_inq_vlist(streamID1);
......@@ -214,61 +213,61 @@ public:
if ((size_t) nts < gridsizemax || operfunc == EOF_TIME)
{
time_space = 1;
grid_space = 0;
timeSpace = 1;
gridSpace = 0;
}
else
{
time_space = 0;
grid_space = 1;
timeSpace = 0;
gridSpace = 1;
}
}
else if (operfunc == EOF_SPATIAL)
{
time_space = 0;
grid_space = 1;
timeSpace = 0;
gridSpace = 1;
}
// reset the requested number of eigen-function to the maximum if neccessary
if (time_space)
if (timeSpace)
{
if (n_eig > nts)
if (numEigen > nts)
{
cdo_warning("Solving in time-space:");
cdo_warning("Number of eigen-functions to write out is bigger than number of time-steps.");
cdo_warning("Setting n_eig to %d.", nts);
cdo_warning("Setting numEigen to %d.", nts);
cdo_warning("If You want to force a solution in grid-space use operator eofspatial");
n_eig = nts;
numEigen = nts;
}
num_eigen_functions = nts;
numEigenFunctions = nts;
}
else if (grid_space)
else if (gridSpace)
{
if (((double) gridsizemax) * gridsizemax > (double) SIZE_MAX) cdo_abort("Grid space too large!");
if ((size_t) n_eig > gridsizemax)
if ((size_t) numEigen > gridsizemax)
{
cdo_warning("Solving in spatial space");
cdo_warning("Number of eigen-functions to write out is bigger than grid size");
cdo_warning("Setting n_eig to %zu", gridsizemax);
cdo_warning("Setting numEigen to %zu", gridsizemax);
cdo_warning("If You want to force a solution in time-space use operator eoftime");
n_eig = gridsizemax;
numEigen = gridsizemax;
}
num_eigen_functions = gridsizemax;
numEigenFunctions = gridsizemax;
}
if (Options::cdoVerbose)
cdo_print("Calculating %d eigenvectors and %zu eigenvalues in %s", n_eig, num_eigen_functions,
(grid_space == 1) ? "grid_space" : "time_space");
cdo_print("Calculating %d eigenvectors and %zu eigenvalues in %s", numEigen, numEigenFunctions,
(gridSpace == 1) ? "grid_space" : "time_space");
weights = Varray<double>(gridsizemax, 1.0);
if (weight_mode == WEIGHT_ON)
if (weightMode == WEIGHT_ON)
{
auto wstatus = gridcell_weights(gridID1, weights);
if (wstatus != 0)
{
weight_mode = WEIGHT_OFF;
weightMode = WEIGHT_OFF;
cdo_warning("Using constant grid cell area weights!");
}
}
......@@ -277,15 +276,15 @@ public:
pack = std::vector<size_t>(gridsizemax);
in = Varray<double>(gridsizemax);
eofdata = std::vector<std::vector<eofdata_t>>(numVars);
eofData2D.resize(numVars);
for (varID = 0; varID < numVars; ++varID)
{
auto nlevs = varList1.vars[varID].nlevels;
eofdata[varID].resize(nlevs);
eofData2D[varID].resize(nlevs);
if (time_space)
for (levelID = 0; levelID < nlevs; ++levelID) eofdata[varID][levelID].data.resize(nts);
if (timeSpace)
for (levelID = 0; levelID < nlevs; ++levelID) eofData2D[varID][levelID].data.resize(nts);
}
if (Options::cdoVerbose) cdo_print("Allocated eigenvalue/eigenvector structures with nts=%d gridsize=%zu", nts, gridsizemax);
......@@ -322,7 +321,7 @@ public:
}
}
if (weight_mode == WEIGHT_ON)
if (weightMode == WEIGHT_ON)
{
sum_w = 0.0;
for (size_t i = 0; i < npack; ++i) sum_w += weights[pack[i]];
......@@ -342,15 +341,16 @@ public:
if (ipack != npack) cdo_abort("Missing values unsupported!");
}
if (grid_space)
auto &eofData = eofData2D[varID][levelID];
if (gridSpace)
{
if (!eofdata[varID][levelID].init)
if (!eofData.init)
{
eofdata[varID][levelID].covar.resize(npack);
for (size_t i = 0; i < npack; ++i) eofdata[varID][levelID].covar[i].resize(npack, 0.0);
eofData.covar.resize(npack);
for (size_t i = 0; i < npack; ++i) eofData.covar[i].resize(npack, 0.0);
}
auto &covar = eofdata[varID][levelID].covar;
auto &covar = eofData.covar;
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(npack, covar, in, pack)
#endif
......@@ -361,21 +361,21 @@ public:
for (size_t jpack = ipack; jpack < npack; ++jpack) covar_i[jpack] += in_i * in[pack[jpack]];
}
}
else if (time_space)
else if (timeSpace)
{
eofdata[varID][levelID].data[tsID].resize(npack);
auto &data = eofdata[varID][levelID].data[tsID];
eofData.data[tsID].resize(npack);
auto &data = eofData.data[tsID];
for (size_t ipack = 0; ipack < npack; ipack++) data[ipack] = in[pack[ipack]];
}
eofdata[varID][levelID].init = true;
eofData.init = true;
}
tsID++;
}
if (grid_space) nts = tsID;
if (gridSpace) nts = tsID;
if (tsID == 1) cdo_abort("File consists of only one timestep!");
......@@ -427,7 +427,7 @@ public:
taxisDefVdatetime(taxisID2, vDateTime);
cdo_def_timestep(streamID2, tsID);
if (tsID < n_eig)
if (tsID < numEigen)
{
taxisDefVdatetime(taxisID3, vDateTime);
cdo_def_timestep(streamID3, tsID);
......@@ -443,21 +443,22 @@ public:
for (levelID = 0; levelID < nlevs; ++levelID)
{
const auto &data = eofdata[varID][levelID].data;
auto &covar = eofdata[varID][levelID].covar;
auto &eofData = eofData2D[varID][levelID];
const auto &data = eofData.data;
auto &covar = eofData.covar;
if (eofdata[varID][levelID].first_call)
if (eofData.first_call)
{
eofdata[varID][levelID].first_call = false;
eofData.first_call = false;
if (Options::cdoVerbose)
cdo_print("Calculating covar matrices for %d levels of var%i (%s)", nlevs, varID + 1, vname);
if (Options::cdoVerbose) cdo_print("processing level %d", levelID + 1);
if (grid_space)
if (gridSpace)
{
if (npack) eofdata[varID][levelID].eig_val.resize(npack);
if (npack) eofData.eigenValues.resize(npack);
for (size_t ipack = 0; ipack < npack; ++ipack)
{
......@@ -475,14 +476,14 @@ public:
}
}
}
else if (time_space)
else if (timeSpace)
{
if (Options::cdoVerbose) cdo_print("allocating covar with %dx%d elements | npack=%zu", nts, nts, npack);
covar.resize(nts);
for (int i = 0; i < nts; ++i) covar[i].resize(nts);
eofdata[varID][levelID].eig_val.resize(nts);
eofData.eigenValues.resize(nts);
#ifdef _OPENMP
#pragma omp parallel for default(none) shared(nts, data, covar, weights, npack, pack, sum_w) schedule(static)
......@@ -503,12 +504,12 @@ public:
}
// Solve the eigen problem
auto &eig_val = eofdata[varID][levelID].eig_val;
if (eigen_mode == JACOBI)
auto &eigenValues = eofData.eigenValues;
if (eigenMode == JACOBI)
// TODO: use return status (>0 okay, -1 did not converge at all)
parallel_eigen_solution_of_symmetric_matrix(covar, eig_val, num_eigen_functions, __func__);
parallel_eigen_solution_of_symmetric_matrix(covar, eigenValues, numEigenFunctions, __func__);
else
eigen_solution_of_symmetric_matrix(covar, eig_val, num_eigen_functions, __func__);
eigen_solution_of_symmetric_matrix(covar, eigenValues, numEigenFunctions, __func__);
// NOW: covar contains the eigenvectors, eig_val the eigenvalues
......@@ -517,24 +518,24 @@ public:
// for ( int i = 0; i < n; i++ ) eig_val[i] *= sum_w;
} // first_call
if (tsID < n_eig)
if (tsID < numEigen)
{
if (grid_space)
if (gridSpace)
scale_eigvec_grid(out, tsID, npack, pack, weights, covar, sum_w);
else if (time_space)
else if (timeSpace)
scale_eigvec_time(out, tsID, nts, npack, pack, weights, covar, data, missval, sum_w);
numMissVals = varray_num_mv(gridsize, out, missval);
cdo_def_field(streamID3, varID, levelID);
cdo_write_field(streamID3, out.data(), numMissVals);
} // loop n_eig
} // loop numEigen
auto eig_val = eofdata[varID][levelID].eig_val.data();
auto eigenValues = eofData.eigenValues.data();
numMissVals = (fp_is_equal(eig_val[tsID], missval)) ? 1 : 0;
numMissVals = (fp_is_equal(eigenValues[tsID], missval)) ? 1 : 0;
cdo_def_field(streamID2, varID, levelID);
cdo_write_field(streamID2, &eig_val[tsID], numMissVals);
cdo_write_field(streamID2, &eigenValues[tsID], numMissVals);
} // loop nlevs
} // loop nvars
}
......
src/eof_mode.cc
+ 1
− 2
View file @ 044c386b
......@@ -34,8 +34,7 @@ get_eigenmode(void)
if (omp_get_max_threads() > 1 && eigen_mode == DANIELSON_LANCZOS)
{
cdo_warning("Requested parallel computation with %i Threads ", omp_get_max_threads());
cdo_warning(" but environmental setting CDO_SVD_MODE causes sequential ");
cdo_warning(" Singular value decomposition");
cdo_warning(" but environmental setting CDO_SVD_MODE causes sequential Singular value decomposition");
}
#endif
......
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