EOFs.cc 18.9 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
18
19
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

/*
   This module contains the following operators:
20

Cedrick Ansorge's avatar
Cedrick Ansorge committed
21
22
     Timeof        eof             EOF in spatial or time space
     Timeof        eofspatial      EOF in spatial space
23
     Timeof        eoftime         EOF in time space
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
*/
25
26
27
28
29
30
/*
 * TODO: 
 * Role of the weights for eofs. Should not be mixed up with division with
 * number of contributing values during summation.
 */

31
32
33
34
#if defined(_OPENMP)
#include <omp.h>
#endif

35
#include <limits.h>  // LONG_MAX
Ralf Mueller's avatar
Ralf Mueller committed
36
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
38
39
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
40
#include "grid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
#include "statistic.h"

Cedrick Ansorge's avatar
Cedrick Ansorge committed
43
44
// NO MISSING VALUE SUPPORT ADDED SO FAR

45
static
46
void scale_eigvec_grid(double *restrict out, int tsID, int npack, const int *restrict pack, const double *restrict weight, double **covar, double sum_w)
47
48
{
  for ( int i = 0; i < npack; ++i )
49
    out[pack[i]] = covar[tsID][i] / sqrt(weight[pack[i]]/sum_w);
50
51
52
}

static
53
void scale_eigvec_time(double *restrict out, int tsID, int nts, int npack, const int *restrict pack, const double *restrict weight,
54
		       double **covar, double **data, double missval, double sum_w)
55
56
{
#if defined(_OPENMP)
57
#pragma omp parallel for default(none) shared(npack, nts, tsID, pack, data, covar, out)
58
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
59
  for ( int i = 0; i < npack; ++i )
60
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
62
      double sum = 0;
      for ( int j = 0; j < nts; ++j )
63
	sum += data[j][i] * covar[tsID][j];
64
65
66
67
      
      out[pack[i]] = sum;
    }
  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
  for ( int j = 0; j < nts; ++j )
69
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
      for ( int i = 0; i < npack; ++i )
71
	out[pack[i]] += data[j][i] * covar[tsID][j];
72
73
74
    }
  */

75
  // Normalizing
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
  double sum = 0;
77
78

#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
#pragma omp parallel for default(none) reduction(+:sum)	\
80
81
  shared(out,weight,pack,npack)
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
  for ( int i = 0; i < npack; ++i )
83
    {
84
      // do not need to account for weights as eigenvectors are non-weighted                               
85
      sum += weight[pack[i]] * out[pack[i]] * out[pack[i]];
86
87
88
89
    }

  if ( sum > 0 )
    {
90
      sum = sqrt(sum/sum_w);
91
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92
#pragma omp parallel for default(none)  shared(npack,pack,sum,out)
93
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94
      for ( int i = 0; i < npack; ++i ) out[pack[i]] /= sum;
95
96
97
98
    }
  else
    {
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
#pragma omp parallel for default(none)  shared(npack,pack,out,missval)
100
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
101
      for ( int i = 0; i < npack; ++i ) out[pack[i]] = missval;
102
103
104
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
105

106
107
108
enum T_EIGEN_MODE get_eigenmode(void)
{
  enum T_EIGEN_MODE eigen_mode = JACOBI;
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

  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");
	}
    }
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

  if ( cdoVerbose ) 
    cdoPrint("Using CDO_SVD_MODE '%s' from %s",
	     eigen_mode==JACOBI?"jacobi":"danielson_lanczos",
	     envstr?"Environment":" default");  

#if defined(_OPENMP)
  if ( omp_get_max_threads() > 1 && eigen_mode == DANIELSON_LANCZOS )  {
    cdoWarning("Requested parallel computation with %i Threads ",omp_get_max_threads());
    cdoWarning("  but environmental setting CDO_SVD_MODE causes sequential ");
    cdoWarning("  Singular value decomposition");
  }
#endif 

  return eigen_mode;
}

141
142
143

