/[escript]/trunk/escript/src/DataTagged.cpp
ViewVC logotype

Contents of /trunk/escript/src/DataTagged.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26