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

Contents of /branches/numpy/escript/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2396 - (show annotations)
Thu Apr 23 23:58:29 2009 UTC (10 years ago) by jfenwick
File size: 30240 byte(s)
Branching to port escript to use numpy rather than numarray

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 int len = other.getNoValues();
176 m_data.resize(len,0.,len);
177 for (int i=0; i<len; i++) {
178 m_data[i]=other.getVectorRO()[i];
179 }
180 }
181
182
183 // Create a new object by copying tags
184 DataTagged::DataTagged(const FunctionSpace& what,
185 const DataTypes::ShapeType& shape,
186 const DataTypes::ValueType& defaultvalue,
187 const DataTagged* tagsource)
188 : parent(what,shape)
189 {
190 // This constructor has not been unit tested yet
191
192 if (defaultvalue.size()!=DataTypes::noValues(shape)) {
193 throw DataException("Programming error - defaultvalue does not match supplied shape.");
194 }
195
196
197 if (!what.canTag())
198 {
199 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
200 }
201
202 if (tagsource!=0)
203 {
204 m_data.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
205
206 DataTagged::DataMapType::const_iterator i;
207 for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
208 addTag(i->first);
209 }
210 }
211 else
212 {
213 m_data.resize(defaultvalue.size());
214 }
215
216 // need to set the default value ....
217 for (int i=0; i<defaultvalue.size(); i++) {
218 m_data[i]=defaultvalue[i];
219 }
220 }
221
222 DataAbstract*
223 DataTagged::deepCopy()
224 {
225 return new DataTagged(*this);
226 }
227
228 DataAbstract*
229 DataTagged::getSlice(const DataTypes::RegionType& region) const
230 {
231 return new DataTagged(*this, region);
232 }
233
234 DataTagged::DataTagged(const DataTagged& other,
235 const DataTypes::RegionType& region)
236 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
237 {
238 // slice constructor
239
240 // get the shape of the slice to copy from other
241 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
242 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
243
244 // allocate enough space in this for all values
245 // (need to add one to allow for the default value)
246 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
247 m_data.resize(len,0.0,len);
248
249 // copy the default value from other to this
250 const DataTypes::ShapeType& otherShape=other.getShape();
251 const DataTypes::ValueType& otherData=other.getVectorRO();
252 DataTypes::copySlice(getVectorRW(),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
253
254 // loop through the tag values copying these
255 DataMapType::const_iterator pos;
256 DataTypes::ValueType::size_type tagOffset=getNoValues();
257 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
258 DataTypes::copySlice(m_data,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
259 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
260 tagOffset+=getNoValues();
261 }
262 }
263
264 void
265 DataTagged::setSlice(const DataAbstract* other,
266 const DataTypes::RegionType& region)
267 {
268
269 // other must be another DataTagged object
270 // Data:setSlice implementation should ensure this
271 const DataTagged* otherTemp=dynamic_cast<const DataTagged*>(other);
272 if (otherTemp==0) {
273 throw DataException("Programming error - casting to DataTagged.");
274 }
275
276 CHECK_FOR_EX_WRITE
277
278 // determine shape of the specified region
279 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
280
281 // modify region specification as needed to match rank of this object
282 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
283
284 // ensure rank/shape of this object is compatible with specified region
285 if (getRank()!=region.size()) {
286 throw DataException("Error - Invalid slice region.");
287 }
288 if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
289 throw DataException (DataTypes::createShapeErrorMessage(
290 "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
291 }
292
293 const DataTypes::ValueType& otherData=otherTemp->getVectorRO();
294 const DataTypes::ShapeType& otherShape=otherTemp->getShape();
295 // copy slice from other default value to this default value
296 DataTypes::copySliceFrom(m_data,getShape(),getDefaultOffset(),otherData,otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
297
298 // loop through tag values in other, adding any which aren't in this, using default value
299 DataMapType::const_iterator pos;
300 for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
301 if (!isCurrentTag(pos->first)) {
302 addTag(pos->first);
303 }
304 }
305
306 // loop through the tag values copying slices from other to this
307 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
308 DataTypes::copySliceFrom(m_data,getShape(),getOffsetForTag(pos->first),otherData, otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
309
310 }
311
312 }
313
314 int
315 DataTagged::getTagNumber(int dpno)
316 {
317 //
318 // Get the number of samples and data-points per sample
319 int numSamples = getNumSamples();
320 int numDataPointsPerSample = getNumDPPSample();
321 int numDataPoints = numSamples * numDataPointsPerSample;
322
323 if (numDataPointsPerSample==0) {
324 throw DataException("DataTagged::getTagNumber error: no data-points associated with this object.");
325 }
326
327 if (dpno<0 || dpno>numDataPoints-1) {
328 throw DataException("DataTagged::getTagNumber error: invalid data-point number supplied.");
329 }
330
331 //
332 // Determine the sample number which corresponds to this data-point number
333 int sampleNo = dpno / numDataPointsPerSample;
334
335 //
336 // Determine the tag number which corresponds to this sample number
337 int tagNo = getFunctionSpace().getTagFromSampleNo(sampleNo);
338
339 //
340 // return the tag number
341 return(tagNo);
342 }
343
344 void
345 DataTagged::setTaggedValue(int tagKey,
346 const DataTypes::ShapeType& pointshape,
347 const ValueType& value,
348 int dataOffset)
349 {
350 if (!DataTypes::checkShape(getShape(), pointshape)) {
351 throw DataException(DataTypes::createShapeErrorMessage(
352 "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
353 }
354 CHECK_FOR_EX_WRITE
355 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
356 if (pos==m_offsetLookup.end()) {
357 // tag couldn't be found so use addTaggedValue
358 addTaggedValue(tagKey,pointshape, value, dataOffset);
359 } else {
360 // copy the values into the data array at the offset determined by m_offsetLookup
361 int offset=pos->second;
362 for (unsigned int i=0; i<getNoValues(); i++) {
363 m_data[offset+i]=value[i+dataOffset];
364 }
365 }
366 }
367
368
369 void
370 DataTagged::addTaggedValues(const TagListType& tagKeys,
371 const ValueBatchType& values,
372 const ShapeType& vShape)
373 {
374 DataTypes::ValueType t(values.size(),0);
375 for (size_t i=0;i<values.size();++i)
376 {
377 t[i]=values[i];
378 }
379 addTaggedValues(tagKeys,t,vShape);
380 }
381
382
383 // Note: The check to see if vShape==our shape is done in the addTaggedValue method
384 void
385 DataTagged::addTaggedValues(const TagListType& tagKeys,
386 const ValueType& values,
387 const ShapeType& vShape)
388 {
389 unsigned int n=getNoValues();
390 unsigned int numVals=values.size()/n;
391 if (values.size()==0) {
392 // copy the current default value for each of the tags
393 TagListType::const_iterator iT;
394 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
395 // the point data view for DataTagged points at the default value
396 addTag(*iT);
397 }
398 } else if (numVals==1 && tagKeys.size()>1) {
399 // assume the one given value will be used for all tag values
400 TagListType::const_iterator iT;
401 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
402 addTaggedValue(*iT, vShape, values,0);
403 }
404 } else {
405 if (tagKeys.size()!=numVals) {
406 stringstream temp;
407 temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
408 << " doesn't match number of values: " << values.size();
409 throw DataException(temp.str());
410 } else {
411 unsigned int i;
412 int offset=0;
413 for (i=0;i<tagKeys.size();i++ ,offset+=n) {
414 addTaggedValue(tagKeys[i],vShape,values,offset);
415 }
416 }
417 }
418 }
419
420
421
422
423 void
424 DataTagged::addTaggedValue(int tagKey,
425 const DataTypes::ShapeType& pointshape,
426 const ValueType& value,
427 int dataOffset)
428 {
429 if (!DataTypes::checkShape(getShape(), pointshape)) {
430 throw DataException(DataTypes::createShapeErrorMessage(
431 "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
432 }
433 CHECK_FOR_EX_WRITE
434 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
435 if (pos!=m_offsetLookup.end()) {
436 // tag already exists so use setTaggedValue
437 setTaggedValue(tagKey,pointshape, value, dataOffset);
438 } else {
439 // save the key and the location of its data in the lookup tab
440 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
441 // add the data given in "value" at the end of m_data
442 // need to make a temp copy of m_data, resize m_data, then copy
443 // all the old values plus the value to be added back into m_data
444 ValueType m_data_temp(m_data);
445 int oldSize=m_data.size();
446 int newSize=m_data.size()+getNoValues();
447 m_data.resize(newSize,0.,newSize);
448 for (int i=0;i<oldSize;i++) {
449 m_data[i]=m_data_temp[i];
450 }
451 for (unsigned int i=0;i<getNoValues();i++) {
452 m_data[oldSize+i]=value[i+dataOffset];
453 }
454 }
455 }
456
457 void
458 DataTagged::addTag(int tagKey)
459 {
460 CHECK_FOR_EX_WRITE
461 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
462 if (pos!=m_offsetLookup.end()) {
463 // tag already exists so use setTaggedValue
464 // setTaggedValue(tagKey,value);
465 } else {
466 // save the key and the location of its data in the lookup tab
467 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data.size()));
468 // add the data given in "value" at the end of m_data
469 // need to make a temp copy of m_data, resize m_data, then copy
470 // all the old values plus the value to be added back into m_data
471 ValueType m_data_temp(m_data);
472 int oldSize=m_data.size();
473 int newSize=m_data.size()+getNoValues();
474 m_data.resize(newSize,0.,newSize);
475 for (int i=0;i<oldSize;i++) {
476 m_data[i]=m_data_temp[i];
477 }
478 for (unsigned int i=0;i<getNoValues();i++) {
479 m_data[oldSize+i]=m_data[m_defaultValueOffset+i];
480 }
481 }
482 }
483
484
485 double*
486 DataTagged::getSampleDataByTag(int tag)
487 {
488 CHECK_FOR_EX_WRITE
489 DataMapType::iterator pos(m_offsetLookup.find(tag));
490 if (pos==m_offsetLookup.end()) {
491 // tag couldn't be found so return the default value
492 return &(m_data[0]);
493 } else {
494 // return the data-point corresponding to the given tag
495 return &(m_data[pos->second]);
496 }
497 }
498
499 string
500 DataTagged::toString() const
501 {
502 using namespace escript::DataTypes;
503 string empty="";
504 stringstream temp;
505 DataMapType::const_iterator i;
506 temp << "Tag(Default)" << endl;
507 temp << pointToString(m_data,getShape(),getDefaultOffset(),empty) << endl;
508 // create a temporary view as the offset will be changed
509 // DataArrayView tempView(getPointDataView().getData(), getPointDataView().getShape());
510 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
511 temp << "Tag(" << i->first << ")" << endl;
512 temp << pointToString(m_data,getShape(),i->second,empty) << endl;
513 // tempView.setOffset(i->second);
514 // temp << tempView.toString() << endl;
515 }
516 return temp.str();
517 }
518
519 DataTypes::ValueType::size_type
520 DataTagged::getPointOffset(int sampleNo,
521 int dataPointNo) const
522 {
523 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
524 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
525 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
526 if (pos!=m_offsetLookup.end()) {
527 offset=pos->second;
528 }
529 return offset;
530 }
531
532 DataTypes::ValueType::size_type
533 DataTagged::getPointOffset(int sampleNo,
534 int dataPointNo)
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::getOffsetForTag(int tag) const
547 {
548 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
549 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
550 if (pos!=m_offsetLookup.end()) {
551 offset=pos->second;
552 }
553 return offset;
554 }
555
556 DataTypes::ValueType::const_reference
557 DataTagged::getDataByTagRO(int tag, DataTypes::ValueType::size_type i) const
558 {
559 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
560 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
561 if (pos!=m_offsetLookup.end()) {
562 offset=pos->second;
563 }
564 return m_data[offset+i];
565 }
566
567 DataTypes::ValueType::reference
568 DataTagged::getDataByTagRW(int tag, DataTypes::ValueType::size_type i)
569 {
570 CHECK_FOR_EX_WRITE
571 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
572 DataTypes::ValueType::size_type offset=m_defaultValueOffset;
573 if (pos!=m_offsetLookup.end()) {
574 offset=pos->second;
575 }
576 return m_data[offset+i];
577 }
578
579 void
580 DataTagged::symmetric(DataAbstract* ev)
581 {
582 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
583 if (temp_ev==0) {
584 throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
585 }
586 const DataTagged::DataMapType& thisLookup=getTagLookup();
587 DataTagged::DataMapType::const_iterator i;
588 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
589 ValueType& evVec=temp_ev->getVectorRW();
590 const ShapeType& evShape=temp_ev->getShape();
591 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
592 temp_ev->addTag(i->first);
593 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
594 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
595 DataMaths::symmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
596 }
597 DataMaths::symmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
598 }
599
600
601 void
602 DataTagged::nonsymmetric(DataAbstract* ev)
603 {
604 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
605 if (temp_ev==0) {
606 throw DataException("Error - DataTagged::nonsymmetric casting to DataTagged failed (probably a programming error).");
607 }
608 const DataTagged::DataMapType& thisLookup=getTagLookup();
609 DataTagged::DataMapType::const_iterator i;
610 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
611 ValueType& evVec=temp_ev->getVectorRW();
612 const ShapeType& evShape=temp_ev->getShape();
613 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
614 temp_ev->addTag(i->first);
615 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
616 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
617 DataMaths::nonsymmetric(m_data,getShape(),offset,evVec, evShape, evoffset);
618 }
619 DataMaths::nonsymmetric(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
620 }
621
622
623 void
624 DataTagged::trace(DataAbstract* ev, int axis_offset)
625 {
626 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
627 if (temp_ev==0) {
628 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
629 }
630 const DataTagged::DataMapType& thisLookup=getTagLookup();
631 DataTagged::DataMapType::const_iterator i;
632 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
633 ValueType& evVec=temp_ev->getVectorRW();
634 const ShapeType& evShape=temp_ev->getShape();
635 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
636 temp_ev->addTag(i->first);
637 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
638 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
639 DataMaths::trace(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
640 }
641 DataMaths::trace(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
642 }
643
644 void
645 DataTagged::transpose(DataAbstract* ev, int axis_offset)
646 {
647 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
648 if (temp_ev==0) {
649 throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
650 }
651 const DataTagged::DataMapType& thisLookup=getTagLookup();
652 DataTagged::DataMapType::const_iterator i;
653 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
654 ValueType& evVec=temp_ev->getVectorRW();
655 const ShapeType& evShape=temp_ev->getShape();
656 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
657 temp_ev->addTag(i->first);
658 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
659 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
660 DataMaths::transpose(m_data,getShape(),offset,evVec, evShape, evoffset, axis_offset);
661 }
662 DataMaths::transpose(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
663 }
664
665 void
666 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
667 {
668 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
669 if (temp_ev==0) {
670 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
671 }
672 const DataTagged::DataMapType& thisLookup=getTagLookup();
673 DataTagged::DataMapType::const_iterator i;
674 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
675 ValueType& evVec=temp_ev->getVectorRW();
676 const ShapeType& evShape=temp_ev->getShape();
677 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
678 temp_ev->addTag(i->first);
679 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
680 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
681 DataMaths::swapaxes(m_data,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
682 }
683 DataMaths::swapaxes(m_data,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
684 }
685
686 void
687 DataTagged::eigenvalues(DataAbstract* ev)
688 {
689 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
690 if (temp_ev==0) {
691 throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (propably a programming error).");
692 }
693 const DataTagged::DataMapType& thisLookup=getTagLookup();
694 DataTagged::DataMapType::const_iterator i;
695 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
696 ValueType& evVec=temp_ev->getVectorRW();
697 const ShapeType& evShape=temp_ev->getShape();
698 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
699 temp_ev->addTag(i->first);
700 // DataArrayView thisView=getDataPointByTag(i->first);
701 // DataArrayView evView=temp_ev->getDataPointByTag(i->first);
702 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
703 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
704 DataMaths::eigenvalues(m_data,getShape(),offset,evVec, evShape, evoffset);
705 }
706 DataMaths::eigenvalues(m_data,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
707 }
708 void
709 DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
710 {
711 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
712 if (temp_ev==0) {
713 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
714 }
715 DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
716 if (temp_V==0) {
717 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (propably a programming error).");
718 }
719 const DataTagged::DataMapType& thisLookup=getTagLookup();
720 DataTagged::DataMapType::const_iterator i;
721 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
722 ValueType& evVec=temp_ev->getVectorRW();
723 const ShapeType& evShape=temp_ev->getShape();
724 ValueType& VVec=temp_V->getVectorRW();
725 const ShapeType& VShape=temp_V->getShape();
726 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
727 temp_ev->addTag(i->first);
728 temp_V->addTag(i->first);
729 /* DataArrayView thisView=getDataPointByTag(i->first);
730 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
731 DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
732 DataTypes::ValueType::size_type offset=getOffsetForTag(i->first);
733 DataTypes::ValueType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
734 DataTypes::ValueType::size_type Voffset=temp_V->getOffsetForTag(i->first);
735 /* DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
736 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
737
738 }
739 DataMaths::eigenvalues_and_eigenvectors(m_data,getShape(),getDefaultOffset(),evVec, evShape,
740 temp_ev->getDefaultOffset(),VVec,VShape,
741 temp_V->getDefaultOffset(), tol);
742
743
744 }
745
746 void
747 DataTagged::setToZero(){
748 CHECK_FOR_EX_WRITE
749 DataTypes::ValueType::size_type n=m_data.size();
750 for (int i=0; i<n ;++i) m_data[i]=0.;
751 }
752
753 void
754 DataTagged::dump(const std::string fileName) const
755 {
756 #ifdef USE_NETCDF
757 const int ldims=DataTypes::maxRank+1;
758 const NcDim* ncdims[ldims];
759 NcVar *var, *tags_var;
760 int rank = getRank();
761 int type= getFunctionSpace().getTypeCode();
762 int ndims =0;
763 long dims[ldims];
764 const double* d_ptr=&(m_data[0]);
765 DataTypes::ShapeType shape = getShape();
766 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
767 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
768 #ifdef PASO_MPI
769 MPI_Status status;
770 #endif
771
772 #ifdef PASO_MPI
773 /* Serialize NetCDF I/O */
774 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, MPI_COMM_WORLD, &status);
775 #endif
776
777 // netCDF error handler
778 NcError err(NcError::verbose_nonfatal);
779 // Create the file.
780 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
781 NcFile dataFile(newFileName, NcFile::Replace);
782 // check if writing was successful
783 if (!dataFile.is_valid())
784 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
785 if (!dataFile.add_att("type_id",1) )
786 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
787 if (!dataFile.add_att("rank",rank) )
788 throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
789 if (!dataFile.add_att("function_space_type",type))
790 throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
791 ndims=rank+1;
792 if ( rank >0 ) {
793 dims[0]=shape[0];
794 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
795 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
796 }
797 if ( rank >1 ) {
798 dims[1]=shape[1];
799 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
800 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
801 }
802 if ( rank >2 ) {
803 dims[2]=shape[2];
804 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
805 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
806 }
807 if ( rank >3 ) {
808 dims[3]=shape[3];
809 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
810 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
811 }
812 const DataTagged::DataMapType& thisLookup=getTagLookup();
813 DataTagged::DataMapType::const_iterator i;
814 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
815 int ntags=1;
816 for (i=thisLookup.begin();i!=thisLookupEnd;i++) ntags++;
817 int* tags =(int*) esysUtils::malloc(ntags*sizeof(int));
818 int c=1;
819 tags[0]=-1;
820 for (i=thisLookup.begin();i!=thisLookupEnd;i++) tags[c++]=i->first;
821 dims[rank]=ntags;
822 if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
823 {
824 esysUtils::free(tags);
825 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
826 }
827 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
828 {
829 esysUtils::free(tags);
830 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
831 }
832 if (! (tags_var->put(tags,dims[rank])) )
833 {
834 esysUtils::free(tags);
835 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
836 }
837 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
838 {
839 esysUtils::free(tags);
840 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
841 }
842 if (! (var->put(d_ptr,dims)) )
843 {
844 esysUtils::free(tags);
845 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
846 }
847 #ifdef PASO_MPI
848 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
849 #endif
850 #else
851 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
852 #endif
853 }
854
855 DataTypes::ValueType&
856 DataTagged::getVectorRW()
857 {
858 CHECK_FOR_EX_WRITE
859 return m_data;
860 }
861
862 const DataTypes::ValueType&
863 DataTagged::getVectorRO() const
864 {
865 return m_data;
866 }
867
868 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26