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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5156 - (show annotations)
Wed Sep 17 00:40:33 2014 UTC (4 years, 7 months ago) by caltinay
File size: 140221 byte(s)
intel 15 doesn't like this.

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