enum T_WEIGHT_MODE get_weightmode(void)
{  
144
  enum T_WEIGHT_MODE weight_mode = WEIGHT_ON;
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165

  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;
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
void *EOFs(void * argument)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
{
Cedrick Ansorge's avatar
Cedrick Ansorge committed
168
  enum {EOF_, EOF_TIME, EOF_SPATIAL};
169

Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
  int nlevs = 0 ;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
171
  int nmiss;
172
  int varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
  int nts = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
  int n = 0;
  int grid_space = 0, time_space = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
  int timer_cov = 0, timer_eig = 0;
177

178
179
180
  int calendar = CALENDAR_STANDARD;
  juldate_t juldate;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
  double missval = 0;
182

Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
  typedef struct {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
185
    bool init;
    bool first_call;
186
    double *eig_val;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
188
    double *covar_array;
    double **covar;
189
    double **data;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
191
  }
  eofdata_t;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192

193
  if ( cdoTimer )
194
195
196
197
198
    {
      timer_cov  = timer_new("Timeof cov");
      timer_eig  = timer_new("Timeof eig");
    }
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
  cdoInitialize(argument);
200

Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
203
204
  cdoOperatorAdd("eof",        EOF_,        0, NULL);
  cdoOperatorAdd("eoftime",    EOF_TIME,    0, NULL);
  cdoOperatorAdd("eofspatial", EOF_SPATIAL, 0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
  // clang-format on
206

Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
  int operatorID = cdoOperatorID();
  int operfunc   = cdoOperatorF1(operatorID);
209
210

  operatorInputArg("Number of eigen functions to write out");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
  int n_eig = parameter2int(operatorArgv()[0]);
212
  
213
  enum T_EIGEN_MODE eigen_mode = get_eigenmode();
214
  enum T_WEIGHT_MODE weight_mode = get_weightmode();
215

216
217
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
  int vlistID1  = pstreamInqVlist(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
  int taxisID1  = vlistInqTaxis(vlistID1);
  int gridID1   = vlistInqVarGrid(vlistID1, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220
  int gridsize  = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
  int nvars     = vlistNvars(vlistID1);
222
  int nrecs;
223

Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
225
226
  int ngrids = vlistNgrids(vlistID1);
  for ( int index = 1; index < ngrids; index++ )
    if ( vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
      cdoAbort("Too many different grids!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228

229
  double *weight = (double *) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
  for ( int i = 0; i < gridsize; ++i ) weight[i] = 1.;
231

232
233
234
235
236
237
238
239
240
  if ( weight_mode == WEIGHT_ON )
    {
      int wstatus = gridWeights(gridID1, weight);
      if ( wstatus != 0  )
	{
	  weight_mode = WEIGHT_OFF;
	  cdoWarning("Using constant grid cell area weights!");
	}
    }
241

Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
  /* eigenvalues */
243

Cedrick Ansorge's avatar
Cedrick Ansorge committed
244
  /* COUNT NUMBER OF TIMESTEPS if EOF_ or EOF_TIME */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
  if ( operfunc == EOF_ || operfunc == EOF_TIME )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
246
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
247
      if ( cdoVerbose ) cdoPrint("Counting timesteps in ifile");
248

Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
250
251
      nts = vlistNtsteps(vlistID1);
      if ( nts == -1 )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
          nts = 0;
253
	  while ( pstreamInqTimestep(streamID1, nts) ) nts++;
254

Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
256
	  if ( cdoVerbose ) cdoPrint("Counted %i timeSteps", nts);
	}
257
258
      else
        if ( cdoVerbose ) cdoPrint("Found %i timeSteps", nts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259

260
      pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
261

262
263
      streamID1 = pstreamOpenRead(cdoStreamName(0));
      vlistID1  = pstreamInqVlist(streamID1);
264
      taxisID1  = vlistInqTaxis(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
265

Uwe Schulzweida's avatar
Uwe Schulzweida committed
266
267
268
269
270
      if ( nts < gridsize || operfunc == EOF_TIME )
	{
	  time_space = 1;
	  grid_space = 0;
	}
271
      else
Cedrick Ansorge's avatar
Cedrick Ansorge committed
272
273
274
275
276
277
278
279
280
        {
          time_space = 0;
          grid_space = 1;
        }
    }
  else if ( operfunc == EOF_SPATIAL )
    {
      time_space = 0;
      grid_space = 1;
281
282
283
284
    }

  /* reset the requested number of eigen-function to the maximum if neccessary */
  if ( time_space )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
285
286
287
288
289
    {
      if ( n_eig > nts )
        {
          cdoWarning("Solving in time-space:");
          cdoWarning("Number of eigen-functions to write out is bigger than number of time-steps.");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
          cdoWarning("Setting n_eig to %d.", nts);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
291
          cdoWarning("If You want to force a solution in grid-space use operator eofspatial");
292
293
294
          n_eig = nts;
        }
      n = nts;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
295
296
    }
  else if ( grid_space )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
297
    {
298
      if ( ((double)gridsize)*gridsize > (double)LONG_MAX ) cdoAbort("Grid space too large!");
299

Cedrick Ansorge's avatar
Cedrick Ansorge committed
300
301
      if ( n_eig > gridsize )
        {
302
303
          cdoWarning("Solving in spatial space");
          cdoWarning("Number of eigen-functions to write out is bigger than grid size");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
304
          cdoWarning("Setting n_eig to %d", gridsize);
305
          cdoWarning("If You want to force a solution in time-space use operator eoftime");
Cedrick Ansorge's avatar
Cedrick Ansorge committed
306
          n_eig = gridsize;
307
308
        }
      n = gridsize;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
309
    }
310

311
  if ( cdoVerbose ) 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312
313
    cdoPrint("Calculating %d eigenvectors and %d eigenvalues in %s",
	     n_eig, n, grid_space==1?"grid_space" : "time_space");
314

315
  /* allocation of temporary fields and output structures */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
  int npack = -1;
317
318
319
  int *pack            = (int *) Malloc(gridsize*sizeof(int));
  double *in           = (double *) Malloc(gridsize*sizeof(double));
  eofdata_t **eofdata  = (eofdata_t **) Malloc(nvars*sizeof(eofdata_t*));
320

Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
322
  for ( varID = 0; varID < nvars; ++varID )
    {
323
324
325
      gridsize = vlistGridsizeMax(vlistID1);
      nlevs    = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
      missval  = vlistInqVarMissval(vlistID1, varID);
326

327
      eofdata[varID] = (eofdata_t *) Malloc(nlevs*sizeof(eofdata_t));
328

Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
      for ( levelID = 0; levelID < nlevs; ++levelID )
330
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
332
	  eofdata[varID][levelID].init = false;
	  eofdata[varID][levelID].first_call = true;
333
	  eofdata[varID][levelID].eig_val = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
335
	  eofdata[varID][levelID].covar_array = NULL;
	  eofdata[varID][levelID].covar = NULL;
336
	  eofdata[varID][levelID].data = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337

338
	  if ( time_space )
339
	    eofdata[varID][levelID].data = (double **) Malloc(nts*sizeof(double *));
340
341
342
        }
    }

343
  if ( cdoVerbose )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
    cdoPrint("Allocated eigenvalue/eigenvector structures with nts=%d gridsize=%d", nts, gridsize);
345

Uwe Schulzweida's avatar
Uwe Schulzweida committed
346
347
  double *covar_array = NULL;
  double **covar = NULL;
348
  double sum_w = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349

Uwe Schulzweida's avatar
Uwe Schulzweida committed
350
  int tsID = 0;
351
352

  /* read the data and create covariance matrices for each var & level */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
  while ( TRUE )
354
    {
355
      nrecs = pstreamInqTimestep(streamID1, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
      if ( nrecs == 0 ) break;
357

358
      for ( int recID = 0; recID < nrecs; recID++ )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
359
        {
360
361
          pstreamInqRecord(streamID1, &varID, &levelID);
          pstreamReadRecord(streamID1, in, &nmiss);
362

Uwe Schulzweida's avatar
Uwe Schulzweida committed
363
	  gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
364
          missval = vlistInqVarMissval(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
	  if ( npack == -1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
367
	    {
	      npack = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368
	      for ( int i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369
		{
Ralf Mueller's avatar
Ralf Mueller committed
370
		  if ( !DBL_IS_EQUAL(weight[i], 0.0) && !DBL_IS_EQUAL(weight[i], missval) &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
371
		       !DBL_IS_EQUAL(in[i], missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
373
374
375
376
                    {
                      pack[npack] = i;
                      npack++;
                    }
                }
377

378
379
380
	      if ( weight_mode == WEIGHT_ON )
		{
		  sum_w = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
381
		  for ( int i = 0; i < npack; i++ ) sum_w += weight[pack[i]];
382
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
	    }
384

Uwe Schulzweida's avatar
Uwe Schulzweida committed
385
386
	  int ipack = 0;
	  for ( int i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
	    {
Ralf Mueller's avatar
Ralf Mueller committed
388
	      if ( !DBL_IS_EQUAL(weight[i], 0.0) && !DBL_IS_EQUAL(weight[i], missval) &&
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389
		   !DBL_IS_EQUAL(in[i], missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
392
                  if ( pack[ipack] != i ) cdoAbort("Missing values unsupported!");
                  ipack++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393
394
		}
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395

Uwe Schulzweida's avatar
Uwe Schulzweida committed
396
	  if ( grid_space )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
397
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398
399
	      if ( !eofdata[varID][levelID].init )
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400
401
                  n = npack;
		  double *covar_array = (double *) Malloc(((size_t)npack)*npack*sizeof(double));
402
		  covar = (double **) Malloc(npack*sizeof(double *));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403
404
405
		  for ( int i = 0; i < npack; ++i ) covar[i] = covar_array + ((size_t)npack)*i;
		  for ( int i = 0; i < npack; ++i )
                    for ( int j = 0; j < npack; ++j ) covar[i][j] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
407
408
409
410
411
412
413

		  eofdata[varID][levelID].covar_array = covar_array;
		  eofdata[varID][levelID].covar       = covar;
		}
	      else
		{
		  covar = eofdata[varID][levelID].covar;
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
414
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
	      for ( int ipack = 0; ipack < npack; ++ipack )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
418
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
		  for ( int jpack = ipack; jpack < npack; ++jpack )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
421
		    covar[ipack][jpack] += in[pack[ipack]] * in[pack[jpack]];
		}
422
	    }
423
424
          else if ( time_space )
	    {
425
	      double *data = (double *) Malloc(npack*sizeof(double));
426
427
	      eofdata[varID][levelID].data[tsID] = data;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
	      for ( int ipack = 0; ipack < npack; ipack++ )
429
		data[ipack] = in[pack[ipack]];
430
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431

Uwe Schulzweida's avatar
Uwe Schulzweida committed
432
	  eofdata[varID][levelID].init = true;
433
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
434

Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
436
      tsID++;
    }
437

Uwe Schulzweida's avatar
Uwe Schulzweida committed
438
439
  if ( grid_space ) nts = tsID;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
440
  if ( tsID == 1 ) cdoAbort("File consists of only one timestep!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441

Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
443
444
  /* write files with eigenvalues (ID3) and eigenvectors (ID2) */

  /* eigenvalues */
445
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446

Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
448
  int vlistID2 = vlistDuplicate(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449
450
451
  taxisDefRdate(taxisID2, 0);
  taxisDefRtime(taxisID2, 0);
  vlistDefTaxis(vlistID2, taxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452

Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
  int gridID2 = gridCreate(GRID_LONLAT, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
454
455
  gridDefXsize(gridID2, 1);
  gridDefYsize(gridID2, 1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
  double xvals = 0, yvals = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
458
  gridDefXvals(gridID2, &xvals);
  gridDefYvals(gridID2, &yvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
459
  for ( int i = 0; i < ngrids; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
461
462
    vlistChangeGridIndex(vlistID2, i, gridID2);

  /*  eigenvectors */
463
  int streamID3 = pstreamOpenWrite(cdoStreamName(2), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464

Uwe Schulzweida's avatar
Uwe Schulzweida committed
465
466
  int vlistID3  = vlistDuplicate(vlistID1);
  int taxisID3  = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
468
469
470
  taxisDefRdate(taxisID3, 0);
  taxisDefRtime(taxisID3, 0);
  vlistDefTaxis(vlistID3, taxisID3);

471
472
  pstreamDefVlist(streamID2, vlistID2);
  pstreamDefVlist(streamID3, vlistID3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473
474
475
476
477

  int vdate = 10101;
  int vtime = 0;
  juldate = juldate_encode(calendar, vdate, vtime);

478
479
480
  double *out = in;
  double *eig_val = NULL;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
481
  int nts_out = nts;
482
  if ( npack < nts ) nts_out = npack;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
484

  for ( tsID = 0; tsID < nts_out; tsID++ )
485
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486
487
      juldate = juldate_add_seconds(60, juldate);
      juldate_decode(calendar, juldate, &vdate, &vtime);
488

Uwe Schulzweida's avatar
Uwe Schulzweida committed
489
490
      taxisDefVdate(taxisID2, vdate);
      taxisDefVtime(taxisID2, vtime);
491
      pstreamDefTimestep(streamID2, tsID);
492

Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
      if ( tsID < n_eig )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
494
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
496
          taxisDefVdate(taxisID3, vdate);
          taxisDefVtime(taxisID3, vtime);
497
          pstreamDefTimestep(streamID3, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
498
        }
499

Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
501
502
503
504
505
      for ( varID = 0; varID < nvars; varID++ )
	{
	  char vname[256];
	  vlistInqVarName(vlistID1, varID, vname);
	  gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
	  nlevs    = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
506

Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
508
	  for ( levelID = 0; levelID < nlevs; levelID++ )
	    {
509
	      double **data = eofdata[varID][levelID].data;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
510

511
	      if ( eofdata[varID][levelID].first_call )
512
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
513
		  eofdata[varID][levelID].first_call = false;
514
515
516

		  if ( cdoVerbose )
		    cdoPrint("Calculating covar matrices for %i levels of var%i (%s)", nlevs, varID, vname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
517
518
		  
		  if ( cdoTimer ) timer_start(timer_cov);
519

Uwe Schulzweida's avatar
Uwe Schulzweida committed
520
		  if ( cdoVerbose ) cdoPrint("processing level %i",levelID);
521
		  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
522
		  if ( grid_space )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
523
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
524
		      if ( npack )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
			{
526
			  eig_val = (double *) Malloc(npack*sizeof(double));
527
			  eofdata[varID][levelID].eig_val = eig_val;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
528
			}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
529
530
531

		      covar = eofdata[varID][levelID].covar;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
		      for ( int ipack = 0; ipack < npack; ++ipack )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
			{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
534
535
536
                          int j;
			  int i = pack[ipack];
			  for ( int jpack = 0; jpack < npack; ++jpack )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
538
539
540
541
542
543
544
545
546
547
548
549
550
			    {
			      if ( jpack < ipack )
				{
				  covar[ipack][jpack] = covar[jpack][ipack];
				}
			      else
				{
				  j = pack[jpack];
				  covar[ipack][jpack] = 
				    covar[ipack][jpack] *   // covariance
				    sqrt(weight[i]) * sqrt(weight[j]) / sum_w /       // weights
				    nts;   // number of data contributing
				}
			    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551
			}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
552
553
		    }
		  else if ( time_space )
554
		    {		      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555
556
		      if ( cdoVerbose )
			cdoPrint("allocating covar with %i x %i elements | npack=%i", nts, nts, npack);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557

558
559
		      covar_array = (double *) Malloc(nts*nts*sizeof(double));
		      covar = (double **) Malloc(nts*sizeof(double *));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560
		      for ( int i = 0; i < nts; ++i ) covar[i] = covar_array + nts*i;
561

562
		      eig_val = (double *) Malloc(nts*sizeof(double));
563
564
565
		      eofdata[varID][levelID].eig_val     = eig_val;
		      eofdata[varID][levelID].covar_array = covar_array;
		      eofdata[varID][levelID].covar       = covar;
566

Uwe Schulzweida's avatar
Uwe Schulzweida committed
567
#if defined(_OPENMP)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
568
#pragma omp parallel for default(shared) schedule(dynamic)
569
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
570
		      for ( int j1 = 0; j1 < nts; j1++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
571
			{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572
573
574
                          double *df1p = data[j1];
			  for ( int j2 = 0; j2 < j1; j2++ ) covar[j1][j2] = covar[j2][j1];
			  for ( int j2 = j1; j2 < nts; j2++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
			    {
576
			      double *df2p = data[j2];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
			      double sum = 0;
			      for ( int i = 0; i < npack; i++ )
579
				sum += weight[pack[i]]*df1p[i]*df2p[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
581
582
			      covar[j1][j2] = sum / sum_w / nts;
			    }
			}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
583
584
585
586
		      
		      if ( cdoVerbose )
			cdoPrint("finished calculation of covar-matrix for var %s", vname);
		    }
587

Uwe Schulzweida's avatar
Uwe Schulzweida committed
588
		  if ( cdoTimer ) timer_stop(timer_cov);
589

Uwe Schulzweida's avatar
Uwe Schulzweida committed
590
591
		  /* SOLVE THE EIGEN PROBLEM */
		  if ( cdoTimer ) timer_start(timer_eig);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592

Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
594
		  if ( eigen_mode == JACOBI ) 
		    // TODO: use return status (>0 okay, -1 did not converge at all) 
595
		    parallel_eigen_solution_of_symmetric_matrix(covar, eig_val, n, __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
596
		  else 
597
		    eigen_solution_of_symmetric_matrix(covar, eig_val, n, __func__);
598

Uwe Schulzweida's avatar
Uwe Schulzweida committed
599
		  if ( cdoTimer ) timer_stop(timer_eig);
600
		  /* NOW: covar contains the eigenvectors, eig_val the eigenvalues */
Cedrick Ansorge's avatar
Cedrick Ansorge committed
601

Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
		  for ( int i = 0; i < gridsize; ++i ) out[i] = missval;
603
	  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604
		  // for ( int i = 0; i < n; i++ ) eig_val[i] *= sum_w;
605
606
607
608
609
610
		} // first_call
	      else
		{
		  covar   = eofdata[varID][levelID].covar;
		  eig_val = eofdata[varID][levelID].eig_val;
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611

612
613
              if ( tsID < n_eig )
                {
614
		  if      ( grid_space ) scale_eigvec_grid(out, tsID, npack, pack, weight, covar, sum_w);
615
		  else if ( time_space ) scale_eigvec_time(out, tsID, nts, npack, pack, weight, covar, data, missval, sum_w);
616
617

                  nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
618
                  for ( int i = 0; i < gridsize; i++ ) if ( DBL_IS_EQUAL(out[i], missval) ) nmiss++;
619

620
621
                  pstreamDefRecord(streamID3, varID, levelID);
                  pstreamWriteRecord(streamID3, out, nmiss);
622
		} // loop n_eig
623
624
625

	      nmiss = 0;
              if ( DBL_IS_EQUAL(eig_val[tsID], missval) ) nmiss = 1;
626
627
              pstreamDefRecord(streamID2, varID, levelID);
              pstreamWriteRecord(streamID2, &eig_val[tsID], nmiss);
628
629
	    } // loop nlevs
	} // loop nvars
Uwe Schulzweida's avatar
Uwe Schulzweida committed
630
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
631

632
  for ( varID = 0; varID < nvars; varID++)
Cedrick Ansorge's avatar
Cedrick Ansorge committed
633
    {
634
      nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
635
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
636
      for( levelID = 0; levelID < nlevs; levelID++ )
637
        { 	  
638
639
640
	  if ( eofdata[varID][levelID].eig_val ) Free(eofdata[varID][levelID].eig_val);
	  if ( eofdata[varID][levelID].covar_array ) Free(eofdata[varID][levelID].covar_array);
	  if ( eofdata[varID][levelID].covar ) Free(eofdata[varID][levelID].covar);
641
642
	  if ( time_space && eofdata[varID][levelID].data )
	    {
643
	      for ( tsID = 0; tsID < nts; tsID++ )
644
645
		if ( eofdata[varID][levelID].data[tsID] ) Free(eofdata[varID][levelID].data[tsID]);
	      Free(eofdata[varID][levelID].data);
646
647
	    }
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
648
      
649
      Free(eofdata[varID]);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
650
    }
651

652
653
654
655
  Free(eofdata);
  Free(in);
  Free(pack);
  Free(weight);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
656

657
658
659
  pstreamClose(streamID3);
  pstreamClose(streamID2);
  pstreamClose(streamID1);
660

Uwe Schulzweida's avatar
Uwe Schulzweida committed
661
662
  vlistDestroy(vlistID2);
  vlistDestroy(vlistID3);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
663
664
665
666
667
668

  gridDestroy(gridID2);

  //  taxisDestroy(taxisID2);
  //  taxisDestroy(taxisID3);

669
670
  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
671
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
}