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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 11 months ago) by caltinay
File size: 30733 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

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 #include "esysUtils/Esys_MPI.h"
22
23 #ifdef USE_NETCDF
24 #include <netcdfcpp.h>
25 #endif
26
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 ESYS_MPI
787 MPI_Status status;
788 #endif
789
790 #ifdef ESYS_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 ESYS_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