diff --git a/src/EOFs.c b/src/EOFs.c
index 664541a43b5993853f424869f7c41a8926e5414a..804faa7abe02d15e33c494f093ab34f2bbf02f22 100644
--- a/src/EOFs.c
+++ b/src/EOFs.c
@@ -28,9 +28,6 @@
* number of contributing values during summation.
*/
-//#define OLD_IMPLEMENTATION
-#define WEIGHTS 1
-
#if defined(_OPENMP)
#include <omp.h>
#endif
@@ -49,12 +46,7 @@ static
void scale_eigvec_grid(double *restrict out, int tsID, int npack, const int *restrict pack, const double *restrict weight, double **covar, double sum_w)
{
for ( int i = 0; i < npack; ++i )
- out[pack[i]] =
-#ifdef OLD_IMPLEMENTATION
- covar[tsID][i] / sqrt(weight[pack[i]]/sum_w);
-#else
- covar[tsID][i];
-#endif
+ out[pack[i]] = covar[tsID][i] / sqrt(weight[pack[i]]/sum_w);
}
static
@@ -90,12 +82,7 @@ void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const
for ( int i = 0; i < npack; ++i )
{
// do not need to account for weights as eigenvectors are non-weighted
-#ifdef OLD_IMPLEMENTATION
- sum += weight[pack[i]] *
-#else
- sum += /*weight[pack[i]] **/
-#endif
- out[pack[i]] * out[pack[i]];
+ sum += weight[pack[i]] * out[pack[i]] * out[pack[i]];
}
if ( sum > 0 )
@@ -118,18 +105,22 @@ void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const
enum T_EIGEN_MODE get_eigenmode(void)
{
- char *envstr = getenv("CDO_SVD_MODE");
-
enum T_EIGEN_MODE eigen_mode = JACOBI;
- if ( envstr && !strncmp(envstr, "danielson_lanczos", 17) )
- eigen_mode = DANIELSON_LANCZOS;
- else if ( envstr && ! strncmp(envstr, "jacobi", 6) )
- eigen_mode = JACOBI;
- else if ( envstr ) {
- cdoWarning("Unknown environmental setting %s for CDO_SVD_MODE. Available options are",envstr);
- cdoWarning(" - 'jacobi' for a one-sided parallelized jacobi algorithm");
- cdoWarning(" - 'danielson_lanzcos' for the D/L algorithm");
- }
+
+ char *envstr = getenv("CDO_SVD_MODE");
+ if ( envstr )
+ {
+ if ( !strncmp(envstr, "danielson_lanczos", 17) )
+ eigen_mode = DANIELSON_LANCZOS;
+ else if ( !strncmp(envstr, "jacobi", 6) )
+ eigen_mode = JACOBI;
+ else
+ {
+ cdoWarning("Unknown environmental setting %s for CDO_SVD_MODE. Available options are", envstr);
+ cdoWarning(" - 'jacobi' for a one-sided parallelized jacobi algorithm");
+ cdoWarning(" - 'danielson_lanzcos' for the D/L algorithm");
+ }
+ }
if ( cdoVerbose )
cdoPrint("Using CDO_SVD_MODE '%s' from %s",
@@ -147,6 +138,31 @@ enum T_EIGEN_MODE get_eigenmode(void)
return eigen_mode;
}
+
+enum T_WEIGHT_MODE get_weightmode(void)
+{
+ enum T_WEIGHT_MODE weight_mode = WEIGHT_OFF;
+
+ char *envstr = getenv("CDO_WEIGHT_MODE");
+ if ( envstr )
+ {
+ if ( !strncmp(envstr, "off", 3) )
+ weight_mode = WEIGHT_OFF;
+ else if ( !strncmp(envstr, "on", 2) )
+ weight_mode = WEIGHT_ON;
+ else
+ cdoWarning("Unknown environmental setting %s for CDO_WEIGHT_MODE. Available options are: on/off", envstr);
+ }
+
+ if ( cdoVerbose )
+ cdoPrint("Using CDO_WEIGHT_MODE '%s' from %s",
+ weight_mode==WEIGHT_OFF?"off":"on",
+ envstr?"Environment":" default");
+
+ return weight_mode;
+}
+
+
void *EOFs(void * argument)
{
enum {EOF_, EOF_TIME, EOF_SPATIAL};
@@ -178,7 +194,7 @@ void *EOFs(void * argument)
}
eofdata_t;
- if ( cdoTimer )
+ if ( cdoTimer )
{
timer_cov = timer_new("Timeof cov");
timer_eig = timer_new("Timeof eig");
@@ -197,6 +213,7 @@ void *EOFs(void * argument)
int n_eig = parameter2int(operatorArgv()[0]);
enum T_EIGEN_MODE eigen_mode = get_eigenmode();
+ enum T_WEIGHT_MODE weight_mode = get_weightmode();
int streamID1 = streamOpenRead(cdoStreamName(0));
int vlistID1 = streamInqVlist(streamID1);
@@ -214,11 +231,17 @@ void *EOFs(void * argument)
}
double *weight = (double *) malloc(gridsize*sizeof(double));
- if ( WEIGHTS )
- gridWeights(gridID1, weight);
- else
- for ( i = 0; i < gridsize; ++i ) weight[i] = 1.;
+ for ( i = 0; i < gridsize; ++i ) weight[i] = 1.;
+ if ( weight_mode == WEIGHT_ON )
+ {
+ int wstatus = gridWeights(gridID1, weight);
+ if ( wstatus != 0 )
+ {
+ weight_mode = WEIGHT_OFF;
+ cdoWarning("Using constant grid cell area weights!");
+ }
+ }
/* eigenvalues */
tsID = 0;
@@ -332,7 +355,7 @@ void *EOFs(void * argument)
int ipack, jpack;
double *covar_array = NULL;
double **covar = NULL;
- double sum_w = 0;
+ double sum_w = 1.;
tsID = 0;
@@ -359,8 +382,11 @@ void *EOFs(void * argument)
pack[npack++] = i;
}
- sum_w = 0;
- for ( i = 0; i < npack; i++ ) sum_w += weight[pack[i]];
+ if ( weight_mode == WEIGHT_ON )
+ {
+ sum_w = 0;
+ for ( i = 0; i < npack; i++ ) sum_w += weight[pack[i]];
+ }
}
ipack = 0;
@@ -527,15 +553,10 @@ void *EOFs(void * argument)
else
{
j = pack[jpack];
-#ifdef OLD_IMPLEMENTATION
covar[ipack][jpack] =
covar[ipack][jpack] * // covariance
sqrt(weight[i]) * sqrt(weight[j]) / sum_w / // weights
nts; // number of data contributing
-#else
- covar[ipack][jpack] =
- covar[ipack][jpack] / nts; // number of data contributing
-#endif
}
}
}
@@ -565,15 +586,9 @@ void *EOFs(void * argument)
sum = 0;
double *df1p = data[j1];
double *df2p = data[j2];
-#ifdef OLD_IMPLEMENTATION
for ( i = 0; i < npack; i++ )
sum += weight[pack[i]]*df1p[i]*df2p[i];
covar[j1][j2] = sum / sum_w / nts;
-#else
- for ( i = 0; i < npack; i++ )
- sum += df1p[i]*df2p[i];
- covar[j1][j2] = sum / nts;
-#endif
}
}
diff --git a/src/cdo_int.h b/src/cdo_int.h
index 4ada68e7666eda66f70f58fecb2b7d84daa7a312..bca11ef781b180f1c28814803bdd1eeff006a1eb 100644
--- a/src/cdo_int.h
+++ b/src/cdo_int.h
@@ -69,7 +69,8 @@ char *strdup(const char *s);
#define SET_DATE(dtstr, date, time) (sprintf(dtstr, "%*d%*d", DATE_LEN-6, date, 6, time))
#define DATE_IS_NEQ(dtstr1, dtstr2, len) (memcmp(dtstr1, dtstr2, len) != 0)
-enum T_EIGEN_MODE {JACOBI, DANIELSON_LANCZOS};
+enum T_WEIGHT_MODE {WEIGHT_OFF, WEIGHT_ON};
+enum T_EIGEN_MODE {JACOBI, DANIELSON_LANCZOS};
#if defined(__xlC__) /* performance problems on IBM */
#ifndef DBL_IS_NAN