← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/io/Checkpoint.cpp

  • Committer: Niclas Jansson
  • Date: 2011-06-10 14:33:43 UTC
  • Revision ID: njansson@csc.kth.se-20110610143343-d21p4am8rghiojfm
Added rudimentary header to binary files

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/io/Checkpoint.cpp
 
1
// Copyright (C) 2009 Niclas Jansson.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2009-09-09
 
5
// Last changed: 2011-04-19
 
6
 
 
7
#include <cstring>
 
8
#include <sstream>
 
9
#include <fstream>
 
10
#include <dolfin/la/Vector.h>
 
11
#include <dolfin/mesh/Mesh.h>
 
12
#include <dolfin/mesh/MeshEditor.h>
 
13
#include <dolfin/mesh/Vertex.h>
 
14
#include <dolfin/function/Function.h>
 
15
#include <dolfin/io/Checkpoint.h>
 
16
 
 
17
#ifdef ENABLE_MPIIO
 
18
#include <mpi.h>
 
19
#endif
 
20
 
 
21
 
 
22
using namespace dolfin;
 
23
//-----------------------------------------------------------------------------
 
24
Checkpoint::Checkpoint() : state(CHECKPOINT), restart_state(OPEN), n(0), 
 
25
                           hdr_initialized(false), disp_initialized(false)
 
26
{
 
27
}
 
28
//-----------------------------------------------------------------------------
 
29
Checkpoint::~Checkpoint()
 
30
{
 
31
}
 
32
//-----------------------------------------------------------------------------
 
33
void Checkpoint::hdr_init(Mesh& mesh, bool static_mesh)
 
