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

Annotation of /branches/arrexp_2137_win/escript/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2212 - (hide annotations)
Wed Jan 14 00:15:00 2009 UTC (10 years, 5 months ago) by jfenwick
File size: 32448 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 jgs 121
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 ksteube 1312
14 ksteube 1811
15 ksteube 1827 #include "Data.h"
16 jgs 474 #include "DataTagged.h"
17 phornby 1628 #include "esysUtils/esys_malloc.h"
18 jgs 496
19 jgs 474 #include "DataConstant.h"
20     #include "DataException.h"
21 gross 1023 #ifdef USE_NETCDF
22 ksteube 1312 #include <netcdfcpp.h>
23 gross 1023 #endif
24 ksteube 1827 #ifdef PASO_MPI
25     #include <mpi.h>
26     #endif
27 jfenwick 1796 #include "DataMaths.h"
28 jgs 82
29 jfenwick 2212
30     #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
31    
32     // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; throw DataException(ss.str());}
33    
34 jgs 82 using namespace std;
35    
36     namespace escript {
37    
38 jgs 121 DataTagged::DataTagged()
39 jfenwick 2005 : parent(FunctionSpace(),DataTypes::scalarShape)
40 jgs 82 {
41 jgs 509 // default constructor
42    
43 jgs 82 // create a scalar default value
44 jgs 151 m_data.resize(1,0.,1);
45 jfenwick 1796 /* DataArrayView temp(m_data,DataTypes::ShapeType());
46     setPointDataView(temp);*/
47 jgs 82 }
48    
49 jfenwick 1796 // DataTagged::DataTagged(const TagListType& tagKeys,
50     // const ValueListType& values,
51     // const DataArrayView& defaultValue,
52     // const FunctionSpace& what)
53     // : DataAbstract(what)
54     // {
55     // // constructor
56     //
57     // // initialise the array of data values
58     // // the default value is always the first item in the values list
59     // int len = defaultValue.noValues();
60     // m_data.resize(len,0.,len);
61     // for (int i=0; i<defaultValue.noValues(); i++) {
62     // m_data[i]=defaultValue.getData(i);
63     // }
64     //
65     // // create the data view
66     // DataArrayView temp(m_data,defaultValue.getShape());
67     // setPointDataView(temp);
68     //
69     // // add remaining tags and values
70     // addTaggedValues(tagKeys,values);
71     // }
72 jgs 509
73 jgs 119 DataTagged::DataTagged(const FunctionSpace& what,
74 jfenwick 1796 const DataTypes::ShapeType &shape,
75 jgs 119 const int tags[],
76 jgs 500 const ValueType& data)
77 jfenwick 2005 : parent(what,shape)
78 jgs 119 {
79 jgs 509 // alternative constructor
80     // not unit_tested tested yet
81 jfenwick 1796 // It is not explicitly unit tested yet, but it is called from DataFactory
82 jgs 509
83 jfenwick 1802 if (!what.canTag())
84     {
85     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
86     }
87 jgs 121 // copy the data
88 jgs 119 m_data=data;
89 jgs 121
90 jgs 119 // create the view of the data
91 jfenwick 1796 // DataArrayView tempView(m_data,shape);
92     // setPointDataView(tempView);
93 jgs 121
94 jfenwick 1796 // we can't rely on the tag array to give us the number of tags so
95     // use the data we have been passed
96     int valsize=DataTypes::noValues(shape);
97     int ntags=data.size()/valsize;
98    
99 jgs 119 // create the tag lookup map
100 jfenwick 1796 // we assume that the first value and first tag are the default value so we skip
101     for (int i=1;i<ntags;++i)
102     {
103     m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
104 jgs 119 }
105     }
106    
107 woo409 757 DataTagged::DataTagged(const FunctionSpace& what,
108 jfenwick 1796 const DataTypes::ShapeType &shape,
109 woo409 757 const TagListType& tags,
110     const ValueType& data)
111 jfenwick 2005 : parent(what,shape)
112 woo409 757 {
113     // alternative constructor
114    
115 jfenwick 1802 if (!what.canTag())
116     {
117     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
118     }
119    
120 woo409 757 // copy the data
121     m_data=data;
122    
123     // create the view of the data
124 jfenwick 1796 // DataArrayView tempView(m_data,shape);
125     // setPointDataView(tempView);
126 woo409 757
127     // create the tag lookup map
128 jfenwick 1796
129     // for (int sampleNo=0; sampleNo<getNumSamples(); sampleNo++) {
130     // m_offsetLookup.insert(DataMapType::value_type(sampleNo,tags[sampleNo]));
131     // }
132    
133     // The above code looks like it will create a map the wrong way around
134    
135     int valsize=DataTypes::noValues(shape);
136     int npoints=(data.size()/valsize)-1;
137     int ntags=tags.size();
138     if (ntags>npoints)
139     { // This throw is not unit tested yet
140     throw DataException("Programming error - Too many tags for the supplied values.");
141 woo409 757 }
142 jfenwick 1796
143     // create the tag lookup map
144     // we assume that the first value is the default value so we skip it (hence the i+1 below)
145     for (int i=0;i<ntags;++i)
146     {
147     m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
148     }
149 woo409 757 }
150    
151    
152 jgs 102 DataTagged::DataTagged(const DataTagged& other)
153 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape()),
154 jfenwick 1946 m_offsetLookup(other.m_offsetLookup),
155     m_data(other.m_data)
156 jgs 82 {
157 jgs 509 // copy constructor
158    
159 jgs 119 // create the data view
160 jfenwick 1796 // DataArrayView temp(m_data,other.getPointDataView().getShape());
161     // setPointDataView(temp);
162 jgs 82 }
163    
164 jgs 102 DataTagged::DataTagged(const DataConstant& other)
165 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape())
166 jgs 82 {
167 jgs 509 // copy constructor
168    
169 jfenwick 1802 if (!other.getFunctionSpace().canTag())
170     {
171     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
172     }
173    
174 jgs 121 // fill the default value with the constant value item from "other"
175 jfenwick 1796 // const DataArrayView& value=other.getPointDataView();
176     int len = other.getNoValues();
177 jgs 151 m_data.resize(len,0.,len);
178 jfenwick 1796 for (int i=0; i<len; i++) {
179     m_data[i]=other.getVector()[i];
180 jgs 121 }
181    
182 jfenwick 1796 // // create the data view
183     // DataArrayView temp(m_data,value.getShape());
184     // setPointDataView(temp);
185 jgs 82 }
186    
187 jfenwick 1796
188     // Create a new object by copying tags
189     DataTagged::DataTagged(const FunctionSpace& what,
190     const DataTypes::ShapeType& shape,
191     const DataTypes::ValueType& defaultvalue,
192     const DataTagged* tagsource)
193 jfenwick 2005 : parent(what,shape)
194 jfenwick 1796 {
195     // This constructor has not been unit tested yet
196    
197     if (defaultvalue.size()!=DataTypes::noValues(shape)) {
198     throw DataException("Programming error - defaultvalue does not match supplied shape.");
199     }
200    
201    
202 jfenwick 1802 if (!what.canTag())
203     {
204     throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
205     }
206    
207 jfenwick 1796 if (tagsource!=0)
208     {
209 jfenwick 1946 m_data.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
210 jfenwick 1796
211     DataTagged::DataMapType::const_iterator i;
212     for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
213     addTag(i->first);
214     }
215     }
216     else
217     {
218     m_data.resize(defaultvalue.size());
219     }
220    
221     // need to set the default value ....
222     for (int i=0; i<defaultvalue.size(); i++) {
223     m_data[i]=defaultvalue[i];
224     }
225     }
226    
227 jfenwick 1799 DataAbstract*
228     DataTagged::deepCopy()
229     {
230     return new DataTagged(*this);
231     }
232 jfenwick 1796
233 jgs 82 DataAbstract*
234 jfenwick 1796 DataTagged::getSlice(const DataTypes::RegionType& region) const
235 jgs 82 {
236 jgs 513 return new DataTagged(*this, region);
237 jgs 82 }
238    
239 jgs 513 DataTagged::DataTagged(const DataTagged& other,
240 jfenwick 1796 const DataTypes::RegionType& region)
241 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
242 jgs 513 {
243     // slice constructor
244    
245     // get the shape of the slice to copy from other
246 jfenwick 1796 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
247     DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
248 jgs 513
249     // allocate enough space in this for all values
250     // (need to add one to allow for the default value)
251 jfenwick 1796 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
252 jgs 513 m_data.resize(len,0.0,len);
253    
254     // create the data view
255 jfenwick 1796 // DataArrayView temp(m_data,regionShape);
256     // setPointDataView(temp);
257 jgs 513
258     // copy the default value from other to this
259 jfenwick 1796 // getDefaultValue().copySlice(other.getDefaultValue(), regionLoopRange);
260     const DataTypes::ShapeType& otherShape=other.getShape();
261     const DataTypes::ValueType& otherData=other.getVector();
262     DataTypes::copySlice(getVector(),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
263 jgs 513
264     // loop through the tag values copying these
265     DataMapType::const_iterator pos;
266 jfenwick 1796 DataTypes::ValueType::size_type tagOffset=getNoValues();
267 jgs 513 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
268 jfenwick 1796 // getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);
269     DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
270 jgs 513 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
271 jfenwick 1796 tagOffset+=getNoValues();
272 jgs 513 }
273     }
274    
275 jgs 82 void
276 jgs 513 DataTagged::setSlice(const DataAbstract* other,
277 jfenwick 1796 const DataTypes::RegionType& region)
278 jgs 82 {
279 jgs 513
280     // other must be another DataTagged object
281     // Data:setSlice implementation should ensure this
282     const DataTagged* otherTemp=dynamic_cast<const DataTagged*>(other);
283     if (otherTemp==0) {
284 jgs 82 throw DataException("Programming error - casting to DataTagged.");
285     }
286 jgs 121
287 jfenwick 2212 CHECK_FOR_EX_WRITE
288    
289 jgs 513 // determine shape of the specified region
290 jfenwick 1796 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
291 jgs 513
292     // modify region specification as needed to match rank of this object
293 jfenwick 1796 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
294 jgs 513
295     // ensure rank/shape of this object is compatible with specified region
296 jfenwick 1796 if (getRank()!=region.size()) {
297 jgs 82 throw DataException("Error - Invalid slice region.");
298     }
299 jfenwick 1796 if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
300     throw DataException (DataTypes::createShapeErrorMessage(
301     "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
302 jgs 82 }
303 jgs 121
304 jfenwick 1796 const DataTypes::ValueType& otherData=otherTemp->getVector();
305     const DataTypes::ShapeType& otherShape=otherTemp->getShape();
306 jgs 513 // copy slice from other default value to this default value
307 jfenwick 1796 // getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);
308     DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
309 jgs 121
310 jgs 541 // loop through tag values in other, adding any which aren't in this, using default value
311     DataMapType::const_iterator pos;
312     for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
313     if (!isCurrentTag(pos->first)) {
314 jfenwick 1796 addTag(pos->first);
315 jgs 541 }
316     }
317    
318 jgs 513 // loop through the tag values copying slices from other to this
319 jgs 121 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
320 jfenwick 1796 // getDataPointByTag(pos->first).copySliceFrom(otherTemp->getDataPointByTag(pos->first), regionLoopRange);
321     DataTypes::copySliceFrom(m_data,getShape(),getOffsetForTag(pos->first),otherData, otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
322    
323 jgs 121 }
324 jgs 513
325 jgs 82 }
326    
327 jgs 149 int
328     DataTagged::getTagNumber(int dpno)
329     {
330     //
331     // Get the number of samples and data-points per sample
332     int numSamples = getNumSamples();
333     int numDataPointsPerSample = getNumDPPSample();
334     int numDataPoints = numSamples * numDataPointsPerSample;
335    
336     if (numDataPointsPerSample==0) {
337     throw DataException("DataTagged::getTagNumber error: no data-points associated with this object.");
338     }
339    
340 jgs 496 if (dpno<0 || dpno>numDataPoints-1) {
341 jgs 149 throw DataException("DataTagged::getTagNumber error: invalid data-point number supplied.");
342     }
343    
344     //
345     // Determine the sample number which corresponds to this data-point number
346     int sampleNo = dpno / numDataPointsPerSample;
347    
348     //
349     // Determine the tag number which corresponds to this sample number
350     int tagNo = getFunctionSpace().getTagFromSampleNo(sampleNo);
351    
352     //
353     // return the tag number
354     return(tagNo);
355     }
356    
357 jgs 82 void
358 jfenwick 1796 DataTagged::setTaggedValue(int tagKey,
359     const DataTypes::ShapeType& pointshape,
360     const ValueType& value,
361     int dataOffset)
362 jgs 509 {
363 jfenwick 1796 if (!DataTypes::checkShape(getShape(), pointshape)) {
364     throw DataException(DataTypes::createShapeErrorMessage(
365     "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
366     }
367 jfenwick 2212 CHECK_FOR_EX_WRITE
368 jfenwick 1796 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
369     if (pos==m_offsetLookup.end()) {
370     // tag couldn't be found so use addTaggedValue
371     addTaggedValue(tagKey,pointshape, value, dataOffset);
372     } else {
373     // copy the values into the data array at the offset determined by m_offsetLookup
374     int offset=pos->second;
375 phornby 2008 for (unsigned int i=0; i<getNoValues(); i++) {
376 jfenwick 1796 m_data[offset+i]=value[i+dataOffset];
377     }
378     }
379 jgs 509 }
380    
381 jfenwick 1796
382 jgs 509 void
383 jfenwick 1796 DataTagged::addTaggedValues(const TagListType& tagKeys,
384     const ValueBatchType& values,
385     const ShapeType& vShape)
386     {
387     DataTypes::ValueType t(values.size(),0);
388 phornby 1834 for (size_t i=0;i<values.size();++i)
389 jfenwick 1796 {
390     t[i]=values[i];
391     }
392     addTaggedValues(tagKeys,t,vShape);
393 jgs 82 }
394    
395 jfenwick 1796
396     // Note: The check to see if vShape==our shape is done in the addTaggedValue method
397 jgs 82 void
398 jgs 509 DataTagged::addTaggedValues(const TagListType& tagKeys,
399 jfenwick 1796 const ValueType& values,
400     const ShapeType& vShape)
401 jgs 509 {
402 jfenwick 1981 unsigned int n=getNoValues();
403     unsigned int numVals=values.size()/n;
404 jgs 509 if (values.size()==0) {
405     // copy the current default value for each of the tags
406     TagListType::const_iterator iT;
407     for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
408     // the point data view for DataTagged points at the default value
409 jfenwick 1796 addTag(*iT);
410 jgs 509 }
411 jfenwick 1796 } else if (numVals==1 && tagKeys.size()>1) {
412 jgs 509 // assume the one given value will be used for all tag values
413     TagListType::const_iterator iT;
414     for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
415 jfenwick 1796 addTaggedValue(*iT, vShape, values,0);
416 jgs 509 }
417     } else {
418 jfenwick 1796 if (tagKeys.size()!=numVals) {
419 jgs 509 stringstream temp;
420     temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
421     << " doesn't match number of values: " << values.size();
422     throw DataException(temp.str());
423     } else {
424 phornby 1455 unsigned int i;
425 jfenwick 1796 int offset=0;
426     for (i=0;i<tagKeys.size();i++ ,offset+=n) {
427     addTaggedValue(tagKeys[i],vShape,values,offset);
428 jgs 509 }
429     }
430     }
431     }
432    
433 jfenwick 1796
434    
435    
436 jgs 509 void
437 jgs 82 DataTagged::addTaggedValue(int tagKey,
438 jfenwick 1796 const DataTypes::ShapeType& pointshape,
439     const ValueType& value,
440     int dataOffset)
441 jgs 82 {
442 jfenwick 1796 if (!DataTypes::checkShape(getShape(), pointshape)) {
443     throw DataException(DataTypes::createShapeErrorMessage(
444     "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
445 jgs 82 }
446 jfenwick 2212 CHECK_FOR_EX_WRITE
447 jgs 121 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
448     if (pos!=m_offsetLookup.end()) {
449     // tag already exists so use setTaggedValue
450 jfenwick 1796 setTaggedValue(tagKey,pointshape, value, dataOffset);
451 jgs 121 } else {
452     // save the key and the location of its data in the lookup tab
453     m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
454     // add the data given in "value" at the end of m_data
455     // need to make a temp copy of m_data, resize m_data, then copy
456     // all the old values plus the value to be added back into m_data
457     ValueType m_data_temp(m_data);
458     int oldSize=m_data.size();
459 jfenwick 1796 int newSize=m_data.size()+getNoValues();
460 jgs 151 m_data.resize(newSize,0.,newSize);
461 jgs 121 for (int i=0;i<oldSize;i++) {
462     m_data[i]=m_data_temp[i];
463     }
464 phornby 2008 for (unsigned int i=0;i<getNoValues();i++) {
465 jfenwick 1796 m_data[oldSize+i]=value[i+dataOffset];
466 jgs 121 }
467     }
468 jgs 82 }
469    
470 jfenwick 1796 void
471     DataTagged::addTag(int tagKey)
472     {
473 jfenwick 2212 CHECK_FOR_EX_WRITE
474 jfenwick 1796 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
475     if (pos!=m_offsetLookup.end()) {
476     // tag already exists so use setTaggedValue
477     // setTaggedValue(tagKey,value);
478     } else {
479     // save the key and the location of its data in the lookup tab
480     m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
481     // add the data given in "value" at the end of m_data
482     // need to make a temp copy of m_data, resize m_data, then copy
483     // all the old values plus the value to be added back into m_data
484     ValueType m_data_temp(m_data);
485     int oldSize=m_data.size();
486     int newSize=m_data.size()+getNoValues();
487     m_data.resize(newSize,0.,newSize);
488     for (int i=0;i<oldSize;i++) {
489     m_data[i]=m_data_temp[i];
490     }
491 phornby 2008 for (unsigned int i=0;i<getNoValues();i++) {
492 jfenwick 1796 m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
493     }
494     }
495     }
496    
497    
498 jgs 82 double*
499     DataTagged::getSampleDataByTag(int tag)
500     {
501 jfenwick 2212 CHECK_FOR_EX_WRITE
502 jgs 82 DataMapType::iterator pos(m_offsetLookup.find(tag));
503     if (pos==m_offsetLookup.end()) {
504     // tag couldn't be found so return the default value
505     return &(m_data[0]);
506     } else {
507     // return the data-point corresponding to the given tag
508     return &(m_data[pos->second]);
509     }
510     }
511    
512     string
513     DataTagged::toString() const
514     {
515 jfenwick 1796 using namespace escript::DataTypes;
516     string empty="";
517 jgs 82 stringstream temp;
518     DataMapType::const_iterator i;
519     temp << "Tag(Default)" << endl;
520 jfenwick 1796 temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
521 jgs 82 // create a temporary view as the offset will be changed
522 jfenwick 1796 // DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
523 jgs 82 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
524     temp << "Tag(" << i->first << ")" << endl;
525 jfenwick 1802 temp << pointToString(m_data,getShape(),i->second,empty) << endl;
526 jfenwick 1796 // tempView.setOffset(i->second);
527     // temp << tempView.toString() << endl;
528 jgs 82 }
529 ksteube 1312 return temp.str();
530 jgs 82 }
531    
532 jfenwick 1796 DataTypes::ValueType::size_type
533 jgs 496 DataTagged::getPointOffset(int sampleNo,
534     int dataPointNo) const
535 jgs 82 {
536 jgs 496 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
537     DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
538 jfenwick 1796 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
539 jgs 82 if (pos!=m_offsetLookup.end()) {
540     offset=pos->second;
541     }
542 jgs 496 return offset;
543 jgs 82 }
544    
545 jfenwick 2005 DataTypes::ValueType::size_type
546     DataTagged::getPointOffset(int sampleNo,
547     int dataPointNo)
548     {
549     int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
550     DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
551     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
552     if (pos!=m_offsetLookup.end()) {
553     offset=pos->second;
554     }
555     return offset;
556     }
557    
558 jfenwick 1796 // DataArrayView
559     // DataTagged::getDataPointByTag(int tag) const
560     // {
561     // DataMapType::const_iterator pos(m_offsetLookup.find(tag));
562     // DataTypes::ValueType::size_type offset=m_defaultValueOffset;
563     // if (pos!=m_offsetLookup.end()) {
564     // offset=pos->second;
565     // }
566     // DataArrayView temp(getPointDataView());
567     // temp.setOffset(offset);
568     // return temp;
569     // }
570     //
571    
572    
573     DataTypes::ValueType::size_type
574     DataTagged::getOffsetForTag(int tag) const
575 jgs 82 {
576 jgs 496 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
577 jfenwick 1796 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
578 jgs 82 if (pos!=m_offsetLookup.end()) {
579     offset=pos->second;
580     }
581 jfenwick 1796 return offset;
582 jgs 82 }
583    
584 jfenwick 1796 DataTypes::ValueType::const_reference
585     DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
586 jgs 82 {
587 jfenwick 1796 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
588     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
589     if (pos!=m_offsetLookup.end()) {
590     offset=pos->second;
591     }
592     return m_data[offset+i];
593 jgs 82 }
594    
595 jfenwick 2212 DataTypes::ValueType::const_reference
596     DataTagged::getDataByTagRO(int tag, DataTypes::ValueType::size_type i) const
597     {
598     DataMapType::const_iterator pos(m_offsetLookup.find(tag));
599     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
600     if (pos!=m_offsetLookup.end()) {
601     offset=pos->second;
602     }
603     return m_data[offset+i];
604     }
605 jfenwick 1796
606 jfenwick 2212
607 jfenwick 1796 DataTypes::ValueType::reference
608     DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
609 jgs 123 {
610 jfenwick 2212 CHECK_FOR_EX_WRITE
611 jfenwick 1796 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
612     DataTypes::ValueType::size_type offset=m_defaultValueOffset;
613     if (pos!=m_offsetLookup.end()) {
614     offset=pos->second;
615     }
616     return m_data[offset+i];
617 jgs 123 }
618    
619 gross 580 void
620 ksteube 775 DataTagged::symmetric(DataAbstract* ev)
621     {
622     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
623     if (temp_ev==0) {
624     throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
625     }
626     const DataTagged::DataMapType& thisLookup=getTagLookup();
627     DataTagged::DataMapType::const_iterator i;
628     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
629 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
630     const ShapeType& evShape=temp_ev->getShape();
631 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
632 jfenwick 1796 temp_ev->addTag(i->first);
633     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
634     // DataArrayView thisView=getDataPointByTag(i->first);
635     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
636     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
637    
638     // DataArrayView::symmetric(thisView,0,evView,0);
639     DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
640 ksteube 775 }
641 jfenwick 1796 // symmetric(m_data,getShape(),getDefaultOffset(),
642     DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
643 ksteube 775 }
644 jfenwick 1796
645    
646 ksteube 775 void
647     DataTagged::nonsymmetric(DataAbstract* ev)
648     {
649     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
650     if (temp_ev==0) {
651     throw DataException("Error - DataTagged::nonsymmetric casting to DataTagged failed (probably a programming error).");
652     }
653     const DataTagged::DataMapType& thisLookup=getTagLookup();
654     DataTagged::DataMapType::const_iterator i;
655     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
656 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
657     const ShapeType& evShape=temp_ev->getShape();
658 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
659 jfenwick 1796 temp_ev->addTag(i->first);
660     /* DataArrayView thisView=getDataPointByTag(i->first);
661     DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
662     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
663     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
664     DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
665 ksteube 775 }
666 jfenwick 1796 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
667 ksteube 775 }
668 jfenwick 1796
669    
670 ksteube 775 void
671 gross 800 DataTagged::trace(DataAbstract* ev, int axis_offset)
672 ksteube 775 {
673     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
674     if (temp_ev==0) {
675 gross 800 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
676 ksteube 775 }
677     const DataTagged::DataMapType& thisLookup=getTagLookup();
678     DataTagged::DataMapType::const_iterator i;
679     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
680 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
681     const ShapeType& evShape=temp_ev->getShape();
682 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
683 jfenwick 1796 temp_ev->addTag(i->first);
684     // DataArrayView thisView=getDataPointByTag(i->first);
685     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
686     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
687     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
688     DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
689 ksteube 775 }
690 jfenwick 1796 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
691 ksteube 775 }
692 gross 800
693 ksteube 775 void
694     DataTagged::transpose(DataAbstract* ev, int axis_offset)
695     {
696     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
697     if (temp_ev==0) {
698     throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
699     }
700     const DataTagged::DataMapType& thisLookup=getTagLookup();
701     DataTagged::DataMapType::const_iterator i;
702     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
703 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
704     const ShapeType& evShape=temp_ev->getShape();
705 ksteube 775 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
706 jfenwick 1796 temp_ev->addTag(i->first);
707     // DataArrayView thisView=getDataPointByTag(i->first);
708     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
709     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
710     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
711     DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
712 ksteube 775 }
713 jfenwick 1796 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
714 ksteube 775 }
715 gross 800
716 ksteube 775 void
717 gross 804 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
718 gross 800 {
719     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
720     if (temp_ev==0) {
721 gross 804 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
722 gross 800 }
723     const DataTagged::DataMapType& thisLookup=getTagLookup();
724     DataTagged::DataMapType::const_iterator i;
725     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
726 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
727     const ShapeType& evShape=temp_ev->getShape();
728 gross 800 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
729 jfenwick 1796 temp_ev->addTag(i->first);
730     /* DataArrayView thisView=getDataPointByTag(i->first);
731     DataArrayView evView=temp_ev->getDataPointByTag(i->first);*/
732     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
733     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
734     DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
735 gross 800 }
736 jfenwick 1796 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
737 gross 800 }
738    
739     void
740 gross 580 DataTagged::eigenvalues(DataAbstract* ev)
741     {
742     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
743     if (temp_ev==0) {
744     throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (propably a programming error).");
745     }
746     const DataTagged::DataMapType& thisLookup=getTagLookup();
747     DataTagged::DataMapType::const_iterator i;
748     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
749 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
750     const ShapeType& evShape=temp_ev->getShape();
751 gross 580 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
752 jfenwick 1796 temp_ev->addTag(i->first);
753     // DataArrayView thisView=getDataPointByTag(i->first);
754     // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
755     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
756     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
757     DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
758 gross 580 }
759 jfenwick 1796 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
760 gross 580 }
761     void
762     DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
763     {
764     DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
765     if (temp_ev==0) {
766     throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
767     }
768     DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
769     if (temp_V==0) {
770     throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
771     }
772     const DataTagged::DataMapType& thisLookup=getTagLookup();
773     DataTagged::DataMapType::const_iterator i;
774     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
775 jfenwick 1796 ValueType& evVec=temp_ev->getVector();
776     const ShapeType& evShape=temp_ev->getShape();
777     ValueType& VVec=temp_V->getVector();
778     const ShapeType& VShape=temp_V->getShape();
779 gross 580 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
780 jfenwick 1796 temp_ev->addTag(i->first);
781     temp_V->addTag(i->first);
782     /* DataArrayView thisView=getDataPointByTag(i->first);
783 gross 580 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
784 jfenwick 1796 DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
785     DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
786     DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
787     DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
788     /* DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
789     DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
790    
791 gross 580 }
792 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
793     temp_ev->getDefaultOffset(),VVec,VShape,
794     temp_V->getDefaultOffset(), tol);
795 jgs 123
796 gross 580
797     }
798    
799 gross 950 void
800 gross 1118 DataTagged::setToZero(){
801 jfenwick 2212 CHECK_FOR_EX_WRITE
802 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
803 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
804     }
805    
806     void
807 gross 950 DataTagged::dump(const std::string fileName) const
808     {
809 gross 1023 #ifdef USE_NETCDF
810 jfenwick 1796 const int ldims=DataTypes::maxRank+1;
811 gross 983 const NcDim* ncdims[ldims];
812     NcVar *var, *tags_var;
813 jfenwick 1796 int rank = getRank();
814 gross 983 int type= getFunctionSpace().getTypeCode();
815     int ndims =0;
816     long dims[ldims];
817 gross 1141 const double* d_ptr=&(m_data[0]);
818 jfenwick 1796 DataTypes::ShapeType shape = getShape();
819 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
820     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
821 ksteube 1827 #ifdef PASO_MPI
822     MPI_Status status;
823     #endif
824 gross 983
825 ksteube 1827 #ifdef PASO_MPI
826     /* Serialize NetCDF I/O */
827     if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
828     #endif
829    
830 gross 983 // netCDF error handler
831     NcError err(NcError::verbose_nonfatal);
832     // Create the file.
833 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
834     NcFile dataFile(newFileName, NcFile::Replace);
835 gross 983 // check if writing was successful
836     if (!dataFile.is_valid())
837     throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
838 gross 1141 if (!dataFile.add_att("type_id",1) )
839 gross 983 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
840     if (!dataFile.add_att("rank",rank) )
841     throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
842     if (!dataFile.add_att("function_space_type",type))
843     throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
844     ndims=rank+1;
845     if ( rank >0 ) {
846     dims[0]=shape[0];
847     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
848 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
849 gross 983 }
850     if ( rank >1 ) {
851     dims[1]=shape[1];
852     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
853 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
854 gross 983 }
855     if ( rank >2 ) {
856     dims[2]=shape[2];
857     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
858 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
859 gross 983 }
860     if ( rank >3 ) {
861     dims[3]=shape[3];
862     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
863 ksteube 1800 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
864 gross 983 }
865     const DataTagged::DataMapType& thisLookup=getTagLookup();
866     DataTagged::DataMapType::const_iterator i;
867     DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
868     int ntags=1;
869     for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
870 phornby 1628 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
871 phornby 1020 int c=1;
872 gross 983 tags[0]=-1;
873     for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
874     dims[rank]=ntags;
875     if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
876 phornby 1020 {
877 phornby 1628 esysUtils::free(tags);
878 phornby 1020 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
879     }
880 gross 983 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
881 phornby 1020 {
882 phornby 1628 esysUtils::free(tags);
883 gross 983 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
884 phornby 1020 }
885 gross 983 if (! (tags_var->put(tags,dims[rank])) )
886 phornby 1020 {
887 phornby 1628 esysUtils::free(tags);
888 gross 983 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
889 phornby 1020 }
890 gross 983 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
891 phornby 1020 {
892 phornby 1628 esysUtils::free(tags);
893 gross 983 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
894 phornby 1020 }
895 gross 1141 if (! (var->put(d_ptr,dims)) )
896 phornby 1020 {
897 phornby 1628 esysUtils::free(tags);
898 gross 983 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
899 phornby 1020 }
900 ksteube 1827 #ifdef PASO_MPI
901     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
902     #endif
903 gross 1023 #else
904     throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
905     #endif
906 gross 950 }
907 jfenwick 1796
908     DataTypes::ValueType&
909     DataTagged::getVector()
910     {
911 jfenwick 2212 CHECK_FOR_EX_WRITE
912 jfenwick 1796 return m_data;
913     }
914    
915     const DataTypes::ValueType&
916     DataTagged::getVector() const
917     {
918 jfenwick 2212 // exclusive write not required for const access
919 jfenwick 1796 return m_data;
920     }
921    
922 jfenwick 2212 const DataTypes::ValueType&
923     DataTagged::getVectorRO() const
924     {
925     return m_data;
926     }
927    
928 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26