/[escript]/branches/subworld2/escriptcore/src/DataTagged.cpp
ViewVC logotype

Contents of /branches/subworld2/escriptcore/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5504 - (show annotations)
Wed Mar 4 22:58:13 2015 UTC (4 years, 1 month ago) by jfenwick
File size: 31186 byte(s)
Again with a more up to date copy


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26