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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3992 - (show annotations)
Wed Sep 26 02:50:48 2012 UTC (6 years, 8 months ago) by caltinay
File size: 135425 byte(s)
Data.cpp: Removed some namespace ambiguities which fixes doxygen warnings, some
cleanup and reindent.
*: doxygen fixes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2012 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 #include "Data.h"
18
19 #include "AbstractContinuousDomain.h"
20 #include "DataConstant.h"
21 #include "DataEmpty.h"
22 #include "DataExpanded.h"
23 #include "DataLazy.h"
24 #include "DataTagged.h"
25 #include "EscriptParams.h"
26 #include "FunctionSpaceException.h"
27 #include "FunctionSpaceFactory.h"
28 #include "UnaryFuncs.h"
29
30 #ifdef IKNOWWHATIMDOING
31 #include "Dodgy.h"
32 #endif
33
34 extern "C" {
35 #include "esysUtils/blocktimer.h"
36 }
37
38 #include <algorithm>
39 #include <fstream>
40 #include <functional>
41 #include <sstream> // so we can throw messages about ranks
42 #include <vector>
43
44 #include <boost/python/dict.hpp>
45 #include <boost/python/extract.hpp>
46 #include <boost/python/long.hpp>
47 #include "WrappedArray.h"
48
49 namespace bp = boost::python;
50 using namespace escript;
51 using namespace std;
52
53 // ensure the current object is not a DataLazy
54 // The idea was that we could add an optional warning whenever a resolve is forced
55 // #define forceResolve() if (isLazy()) {#resolve();}
56
57 #define AUTOLAZYON escriptParams.getAUTOLAZY()
58 #define MAKELAZYOP(X) do {\
59 if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
60 {\
61 DataLazy* c=new DataLazy(borrowDataPtr(),X);\
62 return Data(c);\
63 }\
64 }while(0)
65
66 #define MAKELAZYOPOFF(X,Y) do {\
67 if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
68 {\
69 DataLazy* c=new DataLazy(borrowDataPtr(),X,Y);\
70 return Data(c);\
71 }\
72 }while(0)
73
74 #define MAKELAZYOP2(X,Y,Z) do {\
75 if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
76 {\
77 DataLazy* c=new DataLazy(borrowDataPtr(),X,Y,Z);\
78 return Data(c);\
79 }\
80 }while(0)
81
82 #define MAKELAZYBINSELF(R,X) do {\
83 if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
84 {\
85 DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
86 /* m_data=c->getPtr();*/ set_m_data(c->getPtr());\
87 return (*this);\
88 }\
89 }while(0)
90
91 // like the above but returns a new data rather than *this
92 #define MAKELAZYBIN(R,X) do {\
93 if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
94 {\
95 DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
96 return Data(c);\
97 }\
98 }while(0)
99
100 #define MAKELAZYBIN2(L,R,X) do {\
101 if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
102 {\
103 DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
104 return Data(c);\
105 }\
106 }while(0)
107
108 #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
109
110 namespace
111 {
112
113 template <class ARR>
114 inline
115 bp::tuple
116 pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
117 {
118 bp::list l;
119 unsigned int dim=shape[0];
120 for (size_t i=0;i<dim;++i)
121 {
122 l.append(v[i+offset]);
123 }
124 return bp::tuple(l);
125 }
126
127 template <class ARR>
128 inline
129 bp::tuple
130 pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
131 {
132 unsigned int shape0=shape[0];
133 unsigned int shape1=shape[1];
134 bp::list lj;
135 for (size_t j=0;j<shape0;++j)
136 {
137 bp::list li;
138 for (size_t i=0;i<shape1;++i)
139 {
140 li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
141 }
142 lj.append(bp::tuple(li));
143 }
144 return bp::tuple(lj);
145 }
146
147 template <class ARR>
148 inline
149 bp::tuple
150 pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
151 {
152 unsigned int shape0=shape[0];
153 unsigned int shape1=shape[1];
154 unsigned int shape2=shape[2];
155
156 bp::list lk;
157 for (size_t k=0;k<shape0;++k)
158 {
159 bp::list lj;
160 for (size_t j=0;j<shape1;++j)
161 {
162 bp::list li;
163 for (size_t i=0;i<shape2;++i)
164 {
165 li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
166 }
167 lj.append(bp::tuple(li));
168 }
169 lk.append(bp::tuple(lj));
170 }
171 return bp::tuple(lk);
172 }
173
174 template <class ARR>
175 inline
176 bp::tuple
177 pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
178 {
179 unsigned int shape0=shape[0];
180 unsigned int shape1=shape[1];
181 unsigned int shape2=shape[2];
182 unsigned int shape3=shape[3];
183
184 bp::list ll;
185 for (size_t l=0;l<shape0;++l)
186 {
187 bp::list lk;
188 for (size_t k=0;k<shape1;++k)
189 {
190 bp::list lj;
191 for (size_t j=0;j<shape2;++j)
192 {
193 bp::list li;
194 for (size_t i=0;i<shape3;++i)
195 {
196 li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
197 }
198 lj.append(bp::tuple(li));
199 }
200 lk.append(bp::tuple(lj));
201 }
202 ll.append(bp::tuple(lk));
203 }
204 return bp::tuple(ll);
205 }
206
207
208 // This should be safer once the DataC RO changes have been brought in
209 template <class ARR>
210 bp::tuple
211 pointToTuple( const DataTypes::ShapeType& shape,ARR v)
212 {
213 int rank=shape.size();
214 if (rank==0)
215 {
216 return bp::make_tuple(v[0]);
217 }
218 else if (rank==1)
219 {
220 return pointToTuple1(shape,v,0);
221 }
222 else if (rank==2)
223 {
224 return pointToTuple2(shape,v,0);
225 }
226 else if (rank==3)
227 {
228 return pointToTuple3(shape,v,0);
229 }
230 else if (rank==4)
231 {
232 return pointToTuple4(shape,v,0);
233 }
234 else
235 throw DataException("Unknown rank in pointToTuple.");
236 }
237
238 } // anonymous namespace
239
240 Data::Data()
241 : m_shared(false), m_lazy(false)
242 {
243 //
244 // Default data is type DataEmpty
245 DataAbstract* temp=new DataEmpty();
246 set_m_data(temp->getPtr());
247 m_protected=false;
248 }
249
250 Data::Data(double value,
251 const bp::tuple& shape,
252 const FunctionSpace& what,
253 bool expanded)
254 : m_shared(false), m_lazy(false)
255 {
256 DataTypes::ShapeType dataPointShape;
257 for (int i = 0; i < shape.attr("__len__")(); ++i) {
258 dataPointShape.push_back(bp::extract<const int>(shape[i]));
259 }
260
261 int len = DataTypes::noValues(dataPointShape);
262 DataVector temp_data(len,value,len);
263 initialise(temp_data, dataPointShape, what, expanded);
264 m_protected=false;
265 }
266
267 Data::Data(double value,
268 const DataTypes::ShapeType& dataPointShape,
269 const FunctionSpace& what,
270 bool expanded)
271 : m_shared(false), m_lazy(false)
272 {
273 initialise(value, dataPointShape, what, expanded);
274 m_protected=false;
275 }
276
277 Data::Data(const Data& inData)
278 : m_shared(false), m_lazy(false)
279 {
280 set_m_data(inData.m_data);
281 m_protected=inData.isProtected();
282 }
283
284
285 Data::Data(const Data& inData,
286 const DataTypes::RegionType& region)
287 : m_shared(false), m_lazy(false)
288 {
289 DataAbstract_ptr dat=inData.m_data;
290 if (inData.isLazy())
291 {
292 dat=inData.m_data->resolve();
293 }
294 else
295 {
296 dat=inData.m_data;
297 }
298 //
299 // Create Data which is a slice of another Data
300 DataAbstract* tmp = dat->getSlice(region);
301 set_m_data(DataAbstract_ptr(tmp));
302 m_protected=false;
303
304 }
305
306 Data::Data(const Data& inData,
307 const FunctionSpace& functionspace)
308 : m_shared(false), m_lazy(false)
309 {
310 if (inData.isEmpty())
311 {
312 throw DataException("Error - will not interpolate for instances of DataEmpty.");
313 }
314 if (inData.getFunctionSpace()==functionspace) {
315 set_m_data(inData.m_data);
316 }
317 else
318 {
319
320 if (inData.isConstant()) { // for a constant function, we just need to use the new function space
321 if (!inData.probeInterpolation(functionspace))
322 { // Even though this is constant, we still need to check whether interpolation is allowed
323 throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
324 }
325 // if the data is not lazy, this will just be a cast to DataReady
326 DataReady_ptr dr=inData.m_data->resolve();
327 DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
328 set_m_data(DataAbstract_ptr(dc));
329 } else {
330 Data tmp(0,inData.getDataPointShape(),functionspace,true);
331 // Note: Must use a reference or pointer to a derived object
332 // in order to get polymorphic behaviour. Shouldn't really
333 // be able to create an instance of AbstractDomain but that was done
334 // as a boost:python work around which may no longer be required.
335 /*const AbstractDomain& inDataDomain=inData.getDomain();*/
336 const_Domain_ptr inDataDomain=inData.getDomain();
337 if (inDataDomain==functionspace.getDomain()) {
338 inDataDomain->interpolateOnDomain(tmp,inData);
339 } else {
340 inDataDomain->interpolateACross(tmp,inData);
341 }
342 set_m_data(tmp.m_data);
343 }
344 }
345 m_protected=false;
346 }
347
348 Data::Data(DataAbstract* underlyingdata)
349 : m_shared(false), m_lazy(false)
350 {
351 set_m_data(underlyingdata->getPtr());
352 m_protected=false;
353 }
354
355 Data::Data(DataAbstract_ptr underlyingdata)
356 : m_shared(false), m_lazy(false)
357 {
358 set_m_data(underlyingdata);
359 m_protected=false;
360 }
361
362 Data::Data(const DataTypes::ValueType& value,
363 const DataTypes::ShapeType& shape,
364 const FunctionSpace& what,
365 bool expanded)
366 : m_shared(false), m_lazy(false)
367 {
368 initialise(value,shape,what,expanded);
369 m_protected=false;
370 }
371
372
373 Data::Data(const bp::object& value,
374 const FunctionSpace& what,
375 bool expanded)
376 : m_shared(false), m_lazy(false)
377 {
378 WrappedArray w(value);
379 initialise(w,what,expanded);
380 m_protected=false;
381 }
382
383
384 Data::Data(const WrappedArray& w, const FunctionSpace& what,
385 bool expanded)
386 :m_shared(false), m_lazy(false)
387 {
388 initialise(w,what,expanded);
389 m_protected=false;
390 }
391
392
393 Data::Data(const bp::object& value,
394 const Data& other)
395 : m_shared(false), m_lazy(false)
396 {
397 WrappedArray w(value);
398
399 // extract the shape of the array
400 const DataTypes::ShapeType& tempShape=w.getShape();
401 if (w.getRank()==0) {
402
403 // get the space for the data vector
404 int len1 = DataTypes::noValues(tempShape);
405 DataVector temp_data(len1, 0.0, len1);
406 temp_data.copyFromArray(w,1);
407
408 int len = DataTypes::noValues(other.getDataPointShape());
409
410 DataVector temp2_data(len, temp_data[0], len);
411 DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
412 set_m_data(DataAbstract_ptr(t));
413
414 } else {
415 //
416 // Create a DataConstant with the same sample shape as other
417 DataConstant* t=new DataConstant(w,other.getFunctionSpace());
418 set_m_data(DataAbstract_ptr(t));
419 }
420 m_protected=false;
421 }
422
423
424 Data::~Data()
425 {
426 set_m_data(DataAbstract_ptr());
427 }
428
429
430 // only call in thread safe contexts.
431 // This method should be atomic
432 void Data::set_m_data(DataAbstract_ptr p)
433 {
434 if (m_data.get()!=0) // release old ownership
435 {
436 m_data->removeOwner(this);
437 }
438 if (p.get()!=0)
439 {
440 m_data=p;
441 m_data->addOwner(this);
442 m_shared=m_data->isShared();
443 m_lazy=m_data->isLazy();
444 }
445 }
446
447 void Data::initialise(const WrappedArray& value,
448 const FunctionSpace& what,
449 bool expanded)
450 {
451 //
452 // Construct a Data object of the appropriate type.
453 // Construct the object first as there seems to be a bug which causes
454 // undefined behaviour if an exception is thrown during construction
455 // within the shared_ptr constructor.
456 if (expanded) {
457 DataAbstract* temp=new DataExpanded(value, what);
458 set_m_data(temp->getPtr());
459 } else {
460 DataAbstract* temp=new DataConstant(value, what);
461 set_m_data(temp->getPtr());
462 }
463 }
464
465
466 void
467 Data::initialise(const DataTypes::ValueType& value,
468 const DataTypes::ShapeType& shape,
469 const FunctionSpace& what,
470 bool expanded)
471 {
472 //
473 // Construct a Data object of the appropriate type.
474 // Construct the object first as there seems to be a bug which causes
475 // undefined behaviour if an exception is thrown during construction
476 // within the shared_ptr constructor.
477 if (expanded) {
478 DataAbstract* temp=new DataExpanded(what, shape, value);
479 set_m_data(temp->getPtr());
480 } else {
481 DataAbstract* temp=new DataConstant(what, shape, value);
482 set_m_data(temp->getPtr());
483 }
484 }
485
486 void
487 Data::initialise(const double value,
488 const DataTypes::ShapeType& shape,
489 const FunctionSpace& what,
490 bool expanded)
491 {
492 //
493 // Construct a Data object of the appropriate type.
494 // Construct the object first as there seems to be a bug which causes
495 // undefined behaviour if an exception is thrown during construction
496 // within the shared_ptr constructor.
497 if (expanded) {
498 DataAbstract* temp=new DataExpanded(what, shape, value);
499 set_m_data(temp->getPtr());
500 } else {
501 DataAbstract* temp=new DataConstant(what, shape, value);
502 set_m_data(temp->getPtr());
503 }
504 }
505
506
507 escriptDataC
508 Data::getDataC()
509 {
510 escriptDataC temp;
511 temp.m_dataPtr=(void*)this;
512 return temp;
513 }
514
515 escriptDataC
516 Data::getDataC() const
517 {
518 escriptDataC temp;
519 temp.m_dataPtr=(void*)this;
520 return temp;
521 }
522
523 const bp::tuple
524 Data::getShapeTuple() const
525 {
526 const DataTypes::ShapeType& shape=getDataPointShape();
527 switch(getDataPointRank()) {
528 case 0:
529 return bp::make_tuple();
530 case 1:
531 return bp::make_tuple(bp::long_(shape[0]));
532 case 2:
533 return bp::make_tuple(bp::long_(shape[0]),bp::long_(shape[1]));
534 case 3:
535 return bp::make_tuple(bp::long_(shape[0]),bp::long_(shape[1]),bp::long_(shape[2]));
536 case 4:
537 return bp::make_tuple(bp::long_(shape[0]),bp::long_(shape[1]),bp::long_(shape[2]),bp::long_(shape[3]));
538 default:
539 throw DataException("Error - illegal Data rank.");
540 }
541 }
542
543
544 // The different name is needed because boost has trouble with overloaded functions.
545 // It can't work out what type the function is based soley on its name.
546 // There are ways to fix this involving creating function pointer variables for each form
547 // but there doesn't seem to be a need given that the methods have the same name from the python point of view
548 Data
549 Data::copySelf()
550 {
551 DataAbstract* temp=m_data->deepCopy();
552 return Data(temp);
553 }
554
555 void
556 Data::copy(const Data& other)
557 {
558 DataAbstract* temp=other.m_data->deepCopy();
559 DataAbstract_ptr p=temp->getPtr();
560 set_m_data(p);
561 }
562
563
564 Data
565 Data::delay()
566 {
567 if (!isLazy())
568 {
569 DataLazy* dl=new DataLazy(m_data);
570 return Data(dl);
571 }
572 return *this;
573 }
574
575 void
576 Data::delaySelf()
577 {
578 if (!isLazy())
579 {
580 set_m_data((new DataLazy(m_data))->getPtr());
581 }
582 }
583
584
585 // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
586 // to zero, all the tags from all the DataTags would be in the result.
587 // However since they all have the same value (0) whether they are there or not should not matter.
588 // So I have decided that for all types this method will create a constant 0.
589 // It can be promoted up as required.
590 // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
591 // but we can deal with that if it arises.
592 //
593 void
594 Data::setToZero()
595 {
596 if (isEmpty())
597 {
598 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
599 }
600 if (isLazy())
601 {
602 DataTypes::ValueType v(getNoValues(),0);
603 DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
604 DataLazy* dl=new DataLazy(dc->getPtr());
605 set_m_data(dl->getPtr());
606 }
607 else
608 {
609 exclusiveWrite();
610 m_data->setToZero();
611 }
612 }
613
614
615 void
616 Data::copyWithMask(const Data& other,
617 const Data& mask)
618 {
619 // 1. Interpolate if required so all Datas use the same FS as this
620 // 2. Tag or Expand so that all Data's are the same type
621 // 3. Iterate over the data vectors copying values where mask is >0
622 if (other.isEmpty() || mask.isEmpty())
623 {
624 throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
625 }
626 Data other2(other);
627 Data mask2(mask);
628 other2.resolve();
629 mask2.resolve();
630 this->resolve();
631 FunctionSpace myFS=getFunctionSpace();
632 FunctionSpace oFS=other2.getFunctionSpace();
633 FunctionSpace mFS=mask2.getFunctionSpace();
634 if (oFS!=myFS)
635 {
636 if (other2.probeInterpolation(myFS))
637 {
638 other2=other2.interpolate(myFS);
639 }
640 else
641 {
642 throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
643 }
644 }
645 if (mFS!=myFS)
646 {
647 if (mask2.probeInterpolation(myFS))
648 {
649 mask2=mask2.interpolate(myFS);
650 }
651 else
652 {
653 throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
654 }
655 }
656 // Ensure that all args have the same type
657 if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
658 {
659 this->expand();
660 other2.expand();
661 mask2.expand();
662 }
663 else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
664 {
665 this->tag();
666 other2.tag();
667 mask2.tag();
668 }
669 else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
670 {
671 }
672 else
673 {
674 throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
675 }
676 unsigned int selfrank=getDataPointRank();
677 unsigned int otherrank=other2.getDataPointRank();
678 unsigned int maskrank=mask2.getDataPointRank();
679 if ((selfrank==0) && (otherrank>0 || maskrank>0))
680 {
681 // to get here we must be copying from a large object into a scalar
682 // I am not allowing this.
683 // If you are calling copyWithMask then you are considering keeping some existing values
684 // and so I'm going to assume that you don't want your data objects getting a new shape.
685 throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
686 }
687 exclusiveWrite();
688 // Now we iterate over the elements
689 DataVector& self=getReady()->getVectorRW();;
690 const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
691 const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
692
693 if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
694 {
695 // Not allowing this combination.
696 // it is not clear what the rank of the target should be.
697 // Should it be filled with the scalar (rank stays the same);
698 // or should the target object be reshaped to be a scalar as well.
699 throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
700 }
701 if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
702 {
703 if (mvec[0]>0) // copy whole object if scalar is >0
704 {
705 copy(other);
706 }
707 return;
708 }
709 if (isTagged()) // so all objects involved will also be tagged
710 {
711 // note the !
712 if (!((getDataPointShape()==mask2.getDataPointShape()) &&
713 ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
714 {
715 throw DataException("copyWithMask, shape mismatch.");
716 }
717
718 // We need to consider the possibility that tags are missing or in the wrong order
719 // My guiding assumption here is: All tagged Datas are assumed to have the default value for
720 // all tags which are not explicitly defined
721
722 const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
723 const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
724 DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
725
726 // first, add any tags missing from other or mask
727 const DataTagged::DataMapType& olookup=optr->getTagLookup();
728 const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
729 const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
730 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
731 for (i=olookup.begin();i!=olookup.end();i++)
732 {
733 tptr->addTag(i->first);
734 }
735 for (i=mlookup.begin();i!=mlookup.end();i++) {
736 tptr->addTag(i->first);
737 }
738 // now we know that *this has all the required tags but they aren't guaranteed to be in
739 // the same order
740
741 // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
742 if ((selfrank==otherrank) && (otherrank==maskrank))
743 {
744 for (i=tlookup.begin();i!=tlookup.end();i++)
745 {
746 // get the target offset
747 DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
748 DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
749 DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
750 for (int j=0;j<getDataPointSize();++j)
751 {
752 if (mvec[j+moff]>0)
753 {
754 self[j+toff]=ovec[j+ooff];
755 }
756 }
757 }
758 // now for the default value
759 for (int j=0;j<getDataPointSize();++j)
760 {
761 if (mvec[j+mptr->getDefaultOffset()]>0)
762 {
763 self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
764 }
765 }
766 }
767 else // other is a scalar
768 {
769 for (i=tlookup.begin();i!=tlookup.end();i++)
770 {
771 // get the target offset
772 DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
773 DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
774 DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
775 for (int j=0;j<getDataPointSize();++j)
776 {
777 if (mvec[j+moff]>0)
778 {
779 self[j+toff]=ovec[ooff];
780 }
781 }
782 }
783 // now for the default value
784 for (int j=0;j<getDataPointSize();++j)
785 {
786 if (mvec[j+mptr->getDefaultOffset()]>0)
787 {
788 self[j+tptr->getDefaultOffset()]=ovec[0];
789 }
790 }
791 }
792
793 return; // ugly
794 }
795 // mixed scalar and non-scalar operation
796 if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
797 {
798 size_t num_points=self.size();
799 // OPENMP 3.0 allows unsigned loop vars.
800 #if defined(_OPENMP) && (_OPENMP < 200805)
801 long i;
802 #else
803 size_t i;
804 #endif
805 size_t psize=getDataPointSize();
806 #pragma omp parallel for private(i) schedule(static)
807 for (i=0;i<num_points;++i)
808 {
809 if (mvec[i]>0)
810 {
811 self[i]=ovec[i/psize]; // since this is expanded there is one scalar
812 } // dest point
813 }
814 return; // ugly!
815 }
816 // tagged data is already taken care of so we only need to worry about shapes
817 // special cases with scalars are already dealt with so all we need to worry about is shape
818 if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
819 {
820 ostringstream oss;
821 oss <<"Error - size mismatch in arguments to copyWithMask.";
822 oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
823 oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
824 oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
825 throw DataException(oss.str());
826 }
827 size_t num_points=self.size();
828
829 // OPENMP 3.0 allows unsigned loop vars.
830 #if defined(_OPENMP) && (_OPENMP < 200805)
831 long i;
832 #else
833 size_t i;
834 #endif
835 #pragma omp parallel for private(i) schedule(static)
836 for (i=0;i<num_points;++i)
837 {
838 if (mvec[i]>0)
839 {
840 self[i]=ovec[i];
841 }
842 }
843 }
844
845 bool
846 Data::isExpanded() const
847 {
848 DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
849 return (temp!=0);
850 }
851
852 bool
853 Data::actsExpanded() const
854 {
855 return m_data->actsExpanded();
856 }
857
858
859 bool
860 Data::isTagged() const
861 {
862 DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
863 return (temp!=0);
864 }
865
866 bool
867 Data::isEmpty() const
868 {
869 DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
870 return (temp!=0);
871 }
872
873 bool
874 Data::isConstant() const
875 {
876 DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
877 return (temp!=0);
878 }
879
880 bool
881 Data::isLazy() const
882 {
883 return m_lazy; // not asking m_data because we need to be able to ask this while m_data is changing
884 }
885
886 // at the moment this is synonymous with !isLazy() but that could change
887 bool
888 Data::isReady() const
889 {
890 return (dynamic_cast<DataReady*>(m_data.get())!=0);
891 }
892
893
894 void
895 Data::setProtection()
896 {
897 m_protected=true;
898 }
899
900 bool
901 Data::isProtected() const
902 {
903 return m_protected;
904 }
905
906
907
908 void
909 Data::expand()
910 {
911 if (isConstant()) {
912 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
913 DataAbstract* temp=new DataExpanded(*tempDataConst);
914 set_m_data(temp->getPtr());
915 } else if (isTagged()) {
916 DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
917 DataAbstract* temp=new DataExpanded(*tempDataTag);
918 set_m_data(temp->getPtr());
919 } else if (isExpanded()) {
920 //
921 // do nothing
922 } else if (isEmpty()) {
923 throw DataException("Error - Expansion of DataEmpty not possible.");
924 } else if (isLazy()) {
925 resolve();
926 expand(); // resolve might not give us expanded data
927 } else {
928 throw DataException("Error - Expansion not implemented for this Data type.");
929 }
930 }
931
932 void
933 Data::tag()
934 {
935 if (isConstant()) {
936 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
937 DataAbstract* temp=new DataTagged(*tempDataConst);
938 set_m_data(temp->getPtr());
939 } else if (isTagged()) {
940 // do nothing
941 } else if (isExpanded()) {
942 throw DataException("Error - Creating tag data from DataExpanded not possible.");
943 } else if (isEmpty()) {
944 throw DataException("Error - Creating tag data from DataEmpty not possible.");
945 } else if (isLazy()) {
946 DataAbstract_ptr res=m_data->resolve();
947 if (m_data->isExpanded())
948 {
949 throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
950 }
951 set_m_data(res);
952 tag();
953 } else {
954 throw DataException("Error - Tagging not implemented for this Data type.");
955 }
956 }
957
958 void
959 Data::resolve()
960 {
961 if (isLazy())
962 {
963 set_m_data(m_data->resolve());
964 }
965 }
966
967 void
968 Data::requireWrite()
969 {
970 resolve();
971 exclusiveWrite();
972 }
973
974 Data
975 Data::oneOver() const
976 {
977 MAKELAZYOP(RECIP);
978 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
979 }
980
981 Data
982 Data::wherePositive() const
983 {
984 MAKELAZYOP(GZ);
985 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
986 }
987
988 Data
989 Data::whereNegative() const
990 {
991 MAKELAZYOP(LZ);
992 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
993 }
994
995 Data
996 Data::whereNonNegative() const
997 {
998 MAKELAZYOP(GEZ);
999 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
1000 }
1001
1002 Data
1003 Data::whereNonPositive() const
1004 {
1005 MAKELAZYOP(LEZ);
1006 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1007 }
1008
1009 Data
1010 Data::whereZero(double tol) const
1011 {
1012 // Data dataAbs=abs();
1013 // return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1014 MAKELAZYOPOFF(EZ,tol);
1015 return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1016 }
1017
1018 Data
1019 Data::whereNonZero(double tol) const
1020 {
1021 // Data dataAbs=abs();
1022 // return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1023 MAKELAZYOPOFF(NEZ,tol);
1024 return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1025 }
1026
1027 Data
1028 Data::interpolate(const FunctionSpace& functionspace) const
1029 {
1030 return Data(*this,functionspace);
1031 }
1032
1033 bool
1034 Data::probeInterpolation(const FunctionSpace& functionspace) const
1035 {
1036 return getFunctionSpace().probeInterpolation(functionspace);
1037 }
1038
1039 Data
1040 Data::gradOn(const FunctionSpace& functionspace) const
1041 {
1042 if (isEmpty())
1043 {
1044 throw DataException("Error - operation not permitted on instances of DataEmpty.");
1045 }
1046 double blocktimer_start = blocktimer_time();
1047 if (functionspace.getDomain()!=getDomain())
1048 throw DataException("Error - gradient cannot be calculated on different domains.");
1049 DataTypes::ShapeType grad_shape=getDataPointShape();
1050 grad_shape.push_back(functionspace.getDim());
1051 Data out(0.0,grad_shape,functionspace,true);
1052 getDomain()->setToGradient(out,*this);
1053 blocktimer_increment("grad()", blocktimer_start);
1054 return out;
1055 }
1056
1057 Data
1058 Data::grad() const
1059 {
1060 if (isEmpty())
1061 {
1062 throw DataException("Error - operation not permitted on instances of DataEmpty.");
1063 }
1064 return gradOn(escript::function(*getDomain()));
1065 }
1066
1067 int
1068 Data::getDataPointSize() const
1069 {
1070 return m_data->getNoValues();
1071 }
1072
1073
1074 DataTypes::ValueType::size_type
1075 Data::getLength() const
1076 {
1077 return m_data->getLength();
1078 }
1079
1080
1081 // There is no parallelism here ... elements need to be added in the correct order.
1082 // If we could presize the list and then fill in the elements it might work
1083 // This would need setting elements to be threadsafe.
1084 // Having mulitple C threads calling into one interpreter is aparently a no-no.
1085 const bp::object
1086 Data::toListOfTuples(bool scalarastuple)
1087 {
1088 if (get_MPISize()>1)
1089 {
1090 throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1091 }
1092 unsigned int rank=getDataPointRank();
1093 unsigned int size=getDataPointSize();
1094
1095 int npoints=getNumDataPoints();
1096 expand(); // This will also resolve if required
1097 const DataTypes::ValueType& vec=getReady()->getVectorRO();
1098 bp::list temp;
1099 temp.append(bp::object());
1100 bp::list res(temp*npoints);// presize the list by the "[None] * npoints" trick
1101 if (rank==0)
1102 {
1103 long count;
1104 if (scalarastuple)
1105 {
1106 for (count=0;count<npoints;++count)
1107 {
1108 res[count]=bp::make_tuple(vec[count]);
1109 }
1110 }
1111 else
1112 {
1113 for (count=0;count<npoints;++count)
1114 {
1115 res[count]=vec[count];
1116 }
1117 }
1118 }
1119 else if (rank==1)
1120 {
1121 size_t count;
1122 size_t offset=0;
1123 for (count=0;count<npoints;++count,offset+=size)
1124 {
1125 res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1126 }
1127 }
1128 else if (rank==2)
1129 {
1130 size_t count;
1131 size_t offset=0;
1132 for (count=0;count<npoints;++count,offset+=size)
1133 {
1134 res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1135 }
1136 }
1137 else if (rank==3)
1138 {
1139 size_t count;
1140 size_t offset=0;
1141 for (count=0;count<npoints;++count,offset+=size)
1142 {
1143 res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1144 }
1145 }
1146 else if (rank==4)
1147 {
1148 size_t count;
1149 size_t offset=0;
1150 for (count=0;count<npoints;++count,offset+=size)
1151 {
1152 res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1153 }
1154 }
1155 else
1156 {
1157 throw DataException("Unknown rank in ::toListOfTuples()");
1158 }
1159 return res;
1160 }
1161
1162 const bp::object
1163 Data::getValueOfDataPointAsTuple(int dataPointNo)
1164 {
1165 forceResolve();
1166 if (getNumDataPointsPerSample()>0) {
1167 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1168 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1169 //
1170 // Check a valid sample number has been supplied
1171 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1172 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1173 }
1174
1175 //
1176 // Check a valid data point number has been supplied
1177 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1178 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1179 }
1180 // TODO: global error handling
1181 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1182 return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1183 }
1184 else
1185 {
1186 // The pre-numpy method would return an empty array of the given shape
1187 // I'm going to throw an exception because if we have zero points per sample we have problems
1188 throw DataException("Error - need at least 1 datapoint per sample.");
1189 }
1190 }
1191
1192
1193 void
1194 Data::setValueOfDataPointToPyObject(int dataPointNo, const bp::object& py_object)
1195 {
1196 // this will throw if the value cannot be represented
1197 setValueOfDataPointToArray(dataPointNo,py_object);
1198 }
1199
1200
1201 void
1202 Data::setTupleForGlobalDataPoint(int id, int proc, bp::object v)
1203 {
1204 #ifdef ESYS_MPI
1205 int error=0;
1206 #endif
1207 if( get_MPIRank()==proc )
1208 {
1209 try
1210 {
1211 bp::extract<double> dex(v);
1212 if (dex.check())
1213 {
1214 setValueOfDataPoint(id, dex());
1215 }
1216 else
1217 {
1218 setValueOfDataPointToArray(id, v);
1219 }
1220 }
1221 catch (...)
1222 {
1223 #ifdef ESYS_MPI
1224 error=1;
1225 int e2;
1226 MPI_Allreduce( &error, &e2, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
1227 #endif
1228 // participate in gather
1229 throw;
1230 }
1231 }
1232 #ifdef ESYS_MPI
1233 int e2;
1234 // If we get here, then either we succeeded or it was on another rank
1235 MPI_Allreduce( &error, &e2, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD );
1236 // participate in gather
1237 if (e2)
1238 {
1239 throw DataException("Error in another rank performing setTupleForGlobalDataPoint");
1240 }
1241 #endif
1242 }
1243
1244
1245 void
1246 Data::setValueOfDataPointToArray(int dataPointNo, const bp::object& obj)
1247 {
1248 if (isProtected()) {
1249 throw DataException("Error - attempt to update protected Data object.");
1250 }
1251
1252 WrappedArray w(obj);
1253 //
1254 // check rank
1255 if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1256 throw DataException("Rank of array does not match Data object rank");
1257
1258 //
1259 // check shape of array
1260 for (unsigned int i=0; i<getDataPointRank(); i++) {
1261 if (w.getShape()[i]!=getDataPointShape()[i])
1262 throw DataException("Shape of array does not match Data object rank");
1263 }
1264
1265 exclusiveWrite();
1266
1267 //
1268 // make sure data is expanded:
1269 //
1270 if (!isExpanded()) {
1271 expand();
1272 }
1273 if (getNumDataPointsPerSample()>0) {
1274 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1275 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1276 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1277 } else {
1278 m_data->copyToDataPoint(-1, 0,w);
1279 }
1280 }
1281
1282 void
1283 Data::setValueOfDataPoint(int dataPointNo, const double value)
1284 {
1285 if (isProtected()) {
1286 throw DataException("Error - attempt to update protected Data object.");
1287 }
1288 //
1289 // make sure data is expanded:
1290 exclusiveWrite();
1291 if (!isExpanded()) {
1292 expand();
1293 }
1294 if (getNumDataPointsPerSample()>0) {
1295 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1296 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1297 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1298 } else {
1299 m_data->copyToDataPoint(-1, 0,value);
1300 }
1301 }
1302
1303 const
1304 bp::object
1305 Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1306 {
1307 // This could be lazier than it is now
1308 forceResolve();
1309
1310 // copy datapoint into a buffer
1311 // broadcast buffer to all nodes
1312 // convert buffer to tuple
1313 // return tuple
1314
1315 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1316 size_t length=DataTypes::noValues(dataPointShape);
1317
1318 // added for the MPI communication
1319 double *tmpData = new double[length];
1320
1321 // updated for the MPI case
1322 if( get_MPIRank()==procNo ){
1323 if (getNumDataPointsPerSample()>0) {
1324 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1325 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1326 //
1327 // Check a valid sample number has been supplied
1328 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1329 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
1330 }
1331
1332 //
1333 // Check a valid data point number has been supplied
1334 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1335 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
1336 }
1337 // TODO: global error handling
1338 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1339
1340 memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1341 }
1342 }
1343 #ifdef ESYS_MPI
1344 // broadcast the data to all other processes
1345 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1346 #endif
1347
1348 bp::tuple t=pointToTuple(dataPointShape,tmpData);
1349 delete [] tmpData;
1350 //
1351 // return the loaded array
1352 return t;
1353 }
1354
1355
1356 bp::object
1357 Data::integrateToTuple_const() const
1358 {
1359 if (isLazy())
1360 {
1361 throw DataException("Error - cannot integrate for constant lazy data.");
1362 }
1363 return integrateWorker();
1364 }
1365
1366 bp::object
1367 Data::integrateToTuple()
1368 {
1369 if (isLazy())
1370 {
1371 expand(); // Can't do a non-resolving version of this without changing the domain object
1372 } // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1373 return integrateWorker();
1374
1375 }
1376
1377 bp::object
1378 Data::integrateWorker() const
1379 {
1380 DataTypes::ShapeType shape = getDataPointShape();
1381 int dataPointSize = getDataPointSize();
1382
1383 //
1384 // calculate the integral values
1385 vector<double> integrals(dataPointSize);
1386 vector<double> integrals_local(dataPointSize);
1387 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1388 if (dom==0)
1389 {
1390 throw DataException("Can not integrate over non-continuous domains.");
1391 }
1392 #ifdef ESYS_MPI
1393 dom->setToIntegrals(integrals_local,*this);
1394 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1395 double *tmp = new double[dataPointSize];
1396 double *tmp_local = new double[dataPointSize];
1397 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1398 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1399 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1400 bp::tuple result=pointToTuple(shape,tmp);
1401 delete[] tmp;
1402 delete[] tmp_local;
1403 #else
1404 dom->setToIntegrals(integrals,*this);
1405 /* double *tmp = new double[dataPointSize];
1406 for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1407 bp::tuple result=pointToTuple(shape,integrals);
1408 // delete tmp;
1409 #endif
1410
1411 return result;
1412 }
1413
1414 Data
1415 Data::sin() const
1416 {
1417 MAKELAZYOP(SIN);
1418 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1419 }
1420
1421 Data
1422 Data::cos() const
1423 {
1424 MAKELAZYOP(COS);
1425 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1426 }
1427
1428 Data
1429 Data::tan() const
1430 {
1431 MAKELAZYOP(TAN);
1432 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1433 }
1434
1435 Data
1436 Data::asin() const
1437 {
1438 MAKELAZYOP(ASIN);
1439 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1440 }
1441
1442 Data
1443 Data::acos() const
1444 {
1445 MAKELAZYOP(ACOS);
1446 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1447 }
1448
1449
1450 Data
1451 Data::atan() const
1452 {
1453 MAKELAZYOP(ATAN);
1454 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1455 }
1456
1457 Data
1458 Data::sinh() const
1459 {
1460 MAKELAZYOP(SINH);
1461 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1462 }
1463
1464 Data
1465 Data::cosh() const
1466 {
1467 MAKELAZYOP(COSH);
1468 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1469 }
1470
1471 Data
1472 Data::tanh() const
1473 {
1474 MAKELAZYOP(TANH);
1475 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1476 }
1477
1478
1479 Data
1480 Data::erf() const
1481 {
1482 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1483 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1484 #else
1485 MAKELAZYOP(ERF);
1486 return C_TensorUnaryOperation(*this, ::erf);
1487 #endif
1488 }
1489
1490 Data
1491 Data::asinh() const
1492 {
1493 MAKELAZYOP(ASINH);
1494 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1495 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1496 #else
1497 return C_TensorUnaryOperation(*this, ::asinh);
1498 #endif
1499 }
1500
1501 Data
1502 Data::acosh() const
1503 {
1504 MAKELAZYOP(ACOSH);
1505 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1506 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1507 #else
1508 return C_TensorUnaryOperation(*this, ::acosh);
1509 #endif
1510 }
1511
1512 Data
1513 Data::atanh() const
1514 {
1515 MAKELAZYOP(ATANH);
1516 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1517 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1518 #else
1519 return C_TensorUnaryOperation(*this, ::atanh);
1520 #endif
1521 }
1522
1523 Data
1524 Data::log10() const
1525 {
1526 MAKELAZYOP(LOG10);
1527 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1528 }
1529
1530 Data
1531 Data::log() const
1532 {
1533 MAKELAZYOP(LOG);
1534 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1535 }
1536
1537 Data
1538 Data::sign() const
1539 {
1540 MAKELAZYOP(SIGN);
1541 return C_TensorUnaryOperation(*this, escript::fsign);
1542 }
1543
1544 Data
1545 Data::abs() const
1546 {
1547 MAKELAZYOP(ABS);
1548 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1549 }
1550
1551 Data
1552 Data::neg() const
1553 {
1554 MAKELAZYOP(NEG);
1555 return C_TensorUnaryOperation(*this, negate<double>());
1556 }
1557
1558 Data
1559 Data::pos() const
1560 {
1561 // not doing lazy check here is deliberate.
1562 // since a deep copy of lazy data should be cheap, I'll just let it happen now
1563 Data result;
1564 // perform a deep copy
1565 result.copy(*this);
1566 return result;
1567 }
1568
1569 Data
1570 Data::exp() const
1571 {
1572 MAKELAZYOP(EXP);
1573 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1574 }
1575
1576 Data
1577 Data::sqrt() const
1578 {
1579 MAKELAZYOP(SQRT);
1580 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1581 }
1582
1583 double
1584 Data::Lsup_const() const
1585 {
1586 if (isLazy())
1587 {
1588 throw DataException("Error - cannot compute Lsup for constant lazy data.");
1589 }
1590 return LsupWorker();
1591 }
1592
1593 double
1594 Data::Lsup()
1595 {
1596 if (isLazy())
1597 {
1598 if (!actsExpanded() || CHECK_DO_CRES)
1599 {
1600 resolve();
1601 }
1602 else
1603 {
1604 #ifdef ESYS_MPI
1605 return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1606 #else
1607 return lazyAlgWorker<AbsMax>(0);
1608 #endif
1609 }
1610 }
1611 return LsupWorker();
1612 }
1613
1614 double
1615 Data::sup_const() const
1616 {
1617 if (isLazy())
1618 {
1619 throw DataException("Error - cannot compute sup for constant lazy data.");
1620 }
1621 return supWorker();
1622 }
1623
1624 double
1625 Data::sup()
1626 {
1627 if (isLazy())
1628 {
1629 if (!actsExpanded() || CHECK_DO_CRES)
1630 {
1631 resolve();
1632 }
1633 else
1634 {
1635 #ifdef ESYS_MPI
1636 return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1637 #else
1638 return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1);
1639 #endif
1640 }
1641 }
1642 return supWorker();
1643 }
1644
1645 double
1646 Data::inf_const() const
1647 {
1648 if (isLazy())
1649 {
1650 throw DataException("Error - cannot compute inf for constant lazy data.");
1651 }
1652 return infWorker();
1653 }
1654
1655 double
1656 Data::inf()
1657 {
1658 if (isLazy())
1659 {
1660 if (!actsExpanded() || CHECK_DO_CRES)
1661 {
1662 resolve();
1663 }
1664 else
1665 {
1666 #ifdef ESYS_MPI
1667 return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
1668 #else
1669 return lazyAlgWorker<FMin>(numeric_limits<double>::max());
1670 #endif
1671 }
1672 }
1673 return infWorker();
1674 }
1675
1676 template <class BinaryOp>
1677 double
1678 #ifdef ESYS_MPI
1679 Data::lazyAlgWorker(double init, MPI_Op mpiop_type)
1680 #else
1681 Data::lazyAlgWorker(double init)
1682 #endif
1683 {
1684 if (!isLazy() || !m_data->actsExpanded())
1685 {
1686 throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1687 }
1688 DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1689 EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1690 double val=init;
1691 int i=0;
1692 const size_t numsamples=getNumSamples();
1693 const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1694 BinaryOp operation;
1695 double localValue=0, globalValue;
1696 #pragma omp parallel private(i)
1697 {
1698 double localtot=init;
1699 #pragma omp for schedule(static) private(i)
1700 for (i=0;i<numsamples;++i)
1701 {
1702 size_t roffset=0;
1703 const DataTypes::ValueType* v=dl->resolveSample(i, roffset);
1704 // Now we have the sample, run operation on all points
1705 for (size_t j=0;j<samplesize;++j)
1706 {
1707 localtot=operation(localtot,(*v)[j+roffset]);
1708 }
1709 if (DataMaths::vectorHasNaN(*v,roffset, samplesize))
1710 {
1711 #pragma omp critical
1712 {
1713 localValue=1.0;
1714 }
1715 }
1716 }
1717 #pragma omp critical
1718 val=operation(val,localtot);
1719 }
1720 #ifdef ESYS_MPI
1721 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1722 #else
1723 globalValue=localValue;
1724 #endif
1725 if (globalValue!=0)
1726 {
1727 return makeNaN();
1728 }
1729 #ifdef ESYS_MPI
1730 MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1731 return globalValue;
1732 #else
1733 return val;
1734 #endif
1735 }
1736
1737 // Do not call this on Lazy Data use the proper entry point
1738 double
1739 Data::LsupWorker() const
1740 {
1741 bool haveNaN=getReady()->hasNaN();
1742 double localValue=0;
1743
1744 #ifdef ESYS_MPI
1745 if (haveNaN)
1746 {
1747 localValue=1.0;
1748 }
1749 double globalValue;
1750 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1751 if (globalValue!=0)
1752 {
1753 return makeNaN();
1754 }
1755 #else
1756 if (haveNaN)
1757 {
1758 return makeNaN();
1759 }
1760 #endif
1761
1762 //
1763 // set the initial absolute maximum value to zero
1764
1765 AbsMax abs_max_func;
1766 localValue = algorithm(abs_max_func,0);
1767
1768 #ifdef ESYS_MPI
1769 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1770 return globalValue;
1771 #else
1772 return localValue;
1773 #endif
1774 }
1775
1776 double
1777 Data::supWorker() const
1778 {
1779 bool haveNaN=getReady()->hasNaN();
1780 double localValue=0;
1781
1782 #ifdef ESYS_MPI
1783 if (haveNaN)
1784 {
1785 localValue=1.0;
1786 }
1787 double globalValue;
1788 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1789 if (globalValue!=0)
1790 {
1791 return makeNaN();
1792 }
1793 #else
1794 if (haveNaN)
1795 {
1796 return makeNaN();
1797 }
1798 #endif
1799
1800 //
1801 // set the initial maximum value to min possible double
1802 FMax fmax_func;
1803 localValue = algorithm(fmax_func,numeric_limits<double>::infinity()*-1);
1804 #ifdef ESYS_MPI
1805 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1806 return globalValue;
1807 #else
1808 return localValue;
1809 #endif
1810 }
1811
1812 double
1813 Data::infWorker() const
1814 {
1815 bool haveNaN=getReady()->hasNaN();
1816 double localValue=0;
1817
1818 #ifdef ESYS_MPI
1819 if (haveNaN)
1820 {
1821 localValue=1.0;
1822 }
1823 double globalValue;
1824 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1825 if (globalValue!=0)
1826 {
1827 return makeNaN();
1828 }
1829 #else
1830 if (haveNaN)
1831 {
1832 return makeNaN();
1833 }
1834 #endif
1835 //
1836 // set the initial minimum value to max possible double
1837 FMin fmin_func;
1838 localValue = algorithm(fmin_func,numeric_limits<double>::infinity());
1839 #ifdef ESYS_MPI
1840 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1841 return globalValue;
1842 #else
1843 return localValue;
1844 #endif
1845 }
1846
1847 /* global reduction */
1848
1849
1850 inline Data
1851 Data::minval_nonlazy() const
1852 {
1853 //
1854 // set the initial minimum value to max possible double
1855 FMin fmin_func;
1856 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1857 }
1858
1859
1860 inline Data
1861 Data::maxval_nonlazy() const
1862 {
1863 //
1864 // set the initial maximum value to min possible double
1865 FMax fmax_func;
1866 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1867 }
1868
1869
1870 Data
1871 Data::maxval() const
1872 {
1873 MAKELAZYOP(MAXVAL);
1874 return maxval_nonlazy();
1875 }
1876
1877
1878 Data
1879 Data::minval() const
1880 {
1881 MAKELAZYOP(MINVAL);
1882 return minval_nonlazy();
1883 }
1884
1885
1886 Data
1887 Data::swapaxes(const int axis0, const int axis1) const
1888 {
1889 int axis0_tmp,axis1_tmp;
1890 DataTypes::ShapeType s=getDataPointShape();
1891 DataTypes::ShapeType ev_shape;
1892 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1893 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1894 int rank=getDataPointRank();
1895 if (rank<2) {
1896 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1897 }
1898 if (axis0<0 || axis0>rank-1) {
1899 stringstream e;
1900 e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
1901 throw DataException(e.str());
1902 }
1903 if (axis1<0 || axis1>rank-1) {
1904 stringstream e;
1905 e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
1906 throw DataException(e.str());
1907 }
1908 if (axis0 == axis1) {
1909 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1910 }
1911 MAKELAZYOP2(SWAP,axis0,axis1);
1912 if (axis0 > axis1)
1913 {
1914 axis0_tmp=axis1;
1915 axis1_tmp=axis0;
1916 }
1917 else
1918 {
1919 axis0_tmp=axis0;
1920 axis1_tmp=axis1;
1921 }
1922 for (int i=0; i<rank; i++)
1923 {
1924 if (i == axis0_tmp)
1925 {
1926 ev_shape.push_back(s[axis1_tmp]);
1927 }
1928 else if (i == axis1_tmp)
1929 {
1930 ev_shape.push_back(s[axis0_tmp]);
1931 }
1932 else
1933 {
1934 ev_shape.push_back(s[i]);
1935 }
1936 }
1937 Data ev(0.,ev_shape,getFunctionSpace());
1938 ev.typeMatchRight(*this);
1939 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1940 return ev;
1941 }
1942
1943 Data
1944 Data::symmetric() const
1945 {
1946 // check input
1947 DataTypes::ShapeType s=getDataPointShape();
1948 if (getDataPointRank()==2) {
1949 if(s[0] != s[1])
1950 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1951 }
1952 else if (getDataPointRank()==4) {
1953 if(!(s[0] == s[2] && s[1] == s[3]))
1954 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1955 }
1956 else {
1957 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1958 }
1959 MAKELAZYOP(SYM);
1960 Data ev(0.,getDataPointShape(),getFunctionSpace());
1961 ev.typeMatchRight(*this);
1962 m_data->symmetric(ev.m_data.get());
1963 return ev;
1964 }
1965
1966 Data
1967 Data::nonsymmetric() const
1968 {
1969 MAKELAZYOP(NSYM);
1970 // check input
1971 DataTypes::ShapeType s=getDataPointShape();
1972 if (getDataPointRank()==2) {
1973 if(s[0] != s[1])
1974 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1975 DataTypes::ShapeType ev_shape;
1976 ev_shape.push_back(s[0]);
1977 ev_shape.push_back(s[1]);
1978 Data ev(0.,ev_shape,getFunctionSpace());
1979 ev.typeMatchRight(*this);
1980 m_data->nonsymmetric(ev.m_data.get());
1981 return ev;
1982 }
1983 else if (getDataPointRank()==4) {
1984 if(!(s[0] == s[2] && s[1] == s[3]))
1985 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1986 DataTypes::ShapeType ev_shape;
1987 ev_shape.push_back(s[0]);
1988 ev_shape.push_back(s[1]);
1989 ev_shape.push_back(s[2]);
1990 ev_shape.push_back(s[3]);
1991 Data ev(0.,ev_shape,getFunctionSpace());
1992 ev.typeMatchRight(*this);
1993 m_data->nonsymmetric(ev.m_data.get());
1994 return ev;
1995 }
1996 else {
1997 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1998 }
1999 }
2000
2001 Data
2002 Data::trace(int axis_offset) const
2003 {
2004 MAKELAZYOPOFF(TRACE,axis_offset);
2005 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
2006 {
2007 throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
2008 }
2009 DataTypes::ShapeType s=getDataPointShape();
2010 if (getDataPointRank()==2) {
2011 DataTypes::ShapeType ev_shape;
2012 Data ev(0.,ev_shape,getFunctionSpace());
2013 ev.typeMatchRight(*this);
2014 m_data->trace(ev.m_data.get(), axis_offset);
2015 return ev;
2016 }
2017 if (getDataPointRank()==3) {
2018 DataTypes::ShapeType ev_shape;
2019 if (axis_offset==0) {
2020 int s2=s[2];
2021 ev_shape.push_back(s2);
2022 }
2023 else if (axis_offset==1) {
2024 int s0=s[0];
2025 ev_shape.push_back(s0);
2026 }
2027 Data ev(0.,ev_shape,getFunctionSpace());
2028 ev.typeMatchRight(*this);
2029 m_data->trace(ev.m_data.get(), axis_offset);
2030 return ev;
2031 }
2032 if (getDataPointRank()==4) {
2033 DataTypes::ShapeType ev_shape;
2034 if (axis_offset==0) {
2035 ev_shape.push_back(s[2]);
2036 ev_shape.push_back(s[3]);
2037 }
2038 else if (axis_offset==1) {
2039 ev_shape.push_back(s[0]);
2040 ev_shape.push_back(s[3]);
2041 }
2042 else if (axis_offset==2) {
2043 ev_shape.push_back(s[0]);
2044 ev_shape.push_back(s[1]);
2045 }
2046 Data ev(0.,ev_shape,getFunctionSpace());
2047 ev.typeMatchRight(*this);
2048 m_data->trace(ev.m_data.get(), axis_offset);
2049 return ev;
2050 }
2051 else {
2052 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
2053 }
2054 }
2055
2056 Data
2057 Data::transpose(int axis_offset) const
2058 {
2059 MAKELAZYOPOFF(TRANS,axis_offset);
2060 DataTypes::ShapeType s=getDataPointShape();
2061 DataTypes::ShapeType ev_shape;
2062 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
2063 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
2064 int rank=getDataPointRank();
2065 if (axis_offset<0 || axis_offset>rank) {
2066 stringstream e;
2067 e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
2068 throw DataException(e.str());
2069 }
2070 for (int i=0; i<rank; i++) {
2071 int index = (axis_offset+i)%rank;
2072 ev_shape.push_back(s[index]); // Append to new shape
2073 }
2074 Data ev(0.,ev_shape,getFunctionSpace());
2075 ev.typeMatchRight(*this);
2076 m_data->transpose(ev.m_data.get(), axis_offset);
2077 return ev;
2078 }
2079
2080 Data
2081 Data::eigenvalues() const
2082 {
2083 if (isLazy())
2084 {
2085 Data temp(*this); // to get around the fact that you can't resolve a const Data
2086 temp.resolve();
2087 return temp.eigenvalues();
2088 }
2089 // check input
2090 DataTypes::ShapeType s=getDataPointShape();
2091 if (getDataPointRank()!=2)
2092 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
2093 if(s[0] != s[1])
2094 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
2095 // create return
2096 DataTypes::ShapeType ev_shape(1,s[0]);
2097 Data ev(0.,ev_shape,getFunctionSpace());
2098 ev.typeMatchRight(*this);
2099 m_data->eigenvalues(ev.m_data.get());
2100 return ev;
2101 }
2102
2103 const bp::tuple
2104 Data::eigenvalues_and_eigenvectors(const double tol) const
2105 {
2106 if (isLazy())
2107 {
2108 Data temp(*this); // to get around the fact that you can't resolve a const Data
2109 temp.resolve();
2110 return temp.eigenvalues_and_eigenvectors(tol);
2111 }
2112 DataTypes::ShapeType s=getDataPointShape();
2113 if (getDataPointRank()!=2)
2114 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
2115 if(s[0] != s[1])
2116 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
2117 // create return
2118 DataTypes::ShapeType ev_shape(1,s[0]);
2119 Data ev(0.,ev_shape,getFunctionSpace());
2120 ev.typeMatchRight(*this);
2121 DataTypes::ShapeType V_shape(2,s[0]);
2122 Data V(0.,V_shape,getFunctionSpace());
2123 V.typeMatchRight(*this);
2124 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
2125 return bp::make_tuple(bp::object(ev),bp::object(V));
2126 }
2127
2128 const bp::tuple
2129 Data::minGlobalDataPoint() const
2130 {
2131 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
2132 // abort (for unknown reasons) if there are openmp directives with it in the
2133 // surrounding function
2134
2135 int DataPointNo;
2136 int ProcNo;
2137 calc_minGlobalDataPoint(ProcNo,DataPointNo);
2138 return bp::make_tuple(ProcNo,DataPointNo);
2139 }
2140
2141 void
2142 Data::calc_minGlobalDataPoint(int& ProcNo,
2143 int& DataPointNo) const
2144 {
2145 if (isLazy())
2146 {
2147 Data temp(*this); // to get around the fact that you can't resolve a const Data
2148 temp.resolve();
2149 return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2150 }
2151 int i,j;
2152 int lowi=0,lowj=0;
2153 double min=numeric_limits<double>::max();
2154
2155 Data temp=minval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2156
2157 int numSamples=temp.getNumSamples();
2158 int numDPPSample=temp.getNumDataPointsPerSample();
2159
2160 double local_val, local_min;
2161 #ifdef ESYS_MPI
2162 double next[2];
2163 #endif
2164 int local_lowi=0,local_lowj=0;
2165
2166 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(local_val,local_min)
2167 {
2168 local_min=min;
2169 #pragma omp for private(i,j) schedule(static)
2170 for (i=0; i<numSamples; i++) {
2171 for (j=0; j<numDPPSample; j++) {
2172 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2173 if (local_val<local_min) {
2174 local_min=local_val;
2175 local_lowi=i;
2176 local_lowj=j;
2177 }
2178 }
2179 }
2180 #pragma omp critical
2181 if (local_min<min) { // If we found a smaller value than our sentinel
2182 min=local_min;
2183 lowi=local_lowi;
2184 lowj=local_lowj;
2185 }
2186 } // parallel section
2187
2188 #ifdef ESYS_MPI
2189 // determine the processor on which the minimum occurs
2190 next[0] = min;
2191 next[1] = numSamples;
2192 int lowProc = 0;
2193 double *globalMins = new double[get_MPISize()*2+1];
2194 /*int error =*/ MPI_Gather (next, 2, MPI_DOUBLE, globalMins, 2, MPI_DOUBLE, 0, get_MPIComm() );
2195
2196 if ( get_MPIRank()==0 ) {
2197 for (lowProc=0; lowProc<get_MPISize(); lowProc++)
2198 if (globalMins[lowProc*2+1] > 0) break;
2199 min = globalMins[lowProc*2];
2200 for( i=lowProc+1; i<get_MPISize(); i++ )
2201 if( globalMins[i*2+1]>0 && min>globalMins[i*2] ) {
2202 lowProc = i;
2203 min = globalMins[i*2];
2204 }
2205 }
2206 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2207 DataPointNo = lowj + lowi * numDPPSample;
2208 MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2209 delete [] globalMins;
2210 ProcNo = lowProc;
2211 #else
2212 ProcNo = 0;
2213 DataPointNo = lowj + lowi * numDPPSample;
2214 #endif
2215 }
2216
2217
2218 const bp::tuple
2219 Data::maxGlobalDataPoint() const
2220 {
2221 int DataPointNo;
2222 int ProcNo;
2223 calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2224 return bp::make_tuple(ProcNo,DataPointNo);
2225 }
2226
2227 void
2228 Data::calc_maxGlobalDataPoint(int& ProcNo,
2229 int& DataPointNo) const
2230 {
2231 if (isLazy())
2232 {
2233 Data temp(*this); // to get around the fact that you can't resolve a const Data
2234 temp.resolve();
2235 return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2236 }
2237 int i,j;
2238 int highi=0,highj=0;
2239 //-------------
2240 double max= -numeric_limits<double>::max();
2241
2242 Data temp=maxval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2243
2244 int numSamples=temp.getNumSamples();
2245 int numDPPSample=temp.getNumDataPointsPerSample();
2246
2247 double local_val, local_max;
2248 #ifdef ESYS_MPI
2249 double next[2];
2250 #endif
2251 int local_highi=0,local_highj=0;
2252
2253 #pragma omp parallel firstprivate(local_highi,local_highj) private(local_val,local_max)
2254 {
2255 local_max=max;
2256 #pragma omp for private(i,j) schedule(static)
2257 for (i=0; i<numSamples; i++) {
2258 for (j=0; j<numDPPSample; j++) {
2259 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2260 if (local_val>local_max) {
2261 local_max=local_val;
2262 local_highi=i;
2263 local_highj=j;
2264 }
2265 }
2266 }
2267 #pragma omp critical
2268 if (local_max>max) { // If we found a larger value than our sentinel
2269 max=local_max;
2270 highi=local_highi;
2271 highj=local_highj;
2272 }
2273 }
2274 #ifdef ESYS_MPI
2275 // determine the processor on which the maximum occurs
2276 next[0] = max;
2277 next[1] = numSamples;
2278 int highProc = 0;
2279 double *globalMaxs = new double[get_MPISize()*2+1];
2280 /*int error =*/ MPI_Gather ( next, 2, MPI_DOUBLE, globalMaxs, 2, MPI_DOUBLE, 0, get_MPIComm() );
2281 if( get_MPIRank()==0 ){
2282 for (highProc=0; highProc<get_MPISize(); highProc++)
2283 if (globalMaxs[highProc*2+1] > 0) break;
2284 max = globalMaxs[highProc*2];
2285 for( i=highProc+1; i<get_MPISize(); i++ )
2286 {
2287 if( globalMaxs[i*2+1]>0 && max<globalMaxs[i*2] )
2288 {
2289 highProc = i;
2290 max = globalMaxs[i*2];
2291 }
2292 }
2293 }
2294 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2295 DataPointNo = highj + highi * numDPPSample;
2296 MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2297
2298 delete [] globalMaxs;
2299 ProcNo = highProc;
2300 #else
2301 ProcNo = 0;
2302 DataPointNo = highj + highi * numDPPSample;
2303 #endif
2304 }
2305
2306 void
2307 Data::saveDX(std::string fileName) const
2308 {
2309 if (isEmpty())
2310 {
2311 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2312 }
2313 if (isLazy())
2314 {
2315 Data temp(*this); // to get around the fact that you can't resolve a const Data
2316 temp.resolve();
2317 temp.saveDX(fileName);
2318 return;
2319 }
2320 bp::dict args;
2321 args["data"]=bp::object(this);
2322 getDomain()->saveDX(fileName,args);
2323 }
2324
2325 void
2326 Data::saveVTK(std::string fileName) const
2327 {
2328 if (isEmpty())
2329 {
2330 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2331 }
2332 if (isLazy())
2333 {
2334 Data temp(*this); // to get around the fact that you can't resolve a const Data
2335 temp.resolve();
2336 temp.saveVTK(fileName);
2337 return;
2338 }
2339 bp::dict args;
2340 args["data"]=bp::object(this);
2341 getDomain()->saveVTK(fileName,args,"","");
2342 }
2343
2344
2345
2346 Data&
2347 Data::operator+=(const Data& right)
2348 {
2349 if (isProtected()) {
2350 throw DataException("Error - attempt to update protected Data object.");
2351 }
2352 MAKELAZYBINSELF(right,ADD); // for lazy + is equivalent to +=
2353 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2354 binaryOp(right,plus<double>());
2355 return (*this);
2356 }
2357
2358 Data&
2359 Data::operator+=(const bp::object& right)
2360 {
2361 if (isProtected()) {
2362 throw DataException("Error - attempt to update protected Data object.");
2363 }
2364 Data tmp(right,getFunctionSpace(),false);
2365 (*this)+=tmp;
2366 return *this;
2367 }
2368
2369 // Hmmm, operator= makes a deep copy but the copy constructor does not?
2370 Data&
2371 Data::operator=(const Data& other)
2372 {
2373 m_protected=false; // since any changes should be caught by exclusiveWrite()
2374 set_m_data(other.m_data);
2375 return (*this);
2376 }
2377
2378 Data&
2379 Data::operator-=(const Data& right)
2380 {
2381 if (isProtected()) {
2382 throw DataException("Error - attempt to update protected Data object.");
2383 }
2384 MAKELAZYBINSELF(right,SUB);
2385 exclusiveWrite();
2386 binaryOp(right,minus<double>());
2387 return (*this);
2388 }
2389
2390 Data&
2391 Data::operator-=(const bp::object& right)
2392 {
2393 if (isProtected()) {
2394 throw DataException("Error - attempt to update protected Data object.");
2395 }
2396 Data tmp(right,getFunctionSpace(),false);
2397 (*this)-=tmp;
2398 return (*this);
2399 }
2400
2401 Data&
2402 Data::operator*=(const Data& right)
2403 {
2404 if (isProtected()) {
2405 throw DataException("Error - attempt to update protected Data object.");
2406 }
2407 MAKELAZYBINSELF(right,MUL);
2408 exclusiveWrite();
2409 binaryOp(right,multiplies<double>());
2410 return (*this);
2411 }
2412
2413 Data&
2414 Data::operator*=(const bp::object& right)
2415 {
2416 if (isProtected()) {
2417 throw DataException("Error - attempt to update protected Data object.");
2418 }
2419 Data tmp(right,getFunctionSpace(),false);
2420 (*this)*=tmp;
2421 return (*this);
2422 }
2423
2424 Data&
2425 Data::operator/=(const Data& right)
2426 {
2427 if (isProtected()) {
2428 throw DataException("Error - attempt to update protected Data object.");
2429 }
2430 MAKELAZYBINSELF(right,DIV);
2431 exclusiveWrite();
2432 binaryOp(right,divides<double>());
2433 return (*this);
2434 }
2435
2436 Data&
2437 Data::operator/=(const bp::object& right)
2438 {
2439 if (isProtected()) {
2440 throw DataException("Error - attempt to update protected Data object.");
2441 }
2442 Data tmp(right,getFunctionSpace(),false);
2443 (*this)/=tmp;
2444 return (*this);
2445 }
2446
2447 /* Be careful trying to make this operation lazy.
2448 At time of writing, resolve() and resolveSample() do not throw.
2449 Changing this would mean that any resolve call would need to use MPI (to check for global errors)
2450 */
2451 Data
2452 Data::matrixInverse() const
2453 {
2454 if (isLazy())
2455 {
2456 Data d(*this);
2457 d.resolve();
2458 return d.matrixInverse();
2459 }
2460
2461 Data out(0.,getDataPointShape(),getFunctionSpace());
2462 out.typeMatchRight(*this);
2463 int errcode=m_data->matrixInverse(out.getReadyPtr().get());
2464 #ifdef ESYS_MPI
2465 int globalval=0;
2466 MPI_Allreduce( &errcode, &globalval, 1, MPI_INT, MPI_MAX, get_MPIComm() );
2467 errcode=globalval;
2468 #endif
2469 if (errcode)
2470 {
2471 DataMaths::matrixInverseError(errcode); // throws exceptions
2472 }
2473 return out;
2474 }
2475
2476
2477 Data
2478 Data::rpowO(const bp::object& left) const
2479 {
2480 Data left_d(left,*this);
2481 return left_d.powD(*this);
2482 }
2483
2484 Data
2485 Data::powO(const bp::object& right) const
2486 {
2487 Data tmp(right,getFunctionSpace(),false);
2488 return powD(tmp);
2489 }
2490
2491 Data
2492 Data::powD(const Data& right) const
2493 {
2494 MAKELAZYBIN(right,POW);
2495 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2496 }
2497
2498 //
2499 // NOTE: It is essential to specify the namespace this operator belongs to
2500 Data
2501 escript::operator+(const Data& left, const Data& right)
2502 {
2503 MAKELAZYBIN2(left,right,ADD);
2504 return C_TensorBinaryOperation(left, right, plus<double>());
2505 }
2506
2507 //
2508 // NOTE: It is essential to specify the namespace this operator belongs to
2509 Data
2510 escript::operator-(const Data& left, const Data& right)
2511 {
2512 MAKELAZYBIN2(left,right,SUB);
2513 return C_TensorBinaryOperation(left, right, minus<double>());
2514 }
2515
2516 //
2517 // NOTE: It is essential to specify the namespace this operator belongs to
2518 Data
2519 escript::operator*(const Data& left, const Data& right)
2520 {
2521 MAKELAZYBIN2(left,right,MUL);
2522 return C_TensorBinaryOperation(left, right, multiplies<double>());
2523 }
2524
2525 //
2526 // NOTE: It is essential to specify the namespace this operator belongs to
2527 Data
2528 escript::operator/(const Data& left, const Data& right)
2529 {
2530 MAKELAZYBIN2(left,right,DIV);
2531 return C_TensorBinaryOperation(left, right, divides<double>());
2532 }
2533
2534 //
2535 // NOTE: It is essential to specify the namespace this operator belongs to
2536 Data
2537 escript::operator+(const Data& left, const bp::object& right)
2538 {
2539 Data tmp(right,left.getFunctionSpace(),false);
2540 MAKELAZYBIN2(left,tmp,ADD);
2541 return left+tmp;
2542 }
2543
2544 //
2545 // NOTE: It is essential to specify the namespace this operator belongs to
2546 Data
2547 escript::operator-(const Data& left, const bp::object& right)
2548 {
2549 Data tmp(right,left.getFunctionSpace(),false);
2550 MAKELAZYBIN2(left,tmp,SUB);
2551 return left-tmp;
2552 }
2553
2554 //
2555 // NOTE: It is essential to specify the namespace this operator belongs to
2556 Data
2557 escript::operator*(const Data& left, const bp::object& right)
2558 {
2559 Data tmp(right,left.getFunctionSpace(),false);
2560 MAKELAZYBIN2(left,tmp,MUL);
2561 return left*tmp;
2562 }
2563
2564 //
2565 // NOTE: It is essential to specify the namespace this operator belongs to
2566 Data
2567 escript::operator/(const Data& left, const bp::object& right)
2568 {
2569 Data tmp(right,left.getFunctionSpace(),false);
2570 MAKELAZYBIN2(left,tmp,DIV);
2571 return left/tmp;
2572 }
2573
2574 //
2575 // NOTE: It is essential to specify the namespace this operator belongs to
2576 Data
2577 escript::operator+(const bp::object& left, const Data& right)
2578 {
2579 Data tmp(left,right.getFunctionSpace(),false);
2580 MAKELAZYBIN2(tmp,right,ADD);
2581 return tmp+right;
2582 }
2583
2584 //
2585 // NOTE: It is essential to specify the namespace this operator belongs to
2586 Data
2587 escript::operator-(const bp::object& left, const Data& right)
2588 {
2589 Data tmp(left,right.getFunctionSpace(),false);
2590 MAKELAZYBIN2(tmp,right,SUB);
2591 return tmp-right;
2592 }
2593
2594 //
2595 // NOTE: It is essential to specify the namespace this operator belongs to
2596 Data
2597 escript::operator*(const bp::object& left, const Data& right)
2598 {
2599 Data tmp(left,right.getFunctionSpace(),false);
2600 MAKELAZYBIN2(tmp,right,MUL);
2601 return tmp*right;
2602 }
2603
2604 //
2605 // NOTE: It is essential to specify the namespace this operator belongs to
2606 Data
2607 escript::operator/(const bp::object& left, const Data& right)
2608 {
2609 Data tmp(left,right.getFunctionSpace(),false);
2610 MAKELAZYBIN2(tmp,right,DIV);
2611 return tmp/right;
2612 }
2613
2614
2615 /* TODO */
2616 /* global reduction */
2617 Data
2618 Data::getItem(const bp::object& key) const
2619 {
2620 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2621
2622 if (slice_region.size()!=getDataPointRank()) {
2623 throw DataException("Error - slice size does not match Data rank.");
2624 }
2625
2626 return getSlice(slice_region);
2627 }
2628
2629 /* TODO */
2630 /* global reduction */
2631 Data
2632 Data::getSlice(const DataTypes::RegionType& region) const
2633 {
2634 return Data(*this,region);
2635 }
2636
2637 /* TODO */
2638 /* global reduction */
2639 void
2640 Data::setItemO(const bp::object& key,
2641 const bp::object& value)
2642 {
2643 Data tempData(value,getFunctionSpace());
2644 setItemD(key,tempData);
2645 }
2646
2647 void
2648 Data::setItemD(const bp::object& key,
2649 const Data& value)
2650 {
2651 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2652 if (slice_region.size()!=getDataPointRank()) {
2653 throw DataException("Error - slice size does not match Data rank.");
2654 }
2655 exclusiveWrite();
2656 if (getFunctionSpace()!=value.getFunctionSpace()) {
2657 setSlice(Data(value,getFunctionSpace()),slice_region);
2658 } else {
2659 setSlice(value,slice_region);
2660 }
2661 }
2662
2663 void
2664 Data::setSlice(const Data& value,
2665 const DataTypes::RegionType& region)
2666 {
2667 if (isProtected()) {
2668 throw DataException("Error - attempt to update protected Data object.");
2669 }
2670 forceResolve();
2671 exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2672 Data tempValue(value);
2673 typeMatchLeft(tempValue);
2674 typeMatchRight(tempValue);
2675 getReady()->setSlice(tempValue.m_data.get(),region);
2676 }
2677
2678 void
2679 Data::typeMatchLeft(Data& right) const
2680 {
2681 if (right.isLazy() && !isLazy())
2682 {
2683 right.resolve();
2684 }
2685 if (isExpanded()) {
2686 right.expand();
2687 } else if (isTagged()) {
2688 if (right.isConstant()) {
2689 right.tag();
2690 }
2691 }
2692 }
2693
2694 void
2695 Data::typeMatchRight(const Data& right)
2696 {
2697 if (isLazy() && !right.isLazy())
2698 {
2699 resolve();
2700 }
2701 if (isTagged()) {
2702 if (right.isExpanded()) {
2703 expand();
2704 }
2705 } else if (isConstant()) {
2706 if (right.isExpanded()) {
2707 expand();
2708 } else if (right.isTagged()) {
2709 tag();
2710 }
2711 }
2712 }
2713
2714 // The normal TaggedValue adds the tag if it is not already present
2715 // This form does not. It throws instead.
2716 // This is because the names are maintained by the domain and cannot be added
2717 // without knowing the tag number to map it to.
2718 void
2719 Data::setTaggedValueByName(std::string name,
2720 const bp::object& value)
2721 {
2722 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2723 forceResolve();
2724 exclusiveWrite();
2725 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2726 setTaggedValue(tagKey,value);
2727 }
2728 else
2729 {
2730 throw DataException("Error - unknown tag in setTaggedValueByName.");
2731 }
2732 }
2733
2734 void
2735 Data::setTaggedValue(int tagKey,
2736 const bp::object& value)
2737 {
2738 if (isProtected()) {
2739 throw DataException("Error - attempt to update protected Data object.");
2740 }
2741 //
2742 // Ensure underlying data object is of type DataTagged
2743 forceResolve();
2744 exclusiveWrite();
2745 if (isConstant()) tag();
2746 WrappedArray w(value);
2747
2748 DataVector temp_data2;
2749 temp_data2.copyFromArray(w,1);
2750
2751 m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2752 }
2753
2754
2755 void
2756 Data::setTaggedValueFromCPP(int tagKey,
2757 const DataTypes::ShapeType& pointshape,
2758 const DataTypes::ValueType& value,
2759 int dataOffset)
2760 {
2761 if (isProtected()) {
2762 throw DataException("Error - attempt to update protected Data object.");
2763 }
2764 //
2765 // Ensure underlying data object is of type DataTagged
2766 forceResolve();
2767 if (isConstant()) tag();
2768 exclusiveWrite();
2769 //
2770 // Call DataAbstract::setTaggedValue
2771 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2772 }
2773
2774 int
2775 Data::getTagNumber(int dpno)
2776 {
2777 if (isEmpty())
2778 {
2779 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2780 }
2781 return getFunctionSpace().getTagFromDataPointNo(dpno);
2782 }
2783
2784
2785 ostream& escript::operator<<(ostream& o, const Data& data)
2786 {
2787 o << data.toString();
2788 return o;
2789 }
2790
2791 Data
2792 escript::C_GeneralTensorProduct(Data& arg_0,
2793 Data& arg_1,
2794 int axis_offset,
2795 int transpose)
2796 {
2797 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2798 // SM is the product of the last axis_offset entries in arg_0.getShape().
2799
2800 // deal with any lazy data
2801 // if (arg_0.isLazy()) {arg_0.resolve();}
2802 // if (arg_1.isLazy()) {arg_1.resolve();}
2803 if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2804 {
2805 DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2806 return Data(c);
2807 }
2808
2809 // Interpolate if necessary and find an appropriate function space
2810 Data arg_0_Z, arg_1_Z;
2811 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2812 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2813 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2814 arg_1_Z = Data(arg_1);
2815 }
2816 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2817 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2818 arg_0_Z =Data(arg_0);
2819 }
2820 else {
2821 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2822 }
2823 } else {
2824 arg_0_Z = Data(arg_0);
2825 arg_1_Z = Data(arg_1);
2826 }
2827 // Get rank and shape of inputs
2828 int rank0 = arg_0_Z.getDataPointRank();
2829 int rank1 = arg_1_Z.getDataPointRank();
2830 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2831 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2832
2833 // Prepare for the loops of the product and verify compatibility of shapes
2834 int start0=0, start1=0;
2835 if (transpose == 0) {}
2836 else if (transpose == 1) { start0 = axis_offset; }
2837 else if (transpose == 2) { start1 = rank1-axis_offset; }
2838 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2839
2840
2841 // Adjust the shapes for transpose
2842 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2843 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2844 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2845 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2846
2847 #if 0
2848 // For debugging: show shape after transpose
2849 char tmp[100];
2850 std::string shapeStr;
2851 shapeStr = "(";
2852 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2853 shapeStr += ")";
2854 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2855 shapeStr = "(";
2856 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2857 shapeStr += ")";
2858 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2859 #endif
2860
2861 // Prepare for the loops of the product
2862 int SL=1, SM=1, SR=1;
2863 for (int i=0; i<rank0-axis_offset; i++) {
2864 SL *= tmpShape0[i];
2865 }
2866 for (int i=rank0-axis_offset; i<rank0; i++) {
2867 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2868 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2869 }
2870 SM *= tmpShape0[i];
2871 }
2872 for (int i=axis_offset; i<rank1; i++) {
2873 SR *= tmpShape1[i];
2874 }
2875
2876 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2877 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2878 { // block to limit the scope of out_index
2879 int out_index=0;
2880 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2881 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2882 }
2883
2884 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2885 {
2886 ostringstream os;
2887 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2888 throw DataException(os.str());
2889 }
2890
2891 // Declare output Data object
2892 Data res;
2893
2894 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2895 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2896 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2897 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2898 double *ptr_2 = &(res.getDataAtOffsetRW(0));
2899 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2900 }
2901 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2902
2903 // Prepare the DataConstant input
2904 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2905 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2906
2907 // Borrow DataTagged input from Data object
2908 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2909 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2910
2911 // Prepare a DataTagged output 2
2912 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2913 res.tag();
2914 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2915 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2916
2917 // Prepare offset into DataConstant
2918 int offset_0 = tmp_0->getPointOffset(0,0);
2919 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2920
2921 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2922 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2923
2924 // Compute an MVP for the default
2925 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2926 // Compute an MVP for each tag
2927 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2928 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2929 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2930 tmp_2->addTag(i->first);
2931
2932 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2933 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2934
2935 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2936 }
2937
2938 }
2939 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2940
2941 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2942 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2943 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2944 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2945 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2946 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2947 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2948 int sampleNo_1,dataPointNo_1;
2949 int numSamples_1 = arg_1_Z.getNumSamples();
2950 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2951 int offset_0 = tmp_0->getPointOffset(0,0);
2952 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2953 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2954 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2955 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2956 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2957 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2958 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2959 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2960 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2961 }
2962 }
2963 }
2964 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2965
2966 // Borrow DataTagged input from Data object
2967 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2968 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2969
2970 // Prepare the DataConstant input
2971 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2972 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2973
2974 // Prepare a DataTagged output 2
2975 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2976 res.tag();
2977 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2978 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2979
2980 // Prepare offset into DataConstant
2981 int offset_1 = tmp_1->getPointOffset(0,0);
2982 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2983 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2984 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2985
2986 // Compute an MVP for the default
2987 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2988 // Compute an MVP for each tag
2989 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2990 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2991 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2992
2993 tmp_2->addTag(i->first);
2994 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2995 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2996 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2997 }
2998 }
2999 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
3000
3001 // Borrow DataTagged input from Data object
3002 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3003 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3004
3005 // Borrow DataTagged input from Data object
3006 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3007 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3008
3009 // Prepare a DataTagged output 2
3010 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
3011 res.tag(); // DataTagged output
3012 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3013 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3014
3015 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
3016 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
3017 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3018
3019 // Compute an MVP for the default
3020 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3021 // Merge the tags
3022 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3023 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3024 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3025 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3026 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
3027 }
3028 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3029 tmp_2->addTag(i->first);
3030 }
3031 // Compute an MVP for each tag
3032 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
3033 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
3034 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3035 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
3036 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
3037
3038 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3039 }
3040 }
3041 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
3042
3043 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3044 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3045 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3046 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3047 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3048 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3049 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3050 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3051 int sampleNo_0,dataPointNo_0;
3052 int numSamples_0 = arg_0_Z.getNumSamples();
3053 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3054 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3055 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3056 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
3057 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3058 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3059 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3060 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3061 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3062 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3063 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3064 }
3065 }
3066 }
3067 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
3068
3069 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3070 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3071 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3072 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3073 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3074 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
3075 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3076 int sampleNo_0,dataPointNo_0;
3077 int numSamples_0 = arg_0_Z.getNumSamples();
3078 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3079 int offset_1 = tmp_1->getPointOffset(0,0);
3080 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3081 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3082 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3083 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3084 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3085 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3086 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3087 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3088 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3089 }
3090 }
3091 }
3092 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
3093
3094 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3095 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3096 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3097 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3098 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3099 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3100 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3101 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3102 int sampleNo_0,dataPointNo_0;
3103 int numSamples_0 = arg_0_Z.getNumSamples();
3104 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3105 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3106 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3107 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
3108 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3109 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3110 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3111 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3112 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3113 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3114 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3115 }
3116 }
3117 }
3118 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
3119
3120 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3121 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3122 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3123 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3124 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3125 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3126 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3127 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3128 int sampleNo_0,dataPointNo_0;
3129 int numSamples_0 = arg_0_Z.getNumSamples();
3130 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3131 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3132 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3133 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3134 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3135 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3136 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3137 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3138 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3139 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3140 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3141 }
3142 }
3143 }
3144 else {
3145 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
3146 }
3147
3148 return res;
3149 }
3150
3151 DataAbstract*
3152 Data::borrowData() const
3153 {
3154 return m_data.get();
3155 }
3156
3157 // Not all that happy about returning a non-const from a const
3158 DataAbstract_ptr
3159 Data::borrowDataPtr() const
3160 {
3161 return m_data;
3162 }
3163
3164 // Not all that happy about returning a non-const from a const
3165 DataReady_ptr
3166 Data::borrowReadyPtr() const
3167 {
3168 DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
3169 EsysAssert((dr!=0), "Error - casting to DataReady.");
3170 return dr;
3171 }
3172
3173 std::string
3174 Data::toString() const
3175 {
3176 int localNeedSummary=0;
3177 #ifdef ESYS_MPI
3178 int globalNeedSummary=0;
3179 #endif
3180 if (!m_data->isEmpty() &&
3181 !m_data->isLazy() &&
3182 getLength()>escriptParams.getInt("TOO_MANY_LINES"))
3183 {
3184 localNeedSummary=1;
3185 }
3186
3187 #ifdef ESYS_MPI
3188 MPI_Allreduce( &localNeedSummary, &globalNeedSummary, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3189 localNeedSummary=globalNeedSummary;
3190 #endif
3191
3192 if (localNeedSummary){
3193 stringstream temp;
3194 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
3195 return temp.str();
3196 }
3197 return m_data->toString();
3198 }
3199
3200
3201 // This method is not thread-safe
3202 DataTypes::ValueType::reference
3203 Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
3204 {
3205 checkExclusiveWrite();
3206 return getReady()->getDataAtOffsetRW(i);
3207 }
3208
3209 // This method is not thread-safe
3210 DataTypes::ValueType::const_reference
3211 Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
3212 {
3213 forceResolve();
3214 return getReady()->getDataAtOffsetRO(i);
3215 }
3216
3217
3218 // DataTypes::ValueType::const_reference
3219 // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
3220 // {
3221 // if (isLazy())
3222 // {
3223 // throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
3224 // }
3225 // return getReady()->getDataAtOffsetRO(i);
3226 // }
3227
3228
3229 DataTypes::ValueType::const_reference
3230 Data::getDataPointRO(int sampleNo, int dataPointNo)
3231 {
3232 forceResolve();
3233 if (!isReady())
3234 {
3235 throw DataException("Programmer error -getDataPointRO() not permitted on Lazy Data.");
3236 }
3237 else
3238 {
3239 const DataReady* dr=getReady();
3240 return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
3241 }
3242 }
3243
3244
3245 DataTypes::ValueType::reference
3246 Data::getDataPointRW(int sampleNo, int dataPointNo)
3247 {
3248 checkExclusiveWrite();
3249 DataReady* dr=getReady();
3250 return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
3251 }
3252
3253 Data
3254 Data::interpolateFromTable3DP(bp::object table, double Amin, double Astep,
3255 Data& B, double Bmin, double Bstep,
3256 Data& C, double Cmin, double Cstep,
3257 double undef, bool check_boundaries)
3258 {
3259 WrappedArray t(table);
3260 return interpolateFromTable3D(t, Amin, Astep, undef, B, Bmin, Bstep, C, Cmin, Cstep, check_boundaries);
3261 }
3262
3263 Data
3264 Data::interpolateFromTable2DP(bp::object table, double Amin, double Astep,
3265 Data& B, double Bmin, double Bstep,
3266 double undef, bool check_boundaries)
3267 {
3268 WrappedArray t(table);
3269 return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3270 }
3271
3272 Data
3273 Data::interpolateFromTable1DP(bp::object table, double Amin, double Astep,
3274 double undef,bool check_boundaries)
3275 {
3276 WrappedArray t(table);
3277 return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3278 }
3279
3280
3281 Data
3282 Data::interpolateFromTable1D(const WrappedArray& table, double Amin,
3283 double Astep, double undef, bool check_boundaries)
3284 {
3285 table.convertArray(); // critical! Calling getElt on an unconverted array is not thread safe
3286 int error=0;
3287 if ((getDataPointRank()!=0))
3288 {
3289 throw DataException("Input to 1D interpolation must be scalar");
3290 }
3291 if (table.getRank()!=1)
3292 {
3293 throw DataException("Table for 1D interpolation must be 1D");
3294 }
3295 if (Astep<=0)
3296 {
3297 throw DataException("Astep must be positive&