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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3390 - (show annotations)
Thu Dec 2 00:34:37 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 117932 byte(s)
RandomData added

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