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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 6 months ago) by jfenwick
File size: 30739 byte(s)
Don't panic.
Updating copyright stamps

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26