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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6144 - (show annotations)
Wed Apr 6 05:25:13 2016 UTC (2 years, 7 months ago) by caltinay
File size: 48785 byte(s)
last round of namespacing defines.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include "Data.h"
18 #include "DataConstant.h"
19 #include "DataException.h"
20 #include "DataVectorOps.h"
21 #include "DataTagged.h"
22
23 #include <complex>
24
25 #ifdef ESYS_HAVE_NETCDF
26 #include <netcdfcpp.h>
27 #endif
28
29 #ifdef SLOWSHARECHECK
30 #define CHECK_FOR_EX_WRITE if (isShared()) {throw DataException("Attempt to modify shared object");}
31 #else
32 #define CHECK_FOR_EX_WRITE
33 #endif
34
35 using namespace std;
36
37 namespace escript {
38
39 DataTagged::DataTagged(const FunctionSpace& what,
40 const DataTypes::ShapeType &shape,
41 const int tags[],
42 const DataTypes::RealVectorType& data)
43 : parent(what,shape)
44 {
45 // alternative constructor
46 // not unit_tested tested yet
47 // It is not explicitly unit tested yet, but it is called from DataFactory
48
49 if (!what.canTag())
50 {
51 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
52 }
53 // copy the data
54 m_data_r=data;
55
56 // we can't rely on the tag array to give us the number of tags so
57 // use the data we have been passed
58 int valsize=DataTypes::noValues(shape);
59 int ntags=data.size()/valsize;
60
61 // create the tag lookup map
62 // we assume that the first value and first tag are the default value so we skip
63 for (int i=1;i<ntags;++i)
64 {
65 m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
66 }
67 }
68
69
70 DataTagged::DataTagged(const FunctionSpace& what,
71 const DataTypes::ShapeType &shape,
72 const int tags[],
73 const DataTypes::CplxVectorType& data)
74 : parent(what,shape)
75 {
76 // alternative constructor
77 // not unit_tested tested yet
78 // It is not explicitly unit tested yet, but it is called from DataFactory
79
80 m_iscompl=true;
81 if (!what.canTag())
82 {
83 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
84 }
85 // copy the data
86 m_data_c=data;
87
88 // we can't rely on the tag array to give us the number of tags so
89 // use the data we have been passed
90 int valsize=DataTypes::noValues(shape);
91 int ntags=data.size()/valsize;
92
93 // create the tag lookup map
94 // we assume that the first value and first tag are the default value so we skip
95 for (int i=1;i<ntags;++i)
96 {
97 m_offsetLookup.insert(DataMapType::value_type(tags[i],i*valsize));
98 }
99 }
100
101
102
103 DataTagged::DataTagged(const FunctionSpace& what,
104 const DataTypes::ShapeType &shape,
105 const TagListType& tags,
106 const DataTypes::RealVectorType& data)
107 : parent(what,shape)
108 {
109 // alternative constructor
110
111 if (!what.canTag())
112 {
113 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
114 }
115
116 // copy the data
117 m_data_r=data;
118
119 // The above code looks like it will create a map the wrong way around
120
121 int valsize=DataTypes::noValues(shape);
122 int npoints=(data.size()/valsize)-1;
123 int ntags=tags.size();
124 if (ntags>npoints)
125 { // This throw is not unit tested yet
126 throw DataException("Programming error - Too many tags for the supplied values.");
127 }
128
129 // create the tag lookup map
130 // we assume that the first value is the default value so we skip it (hence the i+1 below)
131 for (int i=0;i<ntags;++i)
132 {
133 m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
134 }
135 }
136
137
138 DataTagged::DataTagged(const FunctionSpace& what,
139 const DataTypes::ShapeType &shape,
140 const TagListType& tags,
141 const DataTypes::CplxVectorType& data)
142 : parent(what,shape)
143 {
144 // alternative constructor
145 m_iscompl=true;
146 if (!what.canTag())
147 {
148 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
149 }
150
151 // copy the data
152 m_data_c=data;
153
154 // The above code looks like it will create a map the wrong way around
155
156 int valsize=DataTypes::noValues(shape);
157 int npoints=(data.size()/valsize)-1;
158 int ntags=tags.size();
159 if (ntags>npoints)
160 { // This throw is not unit tested yet
161 throw DataException("Programming error - Too many tags for the supplied values.");
162 }
163
164 // create the tag lookup map
165 // we assume that the first value is the default value so we skip it (hence the i+1 below)
166 for (int i=0;i<ntags;++i)
167 {
168 m_offsetLookup.insert(DataMapType::value_type(tags[i],(i+1)*valsize));
169 }
170 }
171
172
173
174 DataTagged::DataTagged(const DataTagged& other)
175 : parent(other.getFunctionSpace(),other.getShape()),
176 m_offsetLookup(other.m_offsetLookup),
177 m_data_r(other.m_data_r), m_data_c(other.m_data_c)
178 {
179 // copy constructor
180 m_iscompl=other.m_iscompl;
181 }
182
183 DataTagged::DataTagged(const DataConstant& other)
184 : parent(other.getFunctionSpace(),other.getShape())
185 {
186 // copy constructor
187 m_iscompl=other.isComplex();
188 if (!other.getFunctionSpace().canTag())
189 {
190 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
191 }
192
193 // fill the default value with the constant value item from "other"
194 int len = other.getNoValues();
195 if (m_iscompl)
196 {
197 DataTypes::cplx_t dummy=0;
198 m_data_c.resize(len,0.,len);
199 for (int i=0; i<len; i++) {
200 m_data_c[i]=other.getTypedVectorRO(dummy)[i];
201 }
202 }
203 else
204 {
205 DataTypes::real_t dummy=0;
206 m_data_r.resize(len,0.,len);
207 for (int i=0; i<len; i++) {
208 m_data_r[i]=other.getTypedVectorRO(dummy)[i];
209 }
210 }
211 }
212
213
214 // Create a new object by copying tags
215 DataTagged::DataTagged(const FunctionSpace& what,
216 const DataTypes::ShapeType& shape,
217 const DataTypes::RealVectorType& defaultvalue,
218 const DataTagged* tagsource)
219 : parent(what,shape)
220 {
221 // This constructor has not been unit tested yet
222
223 if (defaultvalue.size()!=DataTypes::noValues(shape)) {
224 throw DataException("Programming error - defaultvalue does not match supplied shape.");
225 }
226
227
228 if (!what.canTag())
229 {
230 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
231 }
232
233 if (tagsource!=0)
234 {
235 m_data_r.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
236
237 DataTagged::DataMapType::const_iterator i;
238 for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
239 addTag(i->first);
240 }
241 }
242 else
243 {
244 m_data_r.resize(defaultvalue.size());
245 }
246
247 // need to set the default value ....
248 for (int i=0; i<defaultvalue.size(); i++) {
249 m_data_r[i]=defaultvalue[i];
250 }
251 }
252
253 // Create a new object by copying tags
254 DataTagged::DataTagged(const FunctionSpace& what,
255 const DataTypes::ShapeType& shape,
256 const DataTypes::CplxVectorType& defaultvalue,
257 const DataTagged* tagsource)
258 : parent(what,shape)
259 {
260 // This constructor has not been unit tested yet
261 m_iscompl=true;
262
263 if (defaultvalue.size()!=DataTypes::noValues(shape)) {
264 throw DataException("Programming error - defaultvalue does not match supplied shape.");
265 }
266
267 if (!what.canTag())
268 {
269 throw DataException("Programming error - DataTag created with a non-taggable FunctionSpace.");
270 }
271
272 if (tagsource!=0)
273 {
274 m_data_r.resize(defaultvalue.size(),0.); // since this is tagged data, we should have blocksize=1
275
276 DataTagged::DataMapType::const_iterator i;
277 for (i=tagsource->getTagLookup().begin();i!=tagsource->getTagLookup().end();i++) {
278 addTag(i->first);
279 }
280 }
281 else
282 {
283 m_data_r.resize(defaultvalue.size());
284 }
285
286 // need to set the default value ....
287 for (int i=0; i<defaultvalue.size(); i++) {
288 m_data_c[i]=defaultvalue[i];
289 }
290 }
291
292
293 DataAbstract*
294 DataTagged::deepCopy() const
295 {
296 return new DataTagged(*this);
297 }
298
299 DataAbstract*
300 DataTagged::getSlice(const DataTypes::RegionType& region) const
301 {
302 return new DataTagged(*this, region);
303 }
304
305 DataTagged::DataTagged(const DataTagged& other,
306 const DataTypes::RegionType& region)
307 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
308 {
309 // slice constructor
310 m_iscompl=other.isComplex();
311
312
313 // get the shape of the slice to copy from other
314 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
315 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
316
317 // allocate enough space in this for all values
318 // (need to add one to allow for the default value)
319 int len = DataTypes::noValues(regionShape)*(other.m_offsetLookup.size()+1);
320 if (m_iscompl)
321 {
322 m_data_c.resize(len,0.0,len);
323 // copy the default value from other to this
324 const DataTypes::ShapeType& otherShape=other.getShape();
325 const DataTypes::CplxVectorType& otherData=other.getTypedVectorRO((DataTypes::cplx_t)0);
326 DataTypes::copySlice(getTypedVectorRW((DataTypes::cplx_t)0),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
327
328 // loop through the tag values copying these
329 DataMapType::const_iterator pos;
330 DataTypes::CplxVectorType::size_type tagOffset=getNoValues();
331 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
332 DataTypes::copySlice(m_data_c,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
333 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
334 tagOffset+=getNoValues();
335 }
336
337
338 }
339 else
340 {
341 m_data_r.resize(len,0.0,len);
342 // copy the default value from other to this
343 const DataTypes::ShapeType& otherShape=other.getShape();
344 const DataTypes::RealVectorType& otherData=other.getTypedVectorRO((DataTypes::real_t)0);
345 DataTypes::copySlice(getTypedVectorRW((DataTypes::real_t)0),getShape(),getDefaultOffset(),otherData,otherShape,other.getDefaultOffset(), regionLoopRange);
346
347 // loop through the tag values copying these
348 DataMapType::const_iterator pos;
349 DataTypes::RealVectorType::size_type tagOffset=getNoValues();
350 for (pos=other.m_offsetLookup.begin();pos!=other.m_offsetLookup.end();pos++){
351 DataTypes::copySlice(m_data_r,getShape(),tagOffset,otherData, otherShape, pos->second, regionLoopRange);
352 m_offsetLookup.insert(DataMapType::value_type(pos->first,tagOffset));
353 tagOffset+=getNoValues();
354 }
355 }
356 }
357
358 void
359 DataTagged::setSlice(const DataAbstract* other,
360 const DataTypes::RegionType& region)
361 {
362
363 // other must be another DataTagged object
364 // Data:setSlice implementation should ensure this
365 const DataTagged* otherTemp=dynamic_cast<const DataTagged*>(other);
366 if (otherTemp==0) {
367 throw DataException("Programming error - casting to DataTagged.");
368 }
369 if (isComplex()!=other->isComplex())
370 {
371 throw DataException("Error - cannot copy between slices of different complexity.");
372 }
373 CHECK_FOR_EX_WRITE
374
375
376
377 // determine shape of the specified region
378 DataTypes::ShapeType regionShape(DataTypes::getResultSliceShape(region));
379
380 // modify region specification as needed to match rank of this object
381 DataTypes::RegionLoopRangeType regionLoopRange=DataTypes::getSliceRegionLoopRange(region);
382
383 // ensure rank/shape of this object is compatible with specified region
384 if (getRank()!=region.size()) {
385 throw DataException("Error - Invalid slice region.");
386 }
387 if (otherTemp->getRank()>0 && !DataTypes::checkShape(other->getShape(),regionShape)) {
388 throw DataException (DataTypes::createShapeErrorMessage(
389 "Error - Couldn't copy slice due to shape mismatch.",regionShape,other->getShape()));
390 }
391
392
393 const DataTypes::ShapeType& otherShape=otherTemp->getShape();
394 if (isComplex()) // from check earlier, other will have the same complexity
395 {
396 // copy slice from other default value to this default value
397 DataTypes::copySliceFrom(m_data_c,getShape(),getDefaultOffset(),otherTemp->getTypedVectorRO((DataTypes::cplx_t)0),
398 otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
399 }
400 else
401 {
402 // copy slice from other default value to this default value
403 DataTypes::copySliceFrom(m_data_r,getShape(),getDefaultOffset(),otherTemp->getTypedVectorRO((DataTypes::real_t)0),
404 otherShape,otherTemp->getDefaultOffset(),regionLoopRange);
405 }
406
407 // loop through tag values in other, adding any which aren't in this, using default value
408 DataMapType::const_iterator pos;
409 for (pos=otherTemp->m_offsetLookup.begin();pos!=otherTemp->m_offsetLookup.end();pos++) {
410 if (!isCurrentTag(pos->first)) {
411 addTag(pos->first);
412 }
413 }
414 if (isComplex())
415 {
416 // loop through the tag values copying slices from other to this
417 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
418 DataTypes::copySliceFrom(m_data_c,getShape(),getOffsetForTag(pos->first),otherTemp->getTypedVectorRO((DataTypes::cplx_t)0),
419 otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
420 }
421 }
422 else
423 {
424 // loop through the tag values copying slices from other to this
425 for (pos=m_offsetLookup.begin();pos!=m_offsetLookup.end();pos++) {
426 DataTypes::copySliceFrom(m_data_r,getShape(),getOffsetForTag(pos->first),otherTemp->getTypedVectorRO((DataTypes::real_t)0),
427 otherShape, otherTemp->getOffsetForTag(pos->first), regionLoopRange);
428 }
429 }
430 }
431
432 int
433 DataTagged::getTagNumber(int dpno)
434 {
435 //
436 // Get the number of samples and data-points per sample
437 int numSamples = getNumSamples();
438 int numDataPointsPerSample = getNumDPPSample();
439 int numDataPoints = numSamples * numDataPointsPerSample;
440
441 if (numDataPointsPerSample==0) {
442 throw DataException("DataTagged::getTagNumber error: no data-points associated with this object.");
443 }
444
445 if (dpno<0 || dpno>numDataPoints-1) {
446 throw DataException("DataTagged::getTagNumber error: invalid data-point number supplied.");
447 }
448
449 //
450 // Determine the sample number which corresponds to this data-point number
451 int sampleNo = dpno / numDataPointsPerSample;
452
453 //
454 // Determine the tag number which corresponds to this sample number
455 int tagNo = getFunctionSpace().getTagFromSampleNo(sampleNo);
456
457 //
458 // return the tag number
459 return(tagNo);
460 }
461
462 void
463 DataTagged::setTaggedValue(int tagKey,
464 const DataTypes::ShapeType& pointshape,
465 const DataTypes::RealVectorType& value,
466 int dataOffset)
467 {
468 if (!DataTypes::checkShape(getShape(), pointshape)) {
469 throw DataException(DataTypes::createShapeErrorMessage(
470 "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
471 }
472 if (isComplex())
473 {
474 throw DataException("Programming Error - attempt to set real value on complex data.");
475 }
476 CHECK_FOR_EX_WRITE
477 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
478 if (pos==m_offsetLookup.end()) {
479 // tag couldn't be found so use addTaggedValue
480 addTaggedValue(tagKey,pointshape, value, dataOffset);
481 } else {
482 // copy the values into the data array at the offset determined by m_offsetLookup
483 int offset=pos->second;
484 for (unsigned int i=0; i<getNoValues(); i++) {
485 m_data_r[offset+i]=value[i+dataOffset];
486 }
487 }
488 }
489
490 void
491 DataTagged::setTaggedValue(int tagKey,
492 const DataTypes::ShapeType& pointshape,
493 const DataTypes::CplxVectorType& value,
494 int dataOffset)
495 {
496 if (!DataTypes::checkShape(getShape(), pointshape)) {
497 throw DataException(DataTypes::createShapeErrorMessage(
498 "Error - Cannot setTaggedValue due to shape mismatch.", pointshape,getShape()));
499 }
500 if (!isComplex())
501 {
502 throw DataException("Programming Error - attempt to set a complex value on real data");
503 }
504 CHECK_FOR_EX_WRITE
505 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
506 if (pos==m_offsetLookup.end()) {
507 // tag couldn't be found so use addTaggedValue
508 addTaggedValue(tagKey,pointshape, value, dataOffset);
509 } else {
510 // copy the values into the data array at the offset determined by m_offsetLookup
511 int offset=pos->second;
512 for (unsigned int i=0; i<getNoValues(); i++) {
513 m_data_c[offset+i]=value[i+dataOffset];
514 }
515 }
516 }
517
518
519 void
520 DataTagged::addTaggedValues(const TagListType& tagKeys,
521 const FloatBatchType& values,
522 const ShapeType& vShape)
523 {
524 DataTypes::RealVectorType t(values.size(),0);
525 for (size_t i=0;i<values.size();++i)
526 {
527 t[i]=values[i];
528 }
529 addTaggedValues(tagKeys,t,vShape);
530 }
531
532
533 // Note: The check to see if vShape==our shape is done in the addTaggedValue method
534 void
535 DataTagged::addTaggedValues(const TagListType& tagKeys,
536 const DataTypes::RealVectorType& values,
537 const ShapeType& vShape)
538 {
539 unsigned int n=getNoValues();
540 unsigned int numVals=values.size()/n;
541 if (values.size()==0) {
542 // copy the current default value for each of the tags
543 TagListType::const_iterator iT;
544 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
545 // the point data view for DataTagged points at the default value
546 addTag(*iT);
547 }
548 } else if (numVals==1 && tagKeys.size()>1) {
549 // assume the one given value will be used for all tag values
550 TagListType::const_iterator iT;
551 for (iT=tagKeys.begin();iT!=tagKeys.end();iT++) {
552 addTaggedValue(*iT, vShape, values,0);
553 }
554 } else {
555 if (tagKeys.size()!=numVals) {
556 stringstream temp;
557 temp << "Error - (addTaggedValue) Number of tags: " << tagKeys.size()
558 << " doesn't match number of values: " << values.size();
559 throw DataException(temp.str());
560 } else {
561 unsigned int i;
562 int offset=0;
563 for (i=0;i<tagKeys.size();i++ ,offset+=n) {
564 addTaggedValue(tagKeys[i],vShape,values,offset);
565 }
566 }
567 }
568 }
569
570
571
572
573 void
574 DataTagged::addTaggedValue(int tagKey,
575 const DataTypes::ShapeType& pointshape,
576 const DataTypes::RealVectorType& value,
577 int dataOffset)
578 {
579 if (!DataTypes::checkShape(getShape(), pointshape)) {
580 throw DataException(DataTypes::createShapeErrorMessage(
581 "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
582 }
583 if (isComplex())
584 {
585 throw DataException("Programming Error - attempt to set a real value on complex data");
586 }
587 CHECK_FOR_EX_WRITE
588 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
589 if (pos!=m_offsetLookup.end()) {
590 // tag already exists so use setTaggedValue
591 setTaggedValue(tagKey,pointshape, value, dataOffset);
592 } else {
593 // save the key and the location of its data in the lookup tab
594 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data_r.size()));
595 // add the data given in "value" at the end of m_data_r
596 // need to make a temp copy of m_data_r, resize m_data_r, then copy
597 // all the old values plus the value to be added back into m_data_r
598 DataTypes::RealVectorType m_data_r_temp(m_data_r);
599 int oldSize=m_data_r.size();
600 int newSize=m_data_r.size()+getNoValues();
601 m_data_r.resize(newSize,0.,newSize);
602 for (int i=0;i<oldSize;i++) {
603 m_data_r[i]=m_data_r_temp[i];
604 }
605 for (unsigned int i=0;i<getNoValues();i++) {
606 m_data_r[oldSize+i]=value[i+dataOffset];
607 }
608 }
609 }
610
611
612 void
613 DataTagged::addTaggedValue(int tagKey,
614 const DataTypes::ShapeType& pointshape,
615 const DataTypes::CplxVectorType& value,
616 int dataOffset)
617 {
618 if (!DataTypes::checkShape(getShape(), pointshape)) {
619 throw DataException(DataTypes::createShapeErrorMessage(
620 "Error - Cannot addTaggedValue due to shape mismatch.", pointshape,getShape()));
621 }
622 if (!isComplex())
623 {
624 throw DataException("Programming error - attempt to set a complex value on real data.");
625 }
626 CHECK_FOR_EX_WRITE
627 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
628 if (pos!=m_offsetLookup.end()) {
629 // tag already exists so use setTaggedValue
630 setTaggedValue(tagKey,pointshape, value, dataOffset);
631 } else {
632 // save the key and the location of its data in the lookup tab
633 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data_c.size()));
634 // add the data given in "value" at the end of m_data_c
635 // need to make a temp copy of m_data_c, resize m_data_c, then copy
636 // all the old values plus the value to be added back into m_data_c
637 DataTypes::CplxVectorType m_data_c_temp(m_data_c);
638 int oldSize=m_data_c.size();
639 int newSize=m_data_c.size()+getNoValues();
640 m_data_c.resize(newSize,0.,newSize);
641 for (int i=0;i<oldSize;i++) {
642 m_data_c[i]=m_data_c_temp[i];
643 }
644 for (unsigned int i=0;i<getNoValues();i++) {
645 m_data_c[oldSize+i]=value[i+dataOffset];
646 }
647 }
648 }
649
650 void
651 DataTagged::addTag(int tagKey)
652 {
653 CHECK_FOR_EX_WRITE
654 DataMapType::iterator pos(m_offsetLookup.find(tagKey));
655 if (pos==m_offsetLookup.end()) {
656 if (isComplex())
657 {
658 // save the key and the location of its data in the lookup tab
659 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data_c.size()));
660 // add the data given in "value" at the end of m_data_c
661 // need to make a temp copy of m_data_c, resize m_data_c, then copy
662 // all the old values plus the value to be added back into m_data_c
663 DataTypes::CplxVectorType m_data_c_temp(m_data_c);
664 int oldSize=m_data_c.size();
665 int newSize=m_data_c.size()+getNoValues();
666 m_data_c.resize(newSize,0.,newSize);
667 for (int i=0;i<oldSize;i++) {
668 m_data_c[i]=m_data_c_temp[i];
669 }
670 for (unsigned int i=0;i<getNoValues();i++) {
671 m_data_c[oldSize+i]=m_data_c[m_defaultValueOffset+i];
672 }
673 }
674 else
675 {
676 // save the key and the location of its data in the lookup tab
677 m_offsetLookup.insert(DataMapType::value_type(tagKey,m_data_r.size()));
678 // add the data given in "value" at the end of m_data_r
679 // need to make a temp copy of m_data_r, resize m_data_r, then copy
680 // all the old values plus the value to be added back into m_data_r
681 DataTypes::RealVectorType m_data_r_temp(m_data_r);
682 int oldSize=m_data_r.size();
683 int newSize=m_data_r.size()+getNoValues();
684 m_data_r.resize(newSize,0.,newSize);
685 for (int i=0;i<oldSize;i++) {
686 m_data_r[i]=m_data_r_temp[i];
687 }
688 for (unsigned int i=0;i<getNoValues();i++) {
689 m_data_r[oldSize+i]=m_data_r[m_defaultValueOffset+i];
690 }
691 }
692 }
693 }
694
695
696 DataTypes::real_t*
697 DataTagged::getSampleDataByTag(int tag, DataTypes::real_t dummy)
698 {
699 CHECK_FOR_EX_WRITE
700 DataMapType::iterator pos(m_offsetLookup.find(tag));
701 if (pos==m_offsetLookup.end()) {
702 // tag couldn't be found so return the default value
703 return &(m_data_r[0]);
704 } else {
705 // return the data-point corresponding to the given tag
706 return &(m_data_r[pos->second]);
707 }
708 }
709
710 DataTypes::cplx_t*
711 DataTagged::getSampleDataByTag(int tag, DataTypes::cplx_t dummy)
712 {
713 CHECK_FOR_EX_WRITE
714 DataMapType::iterator pos(m_offsetLookup.find(tag));
715 if (pos==m_offsetLookup.end()) {
716 // tag couldn't be found so return the default value
717 return &(m_data_c[0]);
718 } else {
719 // return the data-point corresponding to the given tag
720 return &(m_data_c[pos->second]);
721 }
722 }
723
724
725 bool
726 DataTagged::hasNaN() const
727 {
728 bool haveNaN=false;
729 if (isComplex())
730 {
731 #pragma omp parallel for
732 for (DataTypes::CplxVectorType::size_type i=0;i<m_data_c.size();++i)
733 {
734 if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))
735 {
736 #pragma omp critical
737 {
738 haveNaN=true;
739 }
740 }
741 }
742 }
743 else
744 {
745 #pragma omp parallel for
746 for (DataTypes::RealVectorType::size_type i=0;i<m_data_r.size();++i)
747 {
748 if (std::isnan(m_data_r[i]))
749 {
750 #pragma omp critical
751 {
752 haveNaN=true;
753 }
754 }
755 }
756 }
757 return haveNaN;
758 }
759
760 void
761 DataTagged::replaceNaN(double value) {
762 CHECK_FOR_EX_WRITE
763 if (isComplex())
764 {
765 #pragma omp parallel for
766 for (DataTypes::CplxVectorType::size_type i=0;i<m_data_c.size();++i)
767 {
768 if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))
769 {
770 m_data_c[i] = value;
771 }
772 }
773 }
774 else
775 {
776 #pragma omp parallel for
777 for (DataTypes::RealVectorType::size_type i=0;i<m_data_r.size();++i)
778 {
779 if (std::isnan(m_data_r[i]))
780 {
781 m_data_r[i] = value;
782 }
783 }
784 }
785 }
786
787 void
788 DataTagged::replaceNaN(DataTypes::cplx_t value) {
789 CHECK_FOR_EX_WRITE
790 if (isComplex())
791 {
792 #pragma omp parallel for
793 for (DataTypes::CplxVectorType::size_type i=0;i<m_data_c.size();++i)
794 {
795 if (std::isnan(m_data_c[i].real()) || std::isnan(m_data_c[i].imag()))
796 {
797 m_data_c[i] = value;
798 }
799 }
800 }
801 else
802 {
803 complicate();
804 replaceNaN(value);
805 }
806 }
807
808
809 string
810 DataTagged::toString() const
811 {
812 using namespace escript::DataTypes;
813 string empty="";
814 stringstream temp;
815 DataMapType::const_iterator i;
816 temp << "Tag(Default)" << endl;
817
818 if (isComplex())
819 {
820
821 temp << pointToString(m_data_c,getShape(),getDefaultOffset(),empty) << endl;
822 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
823 temp << "Tag(" << i->first << ")" << endl;
824 temp << pointToString(m_data_c,getShape(),i->second,empty) << endl;
825 }
826 }
827 else
828 {
829
830 temp << pointToString(m_data_r,getShape(),getDefaultOffset(),empty) << endl;
831 for (i=m_offsetLookup.begin();i!=m_offsetLookup.end();++i) {
832 temp << "Tag(" << i->first << ")" << endl;
833 temp << pointToString(m_data_r,getShape(),i->second,empty) << endl;
834 }
835 }
836 return temp.str();
837 }
838
839 DataTypes::RealVectorType::size_type
840 DataTagged::getPointOffset(int sampleNo,
841 int dataPointNo) const
842 {
843 int tagKey=getFunctionSpace().getTagFromSampleNo(sampleNo);
844 DataMapType::const_iterator pos(m_offsetLookup.find(tagKey));
845 DataTypes::RealVectorType::size_type offset=m_defaultValueOffset;
846 if (pos!=m_offsetLookup.end()) {
847 offset=pos->second;
848 }
849 return offset;
850 }
851
852 DataTypes::RealVectorType::size_type
853 DataTagged::getOffsetForTag(int tag) const
854 {
855 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
856 DataTypes::RealVectorType::size_type offset=m_defaultValueOffset;
857 if (pos!=m_offsetLookup.end()) {
858 offset=pos->second;
859 }
860 return offset;
861 }
862
863 DataTypes::RealVectorType::const_reference
864 DataTagged::getDataByTagRO(int tag, DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy) const
865 {
866 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
867 DataTypes::RealVectorType::size_type offset=m_defaultValueOffset;
868 if (pos!=m_offsetLookup.end()) {
869 offset=pos->second;
870 }
871 return m_data_r[offset+i];
872 }
873
874 DataTypes::RealVectorType::reference
875 DataTagged::getDataByTagRW(int tag, DataTypes::RealVectorType::size_type i, DataTypes::real_t dummy)
876 {
877 CHECK_FOR_EX_WRITE
878 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
879 DataTypes::RealVectorType::size_type offset=m_defaultValueOffset;
880 if (pos!=m_offsetLookup.end()) {
881 offset=pos->second;
882 }
883 return m_data_r[offset+i];
884 }
885
886 DataTypes::CplxVectorType::const_reference
887 DataTagged::getDataByTagRO(int tag, DataTypes::RealVectorType::size_type i, DataTypes::cplx_t dummy) const
888 {
889 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
890 DataTypes::CplxVectorType::size_type offset=m_defaultValueOffset;
891 if (pos!=m_offsetLookup.end()) {
892 offset=pos->second;
893 }
894 return m_data_c[offset+i];
895 }
896
897 DataTypes::CplxVectorType::reference
898 DataTagged::getDataByTagRW(int tag, DataTypes::RealVectorType::size_type i, DataTypes::cplx_t dummy)
899 {
900 CHECK_FOR_EX_WRITE
901 DataMapType::const_iterator pos(m_offsetLookup.find(tag));
902 DataTypes::CplxVectorType::size_type offset=m_defaultValueOffset;
903 if (pos!=m_offsetLookup.end()) {
904 offset=pos->second;
905 }
906 return m_data_c[offset+i];
907 }
908
909 void
910 DataTagged::symmetric(DataAbstract* ev)
911 {
912 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
913 if (temp_ev==0) {
914 throw DataException("Error - DataTagged::symmetric casting to DataTagged failed (probably a programming error).");
915 }
916 const DataTagged::DataMapType& thisLookup=getTagLookup();
917 DataTagged::DataMapType::const_iterator i;
918 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
919 const ShapeType& evShape=temp_ev->getShape();
920
921 if (isComplex())
922 {
923 DataTypes::CplxVectorType& evVec=temp_ev->getTypedVectorRW((DataTypes::cplx_t)0);
924 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
925 temp_ev->addTag(i->first);
926 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
927 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
928 escript::symmetric(m_data_c,getShape(),offset,evVec, evShape, evoffset);
929 }
930 escript::symmetric(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
931 }
932 else
933 {
934 DataTypes::RealVectorType& evVec=temp_ev->getTypedVectorRW(0.0);
935 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
936 temp_ev->addTag(i->first);
937 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
938 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
939 escript::symmetric(m_data_r,getShape(),offset,evVec, evShape, evoffset);
940 }
941 escript::symmetric(m_data_r,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
942 }
943 }
944
945
946 void
947 DataTagged::antisymmetric(DataAbstract* ev)
948 {
949 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
950 if (temp_ev==0) {
951 throw DataException("Error - DataTagged::antisymmetric casting to DataTagged failed (probably a programming error).");
952 }
953 const DataTagged::DataMapType& thisLookup=getTagLookup();
954 DataTagged::DataMapType::const_iterator i;
955 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
956 const ShapeType& evShape=temp_ev->getShape();
957 if (isComplex())
958 {
959 DataTypes::CplxVectorType& evVec=temp_ev->getTypedVectorRW((DataTypes::cplx_t)0);
960 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
961 temp_ev->addTag(i->first);
962 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
963 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
964 escript::antisymmetric(m_data_c,getShape(),offset,evVec, evShape, evoffset);
965 }
966 escript::antisymmetric(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
967 }
968 else
969 {
970 DataTypes::RealVectorType& evVec=temp_ev->getTypedVectorRW(0.0);
971 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
972 temp_ev->addTag(i->first);
973 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
974 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
975 escript::antisymmetric(m_data_r,getShape(),offset,evVec, evShape, evoffset);
976 }
977 escript::antisymmetric(m_data_r,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
978 }
979 }
980
981 void
982 DataTagged::hermitian(DataAbstract* ev)
983 {
984 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
985 if (temp_ev==0) {
986 throw DataException("Error - DataTagged::hermitian casting to DataTagged failed (probably a programming error).");
987 }
988 if (!isComplex() || !temp_ev->isComplex())
989 {
990 throw DataException("DataTagged::hermitian: do not call this method with real data");
991 }
992
993 const DataTagged::DataMapType& thisLookup=getTagLookup();
994 DataTagged::DataMapType::const_iterator i;
995 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
996 const ShapeType& evShape=temp_ev->getShape();
997
998 DataTypes::CplxVectorType& evVec=temp_ev->getTypedVectorRW((DataTypes::cplx_t)0);
999 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1000 temp_ev->addTag(i->first);
1001 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1002 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1003 escript::hermitian(m_data_c,getShape(),offset,evVec, evShape, evoffset);
1004 }
1005 escript::hermitian(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
1006 }
1007
1008
1009 void
1010 DataTagged::antihermitian(DataAbstract* ev)
1011 {
1012 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1013 if (temp_ev==0) {
1014 throw DataException("Error - DataTagged::antihermitian casting to DataTagged failed (probably a programming error).");
1015 }
1016 if (!isComplex() || !temp_ev->isComplex())
1017 {
1018 throw DataException("DataTagged::antihermitian: do not call this method with real data");
1019 }
1020 const DataTagged::DataMapType& thisLookup=getTagLookup();
1021 DataTagged::DataMapType::const_iterator i;
1022 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1023 const ShapeType& evShape=temp_ev->getShape();
1024 DataTypes::CplxVectorType& evVec=temp_ev->getTypedVectorRW((DataTypes::cplx_t)0);
1025 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1026 temp_ev->addTag(i->first);
1027 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1028 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1029 escript::antihermitian(m_data_c,getShape(),offset,evVec, evShape, evoffset);
1030 }
1031 escript::antihermitian(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset());
1032 }
1033
1034
1035 void
1036 DataTagged::trace(DataAbstract* ev, int axis_offset)
1037 {
1038 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1039 if (temp_ev==0) {
1040 throw DataException("Error - DataTagged::trace casting to DataTagged failed (probably a programming error).");
1041 }
1042 const DataTagged::DataMapType& thisLookup=getTagLookup();
1043 DataTagged::DataMapType::const_iterator i;
1044 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1045 const ShapeType& evShape=temp_ev->getShape();
1046 if (isComplex())
1047 {
1048 DataTypes::CplxVectorType& evVec=temp_ev->getVectorRWC();
1049 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1050 temp_ev->addTag(i->first);
1051 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1052 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1053 escript::trace(m_data_c,getShape(),offset,evVec, evShape, evoffset, axis_offset);
1054 }
1055 escript::trace(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
1056 }
1057 else
1058 {
1059 DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
1060 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1061 temp_ev->addTag(i->first);
1062 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
1063 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1064 escript::trace(m_data_r,getShape(),offset,evVec, evShape, evoffset, axis_offset);
1065 }
1066 escript::trace(m_data_r,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
1067 }
1068 }
1069
1070 void
1071 DataTagged::transpose(DataAbstract* ev, int axis_offset)
1072 {
1073 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1074 if (temp_ev==0) {
1075 throw DataException("Error - DataTagged::transpose casting to DataTagged failed (probably a programming error).");
1076 }
1077 const DataTagged::DataMapType& thisLookup=getTagLookup();
1078 DataTagged::DataMapType::const_iterator i;
1079 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1080 const ShapeType& evShape=temp_ev->getShape();
1081 if (isComplex())
1082 {
1083 DataTypes::CplxVectorType& evVec=temp_ev->getVectorRWC();
1084 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1085 temp_ev->addTag(i->first);
1086 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1087 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1088 escript::transpose(m_data_c,getShape(),offset,evVec, evShape, evoffset, axis_offset);
1089 }
1090 escript::transpose(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
1091 }
1092 else
1093 {
1094 DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
1095 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1096 temp_ev->addTag(i->first);
1097 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
1098 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1099 escript::transpose(m_data_r,getShape(),offset,evVec, evShape, evoffset, axis_offset);
1100 }
1101 escript::transpose(m_data_r,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis_offset);
1102 }
1103 }
1104
1105 void
1106 DataTagged::swapaxes(DataAbstract* ev, int axis0, int axis1)
1107 {
1108 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1109 if (temp_ev==0) {
1110 throw DataException("Error - DataTagged::swapaxes casting to DataTagged failed (probably a programming error).");
1111 }
1112 const DataTagged::DataMapType& thisLookup=getTagLookup();
1113 DataTagged::DataMapType::const_iterator i;
1114 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1115 const ShapeType& evShape=temp_ev->getShape();
1116 if (isComplex())
1117 {
1118 DataTypes::CplxVectorType& evVec=temp_ev->getVectorRWC();
1119 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1120 temp_ev->addTag(i->first);
1121 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1122 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1123 escript::swapaxes(m_data_c,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
1124 }
1125 escript::swapaxes(m_data_c,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
1126 }
1127 else
1128 {
1129 DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
1130 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1131 temp_ev->addTag(i->first);
1132 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
1133 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1134 escript::swapaxes(m_data_r,getShape(),offset,evVec, evShape, evoffset,axis0,axis1);
1135 }
1136 escript::swapaxes(m_data_r,getShape(),getDefaultOffset(),evVec,evShape,temp_ev->getDefaultOffset(),axis0,axis1);
1137 }
1138 }
1139
1140 void
1141 DataTagged::eigenvalues(DataAbstract* ev)
1142 {
1143 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1144 if (temp_ev==0) {
1145 throw DataException("Error - DataTagged::eigenvalues casting to DataTagged failed (probably a programming error).");
1146 }
1147 const DataTagged::DataMapType& thisLookup=getTagLookup();
1148 DataTagged::DataMapType::const_iterator i;
1149 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1150 const ShapeType& evShape=temp_ev->getShape();
1151 if (isComplex())
1152 {
1153 DataTypes::CplxVectorType& evVec=temp_ev->getVectorRWC();
1154
1155 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1156 temp_ev->addTag(i->first);
1157 DataTypes::CplxVectorType::size_type offset=getOffsetForTag(i->first);
1158 DataTypes::CplxVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1159 escript::eigenvalues(m_data_c,getShape(),offset,evVec, evShape, evoffset);
1160 }
1161 escript::eigenvalues(m_data_c,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
1162 }
1163 else
1164 {
1165 DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
1166
1167 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1168 temp_ev->addTag(i->first);
1169 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
1170 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1171 escript::eigenvalues(m_data_r,getShape(),offset,evVec, evShape, evoffset);
1172 }
1173 escript::eigenvalues(m_data_r,getShape(),getDefaultOffset(),evVec, evShape, temp_ev->getDefaultOffset());
1174 }
1175 }
1176 void
1177 DataTagged::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
1178 {
1179 DataTagged* temp_ev=dynamic_cast<DataTagged*>(ev);
1180 if (temp_ev==0) {
1181 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (probably a programming error).");
1182 }
1183 DataTagged* temp_V=dynamic_cast<DataTagged*>(V);
1184 if (temp_V==0) {
1185 throw DataException("Error - DataTagged::eigenvalues_and_eigenvectors casting to DataTagged failed (probably a programming error).");
1186 }
1187 const DataTagged::DataMapType& thisLookup=getTagLookup();
1188 DataTagged::DataMapType::const_iterator i;
1189 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1190 DataTypes::RealVectorType& evVec=temp_ev->getVectorRW();
1191 const ShapeType& evShape=temp_ev->getShape();
1192 DataTypes::RealVectorType& VVec=temp_V->getVectorRW();
1193 const ShapeType& VShape=temp_V->getShape();
1194 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1195 temp_ev->addTag(i->first);
1196 temp_V->addTag(i->first);
1197 /* DataArrayView thisView=getDataPointByTag(i->first);
1198 DataArrayView evView=temp_ev->getDataPointByTag(i->first);
1199 DataArrayView VView=temp_V->getDataPointByTag(i->first);*/
1200 DataTypes::RealVectorType::size_type offset=getOffsetForTag(i->first);
1201 DataTypes::RealVectorType::size_type evoffset=temp_ev->getOffsetForTag(i->first);
1202 DataTypes::RealVectorType::size_type Voffset=temp_V->getOffsetForTag(i->first);
1203 /* DataArrayView::eigenvalues_and_eigenvectors(thisView,0,evView,0,VView,0,tol);*/
1204 escript::eigenvalues_and_eigenvectors(m_data_r,getShape(),offset,evVec, evShape, evoffset,VVec,VShape,Voffset,tol);
1205
1206 }
1207 escript::eigenvalues_and_eigenvectors(m_data_r,getShape(),getDefaultOffset(),evVec, evShape,
1208 temp_ev->getDefaultOffset(),VVec,VShape,
1209 temp_V->getDefaultOffset(), tol);
1210
1211
1212 }
1213
1214 int
1215 DataTagged::matrixInverse(DataAbstract* out) const
1216 {
1217 DataTagged* temp=dynamic_cast<DataTagged*>(out);
1218 if (temp==0)
1219 {
1220 throw DataException("Error - DataTagged::matrixInverse: casting to DataTagged failed (probably a programming error).");
1221 }
1222 if (getRank()!=2)
1223 {
1224 throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
1225 }
1226 const DataTagged::DataMapType& thisLookup=getTagLookup();
1227 DataTagged::DataMapType::const_iterator i;
1228 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1229 DataTypes::RealVectorType& outVec=temp->getVectorRW();
1230 const ShapeType& outShape=temp->getShape();
1231 LapackInverseHelper h(getShape()[0]);
1232 int err=0;
1233 for (i=thisLookup.begin();i!=thisLookupEnd;i++) {
1234 temp->addTag(i->first);
1235 DataTypes::RealVectorType::size_type inoffset=getOffsetForTag(i->first);
1236 DataTypes::RealVectorType::size_type outoffset=temp->getOffsetForTag(i->first);
1237
1238 err=escript::matrix_inverse(m_data_r, getShape(), inoffset, outVec, outShape, outoffset, 1, h);
1239 if (!err) break;
1240 }
1241 if (!err)
1242 {
1243 escript::matrix_inverse(m_data_r, getShape(), getDefaultOffset(), outVec, outShape, temp->getDefaultOffset(), 1, h);
1244 }
1245 return err;
1246 }
1247
1248 void
1249 DataTagged::setToZero(){
1250 CHECK_FOR_EX_WRITE
1251 DataTypes::RealVectorType::size_type n=m_data_r.size();
1252 for (int i=0; i<n ;++i) m_data_r[i]=0.;
1253 }
1254
1255 void
1256 DataTagged::dump(const std::string fileName) const
1257 {
1258 #ifdef ESYS_HAVE_NETCDF
1259 const int ldims=DataTypes::maxRank+1;
1260 const NcDim* ncdims[ldims];
1261 NcVar *var, *tags_var;
1262 int rank = getRank();
1263 int type= getFunctionSpace().getTypeCode();
1264 int ndims =0;
1265 long dims[ldims];
1266 const double* d_ptr=&(m_data_r[0]);
1267 DataTypes::ShapeType shape = getShape();
1268 JMPI mpiInfo(getFunctionSpace().getDomain()->getMPI());
1269 #ifdef ESYS_MPI
1270 const int mpi_iam = mpiInfo->rank;
1271 const int mpi_num = mpiInfo->size;
1272 MPI_Status status;
1273
1274 /* Serialize NetCDF I/O */
1275 if (mpi_iam > 0)
1276 MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81803, mpiInfo->comm, &status);
1277 #endif
1278
1279 // netCDF error handler
1280 NcError err(NcError::verbose_nonfatal);
1281 // Create the file.
1282 const std::string newFileName(mpiInfo->appendRankToFileName(fileName));
1283 NcFile dataFile(newFileName.c_str(), NcFile::Replace);
1284 // check if writing was successful
1285 if (!dataFile.is_valid())
1286 throw DataException("Error - DataTagged:: opening of netCDF file for output failed.");
1287 if (!dataFile.add_att("type_id",1) )
1288 throw DataException("Error - DataTagged:: appending data type to netCDF file failed.");
1289 if (!dataFile.add_att("rank",rank) )
1290 throw DataException("Error - DataTagged:: appending rank attribute to netCDF file failed.");
1291 if (!dataFile.add_att("function_space_type",type))
1292 throw DataException("Error - DataTagged:: appending function space attribute to netCDF file failed.");
1293 ndims=rank+1;
1294 if ( rank >0 ) {
1295 dims[0]=shape[0];
1296 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
1297 throw DataException("Error - DataTagged:: appending ncdimension 0 to netCDF file failed.");
1298 }
1299 if ( rank >1 ) {
1300 dims[1]=shape[1];
1301 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
1302 throw DataException("Error - DataTagged:: appending ncdimension 1 to netCDF file failed.");
1303 }
1304 if ( rank >2 ) {
1305 dims[2]=shape[2];
1306 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
1307 throw DataException("Error - DataTagged:: appending ncdimension 2 to netCDF file failed.");
1308 }
1309 if ( rank >3 ) {
1310 dims[3]=shape[3];
1311 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
1312 throw DataException("Error - DataTagged:: appending ncdimension 3 to netCDF file failed.");
1313 }
1314 const DataTagged::DataMapType& thisLookup=getTagLookup();
1315 DataTagged::DataMapType::const_iterator i;
1316 DataTagged::DataMapType::const_iterator thisLookupEnd=thisLookup.end();
1317 std::vector<int> tags;
1318 tags.push_back(-1);
1319 for (i=thisLookup.begin();i!=thisLookupEnd;i++)
1320 tags.push_back(i->first);
1321 dims[rank]=tags.size();
1322 if (! (ncdims[rank] = dataFile.add_dim("num_tags", dims[rank])) )
1323 {
1324 throw DataException("Error - DataTagged:: appending num_tags to netCDF file failed.");
1325 }
1326 if (! ( tags_var = dataFile.add_var("tags", ncInt, ncdims[rank])) )
1327 {
1328 throw DataException("Error - DataTagged:: appending tags to netCDF file failed.");
1329 }
1330 if (! (tags_var->put(&tags[0], dims[rank])) )
1331 {
1332 throw DataException("Error - DataTagged:: copy tags to netCDF buffer failed.");
1333 }
1334 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
1335 {
1336 throw DataException("Error - DataTagged:: appending variable to netCDF file failed.");
1337 }
1338 if (! (var->put(d_ptr,dims)) )
1339 {
1340 throw DataException("Error - DataTagged:: copy data to netCDF buffer failed.");
1341 }
1342 #ifdef ESYS_MPI
1343 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81803, MPI_COMM_WORLD);
1344 #endif
1345 #else
1346 throw DataException("Error - DataTagged:: dump is not configured with netCDF. Please contact your installation manager.");
1347 #endif
1348 }
1349
1350 DataTypes::RealVectorType&
1351 DataTagged::getVectorRW()
1352 {
1353 CHECK_FOR_EX_WRITE
1354 return m_data_r;
1355 }
1356
1357 const DataTypes::RealVectorType&
1358 DataTagged::getVectorRO() const
1359 {
1360 return m_data_r;
1361 }
1362
1363 DataTypes::CplxVectorType&
1364 DataTagged::getVectorRWC()
1365 {
1366 CHECK_FOR_EX_WRITE
1367 return m_data_c;
1368 }
1369
1370 const DataTypes::CplxVectorType&
1371 DataTagged::getVectorROC() const
1372 {
1373 return m_data_c;
1374 }
1375
1376 DataTypes::RealVectorType&
1377 DataTagged::getTypedVectorRW(DataTypes::real_t dummy)
1378 {
1379 CHECK_FOR_EX_WRITE
1380 return m_data_r;
1381 }
1382
1383 const DataTypes::RealVectorType&
1384 DataTagged::getTypedVectorRO(DataTypes::real_t dummy) const
1385 {
1386 return m_data_r;
1387 }
1388
1389 DataTypes::CplxVectorType&
1390 DataTagged::getTypedVectorRW(DataTypes::cplx_t dummy)
1391 {
1392 CHECK_FOR_EX_WRITE
1393 return m_data_c;
1394 }
1395
1396 const DataTypes::CplxVectorType&
1397 DataTagged::getTypedVectorRO(DataTypes::cplx_t dummy) const
1398 {
1399 return m_data_c;
1400 }
1401
1402 size_t
1403 DataTagged::getTagCount() const
1404 {
1405 return m_offsetLookup.size();
1406 }
1407
1408
1409 void DataTagged::complicate()
1410 {
1411 if (!isComplex())
1412 {
1413 fillComplexFromReal(m_data_r, m_data_c);
1414 this->m_iscompl=true;
1415 m_data_r.resize(0,0,1);
1416 }
1417 }
1418
1419 } // end of namespace
1420

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26