/[escript]/branches/arrexp_2137_win_merge/escript/src/DataExpanded.cpp
ViewVC logotype

Contents of /branches/arrexp_2137_win_merge/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2221 - (show annotations)
Mon Jan 19 06:11:25 2009 UTC (10 years, 5 months ago) by jfenwick
File size: 27263 byte(s)
A second attempt at a thread-safe COW

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "Data.h"
16 #include "DataExpanded.h"
17 #include "DataException.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #ifdef USE_NETCDF
21 #include <netcdfcpp.h>
22 #endif
23 #ifdef PASO_MPI
24 #include <mpi.h>
25 #endif
26
27 #include <boost/python/extract.hpp>
28 #include "DataMaths.h"
29
30 using namespace std;
31 using namespace boost::python;
32 using namespace boost;
33 using namespace escript::DataTypes;
34
35
36 // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
37
38 #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; *((int*)0)=17;throw DataException(ss.str());}
39
40 namespace escript {
41
42 DataExpanded::DataExpanded(const WrappedArray& value,
43 const FunctionSpace& what)
44 : parent(what,value.getShape())
45 {
46 //
47 // initialise the data array for this object
48 initialise(what.getNumSamples(),what.getNumDPPSample());
49 //
50 // copy the given value to every data point
51 copy(value);
52 }
53
54 DataExpanded::DataExpanded(const DataExpanded& other)
55 : parent(other.getFunctionSpace(), other.getShape()),
56 m_data(other.m_data)
57 {
58 }
59
60 DataExpanded::DataExpanded(const DataConstant& other)
61 : parent(other.getFunctionSpace(), other.getShape())
62 {
63 //
64 // initialise the data array for this object
65 initialise(other.getNumSamples(),other.getNumDPPSample());
66 //
67 // DataConstant only has one value, copy this to every data point
68 copy(other);
69 }
70
71 DataExpanded::DataExpanded(const DataTagged& other)
72 : parent(other.getFunctionSpace(), other.getShape())
73 {
74 //
75 // initialise the data array for this object
76 initialise(other.getNumSamples(),other.getNumDPPSample());
77 //
78 // for each data point in this object, extract and copy the corresponding data
79 // value from the given DataTagged object
80 int i,j;
81 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
82 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
83 #pragma omp parallel for private(i,j) schedule(static)
84 for (i=0;i<numRows;i++) {
85 for (j=0;j<numCols;j++) {
86 try {
87 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(),
88 other.getVector(),
89 other.getPointOffset(i,j));
90 }
91 catch (std::exception& e) {
92 cout << e.what() << endl;
93 }
94 }
95 }
96 }
97
98 DataExpanded::DataExpanded(const DataExpanded& other,
99 const DataTypes::RegionType& region)
100 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101 {
102 //
103 // get the shape of the slice
104 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
105 //
106 // initialise this Data object to the shape of the slice
107 initialise(other.getNumSamples(),other.getNumDPPSample());
108 //
109 // copy the data
110 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
111 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
112 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
113 int i,j;
114 #pragma omp parallel for private(i,j) schedule(static)
115 for (i=0;i<numRows;i++) {
116 for (j=0;j<numCols;j++) {
117 try {
118 DataTypes::copySlice(getVector(),getShape(),getPointOffset(i,j),
119 other.getVector(),
120 other.getShape(),
121 other.getPointOffset(i,j),
122 region_loop_range);
123 }
124 catch (std::exception& e) {
125 cout << e.what() << endl;
126 }
127 }
128 }
129 }
130
131 DataExpanded::DataExpanded(const FunctionSpace& what,
132 const DataTypes::ShapeType &shape,
133 const DataTypes::ValueType &data)
134 : parent(what,shape)
135 {
136 EsysAssert(data.size()%getNoValues()==0,
137 "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
138
139 if (data.size()==getNoValues())
140 {
141 ValueType& vec=m_data.getData();
142 //
143 // create the view of the data
144 initialise(what.getNumSamples(),what.getNumDPPSample());
145 // now we copy this value to all elements
146 for (int i=0;i<getLength();)
147 {
148 for (unsigned int j=0;j<getNoValues();++j,++i)
149 {
150 vec[i]=data[j];
151 }
152 }
153 }
154 else
155 {
156 //
157 // copy the data in the correct format
158 m_data.getData()=data;
159 }
160
161
162 }
163
164 DataExpanded::~DataExpanded()
165 {
166 }
167
168 DataAbstract*
169 DataExpanded::deepCopy()
170 {
171 return new DataExpanded(*this);
172 }
173
174
175 DataAbstract*
176 DataExpanded::getSlice(const DataTypes::RegionType& region) const
177 {
178 return new DataExpanded(*this,region);
179 }
180
181 void
182 DataExpanded::setSlice(const DataAbstract* value,
183 const DataTypes::RegionType& region)
184 {
185 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
186 if (tempDataExp==0) {
187 throw DataException("Programming error - casting to DataExpanded.");
188 }
189 CHECK_FOR_EX_WRITE
190 //
191 // get shape of slice
192 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194 //
195 // check shape
196 if (getRank()!=region.size()) {
197 throw DataException("Error - Invalid slice region.");
198 }
199 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
200 throw DataException (DataTypes::createShapeErrorMessage(
201 "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
202 }
203 //
204 // copy the data from the slice into this object
205 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
206 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
207 int i, j;
208 ValueType& vec=getVector();
209 const ShapeType& mshape=getShape();
210 const ValueType& tVec=tempDataExp->getVector();
211 const ShapeType& tShape=tempDataExp->getShape();
212 #pragma omp parallel for private(i,j) schedule(static)
213 for (i=0;i<numRows;i++) {
214 for (j=0;j<numCols;j++) {
215 /* getPointDataView().copySliceFrom(getPointOffset(i,j),
216 tempDataExp->getPointDataView(),
217 tempDataExp->getPointOffset(i,j),
218 region_loop_range);*/
219 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
220 tVec,
221 tShape,
222 tempDataExp->getPointOffset(i,j),
223 region_loop_range);
224
225 }
226 }
227 }
228
229 void
230 DataExpanded::copy(const DataConstant& value)
231 {
232 EsysAssert((checkShape(getShape(), value.getShape())),
233 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
234
235 //
236 // copy a single value to every data point in this object
237 int nRows=m_data.getNumRows();
238 int nCols=m_data.getNumCols();
239 int i,j;
240 #pragma omp parallel for private(i,j) schedule(static)
241 for (i=0;i<nRows;i++) {
242 for (j=0;j<nCols;j++) {
243 // NOTE: An exception may be thown from this call if
244 // DOASSERT is on which of course will play
245 // havoc with the omp threads. Run single threaded
246 // if using DOASSERT.
247 //getPointDataView().copy(getPointOffset(i,j),value);
248 DataTypes::copyPoint(getVector(), getPointOffset(i,j), getNoValues(), value.getVector(), 0);
249 }
250 }
251 }
252
253 void
254 DataExpanded::copy(const WrappedArray& value)
255 {
256 // check the input shape matches this shape
257 if (!DataTypes::checkShape(getShape(),value.getShape())) {
258 throw DataException(DataTypes::createShapeErrorMessage(
259 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
260 value.getShape(),getShape()));
261 }
262 getVector().copyFromArray(value, getNumDPPSample()*getNumSamples());
263 }
264
265
266 void
267 DataExpanded::initialise(int noSamples,
268 int noDataPointsPerSample)
269 {
270 if (noSamples==0) //retain the default empty object
271 {
272 return;
273 }
274 //
275 // resize data array to the required size
276 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
277 }
278
279 string
280 DataExpanded::toString() const
281 {
282 stringstream temp;
283 FunctionSpace fs=getFunctionSpace();
284
285 int offset=0;
286 for (int i=0;i<m_data.getNumRows();i++) {
287 for (int j=0;j<m_data.getNumCols();j++) {
288 offset=getPointOffset(i,j);
289 stringstream suffix;
290 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
291 temp << DataTypes::pointToString(getVector(),getShape(),offset,suffix.str());
292 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
293 temp << endl;
294 }
295 }
296 }
297 string result=temp.str();
298 if (result.empty())
299 {
300 return "(data contains no samples)\n";
301 }
302 return temp.str();
303 }
304
305 DataTypes::ValueType::size_type
306 DataExpanded::getPointOffset(int sampleNo,
307 int dataPointNo) const
308 {
309 return m_data.index(sampleNo,dataPointNo);
310 }
311
312 DataTypes::ValueType::size_type
313 DataExpanded::getPointOffset(int sampleNo,
314 int dataPointNo)
315 {
316 return m_data.index(sampleNo,dataPointNo);
317 }
318
319 DataTypes::ValueType::size_type
320 DataExpanded::getLength() const
321 {
322 return m_data.size();
323 }
324
325
326
327 void
328 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
329 CHECK_FOR_EX_WRITE
330 //
331 // Get the number of samples and data-points per sample.
332 int numSamples = getNumSamples();
333 int numDataPointsPerSample = getNumDPPSample();
334 int dataPointRank = getRank();
335 ShapeType dataPointShape = getShape();
336 if (numSamples*numDataPointsPerSample > 0) {
337 //TODO: global error handling
338 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
339 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
340 }
341 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
342 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
343 }
344 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
345 ValueType& vec=getVector();
346 if (dataPointRank==0) {
347 vec[offset]=value;
348 } else if (dataPointRank==1) {
349 for (int i=0; i<dataPointShape[0]; i++) {
350 vec[offset+i]=value;
351 }
352 } else if (dataPointRank==2) {
353 for (int i=0; i<dataPointShape[0]; i++) {
354 for (int j=0; j<dataPointShape[1]; j++) {
355 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
356 }
357 }
358 } else if (dataPointRank==3) {
359 for (int i=0; i<dataPointShape[0]; i++) {
360 for (int j=0; j<dataPointShape[1]; j++) {
361 for (int k=0; k<dataPointShape[2]; k++) {
362 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
363 }
364 }
365 }
366 } else if (dataPointRank==4) {
367 for (int i=0; i<dataPointShape[0]; i++) {
368 for (int j=0; j<dataPointShape[1]; j++) {
369 for (int k=0; k<dataPointShape[2]; k++) {
370 for (int l=0; l<dataPointShape[3]; l++) {
371 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
372 }
373 }
374 }
375 }
376 }
377 }
378 }
379
380 void
381 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
382 CHECK_FOR_EX_WRITE
383 //
384 // Get the number of samples and data-points per sample.
385 int numSamples = getNumSamples();
386 int numDataPointsPerSample = getNumDPPSample();
387 //
388 // check rank:
389 if (value.getRank()!=getRank())
390 throw DataException("Rank of numarray does not match Data object rank");
391 if (numSamples*numDataPointsPerSample > 0) {
392 //TODO: global error handling
393 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
394 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
395 }
396 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
397 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
398 }
399 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
400 ValueType& vec=getVector();
401 vec.copyFromArrayToOffset(value,offset,1);
402 }
403 }
404
405 void
406 DataExpanded::symmetric(DataAbstract* ev)
407 {
408 int sampleNo,dataPointNo;
409 int numSamples = getNumSamples();
410 int numDataPointsPerSample = getNumDPPSample();
411 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
412 if (temp_ev==0) {
413 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
414 }
415 ValueType& vec=getVector();
416 const ShapeType& shape=getShape();
417 ValueType& evVec=temp_ev->getVector();
418 const ShapeType& evShape=temp_ev->getShape();
419 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
420 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
421 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
422 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
423 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
424 }
425 }
426 }
427 void
428 DataExpanded::nonsymmetric(DataAbstract* ev)
429 {
430 int sampleNo,dataPointNo;
431 int numSamples = getNumSamples();
432 int numDataPointsPerSample = getNumDPPSample();
433 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
434 if (temp_ev==0) {
435 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
436 }
437 ValueType& vec=getVector();
438 const ShapeType& shape=getShape();
439 ValueType& evVec=temp_ev->getVector();
440 const ShapeType& evShape=temp_ev->getShape();
441 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
442 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
443 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
444 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
445 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
446 }
447 }
448 }
449 void
450 DataExpanded::trace(DataAbstract* ev, int axis_offset)
451 {
452 int sampleNo,dataPointNo;
453 int numSamples = getNumSamples();
454 int numDataPointsPerSample = getNumDPPSample();
455 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
456 if (temp_ev==0) {
457 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
458 }
459 ValueType& vec=getVector();
460 const ShapeType& shape=getShape();
461 ValueType& evVec=temp_ev->getVector();
462 const ShapeType& evShape=temp_ev->getShape();
463 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
464 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
465 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
466 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
467 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
468 }
469 }
470 }
471
472 void
473 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
474 {
475 int sampleNo,dataPointNo;
476 int numSamples = getNumSamples();
477 int numDataPointsPerSample = getNumDPPSample();
478 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
479 if (temp_ev==0) {
480 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
481 }
482 ValueType& vec=getVector();
483 const ShapeType& shape=getShape();
484 ValueType& evVec=temp_ev->getVector();
485 const ShapeType& evShape=temp_ev->getShape();
486 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
488 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
489 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
490 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
491 }
492 }
493 }
494
495 void
496 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
497 {
498 int sampleNo,dataPointNo;
499 int numSamples = getNumSamples();
500 int numDataPointsPerSample = getNumDPPSample();
501 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
502 if (temp_ev==0) {
503 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
504 }
505 ValueType& vec=getVector();
506 const ShapeType& shape=getShape();
507 ValueType& evVec=temp_ev->getVector();
508 const ShapeType& evShape=temp_ev->getShape();
509 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
510 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
511 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
512 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
513 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
514 }
515 }
516 }
517 void
518 DataExpanded::eigenvalues(DataAbstract* ev)
519 {
520 int sampleNo,dataPointNo;
521 int numSamples = getNumSamples();
522 int numDataPointsPerSample = getNumDPPSample();
523 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
524 if (temp_ev==0) {
525 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
526 }
527 ValueType& vec=getVector();
528 const ShapeType& shape=getShape();
529 ValueType& evVec=temp_ev->getVector();
530 const ShapeType& evShape=temp_ev->getShape();
531 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
532 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
533 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
534 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
535 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
536 }
537 }
538 }
539 void
540 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
541 {
542 int numSamples = getNumSamples();
543 int numDataPointsPerSample = getNumDPPSample();
544 int sampleNo,dataPointNo;
545 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
546 if (temp_ev==0) {
547 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
548 }
549 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
550 if (temp_V==0) {
551 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
552 }
553 ValueType& vec=getVector();
554 const ShapeType& shape=getShape();
555 ValueType& evVec=temp_ev->getVector();
556 const ShapeType& evShape=temp_ev->getShape();
557 ValueType& VVec=temp_V->getVector();
558 const ShapeType& VShape=temp_V->getShape();
559 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
560 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
561 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
562 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
563 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
564 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
565 }
566 }
567 }
568
569 void
570 DataExpanded::setToZero(){
571 // TODO: Surely there is a more efficient way to do this????
572 // Why is there no memset here? Parallel issues?
573 CHECK_FOR_EX_WRITE
574 int numSamples = getNumSamples();
575 int numDataPointsPerSample = getNumDPPSample();
576 DataTypes::ValueType::size_type n = getNoValues();
577 double* p;
578 int sampleNo,dataPointNo, i;
579 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
580 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
581 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
582 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
583 for (i=0; i<n ;++i) p[i]=0.;
584 }
585 }
586 }
587
588 /* Append MPI rank to file name if multiple MPI processes */
589 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
590 /* Make plenty of room for the mpi_rank number and terminating '\0' */
591 char *newFileName = (char *)malloc(strlen(fileName)+20);
592 strncpy(newFileName, fileName, strlen(fileName)+1);
593 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
594 return(newFileName);
595 }
596
597 void
598 DataExpanded::dump(const std::string fileName) const
599 {
600 #ifdef USE_NETCDF
601 const int ldims=2+DataTypes::maxRank;
602 const NcDim* ncdims[ldims];
603 NcVar *var, *ids;
604 int rank = getRank();
605 int type= getFunctionSpace().getTypeCode();
606 int ndims =0;
607 long dims[ldims];
608 const double* d_ptr=&(m_data[0]);
609 const DataTypes::ShapeType& shape = getShape();
610 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
611 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
612 #ifdef PASO_MPI
613 MPI_Status status;
614 #endif
615
616 #ifdef PASO_MPI
617 /* Serialize NetCDF I/O */
618 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
619 #endif
620
621 // netCDF error handler
622 NcError err(NcError::verbose_nonfatal);
623 // Create the file.
624 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
625 NcFile dataFile(newFileName, NcFile::Replace);
626 // check if writing was successful
627 if (!dataFile.is_valid())
628 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
629 if (!dataFile.add_att("type_id",2) )
630 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
631 if (!dataFile.add_att("rank",rank) )
632 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
633 if (!dataFile.add_att("function_space_type",type))
634 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
635 ndims=rank+2;
636 if ( rank >0 ) {
637 dims[0]=shape[0];
638 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
639 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
640 }
641 if ( rank >1 ) {
642 dims[1]=shape[1];
643 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
644 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
645 }
646 if ( rank >2 ) {
647 dims[2]=shape[2];
648 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
649 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
650 }
651 if ( rank >3 ) {
652 dims[3]=shape[3];
653 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
654 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
655 }
656 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
657 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
658 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
659 dims[rank+1]=getFunctionSpace().getNumSamples();
660 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
661 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
662
663 if (getFunctionSpace().getNumSamples()>0)
664 {
665
666 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
667 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
668 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
669 if (! (ids->put(ids_p,dims[rank+1])) )
670 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
671 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
672 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
673 if (! (var->put(d_ptr,dims)) )
674 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
675 }
676 #ifdef PASO_MPI
677 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
678 #endif
679 #else
680 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
681 #endif
682 }
683
684 void
685 DataExpanded::setTaggedValue(int tagKey,
686 const DataTypes::ShapeType& pointshape,
687 const DataTypes::ValueType& value,
688 int dataOffset)
689 {
690 CHECK_FOR_EX_WRITE
691 int numSamples = getNumSamples();
692 int numDataPointsPerSample = getNumDPPSample();
693 int sampleNo,dataPointNo, i;
694 DataTypes::ValueType::size_type n = getNoValues();
695 double* p;
696 const double* in=&value[0+dataOffset];
697
698 if (value.size() != n) {
699 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
700 }
701
702 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
703 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
704 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
705 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
706 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
707 for (i=0; i<n ;++i) p[i]=in[i];
708 }
709 }
710 }
711 }
712
713
714 void
715 DataExpanded::reorderByReferenceIDs(int *reference_ids)
716 {
717 CHECK_FOR_EX_WRITE
718 int numSamples = getNumSamples();
719 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
720 int sampleNo, sampleNo2,i;
721 double* p,*p2;
722 register double rtmp;
723 FunctionSpace fs=getFunctionSpace();
724
725 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
726 const int id_in=reference_ids[sampleNo];
727 const int id=fs.getReferenceIDOfSample(sampleNo);
728 if (id!=id_in) {
729 bool matched=false;
730 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
731 if (id == reference_ids[sampleNo2]) {
732 p=&(m_data[getPointOffset(sampleNo,0)]);
733 p2=&(m_data[getPointOffset(sampleNo2,0)]);
734 for (i=0; i<n ;i++) {
735 rtmp=p[i];
736 p[i]=p2[i];
737 p2[i]=rtmp;
738 }
739 reference_ids[sampleNo]=id;
740 reference_ids[sampleNo2]=id_in;
741 matched=true;
742 break;
743 }
744 }
745 if (! matched) {
746 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
747 }
748 }
749 }
750 }
751
752 DataTypes::ValueType&
753 DataExpanded::getVector()
754 {
755 CHECK_FOR_EX_WRITE
756 return m_data.getData();
757 }
758
759 const DataTypes::ValueType&
760 DataExpanded::getVector() const
761 {
762 return m_data.getData();
763 }
764
765 const DataTypes::ValueType&
766 DataExpanded::getVectorRO() const
767 {
768 return m_data.getData();
769 }
770
771
772 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26