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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2212 - (show annotations)
Wed Jan 14 00:15:00 2009 UTC (10 years, 3 months ago) by jfenwick
File size: 27954 byte(s)
Executive summary:

This commit adds copy on write checks to operations involving shared data. 

Changes:

new #defines:
~~~~~~~~~~~~~
Data.cpp has ASSIGNMENT_MEANS_DEEPCOPY (defaults to undefined).
Defining this will put the data = operator back to making deep copies instead
of sharing data (now the default.)

Data:
~~~~~
. Added exclusiveWrite method to copy the underlying data if it is shared.
. Some operators which took python objects now call the c++ versions intead of duplicating code.

DataAbstract and offspring:
~~~~~~~~~~~~~~~~~~~~~~~~~~~
. Added method to determine whether the data is currently shared.
. Added getVectorRO to children of DataReady.
. Added getTagRO.

. Operations which modify values in place (or return modifiable pointers) now use
a macro to check for sharing. In the case where a modification attempt is detected, it throws an exception. In the future, I will enable this only for debugging.

. This shold not really have been required but the compiler was not choosing the use the const version as I would have liked. Besides, this makes things explict.

. Moved (and de-inlined) getVector in DataConstant (It was virtual in a parent class).

Unit tests:
~~~~~~~~~~~
Added both python and c++ unit tests to check known cases of sharing and "inplace"
modification operations.

General:
~~~~~~~~
Removed some commented out code.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26