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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3923 - (show annotations)
Mon Jul 9 05:12:12 2012 UTC (6 years, 10 months ago) by jfenwick
File size: 123698 byte(s)
Fixing build error.

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