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

Contents of /branches/arrexp_2137_win/escript/src/DataTagged.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: 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
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 "DataTagged.h"
17 #include "esysUtils/esys_malloc.h"
18
19 #include "DataConstant.h"
20 #include "DataException.h"
21 #ifdef USE_NETCDF
22 #include <netcdfcpp.h>
23 #endif
24 #ifdef PASO_MPI
25 #include <mpi.h>
26 #endif
27 #include "DataMaths.h"
28
29
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 using namespace std;
35
36 namespace escript {
37
38 DataTagged::DataTagged()
39 : parent(FunctionSpace(),DataTypes::scalarShape)
40 {
41 // default constructor
42
43 // create a scalar default value
44 m_data.resize(1,0.,1);
45 /* DataArrayView temp(m_data,DataTypes::ShapeType());
46 setPointDataView(temp);*/
47 }
48
49 // 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
73 DataTagged::DataTagged(const FunctionSpace& what,
74 const DataTypes::ShapeType &shape,
75 const int tags[],
76 const ValueType& data)
77 : parent(what,shape)
78 {
79 // alternative constructor
80 // not unit_tested tested yet
81 // It is not explicitly unit tested yet, but it is called from DataFactory
82
83 if (!what.canTag())
84 {
85 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
86 }
87 // copy the data
88 m_data=data;
89
90 // create the view of the data
91 // DataArrayView tempView(m_data,shape);
92 // setPointDataView(tempView);
93
94 // 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 // create the tag lookup map
100 // 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 }
105 }
106
107 DataTagged::DataTagged(const FunctionSpace& what,
108 const DataTypes::ShapeType &shape,
109 const TagListType& tags,
110 const ValueType& data)
111 : parent(what,shape)
112 {
113 // alternative constructor
114
115 if (!what.canTag())
116 {
117 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
118 }
119
120 // copy the data
121 m_data=data;
122
123 // create the view of the data
124 // DataArrayView tempView(m_data,shape);
125 // setPointDataView(tempView);
126
127 // create the tag lookup map
128
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 }
142
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 }
150
151
152 DataTagged::DataTagged(const DataTagged& other)
153 : parent(other.getFunctionSpace(),other.getShape()),
154 m_offsetLookup(other.m_offsetLookup),
155 m_data(other.m_data)
156 {
157 // copy constructor
158
159 // create the data view
160 // DataArrayView temp(m_data,other.getPointDataView().getShape());
161 // setPointDataView(temp);
162 }
163
164 DataTagged::DataTagged(const DataConstant& other)
165 : parent(other.getFunctionSpace(),other.getShape())
166 {
167 // copy constructor
168
169 if (!other.getFunctionSpace().canTag())
170 {
171 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
172 }
173
174 // fill the default value with the constant value item from "other"
175 // const DataArrayView& value=other.getPointDataView();
176 int len = other.getNoValues();
177 m_data.resize(len,0.,len);
178 for (int i=0; i<len; i++) {
179 m_data[i]=other.getVector()[i];
180 }
181
182 // // create the data view
183 // DataArrayView temp(m_data,value.getShape());
184 // setPointDataView(temp);
185 }
186
187
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 : parent(what,shape)
194 {
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 if (!what.canTag())
203 {
204 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
205 }
206
207 if (tagsource!=0)
208 {
209 m_data.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
210
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 DataAbstract*
228 DataTagged::deepCopy()
229 {
230 return new DataTagged(*this);
231 }
232
233 DataAbstract*
234 DataTagged::getSlice(const DataTypes::RegionType& region) const
235 {
236 return new DataTagged(*this, region);
237 }
238
239 DataTagged::DataTagged(const DataTagged& other,
240 const DataTypes::RegionType& region)
241 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
242 {
243 // slice constructor
244
245 // get the shape of the slice to copy from other
246 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
247 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
248
249 // allocate enough space in this for all values
250 // (need to add one to allow for the default value)
251 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
252 m_data.resize(len,0.0,len);
253
254 // create the data view
255 // DataArrayView temp(m_data,regionShape);
256 // setPointDataView(temp);
257
258 // copy the default value from other to this
259 // 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
264 // loop through the tag values copying these
265 DataMapType::const_iterator pos;
266 DataTypes::ValueType::size_type tagOffset=getNoValues();
267 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
268 // getPointDataView().copySlice(tagOffset,other.getPointDataView(),pos->second,regionLoopRange);
269 DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
270 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
271 tagOffset+=getNoValues();
272 }
273 }
274
275 void
276 DataTagged::setSlice(const DataAbstract* other,
277 const DataTypes::RegionType& region)
278 {
279
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 throw DataException("Programming error - casting to DataTagged.");
285 }
286
287 CHECK_FOR_EX_WRITE
288
289 // determine shape of the specified region
290 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
291
292 // modify region specification as needed to match rank of this object
293 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
294
295 // ensure rank/shape of this object is compatible with specified region
296 if (getRank()!=region.size()) {
297 throw DataException("Error - Invalid slice region.");
298 }
299 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 }
303
304 const DataTypes::ValueType& otherData=otherTemp->getVector();
305 const DataTypes::ShapeType& otherShape=otherTemp->getShape();
306 // copy slice from other default value to this default value
307 // getDefaultValue().copySliceFrom(otherTemp->getDefaultValue(), regionLoopRange);
308 DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
309
310 // 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 addTag(pos->first);
315 }
316 }
317
318 // loop through the tag values copying slices from other to this
319 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
320 // 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 }
324
325 }
326
327 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 if (dpno<0 || dpno>numDataPoints-1) {
341 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 void
358 DataTagged::setTaggedValue(int tagKey,
359 const DataTypes::ShapeType& pointshape,
360 const ValueType& value,
361 int dataOffset)
362 {
363 if (!DataTypes::checkShape(getShape(), pointshape)) {
364 throw DataException(DataTypes::createShapeErrorMessage(
365 "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
366 }
367 CHECK_FOR_EX_WRITE
368 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 for (unsigned int i=0; i<getNoValues(); i++) {
376 m_data[offset+i]=value[i+dataOffset];
377 }
378 }
379 }
380
381
382 void
383 DataTagged::addTaggedValues(const TagListType& tagKeys,
384 const ValueBatchType& values,
385 const ShapeType& vShape)
386 {
387 DataTypes::ValueType t(values.size(),0);
388 for (size_t i=0;i<values.size();++i)
389 {
390 t[i]=values[i];
391 }
392 addTaggedValues(tagKeys,t,vShape);
393 }
394
395
396 // Note: The check to see if vShape==our shape is done in the addTaggedValue method
397 void
398 DataTagged::addTaggedValues(const TagListType& tagKeys,
399 const ValueType& values,
400 const ShapeType& vShape)
401 {
402 unsigned int n=getNoValues();
403 unsigned int numVals=values.size()/n;
404 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 addTag(*iT);
410 }
411 } else if (numVals==1 && tagKeys.size()>1) {
412 // 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 addTaggedValue(*iT, vShape, values,0);
416 }
417 } else {
418 if (tagKeys.size()!=numVals) {
419 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 unsigned int i;
425 int offset=0;
426 for (i=0;i<tagKeys.size();i++ ,offset+=n) {
427 addTaggedValue(tagKeys[i],vShape,values,offset);
428 }
429 }
430 }
431 }
432
433
434
435
436 void
437 DataTagged::addTaggedValue(int tagKey,
438 const DataTypes::ShapeType& pointshape,
439 const ValueType& value,
440 int dataOffset)
441 {
442 if (!DataTypes::checkShape(getShape(), pointshape)) {
443 throw DataException(DataTypes::createShapeErrorMessage(
444 "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
445 }
446 CHECK_FOR_EX_WRITE
447 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
448 if (pos!=m_offsetLookup.end()) {
449 // tag already exists so use setTaggedValue
450 setTaggedValue(tagKey,pointshape, value, dataOffset);
451 } 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 int newSize=m_data.size()+getNoValues();
460 m_data.resize(newSize,0.,newSize);
461 for (int i=0;i<oldSize;i++) {
462 m_data[i]=m_data_temp[i];
463 }
464 for (unsigned int i=0;i<getNoValues();i++) {
465 m_data[oldSize+i]=value[i+dataOffset];
466 }
467 }
468 }
469
470 void
471 DataTagged::addTag(int tagKey)
472 {
473 CHECK_FOR_EX_WRITE
474 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 for (unsigned int i=0;i<getNoValues();i++) {
492 m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
493 }
494 }
495 }
496
497
498 double*
499 DataTagged::getSampleDataByTag(int tag)
500 {
501 CHECK_FOR_EX_WRITE
502 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 using namespace escript::DataTypes;
516 string empty="";
517 stringstream temp;
518 DataMapType::const_iterator i;
519 temp << "Tag(Default)" << endl;
520 temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
521 // create a temporary view as the offset will be changed
522 // DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
523 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
524 temp << "Tag(" << i->first << ")" << endl;
525 temp << pointToString(m_data,getShape(),i->second,empty) << endl;
526 // tempView.setOffset(i->second);
527 // temp << tempView.toString() << endl;
528 }
529 return temp.str();
530 }
531
532 DataTypes::ValueType::size_type
533 DataTagged::getPointOffset(int sampleNo,
534 int dataPointNo) const
535 {
536 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
537 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
538 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
539 if (pos!=m_offsetLookup.end()) {
540 offset=pos->second;
541 }
542 return offset;
543 }
544
545 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 // 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 {
576 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
577 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
578 if (pos!=m_offsetLookup.end()) {
579 offset=pos->second;
580 }
581 return offset;
582 }
583
584 DataTypes::ValueType::const_reference
585 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i) const
586 {
587 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 }
594
595 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
606
607 DataTypes::ValueType::reference
608 DataTagged::getDataByTag(int tag, DataTypes::ValueType::size_type i)
609 {
610 CHECK_FOR_EX_WRITE
611 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 }
618
619 void
620 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 ValueType& evVec=temp_ev->getVector();
630 const ShapeType& evShape=temp_ev->getShape();
631 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
632 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 }
641 // symmetric(m_data,getShape(),getDefaultOffset(),
642 DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
643 }
644
645
646 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 ValueType& evVec=temp_ev->getVector();
657 const ShapeType& evShape=temp_ev->getShape();
658 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
659 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 }
666 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
667 }
668
669
670 void
671 DataTagged::trace(DataAbstract* ev, int axis_offset)
672 {
673 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
674 if (temp_ev==0) {
675 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
676 }
677 const DataTagged::DataMapType& thisLookup=getTagLookup();
678 DataTagged::DataMapType::const_iterator i;
679 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
680 ValueType& evVec=temp_ev->getVector();
681 const ShapeType& evShape=temp_ev->getShape();
682 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
683 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 }
690 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
691 }
692
693 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 ValueType& evVec=temp_ev->getVector();
704 const ShapeType& evShape=temp_ev->getShape();
705 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
706 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 }
713 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
714 }
715
716 void
717 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
718 {
719 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
720 if (temp_ev==0) {
721 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
722 }
723 const DataTagged::DataMapType& thisLookup=getTagLookup();
724 DataTagged::DataMapType::const_iterator i;
725 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
726 ValueType& evVec=temp_ev->getVector();
727 const ShapeType& evShape=temp_ev->getShape();
728 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
729 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 }
736 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
737 }
738
739 void
740 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 ValueType& evVec=temp_ev->getVector();
750 const ShapeType& evShape=temp_ev->getShape();
751 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
752 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 }
759 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
760 }
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 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 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
780 temp_ev->addTag(i->first);
781 temp_V->addTag(i->first);
782 /* DataArrayView thisView=getDataPointByTag(i->first);
783 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
784 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 }
792 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
793 temp_ev->getDefaultOffset(),VVec,VShape,
794 temp_V->getDefaultOffset(), tol);
795
796
797 }
798
799 void
800 DataTagged::setToZero(){
801 CHECK_FOR_EX_WRITE
802 DataTypes::ValueType::size_type n=m_data.size();
803 for (int i=0; i<n ;++i) m_data[i]=0.;
804 }
805
806 void
807 DataTagged::dump(const std::string fileName) const
808 {
809 #ifdef USE_NETCDF
810 const int ldims=DataTypes::maxRank+1;
811 const NcDim* ncdims[ldims];
812 NcVar *var, *tags_var;
813 int rank = getRank();
814 int type= getFunctionSpace().getTypeCode();
815 int ndims =0;
816 long dims[ldims];
817 const double* d_ptr=&(m_data[0]);
818 DataTypes::ShapeType shape = getShape();
819 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
820 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
821 #ifdef PASO_MPI
822 MPI_Status status;
823 #endif
824
825 #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 // netCDF error handler
831 NcError err(NcError::verbose_nonfatal);
832 // Create the file.
833 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
834 NcFile dataFile(newFileName, NcFile::Replace);
835 // check if writing was successful
836 if (!dataFile.is_valid())
837 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
838 if (!dataFile.add_att("type_id",1) )
839 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 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
849 }
850 if ( rank >1 ) {
851 dims[1]=shape[1];
852 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
853 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
854 }
855 if ( rank >2 ) {
856 dims[2]=shape[2];
857 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
858 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
859 }
860 if ( rank >3 ) {
861 dims[3]=shape[3];
862 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
863 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
864 }
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 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
871 int c=1;
872 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 {
877 esysUtils::free(tags);
878 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
879 }
880 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
881 {
882 esysUtils::free(tags);
883 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
884 }
885 if (! (tags_var->put(tags,dims[rank])) )
886 {
887 esysUtils::free(tags);
888 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
889 }
890 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
891 {
892 esysUtils::free(tags);
893 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
894 }
895 if (! (var->put(d_ptr,dims)) )
896 {
897 esysUtils::free(tags);
898 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
899 }
900 #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 #else
904 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
905 #endif
906 }
907
908 DataTypes::ValueType&
909 DataTagged::getVector()
910 {
911 CHECK_FOR_EX_WRITE
912 return m_data;
913 }
914
915 const DataTypes::ValueType&
916 DataTagged::getVector() const
917 {
918 // exclusive write not required for const access
919 return m_data;
920 }
921
922 const DataTypes::ValueType&
923 DataTagged::getVectorRO() const
924 {
925 return m_data;
926 }
927
928 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26