« back to all changes in this revision
Viewing changes to src/read_bp.c
- browse files at revision 15
- compare with another revision
- download diff
- download tarball
- view history from revision 15
- Committer: Package Import Robot
- Author(s): Alastair McKinstry
- Date: 2013-12-09 15:21:31 UTC
- mfrom: (1.1.2)
- Revision ID: package-import@ubuntu.com-20131209152131-jtd4fpmdv3xnunnm
Tags: 1.5.0-1
* New upstream.
* Standards-Version: 3.9.5
* Include latest config.{sub,guess}
* New watch file.
* Create libadios-bin for binaries.
* Standards-Version: 3.9.5
* Include latest config.{sub,guess}
* New watch file.
* Create libadios-bin for binaries.
- files added:
- .pc/debian_paths.patch
- .pc/debian_paths.patch/utils
- .pc/debian_paths.patch/utils/gpp
- .pc/debian_paths.patch/utils/skel
- .pc/debian_paths.patch/utils/skel/bin
- .pc/debian_paths.patch/utils/skel/lib
- .pc/mpi_in_place.patch/src/write
- .pc/python3.patch
- .pc/python3.patch/utils
- .pc/python3.patch/utils/gpp
- .pc/python3.patch/utils/skel
- .pc/python3.patch/utils/skel/bin
- .pc/python3.patch/utils/skel/lib
- .pc/python3.patch/wrappers
- .pc/python3.patch/wrappers/numpy
- .pc/python3.patch/wrappers/numpy/example
- .pc/python3.patch/wrappers/numpy/test
- examples/C/flexpath_arrays
- examples/C/flexpath_arrays/global_range_select
- examples/C/flexpath_arrays/process_select
- examples/C/schema
- examples/coupling
- examples/skel
- examples/skel/genarray3d
- examples/skel/geos5
- examples/skel/grapes
- examples/skel/gts
- examples/skel/s3d
- examples/staging
- examples/staging/stage_write
- examples/staging/stage_write/x.bp.dir
- src/core
- src/nssi
- src/public
- src/read
- src/write
- utils/skel
- utils/skel/bin
- utils/skel/etc
- utils/skel/etc/templates
- utils/skel/lib
- utils/skel/src
- wrappers
- wrappers/java
- wrappers/java/Modules
- wrappers/java/gov
- wrappers/java/gov/ornl
- wrappers/java/gov/ornl/ccs
- wrappers/java/test
- wrappers/matlab
- wrappers/numpy
- wrappers/numpy/Modules
- wrappers/numpy/example
- wrappers/numpy/test
- files removed:
- .pc/errno.patch
- .pc/errno.patch/src
- .pc/nompi_flag.patch
- .pc/nompi_flag.patch/utils
- .pc/nompi_flag.patch/utils/bp2h5
- tools
- tools/matlab
- files modified:
- .pc/applied-patches
added
removed
src/read_bp.c
1
/*
2
* ADIOS is freely available under the terms of the BSD license described
3
* in the COPYING file in the top level directory of this source distribution.
4
*
5
* Copyright (c) 2008 - 2009. UT-BATTELLE, LLC. All rights reserved.
6
*/
7
8
9
/****************************/
10
/* Read method for BP files */
11
/****************************/
12
13
#include <stdlib.h>
14
#include <string.h>
15
#include <math.h>
16
#include "adios.h"
17
#include "bp_utils.h"
18
#include "bp_types.h"
19
#include "adios_types.h"
20
#include "adios_read.h"
21
#include "adios_read_hooks.h"
22
#include "adios_error.h"
23
#include "futils.h"
24
25
#ifdef DMALLOC
26
#include "dmalloc.h"
27
#endif
28
29
MPI_File * get_BP_file_handle(struct BP_file_handle * l, uint32_t file_index)
30
{
31
if (!l)
32
return 0;
33
34
while (l)
35
{
36
if (l->file_index == file_index)
37
return &l->fh;
38
39
l = l->next;
40
}
41
42
return 0;
43
}
44
45
void add_BP_file_handle (struct BP_file_handle ** l, struct BP_file_handle * n)
46
{
47
if (!n)
48
return;
49
50
n->next = *l;
51
*l = n;
52
}
53
54
55
void close_all_BP_files (struct BP_file_handle * l)
56
{
57
struct BP_file_handle * n;
58
59
while (l)
60
{
61
n = l->next;
62
63
MPI_File_close (&l->fh);
64
free (l);
65
66
l = n;
67
}
68
}
69
70
/* Return 0: if file is little endian, 1 if file is big endian
71
* We know if it is different from the current system, so here
72
* we determine the current endianness and report accordingly.
73
*/
74
static int adios_read_bp_get_endianness( uint32_t change_endianness )
75
{
76
int LE = 0;
77
int BE = !LE;
78
int i = 1;
79
char *p = (char *) &i;
80
int current_endianness;
81
if (p[0] == 1) // Lowest address contains the least significant byte
82
current_endianness = LE;
83
else
84
current_endianness = BE;
85
if (change_endianness == adios_flag_yes)
86
return !current_endianness;
87
else
88
return current_endianness;
89
}
90
91
int adios_read_bp_init (MPI_Comm comm) { return 0; }
92
int adios_read_bp_finalize () { return 0; }
93
94
ADIOS_FILE * adios_read_bp_fopen (const char * fname, MPI_Comm comm)
95
{
96
int i, rank;
97
struct BP_FILE * fh;
98
ADIOS_FILE * fp;
99
uint64_t header_size;
100
101
adios_errno = 0;
102
fh = (struct BP_FILE *) malloc (sizeof (struct BP_FILE));
103
if (!fh) {
104
adios_error ( err_no_memory, "Cannot allocate memory for file info.");
105
return NULL;
106
}
107
108
fh->fname = (fname ? strdup (fname) : 0L);
109
fh->sfh = 0;
110
fh->comm = comm;
111
fh->gvar_h = 0;
112
fh->pgs_root = 0;
113
fh->vars_root = 0;
114
fh->attrs_root = 0;
115
fh->b = malloc (sizeof (struct adios_bp_buffer_struct_v1));
116
if (!fh->b) {
117
adios_error ( err_no_memory, "Cannot allocate memory for file info.");
118
return NULL;
119
}
120
fp = (ADIOS_FILE *) malloc (sizeof (ADIOS_FILE));
121
if (!fp) {
122
adios_error ( err_no_memory, "Cannot allocate memory for file info.");
123
return NULL;
124
}
125
126
adios_buffer_struct_init (fh->b);
127
MPI_Comm_rank (comm, &rank);
128
if (bp_read_open (fname, comm, fh))
129
return NULL;
130
131
/* Only rank=0 process reads the footer and it broadcasts to all other processes */
132
if ( rank == 0 ) {
133
if (bp_read_minifooter (fh))
134
return NULL;
135
}
136
MPI_Bcast (&fh->mfooter, sizeof(struct bp_minifooter),MPI_BYTE, 0, comm);
137
138
header_size = fh->mfooter.file_size-fh->mfooter.pgs_index_offset;
139
140
if ( rank != 0) {
141
if (!fh->b->buff) {
142
bp_alloc_aligned (fh->b, header_size);
143
if (!fh->b->buff)
144
return NULL;
145
memset (fh->b->buff, 0, header_size);
146
fh->b->offset = 0;
147
}
148
}
149
MPI_Barrier(comm);
150
MPI_Bcast (fh->b->buff, fh->mfooter.file_size-fh->mfooter.pgs_index_offset, MPI_BYTE, 0 , comm);
151
152
/* Everyone parses the index on its own */
153
bp_parse_pgs (fh);
154
bp_parse_vars (fh);
155
bp_parse_attrs (fh);
156
157
/* fill out ADIOS_FILE struct */
158
fp->fh = (uint64_t) fh;
159
fp->groups_count = fh->gvar_h->group_count;
160
fp->vars_count = fh->mfooter.vars_count;
161
fp->attrs_count = fh->mfooter.attrs_count;
162
fp->tidx_start = fh->tidx_start;
163
fp->ntimesteps = fh->tidx_stop - fh->tidx_start + 1;
164
fp->file_size = fh->mfooter.file_size;
165
fp->version = fh->mfooter.version;
166
fp->endianness = adios_read_bp_get_endianness(fh->mfooter.change_endianness);
167
alloc_namelist (&fp->group_namelist,fp->groups_count);
168
for (i=0;i<fp->groups_count;i++) {
169
if (!fp->group_namelist[i]) {
170
adios_error (err_no_memory, "Could not allocate buffer for %d strings in adios_fopen()", fp->groups_count);
171
adios_read_bp_fclose(fp);
172
return NULL;
173
}
174
else {
175
strcpy(fp->group_namelist[i],fh->gvar_h->namelist[i]);
176
}
177
}
178
return fp;
179
}
180
181
/* This function can be called if user places
182
the wrong sequences of dims for a var
183
*/
184
void adios_read_bp_reset_dimension_order (ADIOS_FILE *fp, int is_fortran)
185
{
186
struct BP_FILE * fh = (struct BP_FILE *)(fp->fh);
187
struct bp_index_pg_struct_v1 ** root = &(fh->pgs_root);
188
struct bp_minifooter * mh = &(fh->mfooter);
189
uint64_t i;
190
191
for (i = 0; i < mh->pgs_count; i++) {
192
is_fortran ? ((*root)->adios_host_language_fortran = adios_flag_yes)
193
: ((*root)->adios_host_language_fortran = adios_flag_no);
194
root = &(*root)->next;
195
}
196
}
197
198
int adios_read_bp_fclose (ADIOS_FILE *fp)
199
{
200
struct BP_FILE * fh = (struct BP_FILE *) fp->fh;
201
struct BP_GROUP_VAR * gh = fh->gvar_h;
202
struct BP_GROUP_ATTR * ah = fh->gattr_h;
203
struct adios_index_var_struct_v1 * vars_root = fh->vars_root, *vr;
204
struct adios_index_attribute_struct_v1 * attrs_root = fh->attrs_root, *ar;
205
struct bp_index_pg_struct_v1 * pgs_root = fh->pgs_root, *pr;
206
int i,j;
207
MPI_File mpi_fh = fh->mpi_fh;
208
209
adios_errno = 0;
210
if (fh->mpi_fh)
211
MPI_File_close (&mpi_fh);
212
213
if (fh->sfh)
214
close_all_BP_files (fh->sfh);
215
216
if (fh->b) {
217
adios_posix_close_internal (fh->b);
218
free(fh->b);
219
}
220
221
/* Free variable structures */
222
/* alloc in bp_utils.c: bp_parse_vars() */
223
while (vars_root) {
224
vr = vars_root;
225
vars_root = vars_root->next;
226
for (j = 0; j < vr->characteristics_count; j++) {
227
// alloc in bp_utils.c:bp_parse_characteristics() <- bp_get_characteristics_data()
228
if (vr->characteristics[j].dims.dims)
229
free (vr->characteristics[j].dims.dims);
230
if (vr->characteristics[j].value)
231
free (vr->characteristics[j].value);
232
// NCSU - Clearing up statistics
233
if (vr->characteristics[j].stats)
234
{
235
uint8_t k = 0, idx = 0;
236
uint8_t i = 0, count = adios_get_stat_set_count(vr->type);
237
238
while (vr->characteristics[j].bitmap >> k)
239
{
240
if ((vr->characteristics[j].bitmap >> k) & 1)
241
{
242
for (i = 0; i < count; i ++)
243
{
244
if (k == adios_statistic_hist)
245
{
246
struct adios_index_characteristics_hist_struct * hist = (struct adios_index_characteristics_hist_struct *) (vr->characteristics [j].stats[i][idx].data);
247
free (hist->breaks);
248
free (hist->frequencies);
249
free (hist);
250
}
251
else
252
free (vr->characteristics[j].stats [i][idx].data);
253
}
254
idx ++;
255
}
256
k ++;
257
}
258
259
for (i = 0; i < count; i ++)
260
free (vr->characteristics[j].stats [i]);
261
262
free (vr->characteristics[j].stats);
263
vr->characteristics[j].stats = 0;
264
}
265
}
266
if (vr->characteristics)
267
free (vr->characteristics);
268
if (vr->group_name)
269
free (vr->group_name);
270
if (vr->var_name)
271
free (vr->var_name);
272
if (vr->var_path)
273
free (vr->var_path);
274
free(vr);
275
}
276
277
/* Free attributes structures */
278
/* alloc in bp_utils.c bp_parse_attrs() */
279
while (attrs_root) {
280
ar = attrs_root;
281
attrs_root = attrs_root->next;
282
for (j = 0; j < ar->characteristics_count; j++) {
283
if (ar->characteristics[j].value)
284
free (ar->characteristics[j].value);
285
}
286
if (ar->characteristics)
287
free (ar->characteristics);
288
if (ar->group_name)
289
free (ar->group_name);
290
if (ar->attr_name)
291
free (ar->attr_name);
292
if (ar->attr_path)
293
free (ar->attr_path);
294
free(ar);
295
}
296
297
298
/* Free process group structures */
299
/* alloc in bp_utils.c bp_parse_pgs() first loop */
300
//printf ("pgs: %d\n", fh->mfooter.pgs_count);
301
while (pgs_root) {
302
pr = pgs_root;
303
pgs_root = pgs_root->next;
304
//printf("%d\tpg pid=%d addr=%x next=%x\n",i, pr->process_id, pr, pr->next);
305
if (pr->group_name)
306
free(pr->group_name);
307
if (pr->time_index_name)
308
free(pr->time_index_name);
309
free(pr);
310
}
311
312
/* Free variable structures in BP_GROUP_VAR */
313
if (gh) {
314
for (j=0;j<2;j++) {
315
for (i=0;i<gh->group_count;i++) {
316
if (gh->time_index[j][i])
317
free(gh->time_index[j][i]);
318
}
319
if (gh->time_index[j])
320
free(gh->time_index[j]);
321
}
322
free (gh->time_index);
323
324
for (i=0;i<gh->group_count;i++) {
325
if (gh->namelist[i])
326
free(gh->namelist[i]);
327
}
328
if (gh->namelist)
329
free (gh->namelist);
330
331
for (i=0;i<fh->mfooter.vars_count;i++) {
332
if (gh->var_namelist[i])
333
free(gh->var_namelist[i]);
334
if (gh->var_offsets[i])
335
free(gh->var_offsets[i]);
336
}
337
if (gh->var_namelist)
338
free (gh->var_namelist);
339
340
if (gh->var_offsets)
341
free(gh->var_offsets);
342
343
if (gh->var_counts_per_group)
344
free(gh->var_counts_per_group);
345
346
if (gh->pg_offsets)
347
free (gh->pg_offsets);
348
349
free (gh);
350
}
351
352
/* Free attribute structures in BP_GROUP_ATTR */
353
if (ah) {
354
for (i = 0; i < fh->mfooter.attrs_count; i++) {
355
if (ah->attr_offsets[i])
356
free(ah->attr_offsets[i]);
357
if (ah->attr_namelist[i])
358
free(ah->attr_namelist[i]);
359
}
360
if (ah->attr_offsets)
361
free(ah->attr_offsets);
362
if (ah->attr_namelist)
363
free(ah->attr_namelist);
364
if (ah->attr_counts_per_group)
365
free(ah->attr_counts_per_group);
366
367
free(ah);
368
}
369
370
if (fh->fname)
371
free (fh->fname);
372
373
if (fh)
374
free (fh);
375
376
free_namelist ((fp->group_namelist),fp->groups_count);
377
free(fp);
378
return 0;
379
}
380
381
382
ADIOS_GROUP * adios_read_bp_gopen (ADIOS_FILE *fp, const char * grpname)
383
{
384
struct BP_FILE * fh = (struct BP_FILE *) fp->fh;
385
int grpid;
386
387
adios_errno = 0;
388
for (grpid=0;grpid<(fh->gvar_h->group_count);grpid++) {
389
if (!strcmp(fh->gvar_h->namelist[grpid], grpname))
390
break;
391
}
392
if (grpid >= fh->gvar_h->group_count) {
393
adios_error ( err_invalid_group, "Invalid group name %s", grpname);
394
return NULL;
395
}
396
return adios_read_bp_gopen_byid(fp, grpid);
397
}
398
399
ADIOS_GROUP * adios_read_bp_gopen_byid (ADIOS_FILE *fp, int grpid)
400
{
401
struct BP_FILE * fh = (struct BP_FILE *) fp->fh;
402
struct BP_GROUP * gh;
403
ADIOS_GROUP * gp;
404
int i, offset;
405
406
adios_errno = 0;
407
if (grpid < 0 || grpid >= fh->gvar_h->group_count) {
408
adios_error ( err_invalid_group, "Invalid group index %d", grpid);
409
return NULL;
410
}
411
412
gh = (struct BP_GROUP *) malloc(sizeof(struct BP_GROUP));
413
if (!gh) {
414
adios_error ( err_no_memory, "Could not allocate memory for group info");
415
return NULL;
416
}
417
418
gp = (ADIOS_GROUP *) malloc(sizeof(ADIOS_GROUP));
419
if (!gp) {
420
adios_error ( err_no_memory, "Could not allocate memory for group info");
421
free(gh);
422
return NULL;
423
}
424
425
/* set offset index of variables (which is a long list of all vars in all groups) in this group */
426
offset = 0;
427
for (i=0; i<grpid; i++)
428
offset += fh->gvar_h->var_counts_per_group[i];
429
430
/* gh->vars_root will point to the list of vars in this group */
431
gh->vars_root = fh->vars_root;
432
for (i=0; i<offset; i++)
433
gh->vars_root = gh->vars_root->next;
434
435
gh->group_id = grpid;
436
gh->vars_offset = offset;
437
gh->vars_count = fh->gvar_h->var_counts_per_group[grpid];
438
439
/* set offset of attributes in this group */
440
offset = 0;
441
for(i=0;i<grpid;i++)
442
offset += fh->gattr_h->attr_counts_per_group[i];
443
444
/* gh->attrs_root will point to the list of vars in this group */
445
gh->attrs_root = fh->attrs_root;
446
for (i=0; i<offset; i++)
447
gh->attrs_root = gh->attrs_root->next;
448
449
gh->attrs_offset = offset;
450
gh->attrs_count = fh->gattr_h->attr_counts_per_group[grpid];
451
452
gh->fh = fh;
453
454
/* fill out ADIOS_GROUP struct */
455
gp->grpid = grpid;
456
gp->gh = (uint64_t) gh;
457
gp->fp = fp;
458
gp->vars_count = gh->vars_count;
459
gp->attrs_count = gh->attrs_count;
460
461
offset = gh->vars_offset;
462
alloc_namelist (&(gp->var_namelist), gp->vars_count);
463
for (i=0;i<gp->vars_count;i++) {
464
if (!gp->var_namelist[i]) {
465
adios_error (err_no_memory, "Could not allocate buffer for %d strings in adios_gopen()", gp->vars_count);
466
adios_read_bp_gclose(gp);
467
return NULL;
468
}
469
else
470
strcpy(gp->var_namelist[i], gh->fh->gvar_h->var_namelist[i+offset]);
471
}
472
473
offset = gh->attrs_offset;
474
alloc_namelist (&(gp->attr_namelist), gp->attrs_count);
475
for (i=0;i<gp->attrs_count;i++) {
476
if (!gp->attr_namelist[i]) {
477
adios_error (err_no_memory, "Could not allocate buffer for %d strings in adios_gopen()", gp->vars_count);
478
adios_read_bp_gclose(gp);
479
return NULL;
480
}
481
else {
482
strcpy(gp->attr_namelist[i], gh->fh->gattr_h->attr_namelist[i+offset]);
483
}
484
}
485
486
return gp;
487
}
488
489
int adios_read_bp_gclose (ADIOS_GROUP *gp)
490
{
491
struct BP_GROUP * gh = (struct BP_GROUP *) gp->gh;
492
493
adios_errno = 0;
494
if (!gh) {
495
adios_error (err_invalid_group_struct, "group handle is NULL!");
496
return err_invalid_group_struct;
497
}
498
else
499
free (gh);
500
501
free_namelist ((gp->var_namelist),gp->vars_count);
502
free_namelist ((gp->attr_namelist),gp->attrs_count);
503
free(gp);
504
return 0;
505
}
506
507
508
509
int adios_read_bp_get_attr (ADIOS_GROUP * gp, const char * attrname, enum ADIOS_DATATYPES * type,
510
int * size, void ** data)
511
{
512
// Find the attribute: full path is stored with a starting /
513
// Like in HDF5, we need to match names given with or without the starting /
514
// startpos is 0 or 1 to indicate if the argument has starting / or not
515
int attrid;
516
int vstartpos = 0, fstartpos = 0;
517
struct BP_GROUP * gh = (struct BP_GROUP *)gp->gh;
518
int offset;
519
520
adios_errno = 0;
521
if (!gp) {
522
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_get_attr()");
523
return adios_errno;
524
}
525
if (!attrname) {
526
adios_error (err_invalid_attrname, "Null pointer passed as attribute name to adios_get_attr()!");
527
return adios_errno;
528
}
529
530
offset = gh->attrs_offset;
531
532
if (attrname[0] == '/')
533
vstartpos = 1;
534
for (attrid=0; attrid<(gp->attrs_count);attrid++) {
535
//if (gp->attr_namelist[attrid][0] == '/')
536
if (gh->fh->gattr_h->attr_namelist[attrid+offset][0] == '/')
537
fstartpos = 1;
538
//if (!strcmp(gp->attr_namelist[attrid]+fstartpos, attrname+vstartpos))
539
if (!strcmp(gh->fh->gattr_h->attr_namelist[attrid+offset]+fstartpos, attrname+vstartpos))
540
break;
541
}
542
if (attrid >= gp->attrs_count) {
543
adios_error ( err_invalid_attrname, "Invalid attribute name %s", attrname);
544
return adios_errno;
545
}
546
547
return adios_read_bp_get_attr_byid(gp, attrid, type, size, data);
548
}
549
550
int adios_read_bp_get_attr_byid (ADIOS_GROUP * gp, int attrid,
551
enum ADIOS_DATATYPES * type, int * size, void ** data)
552
{
553
int i, offset, count;
554
struct BP_GROUP * gh;
555
struct BP_FILE * fh;
556
struct adios_index_attribute_struct_v1 * attr_root;
557
struct adios_index_var_struct_v1 * var_root;
558
int file_is_fortran;
559
560
adios_errno = 0;
561
if (!gp) {
562
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_get_attr()");
563
return adios_errno;
564
}
565
gh = (struct BP_GROUP *) gp->gh;
566
if (!gh) {
567
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
568
return adios_errno;
569
}
570
fh = gh->fh;
571
if (!fh) {
572
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
573
return adios_errno;
574
}
575
if (attrid < 0 || attrid >= gh->attrs_count) {
576
adios_error (err_invalid_attrid, "Invalid attribute id %d (allowed 0..%d)", attrid, gh->attrs_count);
577
return adios_errno;
578
}
579
580
attr_root = gh->attrs_root; /* need to traverse the attribute list of the group */
581
for (i = 0; i < attrid && attr_root; i++)
582
attr_root = attr_root->next;
583
if (i != attrid) {
584
adios_error (err_corrupted_attribute, "Attribute id=%d is valid but was not found in internal data structures!",attrid);
585
return adios_errno;
586
}
587
588
file_is_fortran = (fh->pgs_root->adios_host_language_fortran == adios_flag_yes);
589
590
if (attr_root->characteristics[0].value) {
591
/* Attribute has its own value */
592
*size = bp_get_type_size (attr_root->type, attr_root->characteristics[0].value);
593
*type = attr_root->type;
594
*data = (void *) malloc (*size);
595
if (*data)
596
memcpy(*data, attr_root->characteristics[0].value, *size);
597
}
598
else if (attr_root->characteristics[0].var_id) {
599
/* Attribute is a reference to a variable */
600
/* FIXME: var ids are not unique in BP. If a group of variables are written several
601
times under different path using adios_set_path(), the id of a variable is always
602
the same (should be different). As a temporary fix, we look first for a matching
603
id plus path between an attribute and a variable. If not found, then we look for
604
a match on the ids only.*/
605
var_root = gh->vars_root;
606
while (var_root) {
607
if (var_root->id == attr_root->characteristics[0].var_id &&
608
!strcmp(var_root->var_path, attr_root->attr_path))
609
break;
610
var_root = var_root->next;
611
}
612
if (!var_root) {
613
var_root = gh->vars_root;
614
while (var_root) {
615
if (var_root->id == attr_root->characteristics[0].var_id)
616
break;
617
var_root = var_root->next;
618
}
619
}
620
621
if (!var_root) {
622
adios_error (err_invalid_attribute_reference,
623
"Attribute %s/%s in group %s is a reference to variable ID %d, which is not found",
624
attr_root->attr_path, attr_root->attr_name, attr_root->group_name,
625
attr_root->characteristics[0].var_id);
626
return adios_errno;
627
}
628
629
/* default values in case of error */
630
*data = NULL;
631
*size = 0;
632
*type = attr_root->type;
633
634
/* FIXME: variable and attribute type may not match, then a conversion is needed. */
635
/* Cases:
636
1. attr has no type, var is byte array ==> string
637
2. attr has no type, var is not byte array ==> var type
638
3. attr is string, var is byte array ==> string
639
4. attr type == var type ==> var type
640
5. attr type != var type ==> attr type and conversion needed
641
*/
642
/* Error check: attr cannot reference an array in general */
643
if (var_root->characteristics[0].dims.count > 0) {
644
if ( (var_root->type == adios_byte || var_root->type == adios_unsigned_byte) &&
645
(attr_root->type == adios_unknown || attr_root->type == adios_string) &&
646
(var_root->characteristics[0].dims.count == 1)) {
647
; // this conversions are allowed
648
} else {
649
adios_error (err_invalid_attribute_reference,
650
"Attribute %s/%s in group %s, typeid=%d is a reference to an %d-dimensional array variable "
651
"%s/%s of type %s, which is not supported in ADIOS",
652
attr_root->attr_path, attr_root->attr_name, attr_root->group_name, attr_root->type,
653
var_root->characteristics[0].dims.count,
654
var_root->var_path, var_root->var_name, common_read_type_to_string(var_root->type));
655
return adios_errno;
656
}
657
}
658
659
if ( (attr_root->type == adios_unknown || attr_root->type == adios_string) &&
660
(var_root->type == adios_byte || var_root->type == adios_unsigned_byte) &&
661
(var_root->characteristics[0].dims.count == 1) ) {
662
/* 1D byte arrays are converted to string */
663
/* 1. read in variable */
664
char varname[512];
665
char *tmpdata;
666
uint64_t start, count;
667
int status;
668
start = 0;
669
count = var_root->characteristics[0].dims.dims[0];
670
snprintf(varname, 512, "%s/%s", var_root->var_path, var_root->var_name);
671
tmpdata = (char *) malloc (count+1);
672
if (tmpdata == NULL) {
673
adios_error (err_no_memory,
674
"Cannot allocate memory of %lld bytes for reading in data for attribute %s/%s of group %s.",
675
count, attr_root->attr_path, attr_root->attr_name, attr_root->group_name);
676
return adios_errno;
677
}
678
679
status = adios_read_bp_read_var (gp, varname, &start, &count, tmpdata);
680
681
if (status < 0) {
682
char *msg = strdup(adios_get_last_errmsg());
683
adios_error ((enum ADIOS_ERRCODES) status,
684
"Cannot read data of variable %s/%s for attribute %s/%s of group %s: %s",
685
var_root->var_path, var_root->var_name,
686
attr_root->attr_path, attr_root->attr_name, attr_root->group_name,
687
msg);
688
free(tmpdata);
689
free(msg);
690
return status;
691
}
692
693
*type = adios_string;
694
if (file_is_fortran) {
695
/* Fortran byte array to C string */
696
*data = futils_fstr_to_cstr( tmpdata, (int)count); /* FIXME: supports only 2GB strings... */
697
*size = strlen( (char *)data );
698
free(tmpdata);
699
} else {
700
/* C array to C string */
701
tmpdata[count] = '\0';
702
*size = count+1;
703
*data = tmpdata;
704
}
705
} else {
706
/* other types are inherited */
707
*type = var_root->type;
708
*size = bp_get_type_size (var_root->type, var_root->characteristics[0].value);
709
*data = (void *) malloc (*size);
710
if (*data)
711
memcpy(*data, var_root->characteristics[0].value, *size);
712
}
713
}
714
715
return 0;
716
}
717
718
719
/* Reverse the order in an array in place.
720
use swapping from Fortran/column-major order to ADIOS-read-api/C/row-major order and back
721
*/
722
static void swap_order(int n, uint64_t *array, int *timedim)
723
{
724
int i;
725
uint64_t tmp;
726
for (i=0; i<n/2; i++) {
727
tmp = array[i];
728
array[i] = array[n-1-i];
729
array[n-1-i] = tmp;
730
}
731
if (*timedim > -1)
732
*timedim = (n-1) - *timedim; // swap the time dimension too
733
}
734
735
/* Look up variable id based on variable name.
736
Return index 0..gp->vars_count-1 if found, -1 otherwise
737
*/
738
static int adios_read_bp_find_var(ADIOS_GROUP *gp, const char *varname)
739
{
740
// Find the variable: full path is stored with a starting /
741
// Like in HDF5, we need to match names given with or without the starting /
742
// startpos is 0 or 1 to indicate if the argument has starting / or not
743
int varid;
744
int vstartpos = 0, fstartpos;
745
struct BP_GROUP * gh = (struct BP_GROUP *)gp->gh;
746
int offset;
747
748
adios_errno = 0;
749
if (!gp) {
750
adios_error (err_invalid_group_struct, "Null pointer passed as group");
751
return -1;
752
}
753
if (!varname) {
754
adios_error (err_invalid_varname, "Null pointer passed as variable name!");
755
return -1;
756
}
757
758
/* Search in gp->fh->gvar_h->var_namelist instead of gp->var_namelist, because that
759
one comes back from the user. One idiot (who writes this comment) sorted the
760
list in place after gopen and before inq_var.
761
*/
762
offset = gh->vars_offset;
763
764
if (varname[0] == '/')
765
vstartpos = 1;
766
for (varid=0; varid<(gp->vars_count);varid++) {
767
fstartpos = 0;
768
/* if (gp->var_namelist[varid][0] == '/') */
769
fstartpos = 0;
770
if (gh->fh->gvar_h->var_namelist[varid+offset][0] == '/')
771
fstartpos = 1;
772
/*if (!strcmp(gp->var_namelist[varid]+fstartpos, varname+vstartpos))*/
773
if (!strcmp(gh->fh->gvar_h->var_namelist[varid+offset]+fstartpos, varname+vstartpos))
774
break;
775
}
776
if (varid >= gp->vars_count) {
777
adios_error (err_invalid_varname, "Invalid variable name %s", varname);
778
return -1;
779
}
780
return varid;
781
}
782
783
// NCSU - For custom memory allocation
784
#define CALLOC(var, num, sz, comment)\
785
{\
786
var = calloc (num, sz); \
787
if (!var) {\
788
adios_error_at_line (err_no_memory, __FILE__, __LINE__, "Could not allocate memory for ", comment, " in common_read_get_characteristics"); \
789
return; \
790
}\
791
}
792
793
#define MALLOC(var,sz,comment)\
794
{\
795
var = malloc (sz); \
796
if (!var) {\
797
adios_error_at_line (err_no_memory, __FILE__, __LINE__, "Could not allocate memory for ", comment, " in common_read_get_characteristics"); \
798
return; \
799
}\
800
}\
801
802
// NCSU - Reading the statistics
803
/** Get value and statistics, allocate space for them too */
804
static void adios_read_bp_get_characteristics (struct adios_index_var_struct_v1 * var_root, ADIOS_VARINFO *vi)
805
{
806
int i, j, c, count = 1;
807
int size, sum_size, sum_type;
808
809
vi->value = NULL;
810
811
vi->gmin = vi->gmax = NULL;
812
vi->gavg = NULL;
813
vi->mins = vi->maxs = NULL;
814
vi->avgs = NULL;
815
vi->gstd_dev = NULL;
816
vi->std_devs = NULL;
817
vi->hist = NULL;
818
819
// set value for scalars
820
if (var_root->characteristics [0].value) {
821
size = bp_get_type_size(var_root->type, var_root->characteristics [0].value);
822
vi->value = (void *) malloc (size);
823
824
if (vi->value)
825
memcpy(vi->value, var_root->characteristics [0].value, size);
826
else {
827
adios_error_at_line (err_no_memory, __FILE__, __LINE__, "Could not allocate memory for value in common_read_get_characteristics");
828
return;
829
}
830
} else {
831
vi->value = NULL;
832
}
833
834
int npgs = var_root->characteristics_count, timestep, ntimes = -1;
835
uint64_t gcnt = 0, * cnts;
836
837
double *gsum = NULL, *gsum_square = NULL;
838
double **sums = NULL, **sum_squares = NULL;
839
840
int16_t map[32];
841
memset (map, -1, sizeof(map));
842
843
// Bitmap shows which statistical information has been calculated
844
i = j = 0;
845
while (var_root->characteristics[0].bitmap >> j)
846
{
847
if ((var_root->characteristics[0].bitmap >> j) & 1)
848
map [j] = i ++;
849
j ++;
850
}
851
852
if(vi->timedim >= 0)
853
{
854
ntimes = vi->dims[0];
855
856
if (map[adios_statistic_min] != -1)
857
{
858
MALLOC(vi->mins, ntimes * sizeof(void *), "minimum per timestep");
859
for (i = 0; i < ntimes; i++)
860
vi->mins[i] = 0;
861
}
862
863
if (map[adios_statistic_max] != -1)
864
{
865
MALLOC(vi->maxs, ntimes * sizeof(void *), "maximum per timestep");
866
for (i = 0; i < ntimes; i++)
867
vi->maxs[i] = 0;
868
}
869
870
if (map[adios_statistic_sum] != -1)
871
{
872
MALLOC(sums, ntimes * sizeof(double *), "summation per timestep");
873
MALLOC(vi->avgs, ntimes * sizeof(double *), "average per timestep");
874
875
for (i = 0; i < ntimes; i++)
876
sums[i] = vi->avgs[i] = 0;
877
878
CALLOC(cnts, ntimes, sizeof(uint64_t), "count of elements per timestep");
879
}
880
881
if (map[adios_statistic_sum_square] != -1)
882
{
883
MALLOC(sum_squares, ntimes * sizeof(double *), "summation per timestep");
884
MALLOC(vi->std_devs, ntimes * sizeof(double *), "standard deviation per timestep");
885
886
for (i = 0; i < ntimes; i++)
887
vi->std_devs[i] = sum_squares[i] = 0;
888
}
889
}
890
891
if (map[adios_statistic_hist] != -1 && (var_root->characteristics[0].stats[0][map[adios_statistic_hist]].data))
892
{
893
struct adios_index_characteristics_stat_struct * stats = var_root->characteristics[0].stats[0];
894
struct adios_index_characteristics_hist_struct * hist = stats[map[adios_statistic_hist]].data;
895
int num_breaks = hist->num_breaks;
896
897
MALLOC(vi->hist, sizeof(struct ADIOS_HIST), "histogram");
898
MALLOC(vi->hist->breaks, num_breaks * sizeof(double), "break points of histogram");
899
MALLOC(vi->hist->gfrequencies, (num_breaks + 1) * sizeof(uint32_t), "global frequencies of histogram");
900
901
vi->hist->num_breaks = hist->num_breaks;
902
vi->hist->min = hist->min;
903
vi->hist->max = hist->max;
904
905
memcpy(vi->hist->breaks, hist->breaks, num_breaks * sizeof(double));
906
CALLOC(vi->hist->gfrequencies, (num_breaks + 1), bp_get_type_size(adios_unsigned_integer, ""), "global frequency");
907
908
if (ntimes > 0)
909
{
910
MALLOC(vi->hist->frequenciess, (ntimes * sizeof(int32_t *)), "frequencies for timesteps");
911
for(i = 0; i < ntimes; i++)
912
CALLOC(vi->hist->frequenciess[i], (num_breaks + 1), bp_get_type_size(adios_unsigned_integer, ""), "frequency at timestep");
913
}
914
}
915
916
size = bp_get_type_size (var_root->type, "");
917
sum_size = bp_get_type_size (adios_double, "");
918
919
if (var_root->type == adios_complex || var_root->type == adios_double_complex)
920
{
921
int type;
922
count = 3;
923
timestep = 0;
924
925
if (var_root->type == adios_complex)
926
type = adios_double;
927
else
928
type = adios_long_double;
929
930
// Only a double precision returned for all complex values
931
size = bp_get_type_size (adios_double, "");
932
933
for (i=0; i<var_root->characteristics_count; i++)
934
{
935
if (ntimes > 0)
936
timestep = var_root->characteristics[i].time_index - 1;
937
938
if (!var_root->characteristics[i].stats)
939
continue;
940
941
struct adios_index_characteristics_stat_struct ** stats = var_root->characteristics[i].stats;
942
943
if ((map[adios_statistic_finite] != -1) && (* ((uint8_t *) stats[0][map[adios_statistic_finite]].data) == 0))
944
continue;
945
946
if (map[adios_statistic_min] != -1 && stats[0][map[adios_statistic_min]].data)
947
{
948
double data[3];
949
for (c = 0; c < count; c ++)
950
data[c] = bp_value_to_double(type, stats[c][map[adios_statistic_min]].data);
951
952
if(!vi->gmin) {
953
MALLOC (vi->gmin, count * size, "global minimum")
954
for (c = 0; c < count; c ++)
955
((double * ) vi->gmin)[c] = data[c];
956
957
} else {
958
for (c = 0; c < count; c ++)
959
if (data[c] < ((double *) vi->gmin)[c])
960
((double * ) vi->gmin)[c] = data[c];
961
}
962
963
if (ntimes > 0) {
964
if(!vi->mins[timestep]) {
965
MALLOC (vi->mins[timestep], count * size, "minimum per timestep")
966
for (c = 0; c < count; c ++)
967
((double **) vi->mins)[timestep][c] = data[c];
968
969
} else {
970
for (c = 0; c < count; c ++)
971
if (data[c] < ((double **) vi->mins)[timestep][c])
972
((double **) vi->mins)[timestep][c] = data[c];
973
}
974
}
975
}
976
977
if (map[adios_statistic_max] != -1 && stats[0][map[adios_statistic_max]].data)
978
{
979
double data[3];
980
for (c = 0; c < count; c ++)
981
data[c] = bp_value_to_double(type, stats[c][map[adios_statistic_max]].data);
982
983
if(!vi->gmax) {
984
MALLOC (vi->gmax, count * size, "global minimum")
985
for (c = 0; c < count; c ++)
986
((double * ) vi->gmax)[c] = data[c];
987
988
} else {
989
for (c = 0; c < count; c ++)
990
if (data[c] > ((double *) vi->gmax)[c])
991
((double * ) vi->gmax)[c] = data[c];
992
}
993
994
if (ntimes > 0) {
995
if(!vi->maxs[timestep]) {
996
MALLOC (vi->maxs[timestep], count * size, "minimum per timestep")
997
for (c = 0; c < count; c ++)
998
((double **) vi->maxs)[timestep][c] = data[c];
999
1000
} else {
1001
for (c = 0; c < count; c ++)
1002
if (data[c] > ((double **) vi->maxs)[timestep][c])
1003
((double **) vi->maxs)[timestep][c] = data[c];
1004
}
1005
}
1006
}
1007
1008
if (map[adios_statistic_sum] != -1 && stats[0][map[adios_statistic_sum]].data)
1009
{
1010
double data[3];
1011
for (c = 0; c < count; c ++)
1012
data[c] = bp_value_to_double(type, stats[c][map[adios_statistic_sum]].data);
1013
1014
if(!gsum) {
1015
MALLOC(gsum, count * sum_size, "global summation")
1016
for (c = 0; c < count; c ++)
1017
gsum[c] = data[c];
1018
1019
} else {
1020
for (c = 0; c < count; c ++)
1021
gsum[c] = gsum[c] + data[c];
1022
}
1023
1024
if (ntimes > 0) {
1025
if(!sums[timestep]) {
1026
MALLOC(sums[timestep], count * sum_size, "summation per timestep")
1027
for (c = 0; c < count; c ++)
1028
sums[timestep][c] = data[c];
1029
1030
} else {
1031
for (c = 0; c < count; c ++)
1032
sums[timestep][c] = sums[timestep][c] + data[c];
1033
}
1034
}
1035
}
1036
1037
if (map[adios_statistic_sum_square] != -1 && stats[0][map[adios_statistic_sum_square]].data)
1038
{
1039
double data[3];
1040
for (c = 0; c < count; c ++)
1041
data[c] = bp_value_to_double(type, stats[c][map[adios_statistic_sum_square]].data);
1042
1043
if(!gsum_square) {
1044
MALLOC(gsum_square, count * sum_size, "global summation of squares")
1045
for (c = 0; c < count; c ++)
1046
gsum_square[c] = data[c];
1047
1048
} else {
1049
for (c = 0; c < count; c ++)
1050
gsum_square[c] = gsum_square[c] + data[c];
1051
}
1052
1053
if (ntimes > 0) {
1054
if(!sum_squares[timestep]) {
1055
MALLOC(sum_squares[timestep], count * sum_size, "summation of square per timestep")
1056
for (c = 0; c < count; c ++)
1057
sum_squares[timestep][c] = data[c];
1058
1059
} else {
1060
for (c = 0; c < count; c ++)
1061
sum_squares[timestep][c] = sum_squares[timestep][c] + data[c];
1062
}
1063
}
1064
}
1065
1066
if (map[adios_statistic_cnt] != -1 && stats[0][map[adios_statistic_cnt]].data)
1067
{
1068
if (ntimes > 0)
1069
cnts[timestep] += * ((uint32_t *) stats[0][map[adios_statistic_cnt]].data);
1070
gcnt += * (uint32_t *) stats[0][map[adios_statistic_cnt]].data;
1071
}
1072
}
1073
1074
if(ntimes > 0 && vi->gmin && (map[adios_statistic_sum] != -1) && (map[adios_statistic_sum_square] != -1)) {
1075
// min, max, summation exists only for arrays
1076
// Calculate average / timestep
1077
1078
for(timestep = 0; timestep < ntimes; timestep ++) {
1079
MALLOC(vi->avgs[timestep], count * sum_size, "average per timestep")
1080
for (c = 0; c < count; c ++)
1081
vi->avgs[timestep][c] = sums[timestep][c] / cnts[timestep];
1082
1083
MALLOC(vi->std_devs[timestep], count * sum_size, "standard deviation per timestep")
1084
for (c = 0; c < count; c ++)
1085
vi->std_devs[timestep][c] = sqrt((sum_squares[timestep][c] / cnts[timestep]) - (vi->avgs[timestep][c] * vi->avgs[timestep][c]));
1086
1087
free (sums[timestep]);
1088
free (sum_squares[timestep]);
1089
}
1090
}
1091
1092
// Calculate global average
1093
if(vi->gmin && gsum && (map[adios_statistic_sum] != -1) && (map[adios_statistic_sum_square] != -1)) {
1094
MALLOC(vi->gavg, count * sum_size, "global average")
1095
1096
if(gcnt > 0)
1097
for (c = 0; c < count; c ++)
1098
vi->gavg[c] = gsum[c] / gcnt;
1099
else
1100
for (c = 0; c < count; c ++)
1101
vi->gavg[c] = gsum[c];
1102
1103
MALLOC(vi->gstd_dev, count * sum_size, "global average")
1104
if(vi->gavg && gcnt > 0)
1105
for (c = 0; c < count; c ++)
1106
vi->gstd_dev[c] = sqrt(gsum_square[c] / gcnt - (vi->gavg[c] * vi->gavg[c]));
1107
else
1108
for (c = 0; c < count; c ++)
1109
vi->gstd_dev[c] = 0;
1110
}
1111
}
1112
else
1113
{
1114
timestep = 0;
1115
for (i=0; i<var_root->characteristics_count; i++)
1116
{
1117
if (ntimes > 0)
1118
timestep = var_root->characteristics[i].time_index - 1;
1119
//timestep = i / (npgs / ntimes);
1120
1121
if (!var_root->characteristics[i].stats)
1122
continue;
1123
1124
struct adios_index_characteristics_stat_struct * stats = var_root->characteristics[i].stats[0];
1125
struct adios_index_characteristics_hist_struct * hist = stats[map[adios_statistic_hist]].data;
1126
1127
if (map[adios_statistic_finite] != -1 && (* ((uint8_t *) stats[map[adios_statistic_finite]].data) == 0))
1128
continue;
1129
1130
if (map[adios_statistic_min] != -1 && stats[map[adios_statistic_min]].data)
1131
{
1132
if(!vi->gmin) {
1133
MALLOC (vi->gmin, size, "global minimum")
1134
memcpy(vi->gmin, stats[map[adios_statistic_min]].data, size);
1135
1136
} else if (adios_lt(var_root->type, stats[map[adios_statistic_min]].data, vi->gmin)){
1137
memcpy(vi->gmin, stats[map[adios_statistic_min]].data, size);
1138
}
1139
1140
if (ntimes > 0) {
1141
if(!vi->mins[timestep]) {
1142
MALLOC (vi->mins[timestep], size, "minimum per timestep")
1143
memcpy(vi->mins[timestep], stats[map[adios_statistic_min]].data, size);
1144
1145
} else if (adios_lt(var_root->type, stats[map[adios_statistic_min]].data, vi->mins[timestep])) {
1146
memcpy(vi->mins[timestep], stats[map[adios_statistic_min]].data, size);
1147
}
1148
}
1149
}
1150
1151
if (map[adios_statistic_max] != -1 && stats[map[adios_statistic_max]].data)
1152
{
1153
if(!vi->gmax) {
1154
MALLOC (vi->gmax, size, "global maximum")
1155
memcpy(vi->gmax, stats[map[adios_statistic_max]].data, size);
1156
1157
} else if (adios_lt(var_root->type, vi->gmax, stats[map[adios_statistic_max]].data))
1158
memcpy(vi->gmax, stats[map[adios_statistic_max]].data, size);
1159
1160
if (ntimes > 0) {
1161
if(!vi->maxs[timestep]) {
1162
MALLOC (vi->maxs[timestep], size, "maximum per timestep")
1163
memcpy(vi->maxs[timestep], stats[map[adios_statistic_max]].data, size);
1164
1165
} else if (adios_lt(var_root->type, vi->maxs[timestep], stats[map[adios_statistic_max]].data)) {
1166
memcpy(vi->maxs[timestep], stats[map[adios_statistic_max]].data, size);
1167
}
1168
}
1169
}
1170
1171
if (map[adios_statistic_sum] != -1 && stats[map[adios_statistic_sum]].data)
1172
{
1173
if(!gsum) {
1174
MALLOC(gsum, sum_size, "global summation")
1175
memcpy(gsum, stats[map[adios_statistic_sum]].data, sum_size);
1176
1177
} else {
1178
*gsum = *gsum + * ((double *) stats[map[adios_statistic_sum]].data);
1179
}
1180
1181
if (ntimes > 0) {
1182
if(!sums[timestep]) {
1183
MALLOC(sums[timestep], sum_size, "summation per timestep")
1184
memcpy(sums[timestep], stats[map[adios_statistic_sum]].data, sum_size);
1185
1186
} else {
1187
*sums[timestep] = *sums[timestep] + * ((double *) stats[map[adios_statistic_sum]].data);
1188
}
1189
}
1190
}
1191
1192
if (map[adios_statistic_sum_square] != -1 && stats[map[adios_statistic_sum_square]].data)
1193
{
1194
if(!gsum_square) {
1195
MALLOC(gsum_square, sum_size, "global summation of squares")
1196
memcpy(gsum_square, stats[map[adios_statistic_sum_square]].data, sum_size);
1197
1198
} else {
1199
*gsum_square = *gsum_square + * ((double *) stats[map[adios_statistic_sum_square]].data);
1200
}
1201
1202
if (ntimes > 0) {
1203
if(!sum_squares[timestep]) {
1204
MALLOC(sum_squares[timestep], sum_size, "summation of square per timestep")
1205
memcpy(sum_squares[timestep], stats[map[adios_statistic_sum_square]].data, sum_size);
1206
1207
} else {
1208
*sum_squares[timestep] = *sum_squares[timestep] + * ((double *) stats[map[adios_statistic_sum_square]].data);
1209
}
1210
}
1211
}
1212
1213
if(map[adios_statistic_hist] != -1 && stats[map[adios_statistic_hist]].data)
1214
{
1215
for(j = 0; j <= vi->hist->num_breaks; j++)
1216
{
1217
uint32_t freq = hist->frequencies[j];
1218
vi->hist->gfrequencies[j] += freq;
1219
if (ntimes > 0)
1220
vi->hist->frequenciess[timestep][j] += freq;
1221
}
1222
}
1223
1224
if (map[adios_statistic_cnt] != -1 && stats[map[adios_statistic_cnt]].data)
1225
{
1226
if (ntimes > 0)
1227
cnts[timestep] += * (uint32_t *) stats[map[adios_statistic_cnt]].data;
1228
gcnt += * (uint32_t *) stats[map[adios_statistic_cnt]].data;
1229
}
1230
}
1231
1232
if(ntimes > 0 && vi->gmin && (map[adios_statistic_sum] != -1) && (map[adios_statistic_sum_square] != -1)) {
1233
// min, max, summation exists only for arrays
1234
// Calculate average / timestep
1235
1236
for(timestep = 0; timestep < ntimes; timestep ++) {
1237
MALLOC(vi->avgs[timestep], sum_size, "average per timestep")
1238
*(vi->avgs[timestep]) = *(sums[timestep]) / cnts[timestep];
1239
1240
MALLOC(vi->std_devs[timestep], sum_size, "standard deviation per timestep")
1241
*(vi->std_devs[timestep]) = sqrt(*(sum_squares[timestep]) / cnts[timestep] - ((*(vi->avgs[timestep]) * (*(vi->avgs[timestep])))));
1242
1243
free (sums[timestep]);
1244
free (sum_squares[timestep]);
1245
}
1246
}
1247
1248
// Calculate global average
1249
if(vi->gmin && gsum && (map[adios_statistic_sum] != -1) && (map[adios_statistic_sum_square] != -1)) {
1250
MALLOC(vi->gavg, sum_size, "global average")
1251
if(gcnt > 0)
1252
*vi->gavg = *gsum / gcnt;
1253
else
1254
vi->gavg = gsum;
1255
1256
MALLOC(vi->gstd_dev, sum_size, "global average")
1257
if(vi->gavg && gcnt > 0)
1258
*vi->gstd_dev = sqrt(*gsum_square / gcnt - ((*(vi->gavg)) * (*(vi->gavg))));
1259
else
1260
*vi->gstd_dev = 0;
1261
}
1262
}
1263
1264
if (!vi->value && vi->gmin) {
1265
vi->value = vi->gmin; // arrays have no value but we assign here the minimum
1266
}
1267
1268
if(!vi->gmin) {
1269
vi->gmin = vi->value; // scalars have value but not min
1270
}
1271
if(!vi->gmax) {
1272
vi->gmax = vi->value; // scalars have value but not max
1273
}
1274
1275
if (sums && gsum) {
1276
free (sums);
1277
free (gsum);
1278
}
1279
1280
if (sum_squares && gsum_square) {
1281
free (sum_squares);
1282
free (gsum_square);
1283
}
1284
}
1285
1286
/* get local and global dimensions and offsets from a variable characteristics
1287
return: 1 = it is a global array, 0 = local array
1288
*/
1289
static int adios_read_bp_get_dimensioncharacteristics(struct adios_index_characteristic_struct_v1 *ch,
1290
uint64_t *ldims, uint64_t *gdims, uint64_t *offsets)
1291
{
1292
int is_global=0; // global array or just an array written by one process?
1293
int ndim = ch->dims.count;
1294
int k;
1295
1296
for (k=0; k < ndim; k++) {
1297
ldims[k] = ch->dims.dims[k*3];
1298
gdims[k] = ch->dims.dims[k*3+1];
1299
offsets[k] = ch->dims.dims[k*3+2];
1300
is_global = is_global || gdims[k];
1301
}
1302
return is_global;
1303
}
1304
1305
static void adios_read_bp_get_dimensions (struct adios_index_var_struct_v1 *var_root, int ntsteps, int file_is_fortran,
1306
int *ndim, uint64_t **dims, int *timedim)
1307
{
1308
int i, k;
1309
int is_global; // global array or just an array written by one process?
1310
uint64_t ldims[32];
1311
uint64_t gdims[32];
1312
uint64_t offsets[32];
1313
1314
/* Get dimension information */
1315
*ndim = var_root->characteristics [0].dims.count;
1316
*dims = NULL;
1317
*timedim = -1;
1318
if (*ndim == 0) {
1319
/* done with this scalar variable */
1320
return ;
1321
}
1322
1323
*dims = (uint64_t *) malloc (sizeof(uint64_t) * (*ndim));
1324
memset(*dims,0,sizeof(uint64_t)*(*ndim));
1325
1326
is_global = adios_read_bp_get_dimensioncharacteristics( &(var_root->characteristics[0]),
1327
ldims, gdims, offsets);
1328
1329
if (!is_global) {
1330
/* local array */
1331
for (i=0; i < *ndim; i++) {
1332
/* size of time dimension is always registered as 1 for an array */
1333
if (ldims[i] == 1 && var_root->characteristics_count > 1) {
1334
*timedim = i;
1335
(*dims)[i] = ntsteps;
1336
} else {
1337
(*dims)[i] = ldims[i];
1338
}
1339
gdims[i] = ldims[i];
1340
}
1341
}
1342
else {
1343
/* global array:
1344
time dimension: ldims=1, gdims=0
1345
in C array, it can only be the first dim
1346
in Fortran array, it can only be the last dim
1347
(always the slowest changing dim)
1348
*/
1349
if (gdims[*ndim-1] == 0)
1350
{
1351
if (!file_is_fortran) {
1352
/* first dimension is the time (C array)
1353
* ldims[0] = 1 but gdims does not contain time info and
1354
* gdims[0] is 1st data dimension and
1355
* gdims is shorter by one value than ldims in case of C.
1356
* Therefore, gdims[*ndim-1] = 0 if there is a time dimension.
1357
*/
1358
*timedim = 0;
1359
// error check
1360
if (*ndim > 1 && ldims[0] != 1) {
1361
fprintf(stderr,"ADIOS Error: this is a BP file with C ordering but we didn't find"
1362
"an array to have time dimension in the first dimension. l:g:o = (");
1363
for (i=0; i < *ndim; i++) {
1364
fprintf(stderr,"%llu:%llu:%llu%s", ldims[i], gdims[i], offsets[i], (i<*ndim-1 ? ", " : "") );
1365
}
1366
fprintf(stderr, ")\n");
1367
}
1368
(*dims)[0] = ntsteps;
1369
for (i=1; i < *ndim; i++)
1370
(*dims)[i]=gdims[i-1];
1371
} else {
1372
// last dimension is the time (Fortran array)
1373
*timedim = *ndim - 1;
1374
1375
if (*ndim > 1 && ldims[*ndim-1] != 1) {
1376
fprintf(stderr,"ADIOS Error: this is a BP file with Fortran array ordering but we didn't find"
1377
"an array to have time dimension in the last dimension. l:g:o = (");
1378
for (i=0; i < *ndim; i++) {
1379
fprintf(stderr,"%llu:%llu:%llu%s", ldims[i], gdims[i], offsets[i], (i<*ndim-1 ? ", " : "") );
1380
}
1381
fprintf(stderr, ")\n");
1382
}
1383
for (i=0; i < *ndim-1; i++)
1384
(*dims)[i]=gdims[i];
1385
(*dims)[*timedim] = ntsteps;
1386
}
1387
} else {
1388
// no time dimenstion
1389
for (i=0; i < *ndim; i++)
1390
(*dims)[i]=gdims[i];
1391
}
1392
}
1393
}
1394
1395
ADIOS_VARINFO * adios_read_bp_inq_var (ADIOS_GROUP *gp, const char * varname)
1396
{
1397
int varid = adios_read_bp_find_var(gp, varname);
1398
if (varid < 0)
1399
return NULL;
1400
return adios_read_bp_inq_var_byid(gp, varid);
1401
}
1402
1403
ADIOS_VARINFO * adios_read_bp_inq_var_byid (ADIOS_GROUP *gp, int varid)
1404
{
1405
struct BP_GROUP * gh;
1406
struct BP_FILE * fh;
1407
ADIOS_VARINFO * vi;
1408
int file_is_fortran;
1409
struct adios_index_var_struct_v1 * var_root;
1410
int i,k;
1411
1412
adios_errno = 0;
1413
if (!gp) {
1414
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_inq_var()");
1415
return NULL;
1416
}
1417
gh = (struct BP_GROUP *) gp->gh;
1418
if (!gh) {
1419
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
1420
return NULL;
1421
}
1422
fh = gh->fh;
1423
if (!fh) {
1424
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
1425
return NULL;
1426
}
1427
if (varid < 0 || varid >= gh->vars_count) {
1428
adios_error (err_invalid_varid, "Invalid variable id %d (allowed 0..%d)", varid, gh->vars_count);
1429
return NULL;
1430
}
1431
vi = (ADIOS_VARINFO *) malloc(sizeof(ADIOS_VARINFO));
1432
if (!vi) {
1433
adios_error ( err_no_memory, "Could not allocate memory for variable info");
1434
return NULL;
1435
}
1436
1437
1438
file_is_fortran = (fh->pgs_root->adios_host_language_fortran == adios_flag_yes);
1439
1440
var_root = gh->vars_root; /* first variable of this group. Need to traverse the list */
1441
for (i=0; i<varid && var_root; i++) {
1442
var_root = var_root->next;
1443
}
1444
1445
if (i!=varid) {
1446
adios_error (err_corrupted_variable, "Variable id=%d is valid but was not found in internal data structures!",varid);
1447
return NULL;
1448
}
1449
1450
1451
vi->varid = varid;
1452
1453
vi->type = var_root->type;
1454
if (!var_root->characteristics_count) {
1455
adios_error (err_corrupted_variable, "Variable %s does not have information on dimensions",
1456
gp->var_namelist[varid]);
1457
free(vi);
1458
return NULL;
1459
}
1460
vi->characteristics_count = var_root->characteristics_count;
1461
1462
/* Get value or min/max */
1463
1464
adios_read_bp_get_dimensions (var_root, fh->tidx_stop - fh->tidx_start + 1, file_is_fortran,
1465
&(vi->ndim), &(vi->dims), &(vi->timedim));
1466
1467
if (file_is_fortran != futils_is_called_from_fortran()) {
1468
/* If this is a Fortran written file and this is called from C code/
1469
or this is a C written file and this is called from Fortran code ==>
1470
We need to reverse the order of the dimensions */
1471
swap_order(vi->ndim, vi->dims, &(vi->timedim));
1472
/*printf("File was written from %s and read now from %s, so we swap order of dimensions\n",
1473
(file_is_fortran ? "Fortran" : "C"), (futils_is_called_from_fortran() ? "Fortran" : "C"));*/
1474
}
1475
1476
adios_read_bp_get_characteristics (var_root, vi);
1477
1478
return vi;
1479
}
1480
1481
void adios_read_bp_free_varinfo (ADIOS_VARINFO *vp)
1482
{
1483
if (vp) {
1484
if (vp->dims) free(vp->dims);
1485
if (vp->value) free(vp->value);
1486
if (vp->gmin && vp->gmin != vp->value) free(vp->gmin);
1487
if (vp->gmax && vp->gmax != vp->value) free(vp->gmax);
1488
//if (vp->mins) free(vp->mins);
1489
//if (vp->maxs) free(vp->maxs);
1490
free(vp);
1491
}
1492
}
1493
1494
#define MPI_FILE_READ_OPS \
1495
bp_realloc_aligned(fh->b, slice_size); \
1496
fh->b->offset = 0; \
1497
\
1498
MPI_File_seek (fh->mpi_fh \
1499
,(MPI_Offset)slice_offset \
1500
,MPI_SEEK_SET \
1501
); \
1502
\
1503
MPI_File_read (fh->mpi_fh \
1504
,fh->b->buff \
1505
,slice_size \
1506
,MPI_BYTE \
1507
,&status \
1508
); \
1509
fh->b->offset = 0; \
1510
1511
//We also need to be able to read old .bp which doesn't have 'payload_offset'
1512
#define MPI_FILE_READ_OPS1 \
1513
MPI_File_seek (fh->mpi_fh \
1514
,(MPI_Offset) var_root->characteristics[start_idx + idx].offset \
1515
,MPI_SEEK_SET); \
1516
MPI_File_read (fh->mpi_fh, fh->b->buff, 8, MPI_BYTE, &status); \
1517
tmpcount= *((uint64_t*)fh->b->buff); \
1518
\
1519
bp_realloc_aligned(fh->b, tmpcount + 8); \
1520
fh->b->offset = 0; \
1521
\
1522
MPI_File_seek (fh->mpi_fh \
1523
,(MPI_Offset) (var_root->characteristics[start_idx + idx].offset) \
1524
,MPI_SEEK_SET); \
1525
MPI_File_read (fh->mpi_fh, fh->b->buff, tmpcount + 8, MPI_BYTE, &status); \
1526
fh->b->offset = 0; \
1527
adios_parse_var_data_header_v1 (fh->b, &var_header); \
1528
1529
1530
// To read subfiles
1531
#define MPI_FILE_READ_OPS2 \
1532
bp_realloc_aligned(fh->b, slice_size); \
1533
fh->b->offset = 0; \
1534
\
1535
MPI_File * sfh; \
1536
sfh = get_BP_file_handle (fh->sfh \
1537
,var_root->characteristics[start_idx + idx].file_index \
1538
); \
1539
if (!sfh) \
1540
{ \
1541
int err; \
1542
char * ch, * name_no_path, * name; \
1543
struct BP_file_handle * new_h = \
1544
(struct BP_file_handle *) malloc (sizeof (struct BP_file_handle)); \
1545
new_h->file_index = var_root->characteristics[start_idx + idx].file_index; \
1546
new_h->next = 0; \
1547
if (ch = strrchr (fh->fname, '/')) \
1548
{ \
1549
name_no_path = malloc (strlen (ch + 1) + 1); \
1550
strcpy (name_no_path, ch + 1); \
1551
} \
1552
else \
1553
{ \
1554
name_no_path = malloc (strlen (fh->fname) + 1); \
1555
strcpy (name_no_path, fh->fname); \
1556
} \
1557
\
1558
name = malloc (strlen (fh->fname) + 5 + strlen (name_no_path) + 1 + 10 + 1); \
1559
sprintf (name, "%s.dir/%s.%d", fh->fname, name_no_path, new_h->file_index); \
1560
\
1561
err = MPI_File_open (fh->comm \
1562
,name \
1563
,MPI_MODE_RDONLY \
1564
,(MPI_Info)MPI_INFO_NULL \
1565
,&new_h->fh \
1566
); \
1567
if (err != MPI_SUCCESS) \
1568
{ \
1569
fprintf (stderr, "can not open file %S\n", name); \
1570
return -1; \
1571
} \
1572
\
1573
add_BP_file_handle (&fh->sfh \
1574
,new_h \
1575
); \
1576
sfh = &new_h->fh; \
1577
\
1578
free (name_no_path); \
1579
free (name); \
1580
} \
1581
\
1582
MPI_File_seek (*sfh \
1583
,(MPI_Offset)slice_offset \
1584
,MPI_SEEK_SET \
1585
); \
1586
MPI_File_read (*sfh \
1587
,fh->b->buff \
1588
,slice_size \
1589
,MPI_BYTE \
1590
,&status \
1591
); \
1592
fh->b->offset = 0; \
1593
1594
1595
int64_t adios_read_bp_read_var (ADIOS_GROUP * gp, const char * varname,
1596
const uint64_t * start, const uint64_t * count,
1597
void * data)
1598
{
1599
struct BP_GROUP * gh;
1600
struct BP_FILE * fh;
1601
int varid, has_subfile;
1602
1603
adios_errno = 0;
1604
if (!gp) {
1605
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_read_var()");
1606
return -adios_errno;
1607
}
1608
1609
gh = (struct BP_GROUP *) gp->gh;
1610
if (!gh) {
1611
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
1612
return -adios_errno;
1613
}
1614
1615
fh = gh->fh;
1616
if (!fh) {
1617
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
1618
return -adios_errno;
1619
}
1620
1621
varid = adios_read_bp_find_var(gp, varname);
1622
if (varid < 0 || varid >= gh->vars_count) {
1623
adios_error (err_invalid_varid, "Invalid variable id %d (allowed 0..%d)", varid, gh->vars_count);
1624
return -adios_errno;
1625
}
1626
1627
return adios_read_bp_read_var_byid(gp, varid, start, count, data);
1628
}
1629
1630
/***********************************************
1631
* This routine is to read in data in a 'local *
1632
* array fashion (as opposed to global array) *
1633
* Q. Liu, 11/2010 *
1634
***********************************************/
1635
int64_t adios_read_bp_read_local_var (ADIOS_GROUP * gp, const char * varname,
1636
int vidx, const uint64_t * start,
1637
const uint64_t * count, void * data)
1638
{
1639
struct BP_GROUP * gh;
1640
struct BP_FILE * fh;
1641
struct adios_index_var_struct_v1 * var_root;
1642
struct adios_var_header_struct_v1 var_header;
1643
struct adios_var_payload_struct_v1 var_payload;
1644
int i,j,k, t, varid, start_idx, idx;
1645
int ndim, ndim_notime, has_subfile, file_is_fortran;
1646
uint64_t size, * dims;
1647
uint64_t ldims[32], gdims[32], offsets[32];
1648
uint64_t datasize, nloop, dset_stride,var_stride, total_size=0, items_read;
1649
uint64_t count_notime[32], start_notime[32];
1650
int timedim = -1, temp_timedim, is_global = 0, size_of_type;
1651
uint64_t slice_offset, slice_size, tmpcount = 0;
1652
uint64_t datatimeoffset = 0; // offset in data to write a given timestep
1653
MPI_Status status;
1654
1655
adios_errno = 0;
1656
if (!gp)
1657
{
1658
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_read_var()");
1659
return -adios_errno;
1660
}
1661
1662
gh = (struct BP_GROUP *) gp->gh;
1663
if (!gh)
1664
{
1665
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
1666
return -adios_errno;
1667
}
1668
1669
fh = gh->fh;
1670
if (!fh)
1671
{
1672
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
1673
return -adios_errno;
1674
}
1675
1676
varid = adios_read_bp_find_var(gp, varname);
1677
if (varid < 0 || varid >= gh->vars_count)
1678
{
1679
adios_error (err_invalid_varid, "Invalid variable id %d (allowed 0..%d)", varid, gh->vars_count);
1680
return -adios_errno;
1681
}
1682
1683
/* Check if file is written out by Fortran or C */
1684
file_is_fortran = (fh->pgs_root->adios_host_language_fortran == adios_flag_yes);
1685
1686
/* Check whether we need to handle subfiles */
1687
has_subfile = fh->mfooter.version & ADIOS_VERSION_HAVE_SUBFILE;
1688
1689
var_root = gh->vars_root; /* first variable of this group. Need to traverse the list */
1690
for (i = 0; i< varid && var_root; i++)
1691
{
1692
var_root = var_root->next;
1693
}
1694
1695
if (i != varid)
1696
{
1697
adios_error (err_corrupted_variable,
1698
"Variable id=%d is valid but was not found in internal data structures!",
1699
varid);
1700
return -adios_errno;
1701
}
1702
1703
if (vidx < 0 || vidx >= var_root->characteristics_count)
1704
{
1705
adios_error (err_out_of_bound, "idx=%d is out of bound", vidx);
1706
}
1707
1708
ndim = var_root->characteristics [vidx].dims.count;
1709
1710
/* count_notime/start_notime are working copies of count/start */
1711
for (i = 0; i < ndim; i++)
1712
{
1713
count_notime[i] = count[i];
1714
start_notime[i] = start[i];
1715
}
1716
1717
ndim_notime = ndim;
1718
1719
/* Fortran reader was reported of Fortran dimension order so it gives counts and starts in that order.
1720
We need to swap them here to read correctly in C order */
1721
if (futils_is_called_from_fortran())
1722
{
1723
timedim = -1;
1724
swap_order(ndim_notime, count_notime, &timedim);
1725
swap_order(ndim_notime, start_notime, &timedim);
1726
}
1727
1728
/* items_read = how many data elements are we going to read in total */
1729
items_read = 1;
1730
for (i = 0; i < ndim_notime; i++)
1731
items_read *= count_notime[i];
1732
1733
size_of_type = bp_get_type_size (var_root->type, var_root->characteristics [vidx].value);
1734
1735
/* READ A SCALAR VARIABLE */
1736
if (ndim_notime == 0)
1737
{
1738
slice_size = size_of_type;
1739
start_idx = 0; // OPS macros below need it
1740
idx = vidx; // OPS macros below need it
1741
1742
if (var_root->type == adios_string)
1743
{
1744
// strings are stored without \0 in file
1745
// size_of_type here includes \0 so decrease by one
1746
size_of_type--;
1747
}
1748
1749
/* Old BP files don't have payload_offset characteristic */
1750
if (var_root->characteristics[vidx].payload_offset > 0)
1751
{
1752
slice_offset = var_root->characteristics[vidx].payload_offset;
1753
1754
if (!has_subfile)
1755
{
1756
MPI_FILE_READ_OPS
1757
}
1758
else
1759
{
1760
MPI_FILE_READ_OPS2
1761
}
1762
}
1763
else
1764
{
1765
slice_offset = 0;
1766
MPI_FILE_READ_OPS1
1767
}
1768
1769
memcpy((char *)data + total_size, fh->b->buff + fh->b->offset, size_of_type);
1770
1771
if (fh->mfooter.change_endianness == adios_flag_yes)
1772
change_endianness((char *)data + total_size
1773
,size_of_type
1774
,var_root->type
1775
);
1776
1777
if (var_root->type == adios_string)
1778
{
1779
// add \0 to the end of string
1780
// size_of_type here is the length of string
1781
// FIXME: how would this work for strings written over time?
1782
((char*)data + total_size)[size_of_type] = '\0';
1783
}
1784
1785
total_size += size_of_type;
1786
1787
return total_size;
1788
} /* READ A SCALAR VARIABLE END HERE */
1789
1790
/* READ AN ARRAY VARIABLE */
1791
uint64_t write_offset = 0;
1792
int npg = 0;
1793
tmpcount = 0;
1794
int flag;
1795
datasize = 1;
1796
nloop = 1;
1797
var_stride = 1;
1798
dset_stride = 1;
1799
uint64_t payload_size = size_of_type;
1800
1801
/* To get ldims for the index vidx */
1802
adios_read_bp_get_dimensioncharacteristics( &(var_root->characteristics[vidx]),
1803
ldims, gdims, offsets);
1804
1805
/* Again, a Fortran written file has the dimensions in Fortran order we need to swap here */
1806
/* Only local dims are needed for reading local vars */
1807
if (file_is_fortran)
1808
{
1809
i=-1;
1810
swap_order(ndim, ldims, &(i));
1811
}
1812
1813
/*
1814
printf("ldims = "); for (j = 0; j < ndim; j++) printf("%d ",ldims[j]); printf("\n");
1815
printf("count_notime = "); for (j = 0; j < ndim_notime; j++) printf("%d ",count_notime[j]); printf("\n");
1816
printf("start_notime = "); for (j = 0; j < ndim_notime; j++) printf("%d ",start_notime[j]); printf("\n");
1817
*/
1818
1819
for (j = 0; j < ndim_notime; j++)
1820
{
1821
payload_size *= ldims [j];
1822
1823
if ( (start_notime[j] > ldims[j])
1824
|| (start_notime[j] + count_notime[j] > ldims[j]))
1825
{
1826
adios_error ( err_out_of_bound, "Error: Variable (id=%d) out of bound ("
1827
"the data in dimension %d to read is %llu elements from index %llu"
1828
" but the actual data is [0,%llu])",
1829
varid, j+1, count_notime[j], start_notime[j], ldims[j] - 1);
1830
return -adios_errno;
1831
}
1832
}
1833
1834
/* determined how many (fastest changing) dimensions can we read in in one read */
1835
int break_dim = ndim_notime - 1;
1836
while (break_dim > -1)
1837
{
1838
if (start_notime[break_dim] == 0 && ldims[break_dim] == count_notime[break_dim])
1839
{
1840
datasize *= ldims[break_dim];
1841
}
1842
else
1843
break;
1844
1845
break_dim--;
1846
}
1847
1848
slice_offset = 0;
1849
slice_size = 0;
1850
/* Note: MPI_FILE_READ_OPS - for reading single BP file.
1851
* MPI_FILE_READ_OPS2 - for reading those with subfiles.
1852
* MPI_FILE_READ_OPS1 - for reading old version of BP files
1853
* which don't contain "payload_offset"
1854
* Whenever to use OPS macro, start_idx and idx variable needs to be
1855
* properly set.
1856
*/
1857
1858
start_idx = 0;
1859
idx = vidx;
1860
1861
if (break_dim <= 0)
1862
{
1863
/* The slowest changing dimensions should not be read completely but
1864
we still need to read only one block */
1865
1866
uint64_t size_in_dset = count_notime[0];
1867
uint64_t offset_in_dset = start_notime[0];
1868
1869
slice_size = (break_dim == -1 ? datasize * size_of_type : size_in_dset * datasize * size_of_type);
1870
1871
if (var_root->characteristics[start_idx + idx].payload_offset > 0)
1872
{
1873
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
1874
+ offset_in_dset * datasize * size_of_type;
1875
1876
if (!has_subfile)
1877
{
1878
MPI_FILE_READ_OPS
1879
}
1880
else
1881
{
1882
MPI_FILE_READ_OPS2
1883
}
1884
}
1885
else
1886
{
1887
slice_offset = 0;
1888
MPI_FILE_READ_OPS1
1889
}
1890
1891
memcpy ((char *)data, fh->b->buff + fh->b->offset, slice_size);
1892
if (fh->mfooter.change_endianness == adios_flag_yes)
1893
{
1894
change_endianness((char *)data + write_offset, slice_size, var_root->type);
1895
}
1896
}
1897
else
1898
{
1899
uint64_t stride_offset = 0;
1900
uint64_t * size_in_dset, * offset_in_dset, * offset_in_var;
1901
uint64_t start_in_payload, end_in_payload, s;
1902
uint64_t var_offset;
1903
uint64_t dset_offset;
1904
1905
size_in_dset = (uint64_t *) malloc (8 * ndim_notime);
1906
offset_in_dset = (uint64_t *) malloc (8 * ndim_notime);
1907
offset_in_var = (uint64_t *) malloc (8 * ndim_notime);
1908
1909
if (size_in_dset == 0 || offset_in_dset == 0 || offset_in_var == 0)
1910
{
1911
adios_error (err_no_memory, "Malloc failed in %s at %d\n"
1912
, __FILE__, __LINE__
1913
);
1914
return -adios_errno;
1915
}
1916
1917
for (i = 0; i < ndim_notime ; i++)
1918
{
1919
size_in_dset[i] = count_notime[i];
1920
offset_in_dset[i] = start_notime[i];
1921
offset_in_var[i] = 0;
1922
}
1923
1924
datasize = 1;
1925
var_stride = 1;
1926
for (i = ndim_notime - 1; i >= break_dim; i--)
1927
{
1928
datasize *= size_in_dset[i];
1929
dset_stride *= ldims[i];
1930
var_stride *= count_notime[i];
1931
}
1932
1933
/* Calculate the size of the chunk we are trying to read in */
1934
start_in_payload = 0;
1935
end_in_payload = 0;
1936
s = 1;
1937
for (i = ndim_notime - 1; i >= 0; i--)
1938
{
1939
start_in_payload += s * offset_in_dset[i] * size_of_type;
1940
end_in_payload += s * (offset_in_dset[i] + size_in_dset[i] - 1) * size_of_type;
1941
s *= ldims[i];
1942
}
1943
slice_size = end_in_payload - start_in_payload + 1 * size_of_type;
1944
1945
if (var_root->characteristics[start_idx + idx].payload_offset > 0)
1946
{
1947
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
1948
+ start_in_payload;
1949
if (!has_subfile)
1950
{
1951
MPI_FILE_READ_OPS
1952
}
1953
else
1954
{
1955
MPI_FILE_READ_OPS2
1956
}
1957
1958
for ( i = 0; i < ndim_notime ; i++)
1959
{
1960
offset_in_dset[i] = 0;
1961
}
1962
}
1963
else
1964
{
1965
slice_offset = start_in_payload;
1966
MPI_FILE_READ_OPS1
1967
}
1968
1969
var_offset = 0;
1970
dset_offset = 0;
1971
for (i = 0; i < ndim_notime ; i++)
1972
{
1973
var_offset = offset_in_var[i] + var_offset * count_notime[i];
1974
dset_offset = offset_in_dset[i] + dset_offset * ldims[i];
1975
}
1976
1977
copy_data (data
1978
,fh->b->buff + fh->b->offset
1979
,0
1980
,break_dim
1981
,size_in_dset
1982
,ldims
1983
,count_notime
1984
,var_stride
1985
,dset_stride
1986
,var_offset
1987
,dset_offset
1988
,datasize
1989
,size_of_type
1990
);
1991
1992
free (size_in_dset);
1993
free (offset_in_dset);
1994
free (offset_in_var);
1995
}
1996
1997
total_size += items_read * size_of_type;
1998
1999
return total_size;
2000
}
2001
2002
// The purpose of keeping this function is to be able
2003
// to read in old BP files. Can be deleted later on.
2004
int64_t adios_read_bp_read_var_byid1 (ADIOS_GROUP * gp,
2005
int varid,
2006
const uint64_t * start,
2007
const uint64_t * count,
2008
void * data)
2009
{
2010
struct BP_GROUP * gh;
2011
struct BP_FILE * fh;
2012
int file_is_fortran;
2013
struct adios_index_var_struct_v1 * var_root;
2014
struct adios_var_header_struct_v1 var_header;
2015
struct adios_var_payload_struct_v1 var_payload;
2016
int i,j,k, idx, timestep;
2017
int start_time, stop_time;
2018
int pgoffset, pgcount, next_pgoffset,start_idx, stop_idx;
2019
int ndim, ndim_notime;
2020
uint64_t size;
2021
uint64_t *dims;
2022
uint64_t ldims[32];
2023
uint64_t gdims[32];
2024
uint64_t offsets[32];
2025
uint64_t datasize, nloop, dset_stride,var_stride, total_size=0, items_read;
2026
uint64_t count_notime[32], start_notime[32];
2027
MPI_Status status;
2028
int timedim = -1, temp_timedim, timedim_c;
2029
int rank;
2030
int is_global = 0;
2031
int size_of_type;
2032
uint64_t slice_offset;
2033
uint64_t slice_size;
2034
uint64_t tmpcount = 0;
2035
uint64_t datatimeoffset = 0; // offset in data to write a given timestep
2036
2037
adios_errno = 0;
2038
if (!gp) {
2039
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_read_var()");
2040
return -adios_errno;
2041
}
2042
gh = (struct BP_GROUP *) gp->gh;
2043
if (!gh) {
2044
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
2045
return -adios_errno;
2046
}
2047
fh = gh->fh;
2048
if (!fh) {
2049
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
2050
return -adios_errno;
2051
}
2052
if (varid < 0 || varid >= gh->vars_count) {
2053
adios_error (err_invalid_varid, "Invalid variable id %d (allowed 0..%d)", varid, gh->vars_count);
2054
return -adios_errno;
2055
}
2056
2057
file_is_fortran = (fh->pgs_root->adios_host_language_fortran == adios_flag_yes);
2058
2059
var_root = gh->vars_root; /* first variable of this group. Need to traverse the list */
2060
for (i=0; i<varid && var_root; i++) {
2061
var_root = var_root->next;
2062
}
2063
2064
if (i!=varid) {
2065
adios_error (err_corrupted_variable, "Variable id=%d is valid but was not found in internal data structures!",varid);
2066
return -adios_errno;
2067
}
2068
2069
/* Get dimensions and flip if caller != writer language */
2070
adios_read_bp_get_dimensions (var_root, fh->tidx_stop - fh->tidx_start + 1, file_is_fortran,
2071
&ndim, &dims, &timedim);
2072
2073
/* Here the cases in which .bp written from Fortran and C are considered separately.
2074
1) bp written from Fortran */
2075
if (file_is_fortran)
2076
{
2077
/* Get the timesteps we need to read */
2078
if (timedim > -1)
2079
{
2080
if (timedim != ndim - 1)
2081
{
2082
adios_error (err_no_data_at_timestep,"Variable (id=%d) has wrong time dimension index",
2083
varid);
2084
return -adios_errno;
2085
}
2086
if (futils_is_called_from_fortran())
2087
{
2088
start_time = start[timedim] + fh->tidx_start;
2089
stop_time = start_time + count[timedim] - 1;
2090
}
2091
else
2092
{
2093
start_time = start[0] + fh->tidx_start;
2094
stop_time = start_time + count[0] - 1;
2095
}
2096
}
2097
else
2098
{
2099
/* timeless variable, but we still need to handle the case that
2100
var is not written in the first few timesteps.
2101
This happens in Pixie3D. */
2102
for (i = 0; i < fh->mfooter.time_steps; i++)
2103
{
2104
pgoffset = fh->gvar_h->time_index[0][gh->group_id][i];
2105
if (i < fh->mfooter.time_steps - 1)
2106
next_pgoffset = fh->gvar_h->time_index[0][gh->group_id][i + 1];
2107
else
2108
next_pgoffset = -1;
2109
2110
if (fh->gvar_h->pg_offsets[pgoffset] < var_root->characteristics[0].offset
2111
&& (i == fh->mfooter.time_steps - 1
2112
||fh->gvar_h->pg_offsets[next_pgoffset] > var_root->characteristics[0].offset)
2113
)
2114
{
2115
start_time = fh->tidx_start + i;
2116
stop_time = start_time;
2117
break;
2118
}
2119
}
2120
}
2121
2122
/* flip dims and timedim to C order */
2123
swap_order(ndim, dims, &timedim);
2124
2125
/* Take out the time dimension from start[] and count[] */
2126
/* if we have time dimension */
2127
if (timedim > -1)
2128
{
2129
j = 0;
2130
if (futils_is_called_from_fortran())
2131
temp_timedim = ndim - 1;
2132
else
2133
temp_timedim = 0;
2134
2135
for (i = 0; i < temp_timedim; i++)
2136
{
2137
count_notime[j] = count[i];
2138
start_notime[j] = start[i];
2139
j++;
2140
}
2141
i++; // skip timedim
2142
for (; i < ndim; i++)
2143
{
2144
count_notime[j] = count[i];
2145
start_notime[j] = start[i];
2146
j++;
2147
}
2148
ndim_notime = ndim-1;
2149
}
2150
else
2151
/* if we don't have time dimension */
2152
{
2153
for (i = 0; i < ndim; i++)
2154
{
2155
count_notime[i] = count[i];
2156
start_notime[i] = start[i];
2157
}
2158
ndim_notime = ndim;
2159
}
2160
}
2161
/* 2) .bp written by C */
2162
else
2163
{
2164
/* Get the timesteps we need to read */
2165
if (timedim > -1)
2166
{
2167
/* timedim has to be the 1st dimension. To be extended to handle
2168
the cases timedim at any dimension */
2169
if (timedim != 0)
2170
{
2171
adios_error (err_no_data_at_timestep,"Variable (id=%d) has wrong time dimension",
2172
varid);
2173
return -adios_errno;
2174
}
2175
2176
if (futils_is_called_from_fortran())
2177
{
2178
start_time = start[ndim - 1] + fh->tidx_start;
2179
stop_time = start_time + count[ndim -1] - 1;
2180
}
2181
else
2182
{
2183
start_time = start[0] + fh->tidx_start;
2184
stop_time = start_time + count[0] - 1;
2185
}
2186
2187
start_time = start[timedim] + fh->tidx_start;
2188
stop_time = start_time + count[timedim] - 1;
2189
}
2190
else
2191
{
2192
/* timeless variable */
2193
start_time = fh->tidx_start;
2194
stop_time = fh->tidx_start;
2195
}
2196
2197
/* No need to flip dims, timedim as they are already in C order. */
2198
//swap_order(ndim, dims, &timedim);
2199
2200
/* Take out the time dimension from start[] and count[] */
2201
if (timedim == -1) /* timeless variable */
2202
{
2203
for (i = 0; i < ndim; i++)
2204
{
2205
count_notime[i] = count[i];
2206
start_notime[i] = start[i];
2207
}
2208
ndim_notime = ndim;
2209
}
2210
/* if we have time dimension */
2211
else
2212
{
2213
j = 0;
2214
if (futils_is_called_from_fortran())
2215
temp_timedim = ndim - 1;
2216
else
2217
temp_timedim = 0;
2218
2219
for (i = 0; i < temp_timedim; i++)
2220
{
2221
count_notime[j] = count[i];
2222
start_notime[j] = start[i];
2223
j++;
2224
}
2225
i++; // skip timedim
2226
for (; i < ndim; i++)
2227
{
2228
count_notime[j] = count[i];
2229
start_notime[j] = start[i];
2230
j++;
2231
}
2232
ndim_notime = ndim - 1;
2233
}
2234
}
2235
2236
/* Fortran reader was reported of Fortran dimension order so it gives counts and starts in that order.
2237
We need to swap them here to read correctly in C order */
2238
if ( futils_is_called_from_fortran()) {
2239
swap_order(ndim_notime, count_notime, &timedim);
2240
swap_order(ndim_notime, start_notime, &timedim);
2241
}
2242
2243
/* items_read = how many data elements are we going to read in total (per timestep) */
2244
items_read = 1;
2245
for (i = 0; i < ndim_notime; i++)
2246
items_read *= count_notime[i];
2247
2248
MPI_Comm_rank(gh->fh->comm, &rank);
2249
2250
size_of_type = bp_get_type_size (var_root->type, var_root->characteristics [0].value);
2251
2252
/* For each timestep, do reading separately (they are stored in different sets of process groups */
2253
for (timestep = start_time; timestep <= stop_time; timestep++) {
2254
2255
// pgoffset = the starting offset for the given time step
2256
// pgcount = number of process groups of that time step
2257
pgoffset = fh->gvar_h->time_index[0][gh->group_id][timestep - fh->tidx_start];
2258
pgcount = fh->gvar_h->time_index[1][gh->group_id][timestep - fh->tidx_start];
2259
2260
start_idx = -1;
2261
for (i=0;i<var_root->characteristics_count;i++) {
2262
if ( ( var_root->characteristics[i].offset > fh->gvar_h->pg_offsets[pgoffset])
2263
&& ( (pgoffset + pgcount == fh->mfooter.pgs_count)
2264
||( var_root->characteristics[i].offset < fh->gvar_h->pg_offsets[pgoffset + 1]))
2265
)
2266
{
2267
start_idx = i;
2268
break;
2269
}
2270
}
2271
/*
2272
printf ("var_root->characteristics_count = %d\n", var_root->characteristics_count);
2273
printf ("pg offset 0 = %lld\n", fh->gvar_h->pg_offsets[pgoffset]);
2274
printf ("pg offset 1 = %lld\n", fh->gvar_h->pg_offsets[pgoffset + 1]);
2275
printf ("var offset 3 = %lld\n", var_root->characteristics[3].offset);
2276
printf ("var offset 2 = %lld\n", var_root->characteristics[2].offset);
2277
printf ("var offset 1 = %lld\n", var_root->characteristics[1].offset);
2278
printf ("pgcount = %lld\n", pgcount);
2279
*/
2280
for (i=var_root->characteristics_count-1;i>-1;i--) {
2281
if ( ( var_root->characteristics[i].offset > fh->gvar_h->pg_offsets[pgoffset])
2282
&& ( (pgoffset + pgcount == fh->mfooter.pgs_count)
2283
||( var_root->characteristics[i].offset < fh->gvar_h->pg_offsets[pgoffset + pgcount]))
2284
)
2285
{
2286
stop_idx = i;
2287
break;
2288
}
2289
}
2290
2291
if (start_idx<0) {
2292
adios_error (err_no_data_at_timestep,"Variable (id=%d) has no data at %d time step",
2293
varid, timestep);
2294
return -adios_errno;
2295
}
2296
2297
if (ndim_notime == 0) {
2298
/* READ A SCALAR VARIABLE */
2299
2300
slice_size = size_of_type;
2301
idx = 0; // macros below need it
2302
2303
if (var_root->type == adios_string) {
2304
// strings are stored without \0 in file
2305
// size_of_type here includes \0 so decrease by one
2306
size_of_type--;
2307
}
2308
2309
if (var_root->characteristics[start_idx+idx].payload_offset > 0) {
2310
slice_offset = var_root->characteristics[start_idx+idx].payload_offset;
2311
MPI_FILE_READ_OPS
2312
} else {
2313
slice_offset = 0;
2314
MPI_FILE_READ_OPS1
2315
}
2316
2317
memcpy((char *)data+total_size, fh->b->buff + fh->b->offset, size_of_type);
2318
if (fh->mfooter.change_endianness == adios_flag_yes) {
2319
change_endianness((char *)data+total_size, size_of_type, var_root->type);
2320
}
2321
2322
if (var_root->type == adios_string) {
2323
// add \0 to the end of string
2324
// size_of_type here is the length of string
2325
// FIXME: how would this work for strings written over time?
2326
((char*)data+total_size)[size_of_type] = '\0';
2327
}
2328
2329
total_size += size_of_type;
2330
continue;
2331
}
2332
2333
/* READ AN ARRAY VARIABLE */
2334
//int * idx_table = (int *) malloc (sizeof(int) * pgcount);
2335
//int * idx_table = (int *) malloc (sizeof(int) * (var_root->characteristics_count - start_idx));
2336
int * idx_table = (int *) malloc (sizeof(int) * (stop_idx - start_idx + 1));
2337
2338
uint64_t write_offset = 0;
2339
int npg = 0;
2340
tmpcount = 0;
2341
if (pgcount > var_root->characteristics_count)
2342
pgcount = var_root->characteristics_count;
2343
2344
// loop over the list of pgs to read from one-by-one
2345
for (idx = 0; idx < stop_idx - start_idx + 1; idx++) {
2346
int flag;
2347
datasize = 1;
2348
nloop = 1;
2349
var_stride = 1;
2350
dset_stride = 1;
2351
idx_table[idx] = 1;
2352
uint64_t payload_size = size_of_type;
2353
2354
/* Each pg can have a different sized array, so we need the actual dimensions from it */
2355
is_global = adios_read_bp_get_dimensioncharacteristics( &(var_root->characteristics[start_idx + idx]),
2356
ldims, gdims, offsets);
2357
if (!is_global) {
2358
// we use gdims below, which is 0 for a local array; set to ldims here
2359
for (j = 0; j< ndim; j++) {
2360
gdims[j]=ldims[j];
2361
}
2362
}
2363
2364
/* Again, a Fortran written file has the dimensions in Fortran order we need to swap here */
2365
//if (file_is_fortran != futils_is_called_from_fortran()) {
2366
if (file_is_fortran ) {
2367
i=-1;
2368
swap_order(ndim, gdims, &(i)); // i is dummy
2369
swap_order(ndim, ldims, &(i));
2370
swap_order(ndim, offsets, &(i));
2371
}
2372
2373
/* take out the time dimension */
2374
/* For C, gdims and offset are one size shorter because the timedim part is missing,
2375
so we take it out only for fortran files
2376
*/
2377
if (timedim > -1) {
2378
for (i = timedim; i < ndim-1; i++) {
2379
ldims[i] = ldims[i+1];
2380
if (file_is_fortran) {
2381
gdims[i] = gdims[i+1];
2382
offsets[i] = offsets[i+1];
2383
}
2384
}
2385
}
2386
/*
2387
printf("ldims = "); for (j = 0; j<ndim; j++) printf("%d ",ldims[j]); printf("\n");
2388
printf("gdims = "); for (j = 0; j<ndim; j++) printf("%d ",gdims[j]); printf("\n");
2389
printf("offsets = "); for (j = 0; j<ndim; j++) printf("%d ",offsets[j]); printf("\n");
2390
printf("count_notime = "); for (j = 0; j<ndim_notime; j++) printf("%d ",count_notime[j]); printf("\n");
2391
printf("start_notime = "); for (j = 0; j<ndim_notime; j++) printf("%d ",start_notime[j]); printf("\n");
2392
*/
2393
for (j = 0; j < ndim_notime; j++) {
2394
2395
payload_size *= ldims [j];
2396
2397
if ( (count_notime[j] > gdims[j])
2398
|| (start_notime[j] > gdims[j])
2399
|| (start_notime[j] + count_notime[j] > gdims[j])){
2400
adios_error ( err_out_of_bound, "Error: Variable (id=%d) out of bound ("
2401
"the data in dimension %d to read is %llu elements from index %llu"
2402
" but the actual data is [0,%llu])",
2403
varid, j+1, count_notime[j], start_notime[j], gdims[j] - 1);
2404
return -adios_errno;
2405
}
2406
2407
/* check if there is any data in this pg and this dimension to read in */
2408
flag = (offsets[j] >= start_notime[j]
2409
&& offsets[j] < start_notime[j] + count_notime[j])
2410
|| (offsets[j] < start_notime[j]
2411
&& offsets[j] + ldims[j] > start_notime[j] + count_notime[j])
2412
|| (offsets[j] + ldims[j] > start_notime[j]
2413
&& offsets[j] + ldims[j] <= start_notime[j] + count_notime[j]);
2414
idx_table [idx] = idx_table[idx] && flag;
2415
}
2416
2417
if ( !idx_table[idx] ) {
2418
continue;
2419
}
2420
++npg;
2421
2422
/* determined how many (fastest changing) dimensions can we read in in one read */
2423
int hole_break;
2424
for (i = ndim_notime - 1; i > -1; i--) {
2425
if (offsets[i] == start_notime[i] && ldims[i] == count_notime[i]) {
2426
datasize *= ldims[i];
2427
}
2428
else
2429
break;
2430
}
2431
2432
hole_break = i;
2433
slice_offset = 0;
2434
slice_size = 0;
2435
2436
if (hole_break == -1) {
2437
/* The complete read happens to be exactly one pg, and the entire pg */
2438
/* This means we enter this only once, and npg=1 at the end */
2439
/* This is a rare case. FIXME: cannot eliminate this? */
2440
slice_size = payload_size;
2441
2442
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
2443
slice_offset = var_root->characteristics[start_idx + idx].payload_offset;
2444
MPI_FILE_READ_OPS
2445
} else {
2446
slice_offset = 0;
2447
MPI_FILE_READ_OPS1
2448
}
2449
2450
memcpy( (char *)data, fh->b->buff + fh->b->offset, slice_size);
2451
if (fh->mfooter.change_endianness == adios_flag_yes) {
2452
change_endianness(data, slice_size, var_root->type);
2453
}
2454
}
2455
else if (hole_break == 0)
2456
{
2457
/* The slowest changing dimensions should not be read completely but
2458
we still need to read only one block */
2459
int isize;
2460
uint64_t size_in_dset = 0;
2461
uint64_t offset_in_dset = 0;
2462
uint64_t offset_in_var = 0;
2463
2464
isize = offsets[0] + ldims[0];
2465
if (start_notime[0] >= offsets[0]) {
2466
// head is in
2467
if (start_notime[0]<isize) {
2468
if (start_notime[0] + count_notime[0] > isize)
2469
size_in_dset = isize - start_notime[0];
2470
else
2471
size_in_dset = count_notime[0];
2472
offset_in_dset = start_notime[0] - offsets[0];
2473
offset_in_var = 0;
2474
}
2475
}
2476
else {
2477
// middle is in
2478
if (isize < start_notime[0] + count_notime[0])
2479
size_in_dset = ldims[0];
2480
else
2481
// tail is in
2482
size_in_dset = count_notime[0] + start_notime[0] - offsets[0];
2483
offset_in_dset = 0;
2484
offset_in_var = offsets[0] - start_notime[0];
2485
}
2486
2487
slice_size = size_in_dset * datasize * size_of_type;
2488
write_offset = offset_in_var * datasize * size_of_type;
2489
2490
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
2491
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
2492
+ offset_in_dset * datasize * size_of_type;
2493
MPI_FILE_READ_OPS
2494
} else {
2495
slice_offset = 0;
2496
MPI_FILE_READ_OPS1
2497
}
2498
2499
memcpy ((char *)data + write_offset, fh->b->buff + fh->b->offset, slice_size);
2500
if (fh->mfooter.change_endianness == adios_flag_yes) {
2501
change_endianness((char *)data + write_offset, slice_size, var_root->type);
2502
}
2503
2504
//write_offset += slice_size;
2505
}
2506
else
2507
{
2508
2509
uint64_t stride_offset = 0;
2510
int isize;
2511
uint64_t size_in_dset[10];
2512
uint64_t offset_in_dset[10];
2513
uint64_t offset_in_var[10];
2514
memset(size_in_dset, 0 , 10 * 8);
2515
memset(offset_in_dset, 0 , 10 * 8);
2516
memset(offset_in_var, 0 , 10 * 8);
2517
int hit = 0;
2518
for ( i = 0; i < ndim_notime ; i++) {
2519
isize = offsets[i] + ldims[i];
2520
if (start_notime[i] >= offsets[i]) {
2521
// head is in
2522
if (start_notime[i]<isize) {
2523
if (start_notime[i] + count_notime[i] > isize)
2524
size_in_dset[i] = isize - start_notime[i];
2525
else
2526
size_in_dset[i] = count_notime[i];
2527
offset_in_dset[i] = start_notime[i] - offsets[i];
2528
offset_in_var[i] = 0;
2529
hit = 1 + hit * 10;
2530
}
2531
else
2532
hit = -1;
2533
}
2534
else {
2535
// middle is in
2536
if (isize < start_notime[i] + count_notime[i]) {
2537
size_in_dset[i] = ldims[i];
2538
hit = 2 + hit * 10;
2539
}
2540
else {
2541
// tail is in
2542
size_in_dset[i] = count_notime[i] + start_notime[i] - offsets[i];
2543
hit = 3 + hit * 10;
2544
}
2545
offset_in_dset[i] = 0;
2546
offset_in_var[i] = offsets[i] - start_notime[i];
2547
}
2548
}
2549
2550
datasize = 1;
2551
var_stride = 1;
2552
2553
for ( i = ndim_notime-1; i >= hole_break; i--) {
2554
datasize *= size_in_dset[i];
2555
dset_stride *= ldims[i];
2556
var_stride *= count_notime[i];
2557
}
2558
2559
uint64_t start_in_payload = 0, end_in_payload = 0, s = 1;
2560
for (i = ndim_notime - 1; i > -1; i--) {
2561
start_in_payload += s * offset_in_dset[i] * size_of_type;
2562
end_in_payload += s * (offset_in_dset[i] + size_in_dset[i] - 1) * size_of_type;
2563
s *= ldims[i];
2564
}
2565
2566
slice_size = end_in_payload - start_in_payload + 1 * size_of_type;
2567
2568
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
2569
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
2570
+ start_in_payload;
2571
MPI_FILE_READ_OPS
2572
2573
for ( i = 0; i < ndim_notime ; i++) {
2574
offset_in_dset[i] = 0;
2575
}
2576
} else {
2577
slice_offset = start_in_payload;
2578
MPI_FILE_READ_OPS1
2579
}
2580
2581
uint64_t var_offset = 0;
2582
uint64_t dset_offset = 0;
2583
for ( i = 0; i < hole_break; i++) {
2584
nloop *= size_in_dset[i];
2585
}
2586
2587
for ( i = 0; i < ndim_notime ; i++) {
2588
var_offset = offset_in_var[i] + var_offset * count_notime[i];
2589
dset_offset = offset_in_dset[i] + dset_offset * ldims[i];
2590
}
2591
2592
copy_data (data
2593
,fh->b->buff + fh->b->offset
2594
,0
2595
,hole_break
2596
,size_in_dset
2597
,ldims
2598
,count_notime
2599
,var_stride
2600
,dset_stride
2601
,var_offset
2602
,dset_offset
2603
,datasize
2604
,size_of_type
2605
);
2606
2607
}
2608
} // end for (idx ... loop over pgs
2609
2610
free (idx_table);
2611
2612
total_size += items_read * size_of_type;
2613
// shift target pointer for next read in
2614
data = (char *)data + (items_read * size_of_type);
2615
2616
} // end for (timestep ... loop over timesteps
2617
2618
free (dims);
2619
2620
return total_size;
2621
}
2622
2623
// Search for the start var index.
2624
static int get_var_start_index (struct adios_index_var_struct_v1 * v, int t)
2625
{
2626
int i = 0;
2627
2628
while (i < v->characteristics_count) {
2629
if (v->characteristics[i].time_index == t) {
2630
return i;
2631
}
2632
2633
i++;
2634
}
2635
2636
return -1;
2637
}
2638
2639
// Search for the stop var index
2640
static int get_var_stop_index (struct adios_index_var_struct_v1 * v, int t)
2641
{
2642
int i = v->characteristics_count - 1;
2643
2644
while (i > -1) {
2645
if (v->characteristics[i].time_index == t) {
2646
return i;
2647
}
2648
2649
i--;
2650
}
2651
2652
return -1;
2653
}
2654
2655
int64_t adios_read_bp_read_var_byid2 (ADIOS_GROUP * gp,
2656
int varid,
2657
const uint64_t * start,
2658
const uint64_t * count,
2659
void * data)
2660
{
2661
struct BP_GROUP * gh;
2662
struct BP_FILE * fh;
2663
struct adios_index_var_struct_v1 * var_root;
2664
struct adios_var_header_struct_v1 var_header;
2665
struct adios_var_payload_struct_v1 var_payload;
2666
int i,j,k, idx, t;
2667
int start_time, stop_time;
2668
int start_idx, stop_idx;
2669
int ndim, ndim_notime, has_subfile, file_is_fortran;
2670
uint64_t size, * dims;
2671
uint64_t ldims[32], gdims[32], offsets[32];
2672
uint64_t datasize, nloop, dset_stride,var_stride, total_size=0, items_read;
2673
uint64_t count_notime[32], start_notime[32];
2674
int timedim = -1, temp_timedim, is_global = 0, size_of_type;
2675
uint64_t slice_offset, slice_size, tmpcount = 0;
2676
uint64_t datatimeoffset = 0; // offset in data to write a given timestep
2677
MPI_Status status;
2678
2679
gh = (struct BP_GROUP *) gp->gh;
2680
fh = gh->fh;
2681
2682
file_is_fortran = (fh->pgs_root->adios_host_language_fortran == adios_flag_yes);
2683
has_subfile = fh->mfooter.version & ADIOS_VERSION_HAVE_SUBFILE;
2684
2685
var_root = gh->vars_root; /* first variable of this group. Need to traverse the list */
2686
for (i = 0; i< varid && var_root; i++) {
2687
var_root = var_root->next;
2688
}
2689
2690
if (i!=varid) {
2691
adios_error (err_corrupted_variable,
2692
"Variable id=%d is valid but was not found in internal data structures!",
2693
varid);
2694
return -adios_errno;
2695
}
2696
2697
/* Get dimensions and flip if caller != writer language */
2698
adios_read_bp_get_dimensions (var_root, fh->tidx_stop - fh->tidx_start + 1, file_is_fortran,
2699
&ndim, &dims, &timedim);
2700
2701
/* In a Fortran written files, dimensions are in reversed order for C */
2702
//if ( file_is_fortran != futils_is_called_from_fortran() )
2703
if (file_is_fortran)
2704
swap_order(ndim, dims, &timedim);
2705
2706
/* Take out the time dimension from start[] and count[] */
2707
if (timedim == -1) {
2708
/* For timeless var, we still search from fh->tidx_start to fh->tidx_stop
2709
to handle the situation that some variables are dumped out in selected timesteps
2710
*/
2711
start_time = fh->tidx_start;
2712
stop_time = fh->tidx_stop;
2713
2714
for (i = 0; i < ndim; i++) {
2715
count_notime[i] = count[i];
2716
start_notime[i] = start[i];
2717
}
2718
ndim_notime = ndim;
2719
} else {
2720
j = 0;
2721
if (futils_is_called_from_fortran())
2722
temp_timedim = ndim - 1;
2723
else
2724
temp_timedim = 0;
2725
2726
start_time = start[temp_timedim] + fh->tidx_start;
2727
stop_time = start_time + count[temp_timedim] - 1;
2728
2729
for (i = 0; i < temp_timedim; i++) {
2730
count_notime[j] = count[i];
2731
start_notime[j] = start[i];
2732
j++;
2733
}
2734
i++; // skip timedim
2735
for (; i < ndim; i++) {
2736
count_notime[j] = count[i];
2737
start_notime[j] = start[i];
2738
j++;
2739
}
2740
ndim_notime = ndim-1;
2741
}
2742
2743
/* Fortran reader was reported of Fortran dimension order so it gives counts and starts in that order.
2744
We need to swap them here to read correctly in C order */
2745
if ( futils_is_called_from_fortran()) {
2746
swap_order(ndim_notime, count_notime, &timedim);
2747
swap_order(ndim_notime, start_notime, &timedim);
2748
}
2749
2750
/* items_read = how many data elements are we going to read in total (per timestep) */
2751
items_read = 1;
2752
for (i = 0; i < ndim_notime; i++)
2753
items_read *= count_notime[i];
2754
2755
size_of_type = bp_get_type_size (var_root->type, var_root->characteristics [0].value);
2756
2757
/* For each timestep, do reading separately (they are stored in different sets of process groups */
2758
for (t = start_time; t <= stop_time; t++) {
2759
start_idx = get_var_start_index(var_root, t);
2760
stop_idx = get_var_stop_index(var_root, t);
2761
2762
if (start_idx < 0 || stop_idx < 0) {
2763
adios_error (err_no_data_at_timestep,"Variable (id=%d) has no data at %d time step",
2764
varid, t);
2765
// return -adios_errno;
2766
continue;
2767
}
2768
2769
if (ndim_notime == 0) {
2770
/* READ A SCALAR VARIABLE */
2771
slice_size = size_of_type;
2772
idx = 0; // macros below need it
2773
2774
if (var_root->type == adios_string) {
2775
// strings are stored without \0 in file
2776
// size_of_type here includes \0 so decrease by one
2777
size_of_type--;
2778
}
2779
2780
if (var_root->characteristics[start_idx+idx].payload_offset > 0) {
2781
slice_offset = var_root->characteristics[start_idx+idx].payload_offset;
2782
if (!has_subfile) {
2783
MPI_FILE_READ_OPS
2784
} else {
2785
MPI_FILE_READ_OPS2
2786
}
2787
} else {
2788
slice_offset = 0;
2789
MPI_FILE_READ_OPS1
2790
}
2791
2792
memcpy((char *)data+total_size, fh->b->buff + fh->b->offset, size_of_type);
2793
if (fh->mfooter.change_endianness == adios_flag_yes) {
2794
change_endianness((char *)data+total_size, size_of_type, var_root->type);
2795
}
2796
2797
if (var_root->type == adios_string) {
2798
// add \0 to the end of string
2799
// size_of_type here is the length of string
2800
// FIXME: how would this work for strings written over time?
2801
((char*)data+total_size)[size_of_type] = '\0';
2802
}
2803
2804
total_size += size_of_type;
2805
2806
if (timedim == -1)
2807
break;
2808
else
2809
continue;
2810
}
2811
2812
/* READ AN ARRAY VARIABLE */
2813
int * idx_table = (int *) malloc (sizeof(int) * (stop_idx - start_idx + 1));
2814
2815
uint64_t write_offset = 0;
2816
int npg = 0;
2817
tmpcount = 0;
2818
// loop over the list of pgs to read from one-by-one
2819
for (idx = 0; idx < stop_idx - start_idx + 1; idx++) {
2820
int flag;
2821
datasize = 1;
2822
nloop = 1;
2823
var_stride = 1;
2824
dset_stride = 1;
2825
idx_table[idx] = 1;
2826
uint64_t payload_size = size_of_type;
2827
2828
/* Each pg can have a different sized array, so we need the actual dimensions from it */
2829
is_global = adios_read_bp_get_dimensioncharacteristics( &(var_root->characteristics[start_idx + idx]),
2830
ldims, gdims, offsets);
2831
if (!is_global) {
2832
// we use gdims below, which is 0 for a local array; set to ldims here
2833
for (j = 0; j< ndim; j++) {
2834
gdims[j]=ldims[j];
2835
}
2836
// we need to read only the first PG, not all, so let's prevent a second loop
2837
stop_idx = start_idx;
2838
}
2839
2840
/* Again, a Fortran written file has the dimensions in Fortran order we need to swap here */
2841
//if (file_is_fortran != futils_is_called_from_fortran()) {
2842
if (file_is_fortran) {
2843
i=-1;
2844
swap_order(ndim, gdims, &(i)); // i is dummy
2845
swap_order(ndim, ldims, &(i));
2846
swap_order(ndim, offsets, &(i));
2847
}
2848
2849
/* take out the time dimension */
2850
/* For C, gdims and offset are one size shorter because the timedim part is missing,
2851
so we take it out only for fortran files
2852
*/
2853
if (timedim > -1) {
2854
for (i = timedim; i < ndim-1; i++) {
2855
ldims[i] = ldims[i+1];
2856
if (file_is_fortran) {
2857
gdims[i] = gdims[i+1];
2858
offsets[i] = offsets[i+1];
2859
}
2860
}
2861
}
2862
2863
/*
2864
printf("ldims = "); for (j = 0; j<ndim; j++) printf("%d ",ldims[j]); printf("\n");
2865
printf("gdims = "); for (j = 0; j<ndim; j++) printf("%d ",gdims[j]); printf("\n");
2866
printf("offsets = "); for (j = 0; j<ndim; j++) printf("%d ",offsets[j]); printf("\n");
2867
printf("count_notime = "); for (j = 0; j<ndim_notime; j++) printf("%d ",count_notime[j]); printf("\n");
2868
printf("start_notime = "); for (j = 0; j<ndim_notime; j++) printf("%d ",start_notime[j]); printf("\n");
2869
*/
2870
2871
for (j = 0; j < ndim_notime; j++) {
2872
2873
payload_size *= ldims [j];
2874
2875
if ( (count_notime[j] > gdims[j])
2876
|| (start_notime[j] > gdims[j])
2877
|| (start_notime[j] + count_notime[j] > gdims[j])){
2878
adios_error ( err_out_of_bound, "Error: Variable (id=%d) out of bound ("
2879
"the data in dimension %d to read is %llu elements from index %llu"
2880
" but the actual data is [0,%llu])",
2881
varid, j+1, count_notime[j], start_notime[j], gdims[j] - 1);
2882
return -adios_errno;
2883
}
2884
2885
/* check if there is any data in this pg and this dimension to read in */
2886
flag = (offsets[j] >= start_notime[j]
2887
&& offsets[j] < start_notime[j] + count_notime[j])
2888
|| (offsets[j] < start_notime[j]
2889
&& offsets[j] + ldims[j] > start_notime[j] + count_notime[j])
2890
|| (offsets[j] + ldims[j] > start_notime[j]
2891
&& offsets[j] + ldims[j] <= start_notime[j] + count_notime[j]);
2892
idx_table [idx] = idx_table[idx] && flag;
2893
}
2894
2895
if ( !idx_table[idx] ) {
2896
continue;
2897
}
2898
++npg;
2899
2900
/* determined how many (fastest changing) dimensions can we read in in one read */
2901
int hole_break;
2902
for (i = ndim_notime - 1; i > -1; i--) {
2903
if (offsets[i] == start_notime[i] && ldims[i] == count_notime[i]) {
2904
datasize *= ldims[i];
2905
}
2906
else
2907
break;
2908
}
2909
2910
hole_break = i;
2911
slice_offset = 0;
2912
slice_size = 0;
2913
if (hole_break == -1) {
2914
/* The complete read happens to be exactly one pg, and the entire pg */
2915
/* This means we enter this only once, and npg=1 at the end */
2916
/* This is a rare case. FIXME: cannot eliminate this? */
2917
slice_size = payload_size;
2918
2919
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
2920
slice_offset = var_root->characteristics[start_idx + idx].payload_offset;
2921
if (!has_subfile) {
2922
MPI_FILE_READ_OPS
2923
} else {
2924
MPI_FILE_READ_OPS2
2925
}
2926
} else {
2927
slice_offset = 0;
2928
MPI_FILE_READ_OPS1
2929
}
2930
2931
memcpy( (char *)data, fh->b->buff + fh->b->offset, slice_size);
2932
if (fh->mfooter.change_endianness == adios_flag_yes) {
2933
change_endianness(data, slice_size, var_root->type);
2934
}
2935
}
2936
else if (hole_break == 0)
2937
{
2938
/* The slowest changing dimensions should not be read completely but
2939
we still need to read only one block */
2940
int isize;
2941
uint64_t size_in_dset = 0;
2942
uint64_t offset_in_dset = 0;
2943
uint64_t offset_in_var = 0;
2944
2945
isize = offsets[0] + ldims[0];
2946
if (start_notime[0] >= offsets[0]) {
2947
// head is in
2948
if (start_notime[0]<isize) {
2949
if (start_notime[0] + count_notime[0] > isize)
2950
size_in_dset = isize - start_notime[0];
2951
else
2952
size_in_dset = count_notime[0];
2953
offset_in_dset = start_notime[0] - offsets[0];
2954
offset_in_var = 0;
2955
}
2956
}
2957
else {
2958
// middle is in
2959
if (isize < start_notime[0] + count_notime[0])
2960
size_in_dset = ldims[0];
2961
else
2962
// tail is in
2963
size_in_dset = count_notime[0] + start_notime[0] - offsets[0];
2964
offset_in_dset = 0;
2965
offset_in_var = offsets[0] - start_notime[0];
2966
}
2967
2968
slice_size = size_in_dset * datasize * size_of_type;
2969
write_offset = offset_in_var * datasize * size_of_type;
2970
2971
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
2972
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
2973
+ offset_in_dset * datasize * size_of_type;
2974
if (!has_subfile) {
2975
MPI_FILE_READ_OPS
2976
} else {
2977
MPI_FILE_READ_OPS2
2978
}
2979
2980
} else {
2981
slice_offset = 0;
2982
MPI_FILE_READ_OPS1
2983
}
2984
2985
memcpy ((char *)data + write_offset, fh->b->buff + fh->b->offset, slice_size);
2986
if (fh->mfooter.change_endianness == adios_flag_yes) {
2987
change_endianness((char *)data + write_offset, slice_size, var_root->type);
2988
}
2989
2990
//write_offset += slice_size;
2991
}
2992
else
2993
{
2994
2995
uint64_t stride_offset = 0;
2996
int isize;
2997
uint64_t size_in_dset[10];
2998
uint64_t offset_in_dset[10];
2999
uint64_t offset_in_var[10];
3000
memset(size_in_dset, 0 , 10 * 8);
3001
memset(offset_in_dset, 0 , 10 * 8);
3002
memset(offset_in_var, 0 , 10 * 8);
3003
int hit = 0;
3004
for ( i = 0; i < ndim_notime ; i++) {
3005
isize = offsets[i] + ldims[i];
3006
if (start_notime[i] >= offsets[i]) {
3007
// head is in
3008
if (start_notime[i]<isize) {
3009
if (start_notime[i] + count_notime[i] > isize)
3010
size_in_dset[i] = isize - start_notime[i];
3011
else
3012
size_in_dset[i] = count_notime[i];
3013
offset_in_dset[i] = start_notime[i] - offsets[i];
3014
offset_in_var[i] = 0;
3015
hit = 1 + hit * 10;
3016
}
3017
else
3018
hit = -1;
3019
}
3020
else {
3021
// middle is in
3022
if (isize < start_notime[i] + count_notime[i]) {
3023
size_in_dset[i] = ldims[i];
3024
hit = 2 + hit * 10;
3025
}
3026
else {
3027
// tail is in
3028
size_in_dset[i] = count_notime[i] + start_notime[i] - offsets[i];
3029
hit = 3 + hit * 10;
3030
}
3031
offset_in_dset[i] = 0;
3032
offset_in_var[i] = offsets[i] - start_notime[i];
3033
}
3034
}
3035
3036
datasize = 1;
3037
var_stride = 1;
3038
3039
for ( i = ndim_notime-1; i >= hole_break; i--) {
3040
datasize *= size_in_dset[i];
3041
dset_stride *= ldims[i];
3042
var_stride *= count_notime[i];
3043
}
3044
3045
uint64_t start_in_payload = 0, end_in_payload = 0, s = 1;
3046
for (i = ndim_notime - 1; i > -1; i--) {
3047
start_in_payload += s * offset_in_dset[i] * size_of_type;
3048
end_in_payload += s * (offset_in_dset[i] + size_in_dset[i] - 1) * size_of_type;
3049
s *= ldims[i];
3050
}
3051
3052
slice_size = end_in_payload - start_in_payload + 1 * size_of_type;
3053
3054
if (var_root->characteristics[start_idx + idx].payload_offset > 0) {
3055
slice_offset = var_root->characteristics[start_idx + idx].payload_offset
3056
+ start_in_payload;
3057
if (!has_subfile) {
3058
MPI_FILE_READ_OPS
3059
} else {
3060
MPI_FILE_READ_OPS2
3061
}
3062
3063
for ( i = 0; i < ndim_notime ; i++) {
3064
offset_in_dset[i] = 0;
3065
}
3066
} else {
3067
slice_offset = start_in_payload;
3068
MPI_FILE_READ_OPS1
3069
}
3070
3071
uint64_t var_offset = 0;
3072
uint64_t dset_offset = 0;
3073
for ( i = 0; i < hole_break; i++) {
3074
nloop *= size_in_dset[i];
3075
}
3076
3077
for ( i = 0; i < ndim_notime ; i++) {
3078
var_offset = offset_in_var[i] + var_offset * count_notime[i];
3079
dset_offset = offset_in_dset[i] + dset_offset * ldims[i];
3080
}
3081
3082
copy_data (data
3083
,fh->b->buff + fh->b->offset
3084
,0
3085
,hole_break
3086
,size_in_dset
3087
,ldims
3088
,count_notime
3089
,var_stride
3090
,dset_stride
3091
,var_offset
3092
,dset_offset
3093
,datasize
3094
,size_of_type
3095
);
3096
}
3097
} // end for (idx ... loop over pgs
3098
3099
free (idx_table);
3100
3101
total_size += items_read * size_of_type;
3102
// shift target pointer for next read in
3103
data = (char *)data + (items_read * size_of_type);
3104
3105
if (timedim == -1)
3106
break;
3107
} // end for (timestep ... loop over timesteps
3108
3109
free (dims);
3110
3111
return total_size;
3112
}
3113
3114
int64_t adios_read_bp_read_var_byid (ADIOS_GROUP * gp,
3115
int varid,
3116
const uint64_t * start,
3117
const uint64_t * count,
3118
void * data)
3119
{
3120
struct BP_GROUP * gh;
3121
struct BP_FILE * fh;
3122
int has_time_index_characteristic;
3123
3124
adios_errno = 0;
3125
if (!gp) {
3126
adios_error (err_invalid_group_struct, "Null pointer passed as group to adios_read_var()");
3127
return -adios_errno;
3128
}
3129
3130
gh = (struct BP_GROUP *) gp->gh;
3131
if (!gh) {
3132
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh group handle is NULL!");
3133
return -adios_errno;
3134
}
3135
3136
fh = gh->fh;
3137
if (!fh) {
3138
adios_error (err_invalid_group_struct, "Invalid ADIOS_GROUP struct: .gh->fh file handle is NULL!");
3139
return -adios_errno;
3140
}
3141
3142
has_time_index_characteristic = fh->mfooter.version & ADIOS_VERSION_HAVE_TIME_INDEX_CHARACTERISTIC;
3143
if (!has_time_index_characteristic) {
3144
// read older file format. Can be deleted later on.
3145
return adios_read_bp_read_var_byid1(gp, varid, start, count, data);
3146
} else {
3147
return adios_read_bp_read_var_byid2(gp, varid, start, count, data);
3148
}
3149
}
3150
3151
// NCSU - Timer series analysis, correlation
3152
double adios_stat_cor (ADIOS_VARINFO * vix, ADIOS_VARINFO * viy, char * characteristic, uint32_t time_start, uint32_t time_end, uint32_t lag)
3153
{
3154
int i,j;
3155
3156
double avg_x = 0.0, avg_y = 0.0, avg_lag = 0.0;
3157
double var_x = 0.0, var_y = 0.0, var_lag = 0.0;
3158
double cov = 0;
3159
3160
if (vix == NULL)
3161
{
3162
fprintf(stderr, "Variable not defined\n");
3163
return 0;
3164
}
3165
3166
// If the vix and viy are not time series objects, return.
3167
if ((vix->timedim < 0) && (viy->timedim < 0))
3168
{
3169
fprintf(stderr, "Covariance must involve timeseries data\n");
3170
return 0;
3171
}
3172
3173
uint32_t min = vix->dims[0] - 1;
3174
if (viy && (min > viy->dims[0] - 1))
3175
min = viy->dims[0] - 1;
3176
3177
if(time_start == 0 && time_end == 0)
3178
{ //global covariance
3179
if(viy == NULL) {
3180
fprintf(stderr, "Must have two variables for global covariance\n");
3181
return 0;
3182
}
3183
3184
// Assign vix to viy, and calculate covariance
3185
viy = vix;
3186
time_start = 0;
3187
time_end = min;
3188
}
3189
// Check the bounds of time
3190
if ( (time_start >= 0) && (time_start <= min)
3191
&& (time_end >= 0) && (time_end <= min)
3192
&& (time_start <= time_end))
3193
{
3194
if(viy == NULL) //user must want to run covariance against itself
3195
{
3196
if(! (time_end+lag) > min)
3197
{
3198
fprintf(stderr, "Must leave enough timesteps for lag\n");
3199
return 0;
3200
}
3201
3202
if (strcmp(characteristic, "average") == 0 || strcmp(characteristic, "avg") == 0)
3203
{
3204
for (i = time_start; i <= time_end; i ++)
3205
{
3206
avg_x += bp_value_to_double (adios_double, vix->avgs[i]) / (time_end - time_start + 1);
3207
avg_lag += bp_value_to_double (adios_double, vix->avgs[i + lag]) / (time_end - time_start + 1);
3208
}
3209
3210
for (i = time_start; i <= time_end; i ++)
3211
{
3212
double val_x = bp_value_to_double (adios_double, vix->avgs[i]);
3213
double val_lag = bp_value_to_double (adios_double, vix->avgs[i + lag]);
3214
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3215
var_lag += (val_lag - avg_lag) * (val_lag - avg_lag) / (time_end - time_start + 1);
3216
cov += (val_x - avg_x) * (val_lag - avg_lag) / (time_end - time_start + 1);
3217
}
3218
}
3219
else if (strcmp(characteristic, "standard deviation") == 0 || strcmp(characteristic, "std_dev") == 0)
3220
{
3221
for (i = time_start; i <= time_end; i ++)
3222
{
3223
avg_x += bp_value_to_double (adios_double, vix->std_devs[i]) / (time_end - time_start + 1);
3224
avg_lag += bp_value_to_double (adios_double, vix->std_devs[i + lag]) / (time_end - time_start + 1);
3225
}
3226
3227
for (i = time_start; i <= time_end; i ++)
3228
{
3229
double val_x = bp_value_to_double (adios_double, vix->std_devs[i]);
3230
double val_lag = bp_value_to_double (adios_double, vix->std_devs[i + lag]);
3231
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3232
var_lag += (val_lag - avg_lag) * (val_lag - avg_lag) / (time_end - time_start + 1);
3233
cov += (val_x - avg_x) * (val_lag - avg_lag) / (time_end - time_start + 1);
3234
}
3235
}
3236
else if (strcmp(characteristic, "minimum") == 0 || strcmp(characteristic, "min") == 0)
3237
{
3238
for (i = time_start; i <= time_end; i ++)
3239
{
3240
avg_x += bp_value_to_double (vix->type, vix->mins[i]) / (time_end - time_start + 1);
3241
avg_lag += bp_value_to_double (vix->type, vix->mins[i + lag]) / (time_end - time_start + 1);
3242
}
3243
3244
for (i = time_start; i <= time_end; i ++)
3245
{
3246
double val_x = bp_value_to_double (vix->type, vix->mins[i]);
3247
double val_lag = bp_value_to_double (vix->type, vix->mins[i + lag]);
3248
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3249
var_lag += (val_lag - avg_lag) * (val_lag - avg_lag) / (time_end - time_start + 1);
3250
cov += (val_x - avg_x) * (val_lag - avg_lag) / (time_end - time_start + 1);
3251
}
3252
}
3253
else if (strcmp(characteristic, "maximum") == 0 || strcmp(characteristic, "max") == 0)
3254
{
3255
for (i = time_start; i <= time_end; i ++)
3256
{
3257
avg_x += bp_value_to_double (vix->type, vix->maxs[i]) / (time_end - time_start + 1);
3258
avg_lag += bp_value_to_double (vix->type, vix->maxs[i]) / (time_end - time_start + 1);
3259
}
3260
3261
for (i = time_start; i <= time_end; i ++)
3262
{
3263
double val_x = bp_value_to_double (vix->type, vix->maxs[i]);
3264
double val_lag = bp_value_to_double (vix->type, vix->maxs[i + lag]);
3265
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3266
var_lag += (val_lag - avg_lag) * (val_lag - avg_lag) / (time_end - time_start + 1);
3267
cov += (val_x - avg_x) * (val_lag - avg_lag) / (time_end - time_start + 1);
3268
}
3269
}
3270
else
3271
{
3272
fprintf (stderr, "Unknown characteristic\n");
3273
return 0;
3274
}
3275
return cov / (sqrt (var_x) * sqrt (var_lag));
3276
}
3277
else
3278
{
3279
if (strcmp(characteristic, "average") == 0 || strcmp(characteristic, "avg") == 0)
3280
{
3281
for (i = time_start; i <= time_end; i ++)
3282
{
3283
avg_x += bp_value_to_double(adios_double, vix->avgs[i]) / (time_end - time_start + 1);
3284
avg_y += bp_value_to_double(adios_double, viy->avgs[i]) / (time_end - time_start + 1);
3285
}
3286
for (i = time_start; i <= time_end; i ++)
3287
{
3288
double val_x = bp_value_to_double (adios_double, vix->avgs[i]);
3289
double val_y = bp_value_to_double (adios_double, viy->avgs[i]);
3290
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3291
var_y += (val_y - avg_y) * (val_y - avg_y) / (time_end - time_start + 1);
3292
cov += (val_x - avg_x) * (val_y - avg_y) / (time_end - time_start + 1);
3293
}
3294
}
3295
else if (strcmp(characteristic, "standard deviation") == 0 || strcmp(characteristic, "std_dev") == 0)
3296
{
3297
for (i = time_start; i <= time_end; i ++)
3298
{
3299
avg_x += bp_value_to_double(adios_double, vix->std_devs[i]) / (time_end - time_start + 1);
3300
avg_y += bp_value_to_double(adios_double, viy->std_devs[i]) / (time_end - time_start + 1);
3301
}
3302
for (i = time_start; i <= time_end; i ++)
3303
{
3304
double val_x = bp_value_to_double (adios_double, vix->std_devs[i]);
3305
double val_y = bp_value_to_double (adios_double, viy->std_devs[i]);
3306
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3307
var_y += (val_y - avg_y) * (val_y - avg_y) / (time_end - time_start + 1);
3308
cov += (val_x - avg_x) * (val_y - avg_y) / (time_end - time_start + 1);
3309
}
3310
}
3311
else if (strcmp(characteristic, "minimum") == 0 || strcmp(characteristic, "min") == 0)
3312
{
3313
for (i = time_start; i <= time_end; i ++)
3314
{
3315
avg_x += bp_value_to_double(vix->type, vix->mins[i]) / (time_end - time_start + 1);
3316
avg_y += bp_value_to_double(viy->type, viy->mins[i]) / (time_end - time_start + 1);
3317
}
3318
for (i = time_start; i <= time_end; i ++)
3319
{
3320
double val_x = bp_value_to_double (vix->type, vix->mins[i]);
3321
double val_y = bp_value_to_double (viy->type, viy->mins[i]);
3322
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3323
var_y += (val_y - avg_y) * (val_y - avg_y) / (time_end - time_start + 1);
3324
cov += (val_x - avg_x) * (val_y - avg_y) / (time_end - time_start + 1);
3325
}
3326
}
3327
else if (strcmp(characteristic, "maximum") == 0 || strcmp(characteristic, "max") == 0)
3328
{
3329
for (i = time_start; i <= time_end; i ++)
3330
{
3331
avg_x += bp_value_to_double(vix->type, vix->maxs[i]) / (time_end - time_start + 1);
3332
avg_y += bp_value_to_double(vix->type, viy->maxs[i]) / (time_end - time_start + 1);
3333
}
3334
for (i = time_start; i <= time_end; i ++)
3335
{
3336
double val_x = bp_value_to_double (vix->type, vix->maxs[i]);
3337
double val_y = bp_value_to_double (viy->type, viy->maxs[i]);
3338
var_x += (val_x - avg_x) * (val_x - avg_x) / (time_end - time_start + 1);
3339
var_y += (val_y - avg_y) * (val_y - avg_y) / (time_end - time_start + 1);
3340
cov += (val_x - avg_x) * (val_y - avg_y) / (time_end - time_start + 1);
3341
}
3342
}
3343
else
3344
{
3345
fprintf (stderr, "Unknown characteristic\n");
3346
return 0;
3347
}
3348
return cov / (sqrt (var_x) * sqrt (var_y));
3349
}
3350
}
3351
else
3352
{
3353
fprintf (stderr, "Time values out of bounds\n");
3354
return 0;
3355
}
3356
}
3357
3358
// NCSU - Time series analysis, covariance
3359
//covariance(x,y) = sum(i=1,..N) [(x_1 - x_mean)(y_i - y_mean)]/N
3360
double adios_stat_cov (ADIOS_VARINFO * vix, ADIOS_VARINFO * viy, char * characteristic, uint32_t time_start, uint32_t time_end, uint32_t lag)
3361
{
3362
int i,j;
3363
3364
double avg_x = 0.0, avg_y = 0.0, avg_lag = 0.0;
3365
double cov = 0;
3366
3367
if (vix == NULL)
3368
{
3369
fprintf(stderr, "Variable not defined\n");
3370
return 0;
3371
}
3372
3373
// If the vix and viy are not time series objects, return.
3374
if ((vix->timedim < 0) && (viy->timedim < 0))
3375
{
3376
fprintf(stderr, "Covariance must involve timeseries data\n");
3377
return 0;
3378
}
3379
3380
uint32_t min = vix->dims[0] - 1;
3381
if (viy && (min > viy->dims[0] - 1))
3382
min = viy->dims[0] - 1;
3383
3384
if(time_start == 0 && time_end == 0)
3385
{ //global covariance
3386
if(viy == NULL) {
3387
fprintf(stderr, "Must have two variables for global covariance\n");
3388
return 0;
3389
}
3390
3391
// Assign vix to viy, and calculate covariance
3392
viy = vix;
3393
time_start = 0;
3394
time_end = min;
3395
}
3396
// Check the bounds of time
3397
if ( (time_start >= 0) && (time_start <= min)
3398
&& (time_end >= 0) && (time_end <= min)
3399
&& (time_start <= time_end))
3400
{
3401
if(viy == NULL) //user must want to run covariance against itself
3402
{
3403
if(! (time_end+lag) > min)
3404
{
3405
fprintf(stderr, "Must leave enough timesteps for lag\n");
3406
return 0;
3407
}
3408
3409
if (strcmp(characteristic, "average") == 0 || strcmp(characteristic, "avg") == 0)
3410
{
3411
for (i = time_start; i <= time_end; i ++)
3412
{
3413
avg_x += bp_value_to_double (adios_double, vix->avgs[i]) / (time_end - time_start + 1);
3414
avg_lag += bp_value_to_double (adios_double, vix->avgs[i + lag]) / (time_end - time_start + 1);
3415
}
3416
3417
for (i = time_start; i <= time_end; i ++)
3418
cov += (bp_value_to_double (adios_double, vix->avgs[i]) - avg_x) * (bp_value_to_double (adios_double, vix->avgs[i+lag]) - avg_lag) / (time_end - time_start + 1);
3419
}
3420
else if (strcmp(characteristic, "standard deviation") == 0 || strcmp(characteristic, "std_dev") == 0)
3421
{
3422
for (i = time_start; i <= time_end; i ++)
3423
{
3424
avg_x += bp_value_to_double (adios_double, vix->std_devs[i]) / (time_end - time_start + 1);
3425
avg_lag += bp_value_to_double (adios_double, vix->std_devs[i + lag]) / (time_end - time_start + 1);
3426
}
3427
3428
for (i = time_start; i <= time_end; i ++)
3429
cov += (bp_value_to_double (adios_double, vix->std_devs[i]) - avg_x) * (bp_value_to_double (adios_double, vix->std_devs[i+lag]) - avg_lag) / (time_end - time_start + 1);
3430
}
3431
else if (strcmp(characteristic, "minimum") == 0 || strcmp(characteristic, "min") == 0)
3432
{
3433
for (i = time_start; i <= time_end; i ++)
3434
{
3435
avg_x += bp_value_to_double (vix->type, vix->mins[i]) / (time_end - time_start + 1);
3436
avg_lag += bp_value_to_double (vix->type, vix->mins[i + lag]) / (time_end - time_start + 1);
3437
}
3438
3439
for (i = time_start; i <= time_end; i ++)
3440
cov += (bp_value_to_double (vix->type, vix->mins[i]) - avg_x) * (bp_value_to_double (vix->type, vix->mins[i+lag]) - avg_lag) / (time_end - time_start + 1);
3441
}
3442
else if (strcmp(characteristic, "maximum") == 0 || strcmp(characteristic, "max") == 0)
3443
{
3444
for (i = time_start; i <= time_end; i ++)
3445
{
3446
avg_x += bp_value_to_double (vix->type, vix->maxs[i]) / (time_end - time_start + 1);
3447
avg_lag += bp_value_to_double (vix->type, vix->maxs[i + lag]) / (time_end - time_start + 1);
3448
}
3449
3450
for (i = time_start; i <= time_end; i ++)
3451
cov += (bp_value_to_double (vix->type, vix->maxs[i]) - avg_x) * (bp_value_to_double (vix->type, vix->maxs[i+lag]) - avg_lag) / (time_end - time_start + 1);
3452
}
3453
else
3454
{
3455
fprintf (stderr, "Unknown characteristic\n");
3456
return 0;
3457
}
3458
}
3459
else
3460
{
3461
if (strcmp(characteristic, "average") == 0 || strcmp(characteristic, "avg") == 0)
3462
{
3463
for (i = time_start; i <= time_end; i ++)
3464
{
3465
avg_x += bp_value_to_double(adios_double, vix->avgs[i]) / (time_end - time_start + 1);
3466
avg_y += bp_value_to_double(adios_double, viy->avgs[i]) / (time_end - time_start + 1);
3467
}
3468
for (i = time_start; i <= time_end; i ++)
3469
{
3470
cov += (bp_value_to_double(adios_double, vix->avgs[i]) - avg_x) * (bp_value_to_double(adios_double, viy->avgs[i]) - avg_y) / (time_end - time_start + 1);
3471
}
3472
}
3473
else if (strcmp(characteristic, "standard deviation") == 0 || strcmp(characteristic, "std_dev") == 0)
3474
{
3475
for (i = time_start; i <= time_end; i ++)
3476
{
3477
avg_x += bp_value_to_double(adios_double, vix->std_devs[i]) / (time_end - time_start + 1);
3478
avg_y += bp_value_to_double(adios_double, viy->std_devs[i]) / (time_end - time_start + 1);
3479
}
3480
for (i = time_start; i <= time_end; i ++)
3481
{
3482
cov += (bp_value_to_double(adios_double, vix->std_devs[i]) - avg_x) * (bp_value_to_double(adios_double, viy->std_devs[i]) - avg_y) / (time_end - time_start + 1);
3483
}
3484
}
3485
else if (strcmp(characteristic, "minimum") == 0 || strcmp(characteristic, "min") == 0)
3486
{
3487
for (i = time_start; i <= time_end; i ++)
3488
{
3489
avg_x += bp_value_to_double(vix->type, vix->mins[i]) / (time_end - time_start + 1);
3490
avg_y += bp_value_to_double(viy->type, viy->mins[i]) / (time_end - time_start + 1);
3491
}
3492
for (i = time_start; i <= time_end; i ++)
3493
{
3494
cov += (bp_value_to_double(vix->type, vix->mins[i]) - avg_x) * (bp_value_to_double(viy->type, viy->mins[i]) - avg_y) / (time_end - time_start + 1);
3495
}
3496
}
3497
else if (strcmp(characteristic, "maximum") == 0 || strcmp(characteristic, "max") == 0)
3498
{
3499
for (i = time_start; i <= time_end; i ++)
3500
{
3501
avg_x += bp_value_to_double(vix->type, vix->maxs[i]) / (time_end - time_start + 1);
3502
avg_y += bp_value_to_double(vix->type, viy->maxs[i]) / (time_end - time_start + 1);
3503
}
3504
for (i = time_start; i <= time_end; i ++)
3505
{
3506
cov += (bp_value_to_double(vix->type, vix->maxs[i]) - avg_x) * (bp_value_to_double(viy->type, viy->maxs[i]) - avg_y) / (time_end - time_start + 1);
3507
}
3508
}
3509
else
3510
{
3511
fprintf (stderr, "Unknown characteristic\n");
3512
return 0;
3513
}
3514
}
3515
}
3516
else
3517
{
3518
fprintf (stderr, "Time values out of bounds\n");
3519
}
3520
return cov;
3521
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1