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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2742 - (hide annotations)
Thu Nov 12 06:03:37 2009 UTC (10 years ago) by jfenwick
File size: 31518 byte(s)
Merging changes from the lapack branch.

The inverse() operation has been moved into c++. [No lazy support for this operation yet.]
Optional Lapack support has been added for matrices larger than 3x3. 
service0 is set to use mkl_lapack.



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26