34
{
 
35
 
 
36
  if(!hdr_initialized || !static_mesh)
 
37
  {    
 
38
    hdr.num_coords = mesh.numVertices() * mesh.geometry().dim();
 
39
    hdr.num_entities = mesh.type().numEntities(0);
 
40
    hdr.num_centities = mesh.numCells() * hdr.num_entities;
 
41
    hdr.type = mesh.type().cellType();
 
42
    hdr.tdim = mesh.topology().dim();
 
43
    hdr.gdim = mesh.geometry().dim();
 
44
    hdr.num_vertices = mesh.numVertices();
 
45
    hdr.num_cells = mesh.numCells();  
 
46
    hdr.num_ghosts = mesh.distdata().num_ghost(0);
 
47
    hdr.num_shared = mesh.distdata().num_shared(0);
 
48
  
 
49
    
 
50
#ifdef ENABLE_MPIIO
 
51
    uint local_data[5] = {hdr.num_coords, 
 
52
                          hdr.num_centities,
 
53
                          hdr.num_vertices,
 
54
                          (2*hdr.num_ghosts),
 
55
                          hdr.num_shared};
 
56
 
 
57
    memset(&hdr.offsets[0], 0, 5 * sizeof(uint));
 
58
#if ( MPI_VERSION > 1 )
 
59
    MPI_Exscan(&local_data[0], &hdr.offsets[0], 5,
 
60
               MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
61
#else
 
62
    MPI_Scan(&local_data[0], &hdr.offsets[0], 5, 
 
63
             MPI_UNSIGNED, MPI_SUM, MPI::DOLFIN_COMM);
 
64
    
 
65
    hdr.offsets[0] -= local_data[0];
 
66
    hdr.offsets[1] -= local_data[1];
 
67
    hdr.offsets[2] -= local_data[2];
 
68
    hdr.offsets[3] -= local_data[3];
 
69
    hdr.offsets[4] -= local_data[4];
 
70
#endif
 
71
    
 
72
    if (!disp_initialized || !static_mesh) 
 
73
    {
 
74
      memset(&hdr.disp[0], 0, 5 * sizeof(uint));
 
75
      MPI_Allreduce(&local_data[0], &hdr.disp[0], 5, MPI_UNSIGNED,
 
76
                    MPI_SUM, MPI::DOLFIN_COMM);
 
77
      disp_initialized = true;
 
78
    }
 
79
 
 
80
#endif
 
81
  }
 
82
  
 
83
  hdr_initialized = true;
 
84
 
 
85
}
 
86
//-----------------------------------------------------------------------------
 
87
void Checkpoint::write(std::string fname, uint id, real t, Mesh& mesh, 
 
88
                       std::vector<Function *> func,
 
89
                       std::vector<Vector *> vec, bool static_mesh)
 
90
{
 
91
 
 
92
  message("Writing checkpoint (%s%d) at time %g", fname.c_str(), n%2, t);
 
93
  std::ostringstream _fname;
 
94
#ifndef ENABLE_MPIIO
 
95
  if( MPI::numProcesses() > 1) 
 
96
    _fname << fname << (n++)%2 << "_" <<  MPI::processNumber() << ".chkp";
 
97
  else
 
98
    _fname << fname << (n++)%2 << ".chkp";
 
99
 
 
100
  std::ofstream out(_fname.str().c_str(), std::ofstream::binary);
 
101
 
 
102
  out.write((char *) &id, sizeof(uint));
 
103
  out.write((char *) &t, sizeof(real));
 
104
 
 
105
#else
 
106
  _fname << fname << (n++)%2 << ".chkp";
 
107
    
 
108
  MPI_File out;
 
109
  MPI_File_open(dolfin::MPI::DOLFIN_COMM, (char *) _fname.str().c_str(),
 
110
                MPI_MODE_WRONLY | MPI_MODE_CREATE, MPI_INFO_NULL, &out);
 
111
  MPI_File_set_size(out, 0);
 
112
  
 
113
  byte_offset = 0;
 
114
  MPI_File_write_all(out, &id, 1, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
115
  MPI_File_write_all(out, &t, 1, MPI_DOUBLE, MPI_STATUS_IGNORE);
 
116
 
 
117
  byte_offset += sizeof(uint);
 
118
  byte_offset += sizeof(real);
 
119
#endif
 
120
 
 
121
  hdr_init(mesh, static_mesh);
 
122
  write(mesh, out);
 
123
  write(func, out);
 
124
  write(vec, out);
 
125
  
 
126
#ifdef ENABLE_MPIIO
 
127
  MPI_File_close(&out);
 
128
#else
 
129
  out.close();
 
130
#endif
 
131
 
 
132
 
 
133
  state = RESTART;
 
134
 
 
135
}
 
136
//-----------------------------------------------------------------------------
 
137
void Checkpoint::restart(std::string fname)
 
138
{
 
139
 
 
140
  if (restart_state != OPEN)
 
141
    error("Shut her down, Scotty, she's sucking mud again!");
 
142
 
 
143
  std::ostringstream _fname;
 
144
#ifndef ENABLE_MPIIO
 
145
  if( MPI::numProcesses() > 1) 
 
146
    _fname << fname << "_" << MPI::processNumber() << ".chkp";
 
147
  else
 
148
    _fname << fname << ".chkp";
 
149
#else
 
150
  _fname << fname << ".chkp";
 
151
#endif
 
152
 
 
153
#ifdef ENABLE_MPIIO
 
154
  MPI_File_open(dolfin::MPI::DOLFIN_COMM, (char *)  _fname.str().c_str(),
 
155
                MPI_MODE_RDONLY, MPI_INFO_NULL, &in);
 
156
  MPI_File_read_all(in, &_id, 1, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
157
  MPI_File_read_all(in, &_t, 1, MPI_DOUBLE, MPI_STATUS_IGNORE);
 
158
  byte_offset = sizeof(uint) + sizeof(real);
 
159
#else
 
160
  in.open(_fname.str().c_str(), std::ifstream::binary);
 
161
  in.read((char *) &_id, sizeof(uint));
 
162
  in.read((char *) &_t, sizeof(real));
 
163
#endif
 
164
  message("Restarting from time %g checkpoint id %d", _t, _id);
 
165
 
 
166
  state = RESTART;
 
167
  restart_state = MESH;
 
168
      
 
169
}
 
170
//-----------------------------------------------------------------------------
 
171
void Checkpoint::load(Mesh& mesh)
 
172
{ 
 
173
  if (restart_state != MESH)
 
174
    error("Shut her down, Scotty, she's sucking mud again!");
 
175
 
 
176
#ifdef ENABLE_MPIIO
 
177
  uint pe_size = MPI::numProcesses();
 
178
  uint pe_rank = MPI::processNumber();
 
179
 
 
180
  MPI_File_read_at_all(in, byte_offset + pe_rank * sizeof(chkp_mesh_hdr),
 
181
                       &hdr, sizeof(chkp_mesh_hdr), MPI_BYTE, MPI_STATUS_IGNORE);
 
182
  byte_offset += pe_size * sizeof(chkp_mesh_hdr);
 
183
#else
 
184
  in.read((char *)&hdr, sizeof(chkp_mesh_hdr));
 
185
#endif  
 
186
  real *coords = new real[hdr.num_coords];
 
187
 
 
188
#ifdef ENABLE_MPIIO
 
189
  MPI_File_read_at_all(in, byte_offset + hdr.offsets[0] * sizeof(real),
 
190
                       coords, hdr.num_coords, MPI_DOUBLE, MPI_STATUS_IGNORE);
 
191
  byte_offset += hdr.disp[0] * sizeof(real);
 
192
#else
 
193
  in.read((char *)coords, (hdr.num_coords) * sizeof(real));
 
194
#endif
 
195
 
 
196
  Mesh _mesh;
 
197
  MeshEditor editor;  
 
198
  editor.open(_mesh, hdr.type, hdr.tdim, hdr.gdim);
 
199
  editor.initVertices(hdr.num_vertices);
 
200
 
 
201
  uint vi = 0;
 
202
  for(uint i = 0 ; i < hdr.num_coords; i += hdr.gdim)
 
203
  {
 
204
    switch(hdr.gdim)
 
205
    {      
 
206
    case 2:
 
207
      editor.addVertex(vi++, coords[i], coords[i+1]); break;
 
208
    case 3:
 
209
      editor.addVertex(vi++, coords[i], coords[i+1], coords[i+2]); break;
 
210
    }
 
211
  }
 
212
  
 
213
  delete[] coords;
 
214
 
 
215
  editor.initCells(hdr.num_cells);
 
216
 
 
217
  uint *cells = new uint[hdr.num_centities];
 
218
#ifdef ENABLE_MPIIO
 
219
  MPI_File_read_at_all(in, byte_offset + hdr.offsets[1] * sizeof(uint), 
 
220
                       cells, hdr.num_centities, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
221
  byte_offset += hdr.disp[1] * sizeof(uint);
 
222
#else
 
223
  in.read((char *)cells, (hdr.num_centities) * sizeof(uint));
 
224
#endif
 
225
 
 
226
  Array<uint> v;
 
227
  uint ci = 0;
 
228
  for (uint i = 0; i < hdr.num_centities; i += hdr.num_entities)
 
229
  {
 
230
    v.clear();
 
231
    for (uint j = 0; j < hdr.num_entities; j++)
 
232
      v.push_back(cells[i + j]);
 
233
    editor.addCell(ci++, v);    
 
234
  }  
 
235
  editor.close();
 
236
  delete[] cells;
 
237
 
 
238
  if (MPI::numProcesses() > 1) 
 
239
  {
 
240
    uint *mapping = new uint[_mesh.numVertices()];
 
241
#ifdef ENABLE_MPIIO
 
242
    MPI_File_read_at_all(in, byte_offset + hdr.offsets[2] * sizeof(uint),
 
243
                         mapping, hdr.num_vertices, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
244
    byte_offset += hdr.disp[2] * sizeof(uint);
 
245
#else
 
246
    in.read((char *)mapping, hdr.num_vertices * sizeof(uint));
 
247
#endif
 
248
    for (VertexIterator v(_mesh); !v.end(); ++v)
 
249
      _mesh.distdata().set_map(v->index(), mapping[v->index()], 0);
 
250
    delete[] mapping;
 
251
 
 
252
    uint *ghosts = new uint[2 * hdr.num_ghosts];
 
253
#ifdef ENABLE_MPIIO
 
254
    MPI_File_read_at_all(in, byte_offset + hdr.offsets[3] * sizeof(uint),
 
255
                         ghosts, 2*hdr.num_ghosts, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
256
    byte_offset += hdr.disp[3] * sizeof(uint);
 
257
#else
 
258
    in.read((char *)ghosts, 2*hdr.num_ghosts * sizeof(uint));
 
259
#endif
 
260
    for (uint i = 0; i < 2 * hdr.num_ghosts; i += 2)
 
261
    {
 
262
      _mesh.distdata().set_ghost(ghosts[i], 0);
 
263
      _mesh.distdata().set_ghost_owner(ghosts[i], ghosts[i+1], 0);
 
264
    }
 
265
    delete[] ghosts;
 
266
    
 
267
    uint *shared = new uint[hdr.num_shared];
 
268
#ifdef ENABLE_MPIIO
 
269
    MPI_File_read_at_all(in, byte_offset + hdr.offsets[4] * sizeof(uint),
 
270
                         shared, hdr.num_shared, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
271
    byte_offset += hdr.disp[4] * sizeof(uint);
 
272
#else
 
273
    in.read((char *)shared, hdr.num_shared * sizeof(uint));
 
274
#endif
 
275
    for(uint i = 0; i < hdr.num_shared; i++)
 
276
      _mesh.distdata().set_shared(shared[i], 0);
 
277
    delete[] shared;    
 
278
  }
 
279
 
 
280
  mesh = _mesh;
 
281
  mesh.distdata().invalid_numbering();
 
282
  mesh.renumber();
 
283
 
 
284
  restart_state = FUNC;
 
285
 
 
286
}
 
287
//-----------------------------------------------------------------------------
 
288
void Checkpoint::load(std::vector<Function *> func)
 
289
{
 
290
  if (restart_state != FUNC)
 
291
    error("Shut her down, Scotty, she's sucking mud again!");
 
292
 
 
293
  std::vector<Function *>::iterator it;
 
294
  uint local_size;
 
295
  uint pe_rank = MPI::processNumber();
 
296
  uint pe_size = MPI::numProcesses();
 
297
  uint vector_offset[3];
 
298
  // FIXME store max(local_size) 
 
299
  for (it = func.begin(); it != func.end(); ++it)
 
300
  {
 
301
#ifdef ENABLE_MPIIO
 
302
    MPI_File_read_at_all(in, byte_offset + pe_rank * 3 * sizeof(uint),
 
303
                         &vector_offset[0], 3, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
304
    byte_offset += pe_size * 3 * sizeof(uint);
 
305
    local_size = vector_offset[1];   
 
306
#else
 
307
    in.read((char *)&local_size, sizeof(uint));
 
308
#endif
 
309
    real *values = new real[local_size];
 
310
#ifdef ENABLE_MPIIO
 
311
    MPI_File_read_at_all(in, byte_offset + vector_offset[0] * sizeof(real),
 
312
                         values, vector_offset[1], MPI_DOUBLE, MPI_STATUS_IGNORE);
 
313
    byte_offset += vector_offset[2] * sizeof(real);
 
314
#else
 
315
    in.read((char *)values, local_size * sizeof(real));    
 
316
#endif
 
317
 
 
318
    (*it)->vector().set(values);
 
319
    (*it)->vector().apply();
 
320
 
 
321
    delete[] values;
 
322
  }
 
323
 
 
324
  restart_state = VEC;
 
325
}
 
326
//-----------------------------------------------------------------------------
 
327
void Checkpoint::load(std::vector<Vector *> vec)
 
328
{
 
329
  if (restart_state != VEC)
 
330
    error("Shut her down, Scotty, she's sucking mud again!");
 
331
  
 
332
  std::vector<Vector *>::iterator it;
 
333
  uint local_size;
 
334
  uint pe_rank = MPI::processNumber();
 
335
  uint pe_size = MPI::numProcesses();
 
336
  uint vector_offset[3];
 
337
  for (it = vec.begin(); it != vec.end(); ++it)
 
338
  {
 
339
#ifdef ENABLE_MPIIO
 
340
    MPI_File_read_at_all(in, byte_offset + pe_rank * 3 * sizeof(uint),
 
341
                         &vector_offset[0], 3, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
342
    byte_offset += pe_size * 3 * sizeof(uint);
 
343
    local_size = vector_offset[1];
 
344
#else
 
345
    in.read((char *)&local_size, sizeof(uint));
 
346
#endif
 
347
    real *values = new real[local_size];
 
348
#ifdef ENABLE_MPIIO
 
349
    MPI_File_read_at_all(in, byte_offset + vector_offset[0] * sizeof(real),
 
350
                         values, vector_offset[1], MPI_DOUBLE, MPI_STATUS_IGNORE);
 
351
    byte_offset += vector_offset[3] * sizeof(double);
 
352
#else
 
353
    in.read((char *)values, local_size * sizeof(real));
 
354
#endif
 
355
    (*it)->set(values);
 
356
    (*it)->apply();
 
357
    delete[] values;
 
358
  }
 
359
  
 
360
#ifdef ENABLE_MPIIO
 
361
  MPI_File_close(&in);
 
362
#else
 
363
  in.close();
 
364
#endif
 
365
}
 
366
//-----------------------------------------------------------------------------
 
367
void Checkpoint::write(Mesh& mesh, chkp_outstream& out)
 
368
{
 
369
 
 
370
#ifdef ENABLE_MPIIO
 
371
  uint pe_size = MPI::numProcesses();
 
372
  uint pe_rank = MPI::processNumber();
 
373
  
 
374
  MPI_File_write_at_all(out, byte_offset + pe_rank * sizeof(chkp_mesh_hdr),
 
375
                        &hdr, sizeof(chkp_mesh_hdr), MPI_BYTE, MPI_STATUS_IGNORE);
 
376
  byte_offset += pe_size * sizeof(chkp_mesh_hdr);
 
377
  
 
378
  MPI_File_write_at_all(out, byte_offset + hdr.offsets[0] * sizeof(real),
 
379
                        mesh.coordinates(), hdr.num_coords, MPI_DOUBLE, MPI_STATUS_IGNORE);
 
380
  byte_offset += hdr.disp[0] * sizeof(real);
 
381
  
 
382
  MPI_File_write_at_all(out, byte_offset + hdr.offsets[1] * sizeof(uint),
 
383
                        mesh.cells(), hdr.num_centities, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
384
  byte_offset += hdr.disp[1] * sizeof(uint);
 
385
 
 
386
#else
 
387
  out.write((char *)&hdr, sizeof(chkp_mesh_hdr));
 
388
  out.write((char *)mesh.coordinates(), hdr.num_coords * sizeof(real));
 
389
  out.write((char *)mesh.cells(), hdr.num_centities * sizeof(uint));
 
390
#endif
 
391
 
 
392
  if (MPI::numProcesses() > 1) 
 
393
  {
 
394
    uint *mapping = new uint[mesh.numVertices()];
 
395
    for (VertexIterator v(mesh); !v.end(); ++v)
 
396
      mapping[v->index()] = mesh.distdata().get_global(*v);    
 
397
#ifdef ENABLE_MPIIO
 
398
    MPI_File_write_at_all(out, byte_offset + hdr.offsets[2] * sizeof(uint),
 
399
                          mapping, hdr.num_vertices, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
400
    byte_offset += hdr.disp[2] * sizeof(uint);   
 
401
#else
 
402
    out.write((char *)mapping, hdr.num_vertices * sizeof(uint));
 
403
#endif
 
404
    delete[] mapping;
 
405
 
 
406
    uint *ghosts = new uint[2 * hdr.num_ghosts];
 
407
    uint *gp = &ghosts[0];
 
408
    for (MeshGhostIterator g(mesh.distdata(), 0); !g.end(); ++g)
 
409
    {
 
410
      *gp++ = g.index();
 
411
      *gp++ = g.owner();
 
412
    }
 
413
#ifdef ENABLE_MPIIO
 
414
    MPI_File_write_at_all(out, byte_offset + hdr.offsets[3] * sizeof(uint),
 
415
                          ghosts, 2 * hdr.num_ghosts, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
416
    byte_offset += hdr.disp[3] * sizeof(uint);
 
417
#else
 
418
    out.write((char *)ghosts, 2 * hdr.num_ghosts * sizeof(uint));
 
419
#endif
 
420
    delete[] ghosts;
 
421
 
 
422
    uint *shared = new uint[hdr.num_shared];
 
423
    uint *sp = &shared[0];
 
424
    for (MeshSharedIterator s(mesh.distdata(), 0); !s.end(); ++s)
 
425
      *sp++ = s.index();
 
426
#ifdef ENABLE_MPIIO
 
427
    MPI_File_write_at_all(out, byte_offset + hdr.offsets[4] * sizeof(uint),
 
428
                          shared, hdr.num_shared, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
429
    byte_offset += hdr.disp[4] * sizeof(uint);
 
430
#else
 
431
    out.write((char *)shared, hdr.num_shared * sizeof(uint));
 
432
#endif
 
433
    delete[] shared;
 
434
 
 
435
  }
 
436
 
 
437
}
 
438
//-----------------------------------------------------------------------------
 
439
void Checkpoint::write(std::vector<Function *> func, chkp_outstream& out)
 
440
{
 
441
  std::vector<Function *>::iterator it;
 
442
 
 
443
  uint max_size = 0;
 
444
  for (it = func.begin(); it != func.end(); ++it)
 
445
  {
 
446
    if((*it)->type() != Function::discrete)
 
447
      error("Checkpoint restart only implemented for discrete functions");    
 
448
    max_size = std::max(max_size, (*it)->vector().local_size());
 
449
  }
 
450
  
 
451
  real *values = new real[max_size];
 
452
#ifdef ENABLE_MPIIO
 
453
  uint vector_offset[3];
 
454
  uint pe_size = MPI::numProcesses();
 
455
  uint pe_rank = MPI::processNumber();
 
456
#endif
 
457
 
 
458
  for (it = func.begin(); it != func.end(); ++it)
 
459
  {
 
460
    uint local_size = (*it)->vector().local_size();
 
461
    (*it)->vector().get(values);
 
462
 
 
463
#ifdef ENABLE_MPIIO
 
464
    vector_offset[0] = (*it)->vector().offset();
 
465
    vector_offset[1] = (*it)->vector().local_size();
 
466
    vector_offset[2] = (*it)->vector().size();
 
467
    
 
468
    MPI_File_write_at_all(out, byte_offset + pe_rank * 3 * sizeof(uint),
 
469
                          &vector_offset[0], 3, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
470
    byte_offset += pe_size * 3 * sizeof(uint);
 
471
 
 
472
    MPI_File_write_at_all(out, byte_offset + vector_offset[0] * sizeof(real),
 
473
                          values, vector_offset[1], MPI_DOUBLE, MPI_STATUS_IGNORE);
 
474
    byte_offset += (*it)->vector().size() * sizeof(real);
 
475
#else
 
476
    out.write((char *)&local_size, sizeof(uint));
 
477
    out.write((char *)values, (*it)->vector().local_size() * sizeof(real));
 
478
#endif
 
479
  }
 
480
  delete[] values;
 
481
  
 
482
}
 
483
//-----------------------------------------------------------------------------
 
484
void Checkpoint::write(std::vector<Vector *> vec, chkp_outstream& out)
 
485
{
 
486
  std::vector<Vector *>::iterator it;
 
487
 
 
488
  uint max_size = 0;
 
489
  for (it = vec.begin(); it != vec.end(); ++it)
 
490
    max_size = std::max(max_size, (*it)->local_size());
 
491
  
 
492
  real *values = new real[max_size];
 
493
#ifdef ENABLE_MPIIO
 
494
  uint vector_offset[3];
 
495
  uint pe_size = MPI::numProcesses();
 
496
  uint pe_rank = MPI::processNumber();
 
497
#endif
 
498
  for (it = vec.begin(); it != vec.end(); ++it)
 
499
  {
 
500
 
 
501
    (*it)->get(values);
 
502
 
 
503
#ifdef ENABLE_MPIIO
 
504
    vector_offset[0] = (*it)->offset();
 
505
    vector_offset[1] = (*it)->local_size();
 
506
    vector_offset[2] = (*it)->size();
 
507
    
 
508
    MPI_File_write_at_all(out, byte_offset + pe_rank * 3 * sizeof(uint), 
 
509
                          &vector_offset[0], 3, MPI_UNSIGNED, MPI_STATUS_IGNORE);
 
510
    byte_offset += pe_size * 3 * sizeof(uint);
 
511
 
 
512
    MPI_File_write_at_all(out, byte_offset + vector_offset[0] * sizeof(real),
 
513
                          values, vector_offset[1], MPI_DOUBLE, MPI_STATUS_IGNORE);
 
514
    byte_offset += vector_offset[2] * sizeof(real);
 
515
#else
 
516
    uint local_size = (*it)->local_size();
 
517
    out.write((char *)&local_size, sizeof(uint));
 
518
    out.write((char *)values, (*it)->local_size() * sizeof(real));
 
519
#endif
 
520
  }
 
521
  delete[] values;
 
522
  
 
523
}
 
524
//-----------------------------------------------------------------------------