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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 5 months ago) by caltinay
File size: 138848 byte(s)
Assorted spelling fixes.

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