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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3786 - (show annotations)
Wed Jan 25 04:46:39 2012 UTC (7 years, 8 months ago) by jfenwick
File size: 123107 byte(s)
Made 2D interpolation call consistant with 1D and 3D.
Doco updates later

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "Data.h"
16
17 #include "DataExpanded.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #include "DataEmpty.h"
21 #include "DataLazy.h"
22 #include "FunctionSpaceFactory.h"
23 #include "AbstractContinuousDomain.h"
24 #include "UnaryFuncs.h"
25 #include "FunctionSpaceException.h"
26 #include "EscriptParams.h"
27
28 #ifdef IKNOWWHATIMDOING
29 #include "Dodgy.h"
30 #endif
31
32 extern "C" {
33 #include "esysUtils/blocktimer.h"
34 }
35
36 #include <fstream>
37 #include <algorithm>
38 #include <vector>
39 #include <functional>
40 #include <sstream> // so we can throw messages about ranks
41
42 #include <boost/python/dict.hpp>
43 #include <boost/python/extract.hpp>
44 #include <boost/python/long.hpp>
45 #include "WrappedArray.h"
46
47 using namespace std;
48 using namespace boost::python;
49 using namespace boost;
50 using namespace escript;
51
52 // ensure the current object is not a DataLazy
53 // The idea was that we could add an optional warning whenever a resolve is forced
54 // #define forceResolve() if (isLazy()) {#resolve();}
55
56 #define AUTOLAZYON escriptParams.getAUTOLAZY()
57 #define MAKELAZYOP(X) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
58 {\
59 DataLazy* c=new DataLazy(borrowDataPtr(),X);\
60 return Data(c);\
61 }
62 #define MAKELAZYOPOFF(X,Y) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
63 {\
64 DataLazy* c=new DataLazy(borrowDataPtr(),X,Y);\
65 return Data(c);\
66 }
67
68 #define MAKELAZYOP2(X,Y,Z) if (isLazy() || (AUTOLAZYON && m_data->isExpanded())) \
69 {\
70 DataLazy* c=new DataLazy(borrowDataPtr(),X,Y,Z);\
71 return Data(c);\
72 }
73
74 #define MAKELAZYBINSELF(R,X) if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
75 {\
76 DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
77 /* m_data=c->getPtr();*/ set_m_data(c->getPtr());\
78 return (*this);\
79 }
80
81 // like the above but returns a new data rather than *this
82 #define MAKELAZYBIN(R,X) if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
83 {\
84 DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
85 return Data(c);\
86 }
87
88 #define MAKELAZYBIN2(L,R,X) if (L.isLazy() || R.isLazy() || (AUTOLAZYON && (L.isExpanded() || R.isExpanded()))) \
89 {\
90 DataLazy* c=new DataLazy(L.borrowDataPtr(),R.borrowDataPtr(),X);\
91 return Data(c);\
92 }
93
94 #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
95
96 namespace
97 {
98
99 template <class ARR>
100 inline
101 boost::python::tuple
102 pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
103 {
104 using namespace boost::python;
105 using boost::python::tuple;
106 using boost::python::list;
107
108 list l;
109 unsigned int dim=shape[0];
110 for (size_t i=0;i<dim;++i)
111 {
112 l.append(v[i+offset]);
113 }
114 return tuple(l);
115 }
116
117 template <class ARR>
118 inline
119 boost::python::tuple
120 pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
121 {
122 using namespace boost::python;
123 using boost::python::tuple;
124 using boost::python::list;
125
126 unsigned int shape0=shape[0];
127 unsigned int shape1=shape[1];
128 list lj;
129 for (size_t j=0;j<shape0;++j)
130 {
131 list li;
132 for (size_t i=0;i<shape1;++i)
133 {
134 li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
135 }
136 lj.append(tuple(li));
137 }
138 return tuple(lj);
139 }
140
141 template <class ARR>
142 inline
143 boost::python::tuple
144 pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
145 {
146 using namespace boost::python;
147 using boost::python::tuple;
148 using boost::python::list;
149
150 unsigned int shape0=shape[0];
151 unsigned int shape1=shape[1];
152 unsigned int shape2=shape[2];
153
154 list lk;
155 for (size_t k=0;k<shape0;++k)
156 {
157 list lj;
158 for (size_t j=0;j<shape1;++j)
159 {
160 list li;
161 for (size_t i=0;i<shape2;++i)
162 {
163 li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
164 }
165 lj.append(tuple(li));
166 }
167 lk.append(tuple(lj));
168 }
169 return tuple(lk);
170 }
171
172 template <class ARR>
173 inline
174 boost::python::tuple
175 pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
176 {
177 using namespace boost::python;
178 using boost::python::tuple;
179 using boost::python::list;
180
181 unsigned int shape0=shape[0];
182 unsigned int shape1=shape[1];
183 unsigned int shape2=shape[2];
184 unsigned int shape3=shape[3];
185
186 list ll;
187 for (size_t l=0;l<shape0;++l)
188 {
189 list lk;
190 for (size_t k=0;k<shape1;++k)
191 {
192 list lj;
193 for (size_t j=0;j<shape2;++j)
194 {
195 list li;
196 for (size_t i=0;i<shape3;++i)
197 {
198 li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
199 }
200 lj.append(tuple(li));
201 }
202 lk.append(tuple(lj));
203 }
204 ll.append(tuple(lk));
205 }
206 return tuple(ll);
207 }
208
209
210 // This should be safer once the DataC RO changes have been brought in
211 template <class ARR>
212 boost::python::tuple
213 pointToTuple( const DataTypes::ShapeType& shape,ARR v)
214 {
215 using namespace boost::python;
216 using boost::python::list;
217 int rank=shape.size();
218 if (rank==0)
219 {
220 return make_tuple(v[0]);
221 }
222 else if (rank==1)
223 {
224 return pointToTuple1(shape,v,0);
225 }
226 else if (rank==2)
227 {
228 return pointToTuple2(shape,v,0);
229 }
230 else if (rank==3)
231 {
232 return pointToTuple3(shape,v,0);
233 }
234 else if (rank==4)
235 {
236 return pointToTuple4(shape,v,0);
237 }
238 else
239 throw DataException("Unknown rank in pointToTuple.");
240 }
241
242 } // anonymous namespace
243
244 Data::Data()
245 : m_shared(false), m_lazy(false)
246 {
247 //
248 // Default data is type DataEmpty
249 DataAbstract* temp=new DataEmpty();
250 // m_data=temp->getPtr();
251 set_m_data(temp->getPtr());
252 m_protected=false;
253 }
254
255 Data::Data(double value,
256 const tuple& shape,
257 const FunctionSpace& what,
258 bool expanded)
259 : m_shared(false), m_lazy(false)
260 {
261 DataTypes::ShapeType dataPointShape;
262 for (int i = 0; i < shape.attr("__len__")(); ++i) {
263 dataPointShape.push_back(extract<const int>(shape[i]));
264 }
265
266 int len = DataTypes::noValues(dataPointShape);
267 DataVector temp_data(len,value,len);
268 initialise(temp_data, dataPointShape, what, expanded);
269 m_protected=false;
270 }
271
272 Data::Data(double value,
273 const DataTypes::ShapeType& dataPointShape,
274 const FunctionSpace& what,
275 bool expanded)
276 : m_shared(false), m_lazy(false)
277 {
278 // int len = DataTypes::noValues(dataPointShape);
279 // DataVector temp_data(len,value,len);
280 initialise(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 if( get_MPIRank()==proc )
1236 {
1237 boost::python::extract<double> dex(v);
1238 if (dex.check())
1239 {
1240 setValueOfDataPoint(id, dex());
1241 }
1242 else
1243 {
1244 setValueOfDataPointToArray(id, v);
1245 }
1246 }
1247 }
1248
1249
1250 void
1251 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1252 {
1253 if (isProtected()) {
1254 throw DataException("Error - attempt to update protected Data object.");
1255 }
1256
1257 WrappedArray w(obj);
1258 //
1259 // check rank
1260 if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1261 throw DataException("Rank of array does not match Data object rank");
1262
1263 //
1264 // check shape of array
1265 for (unsigned int i=0; i<getDataPointRank(); i++) {
1266 if (w.getShape()[i]!=getDataPointShape()[i])
1267 throw DataException("Shape of array does not match Data object rank");
1268 }
1269
1270 exclusiveWrite();
1271
1272 //
1273 // make sure data is expanded:
1274 //
1275 if (!isExpanded()) {
1276 expand();
1277 }
1278 if (getNumDataPointsPerSample()>0) {
1279 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1280 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1281 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1282 } else {
1283 m_data->copyToDataPoint(-1, 0,w);
1284 }
1285 }
1286
1287 void
1288 Data::setValueOfDataPoint(int dataPointNo, const double value)
1289 {
1290 if (isProtected()) {
1291 throw DataException("Error - attempt to update protected Data object.");
1292 }
1293 //
1294 // make sure data is expanded:
1295 exclusiveWrite();
1296 if (!isExpanded()) {
1297 expand();
1298 }
1299 if (getNumDataPointsPerSample()>0) {
1300 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1301 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1302 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1303 } else {
1304 m_data->copyToDataPoint(-1, 0,value);
1305 }
1306 }
1307
1308 const
1309 boost::python::object
1310 Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1311 {
1312 // This could be lazier than it is now
1313 forceResolve();
1314
1315 // copy datapoint into a buffer
1316 // broadcast buffer to all nodes
1317 // convert buffer to tuple
1318 // return tuple
1319
1320 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1321 size_t length=DataTypes::noValues(dataPointShape);
1322
1323 // added for the MPI communication
1324 double *tmpData = new double[length];
1325
1326 // updated for the MPI case
1327 if( get_MPIRank()==procNo ){
1328 if (getNumDataPointsPerSample()>0) {
1329 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1330 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1331 //
1332 // Check a valid sample number has been supplied
1333 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1334 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
1335 }
1336
1337 //
1338 // Check a valid data point number has been supplied
1339 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1340 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
1341 }
1342 // TODO: global error handling
1343 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1344
1345 memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1346 }
1347 }
1348 #ifdef ESYS_MPI
1349 // broadcast the data to all other processes
1350 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1351 #endif
1352
1353 boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1354 delete [] tmpData;
1355 //
1356 // return the loaded array
1357 return t;
1358
1359 }
1360
1361
1362 boost::python::object
1363 Data::integrateToTuple_const() const
1364 {
1365 if (isLazy())
1366 {
1367 throw DataException("Error - cannot integrate for constant lazy data.");
1368 }
1369 return integrateWorker();
1370 }
1371
1372 boost::python::object
1373 Data::integrateToTuple()
1374 {
1375 if (isLazy())
1376 {
1377 expand(); // Can't do a non-resolving version of this without changing the domain object
1378 } // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1379 return integrateWorker();
1380
1381 }
1382
1383 boost::python::object
1384 Data::integrateWorker() const
1385 {
1386 DataTypes::ShapeType shape = getDataPointShape();
1387 int dataPointSize = getDataPointSize();
1388
1389 //
1390 // calculate the integral values
1391 vector<double> integrals(dataPointSize);
1392 vector<double> integrals_local(dataPointSize);
1393 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1394 if (dom==0)
1395 {
1396 throw DataException("Can not integrate over non-continuous domains.");
1397 }
1398 #ifdef ESYS_MPI
1399 dom->setToIntegrals(integrals_local,*this);
1400 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1401 double *tmp = new double[dataPointSize];
1402 double *tmp_local = new double[dataPointSize];
1403 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1404 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1405 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1406 tuple result=pointToTuple(shape,tmp);
1407 delete[] tmp;
1408 delete[] tmp_local;
1409 #else
1410 dom->setToIntegrals(integrals,*this);
1411 /* double *tmp = new double[dataPointSize];
1412 for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1413 tuple result=pointToTuple(shape,integrals);
1414 // delete tmp;
1415 #endif
1416
1417
1418 return result;
1419 }
1420
1421 Data
1422 Data::sin() const
1423 {
1424 MAKELAZYOP(SIN)
1425 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1426 }
1427
1428 Data
1429 Data::cos() const
1430 {
1431 MAKELAZYOP(COS)
1432 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1433 }
1434
1435 Data
1436 Data::tan() const
1437 {
1438 MAKELAZYOP(TAN)
1439 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1440 }
1441
1442 Data
1443 Data::asin() const
1444 {
1445 MAKELAZYOP(ASIN)
1446 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1447 }
1448
1449 Data
1450 Data::acos() const
1451 {
1452 MAKELAZYOP(ACOS)
1453 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1454 }
1455
1456
1457 Data
1458 Data::atan() const
1459 {
1460 MAKELAZYOP(ATAN)
1461 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1462 }
1463
1464 Data
1465 Data::sinh() const
1466 {
1467 MAKELAZYOP(SINH)
1468 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1469 }
1470
1471 Data
1472 Data::cosh() const
1473 {
1474 MAKELAZYOP(COSH)
1475 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1476 }
1477
1478 Data
1479 Data::tanh() const
1480 {
1481 MAKELAZYOP(TANH)
1482 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1483 }
1484
1485
1486 Data
1487 Data::erf() const
1488 {
1489 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1490 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1491 #else
1492 MAKELAZYOP(ERF)
1493 return C_TensorUnaryOperation(*this, ::erf);
1494 #endif
1495 }
1496
1497 Data
1498 Data::asinh() const
1499 {
1500 MAKELAZYOP(ASINH)
1501 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1502 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1503 #else
1504 return C_TensorUnaryOperation(*this, ::asinh);
1505 #endif
1506 }
1507
1508 Data
1509 Data::acosh() const
1510 {
1511 MAKELAZYOP(ACOSH)
1512 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1513 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1514 #else
1515 return C_TensorUnaryOperation(*this, ::acosh);
1516 #endif
1517 }
1518
1519 Data
1520 Data::atanh() const
1521 {
1522 MAKELAZYOP(ATANH)
1523 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1524 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1525 #else
1526 return C_TensorUnaryOperation(*this, ::atanh);
1527 #endif
1528 }
1529
1530 Data
1531 Data::log10() const
1532 {
1533 MAKELAZYOP(LOG10)
1534 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1535 }
1536
1537 Data
1538 Data::log() const
1539 {
1540 MAKELAZYOP(LOG)
1541 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1542 }
1543
1544 Data
1545 Data::sign() const
1546 {
1547 MAKELAZYOP(SIGN)
1548 return C_TensorUnaryOperation(*this, escript::fsign);
1549 }
1550
1551 Data
1552 Data::abs() const
1553 {
1554 MAKELAZYOP(ABS)
1555 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1556 }
1557
1558 Data
1559 Data::neg() const
1560 {
1561 MAKELAZYOP(NEG)
1562 return C_TensorUnaryOperation(*this, negate<double>());
1563 }
1564
1565 Data
1566 Data::pos() const
1567 {
1568 // not doing lazy check here is deliberate.
1569 // since a deep copy of lazy data should be cheap, I'll just let it happen now
1570 Data result;
1571 // perform a deep copy
1572 result.copy(*this);
1573 return result;
1574 }
1575
1576 Data
1577 Data::exp() const
1578 {
1579 MAKELAZYOP(EXP)
1580 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1581 }
1582
1583 Data
1584 Data::sqrt() const
1585 {
1586 MAKELAZYOP(SQRT)
1587 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1588 }
1589
1590 double
1591 Data::Lsup_const() const
1592 {
1593 if (isLazy())
1594 {
1595 throw DataException("Error - cannot compute Lsup for constant lazy data.");
1596 }
1597 return LsupWorker();
1598 }
1599
1600 double
1601 Data::Lsup()
1602 {
1603 if (isLazy())
1604 {
1605 if (!actsExpanded() || CHECK_DO_CRES)
1606 {
1607 resolve();
1608 }
1609 else
1610 {
1611 #ifdef ESYS_MPI
1612 return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1613 #else
1614 return lazyAlgWorker<AbsMax>(0);
1615 #endif
1616 }
1617 }
1618 return LsupWorker();
1619 }
1620
1621 double
1622 Data::sup_const() const
1623 {
1624 if (isLazy())
1625 {
1626 throw DataException("Error - cannot compute sup for constant lazy data.");
1627 }
1628 return supWorker();
1629 }
1630
1631 double
1632 Data::sup()
1633 {
1634 if (isLazy())
1635 {
1636 if (!actsExpanded() || CHECK_DO_CRES)
1637 {
1638 resolve();
1639 }
1640 else
1641 {
1642 #ifdef ESYS_MPI
1643 return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1644 #else
1645 return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1);
1646 #endif
1647 }
1648 }
1649 return supWorker();
1650 }
1651
1652 double
1653 Data::inf_const() const
1654 {
1655 if (isLazy())
1656 {
1657 throw DataException("Error - cannot compute inf for constant lazy data.");
1658 }
1659 return infWorker();
1660 }
1661
1662 double
1663 Data::inf()
1664 {
1665 if (isLazy())
1666 {
1667 if (!actsExpanded() || CHECK_DO_CRES)
1668 {
1669 resolve();
1670 }
1671 else
1672 {
1673 #ifdef ESYS_MPI
1674 return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
1675 #else
1676 return lazyAlgWorker<FMin>(numeric_limits<double>::max());
1677 #endif
1678 }
1679 }
1680 return infWorker();
1681 }
1682
1683 template <class BinaryOp>
1684 double
1685 #ifdef ESYS_MPI
1686 Data::lazyAlgWorker(double init, MPI_Op mpiop_type)
1687 #else
1688 Data::lazyAlgWorker(double init)
1689 #endif
1690 {
1691 if (!isLazy() || !m_data->actsExpanded())
1692 {
1693 throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1694 }
1695 DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1696 EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1697 double val=init;
1698 int i=0;
1699 const size_t numsamples=getNumSamples();
1700 const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1701 BinaryOp operation;
1702 double localValue=0, globalValue;
1703 #pragma omp parallel private(i)
1704 {
1705 double localtot=init;
1706 #pragma omp for schedule(static) private(i)
1707 for (i=0;i<numsamples;++i)
1708 {
1709 size_t roffset=0;
1710 const DataTypes::ValueType* v=dl->resolveSample(i, roffset);
1711 // Now we have the sample, run operation on all points
1712 for (size_t j=0;j<samplesize;++j)
1713 {
1714 localtot=operation(localtot,(*v)[j+roffset]);
1715 }
1716 if (DataMaths::vectorHasNaN(*v,roffset, samplesize))
1717 {
1718 #pragma omp critical
1719 {
1720 localValue=1.0;
1721 }
1722 }
1723 }
1724 #pragma omp critical
1725 val=operation(val,localtot);
1726 }
1727 #ifdef ESYS_MPI
1728 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1729 #else
1730 globalValue=localValue;
1731 #endif
1732 if (globalValue!=0)
1733 {
1734 return makeNaN();
1735 }
1736 #ifdef ESYS_MPI
1737 MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1738 return globalValue;
1739 #else
1740 return val;
1741 #endif
1742 }
1743
1744 // Do not call this on Lazy Data use the proper entry point
1745 double
1746 Data::LsupWorker() const
1747 {
1748 bool haveNaN=getReady()->hasNaN();
1749 double localValue=0;
1750
1751 #ifdef ESYS_MPI
1752 if (haveNaN)
1753 {
1754 localValue=1.0;
1755 }
1756 double globalValue;
1757 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1758 if (globalValue!=0)
1759 {
1760 return makeNaN();
1761 }
1762 #else
1763 if (haveNaN)
1764 {
1765 return makeNaN();
1766 }
1767 #endif
1768
1769 //
1770 // set the initial absolute maximum value to zero
1771
1772 AbsMax abs_max_func;
1773 localValue = algorithm(abs_max_func,0);
1774
1775 #ifdef ESYS_MPI
1776
1777 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1778 return globalValue;
1779 #else
1780 return localValue;
1781 #endif
1782 }
1783
1784 double
1785 Data::supWorker() const
1786 {
1787 bool haveNaN=getReady()->hasNaN();
1788 double localValue=0;
1789
1790 #ifdef ESYS_MPI
1791 if (haveNaN)
1792 {
1793 localValue=1.0;
1794 }
1795 double globalValue;
1796 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1797 if (globalValue!=0)
1798 {
1799 return makeNaN();
1800 }
1801 #else
1802 if (haveNaN)
1803 {
1804 return makeNaN();
1805 }
1806 #endif
1807
1808 //
1809 // set the initial maximum value to min possible double
1810 FMax fmax_func;
1811 localValue = algorithm(fmax_func,numeric_limits<double>::infinity()*-1);
1812 #ifdef ESYS_MPI
1813 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1814 return globalValue;
1815 #else
1816 return localValue;
1817 #endif
1818 }
1819
1820 double
1821 Data::infWorker() const
1822 {
1823 bool haveNaN=getReady()->hasNaN();
1824 double localValue=0;
1825
1826 #ifdef ESYS_MPI
1827 if (haveNaN)
1828 {
1829 localValue=1.0;
1830 }
1831 double globalValue;
1832 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1833 if (globalValue!=0)
1834 {
1835 return makeNaN();
1836 }
1837 #else
1838 if (haveNaN)
1839 {
1840 return makeNaN();
1841 }
1842 #endif
1843 //
1844 // set the initial minimum value to max possible double
1845 FMin fmin_func;
1846 localValue = algorithm(fmin_func,numeric_limits<double>::infinity());
1847 #ifdef ESYS_MPI
1848 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1849 return globalValue;
1850 #else
1851 return localValue;
1852 #endif
1853 }
1854
1855 /* global reduction */
1856
1857
1858 inline Data
1859 Data::minval_nonlazy() const
1860 {
1861 //
1862 // set the initial minimum value to max possible double
1863 FMin fmin_func;
1864 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1865 }
1866
1867
1868 inline Data
1869 Data::maxval_nonlazy() const
1870 {
1871 //
1872 // set the initial maximum value to min possible double
1873 FMax fmax_func;
1874 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1875 }
1876
1877
1878
1879 Data
1880 Data::maxval() const
1881 {
1882 MAKELAZYOP(MAXVAL)
1883 return maxval_nonlazy();
1884 }
1885
1886
1887 Data
1888 Data::minval() const
1889 {
1890 MAKELAZYOP(MINVAL)
1891 return minval_nonlazy();
1892 }
1893
1894
1895 Data
1896 Data::swapaxes(const int axis0, const int axis1) const
1897 {
1898 int axis0_tmp,axis1_tmp;
1899 DataTypes::ShapeType s=getDataPointShape();
1900 DataTypes::ShapeType ev_shape;
1901 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1902 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1903 int rank=getDataPointRank();
1904 if (rank<2) {
1905 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1906 }
1907 if (axis0<0 || axis0>rank-1) {
1908 stringstream e;
1909 e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
1910 throw DataException(e.str());
1911 }
1912 if (axis1<0 || axis1>rank-1) {
1913 stringstream e;
1914 e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
1915 throw DataException(e.str());
1916 }
1917 if (axis0 == axis1) {
1918 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1919 }
1920 MAKELAZYOP2(SWAP,axis0,axis1)
1921 if (axis0 > axis1)
1922 {
1923 axis0_tmp=axis1;
1924 axis1_tmp=axis0;
1925 }
1926 else
1927 {
1928 axis0_tmp=axis0;
1929 axis1_tmp=axis1;
1930 }
1931 for (int i=0; i<rank; i++)
1932 {
1933 if (i == axis0_tmp)
1934 {
1935 ev_shape.push_back(s[axis1_tmp]);
1936 }
1937 else if (i == axis1_tmp)
1938 {
1939 ev_shape.push_back(s[axis0_tmp]);
1940 }
1941 else
1942 {
1943 ev_shape.push_back(s[i]);
1944 }
1945 }
1946 Data ev(0.,ev_shape,getFunctionSpace());
1947 ev.typeMatchRight(*this);
1948 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1949 return ev;
1950 }
1951
1952 Data
1953 Data::symmetric() const
1954 {
1955 // check input
1956 DataTypes::ShapeType s=getDataPointShape();
1957 if (getDataPointRank()==2) {
1958 if(s[0] != s[1])
1959 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1960 }
1961 else if (getDataPointRank()==4) {
1962 if(!(s[0] == s[2] && s[1] == s[3]))
1963 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1964 }
1965 else {
1966 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1967 }
1968 MAKELAZYOP(SYM)
1969 Data ev(0.,getDataPointShape(),getFunctionSpace());
1970 ev.typeMatchRight(*this);
1971 m_data->symmetric(ev.m_data.get());
1972 return ev;
1973 }
1974
1975 Data
1976 Data::nonsymmetric() const
1977 {
1978 MAKELAZYOP(NSYM)
1979 // check input
1980 DataTypes::ShapeType s=getDataPointShape();
1981 if (getDataPointRank()==2) {
1982 if(s[0] != s[1])
1983 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1984 DataTypes::ShapeType ev_shape;
1985 ev_shape.push_back(s[0]);
1986 ev_shape.push_back(s[1]);
1987 Data ev(0.,ev_shape,getFunctionSpace());
1988 ev.typeMatchRight(*this);
1989 m_data->nonsymmetric(ev.m_data.get());
1990 return ev;
1991 }
1992 else if (getDataPointRank()==4) {
1993 if(!(s[0] == s[2] && s[1] == s[3]))
1994 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1995 DataTypes::ShapeType ev_shape;
1996 ev_shape.push_back(s[0]);
1997 ev_shape.push_back(s[1]);
1998 ev_shape.push_back(s[2]);
1999 ev_shape.push_back(s[3]);
2000 Data ev(0.,ev_shape,getFunctionSpace());
2001 ev.typeMatchRight(*this);
2002 m_data->nonsymmetric(ev.m_data.get());
2003 return ev;
2004 }
2005 else {
2006 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
2007 }
2008 }
2009
2010 Data
2011 Data::trace(int axis_offset) const
2012 {
2013 MAKELAZYOPOFF(TRACE,axis_offset)
2014 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
2015 {
2016 throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
2017 }
2018 DataTypes::ShapeType s=getDataPointShape();
2019 if (getDataPointRank()==2) {
2020 DataTypes::ShapeType ev_shape;
2021 Data ev(0.,ev_shape,getFunctionSpace());
2022 ev.typeMatchRight(*this);
2023 m_data->trace(ev.m_data.get(), axis_offset);
2024 return ev;
2025 }
2026 if (getDataPointRank()==3) {
2027 DataTypes::ShapeType ev_shape;
2028 if (axis_offset==0) {
2029 int s2=s[2];
2030 ev_shape.push_back(s2);
2031 }
2032 else if (axis_offset==1) {
2033 int s0=s[0];
2034 ev_shape.push_back(s0);
2035 }
2036 Data ev(0.,ev_shape,getFunctionSpace());
2037 ev.typeMatchRight(*this);
2038 m_data->trace(ev.m_data.get(), axis_offset);
2039 return ev;
2040 }
2041 if (getDataPointRank()==4) {
2042 DataTypes::ShapeType ev_shape;
2043 if (axis_offset==0) {
2044 ev_shape.push_back(s[2]);
2045 ev_shape.push_back(s[3]);
2046 }
2047 else if (axis_offset==1) {
2048 ev_shape.push_back(s[0]);
2049 ev_shape.push_back(s[3]);
2050 }
2051 else if (axis_offset==2) {
2052 ev_shape.push_back(s[0]);
2053 ev_shape.push_back(s[1]);
2054 }
2055 Data ev(0.,ev_shape,getFunctionSpace());
2056 ev.typeMatchRight(*this);
2057 m_data->trace(ev.m_data.get(), axis_offset);
2058 return ev;
2059 }
2060 else {
2061 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
2062 }
2063 }
2064
2065 Data
2066 Data::transpose(int axis_offset) const
2067 {
2068 MAKELAZYOPOFF(TRANS,axis_offset)
2069 DataTypes::ShapeType s=getDataPointShape();
2070 DataTypes::ShapeType ev_shape;
2071 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
2072 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
2073 int rank=getDataPointRank();
2074 if (axis_offset<0 || axis_offset>rank) {
2075 stringstream e;
2076 e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
2077 throw DataException(e.str());
2078 }
2079 for (int i=0; i<rank; i++) {
2080
2081 int index = (axis_offset+i)%rank;
2082 ev_shape.push_back(s[index]); // Append to new shape
2083 }
2084 Data ev(0.,ev_shape,getFunctionSpace());
2085 ev.typeMatchRight(*this);
2086 m_data->transpose(ev.m_data.get(), axis_offset);
2087 return ev;
2088 }
2089
2090 Data
2091 Data::eigenvalues() const
2092 {
2093 if (isLazy())
2094 {
2095 Data temp(*this); // to get around the fact that you can't resolve a const Data
2096 temp.resolve();
2097 return temp.eigenvalues();
2098 }
2099 // check input
2100 DataTypes::ShapeType s=getDataPointShape();
2101 if (getDataPointRank()!=2)
2102 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
2103 if(s[0] != s[1])
2104 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
2105 // create return
2106 DataTypes::ShapeType ev_shape(1,s[0]);
2107 Data ev(0.,ev_shape,getFunctionSpace());
2108 ev.typeMatchRight(*this);
2109 m_data->eigenvalues(ev.m_data.get());
2110 return ev;
2111 }
2112
2113 const boost::python::tuple
2114 Data::eigenvalues_and_eigenvectors(const double tol) const
2115 {
2116 if (isLazy())
2117 {
2118 Data temp(*this); // to get around the fact that you can't resolve a const Data
2119 temp.resolve();
2120 return temp.eigenvalues_and_eigenvectors(tol);
2121 }
2122 DataTypes::ShapeType s=getDataPointShape();
2123 if (getDataPointRank()!=2)
2124 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
2125 if(s[0] != s[1])
2126 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
2127 // create return
2128 DataTypes::ShapeType ev_shape(1,s[0]);
2129 Data ev(0.,ev_shape,getFunctionSpace());
2130 ev.typeMatchRight(*this);
2131 DataTypes::ShapeType V_shape(2,s[0]);
2132 Data V(0.,V_shape,getFunctionSpace());
2133 V.typeMatchRight(*this);
2134 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
2135 return make_tuple(boost::python::object(ev),boost::python::object(V));
2136 }
2137
2138 const boost::python::tuple
2139 Data::minGlobalDataPoint() const
2140 {
2141 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
2142 // abort (for unknown reasons) if there are openmp directives with it in the
2143 // surrounding function
2144
2145 int DataPointNo;
2146 int ProcNo;
2147 calc_minGlobalDataPoint(ProcNo,DataPointNo);
2148 return make_tuple(ProcNo,DataPointNo);
2149 }
2150
2151 void
2152 Data::calc_minGlobalDataPoint(int& ProcNo,
2153 int& DataPointNo) const
2154 {
2155 if (isLazy())
2156 {
2157 Data temp(*this); // to get around the fact that you can't resolve a const Data
2158 temp.resolve();
2159 return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2160 }
2161 int i,j;
2162 int lowi=0,lowj=0;
2163 double min=numeric_limits<double>::max();
2164
2165 Data temp=minval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2166
2167 int numSamples=temp.getNumSamples();
2168 int numDPPSample=temp.getNumDataPointsPerSample();
2169
2170 double local_val, local_min;
2171 #ifdef ESYS_MPI
2172 double next[2];
2173 #endif
2174 int local_lowi=0,local_lowj=0;
2175
2176 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(local_val,local_min)
2177 {
2178 local_min=min;
2179 #pragma omp for private(i,j) schedule(static)
2180 for (i=0; i<numSamples; i++) {
2181 for (j=0; j<numDPPSample; j++) {
2182 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2183 if (local_val<local_min) {
2184 local_min=local_val;
2185 local_lowi=i;
2186 local_lowj=j;
2187 }
2188 }
2189 }
2190 #pragma omp critical
2191 if (local_min<min) { // If we found a smaller value than our sentinel
2192 min=local_min;
2193 lowi=local_lowi;
2194 lowj=local_lowj;
2195 }
2196 }
2197
2198 #ifdef ESYS_MPI
2199 // determine the processor on which the minimum occurs
2200 next[0] = min;
2201 next[1] = numSamples;
2202 int lowProc = 0;
2203 double *globalMins = new double[get_MPISize()*2+1];
2204 /*int error =*/ MPI_Gather (next, 2, MPI_DOUBLE, globalMins, 2, MPI_DOUBLE, 0, get_MPIComm() );
2205
2206 if( get_MPIRank()==0 ){
2207 for (lowProc=0; lowProc<get_MPISize(); lowProc++)
2208 if (globalMins[lowProc*2+1] > 0) break;
2209 min = globalMins[lowProc*2];
2210 for( i=lowProc+1; i<get_MPISize(); i++ )
2211 if( globalMins[i*2+1]>0 && min>globalMins[i*2] ){
2212 lowProc = i;
2213 min = globalMins[i*2];
2214 }
2215 }
2216 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2217 DataPointNo = lowj + lowi * numDPPSample;
2218 MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2219 delete [] globalMins;
2220 ProcNo = lowProc;
2221 #else
2222 ProcNo = 0;
2223 DataPointNo = lowj + lowi * numDPPSample;
2224 #endif
2225 }
2226
2227
2228 const boost::python::tuple
2229 Data::maxGlobalDataPoint() const
2230 {
2231 int DataPointNo;
2232 int ProcNo;
2233 calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2234 return make_tuple(ProcNo,DataPointNo);
2235 }
2236
2237 void
2238 Data::calc_maxGlobalDataPoint(int& ProcNo,
2239 int& DataPointNo) const
2240 {
2241 if (isLazy())
2242 {
2243 Data temp(*this); // to get around the fact that you can't resolve a const Data
2244 temp.resolve();
2245 return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2246 }
2247 int i,j;
2248 int highi=0,highj=0;
2249 //-------------
2250 double max= -numeric_limits<double>::max();
2251
2252 Data temp=maxval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2253
2254 int numSamples=temp.getNumSamples();
2255 int numDPPSample=temp.getNumDataPointsPerSample();
2256
2257 double local_val, local_max;
2258 #ifdef ESYS_MPI
2259 double next[2];
2260 #endif
2261 int local_highi=0,local_highj=0;
2262
2263 #pragma omp parallel firstprivate(local_highi,local_highj) private(local_val,local_max)
2264 {
2265 local_max=max;
2266 #pragma omp for private(i,j) schedule(static)
2267 for (i=0; i<numSamples; i++) {
2268 for (j=0; j<numDPPSample; j++) {
2269 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2270 if (local_val>local_max) {
2271 local_max=local_val;
2272 local_highi=i;
2273 local_highj=j;
2274 }
2275 }
2276 }
2277 #pragma omp critical
2278 if (local_max>max) { // If we found a larger value than our sentinel
2279 max=local_max;
2280 highi=local_highi;
2281 highj=local_highj;
2282 }
2283 }
2284 #ifdef ESYS_MPI
2285 // determine the processor on which the maximum occurs
2286 next[0] = max;
2287 next[1] = numSamples;
2288 int highProc = 0;
2289 double *globalMaxs = new double[get_MPISize()*2+1];
2290 /*int error =*/ MPI_Gather ( next, 2, MPI_DOUBLE, globalMaxs, 2, MPI_DOUBLE, 0, get_MPIComm() );
2291 if( get_MPIRank()==0 ){
2292 for (highProc=0; highProc<get_MPISize(); highProc++)
2293 if (globalMaxs[highProc*2+1] > 0) break;
2294 max = globalMaxs[highProc*2];
2295 for( i=highProc+1; i<get_MPISize(); i++ )
2296 {
2297 if( globalMaxs[i*2+1]>0 && max<globalMaxs[i*2] )
2298 {
2299 highProc = i;
2300 max = globalMaxs[i*2];
2301 }
2302 }
2303 }
2304 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2305 DataPointNo = highj + highi * numDPPSample;
2306 MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2307
2308 delete [] globalMaxs;
2309 ProcNo = highProc;
2310 #else
2311 ProcNo = 0;
2312 DataPointNo = highj + highi * numDPPSample;
2313 #endif
2314 }
2315
2316 void
2317 Data::saveDX(std::string fileName) const
2318 {
2319 if (isEmpty())
2320 {
2321 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2322 }
2323 if (isLazy())
2324 {
2325 Data temp(*this); // to get around the fact that you can't resolve a const Data
2326 temp.resolve();
2327 temp.saveDX(fileName);
2328 return;
2329 }
2330 boost::python::dict args;
2331 args["data"]=boost::python::object(this);
2332 getDomain()->saveDX(fileName,args);
2333 return;
2334 }
2335
2336 void
2337 Data::saveVTK(std::string fileName) const
2338 {
2339 if (isEmpty())
2340 {
2341 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2342 }
2343 if (isLazy())
2344 {
2345 Data temp(*this); // to get around the fact that you can't resolve a const Data
2346 temp.resolve();
2347 temp.saveVTK(fileName);
2348 return;
2349 }
2350 boost::python::dict args;
2351 args["data"]=boost::python::object(this);
2352 getDomain()->saveVTK(fileName,args,"","");
2353 return;
2354 }
2355
2356
2357
2358 Data&
2359 Data::operator+=(const Data& right)
2360 {
2361 if (isProtected()) {
2362 throw DataException("Error - attempt to update protected Data object.");
2363 }
2364 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
2365 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2366 binaryOp(right,plus<double>());
2367 return (*this);
2368 }
2369
2370 Data&
2371 Data::operator+=(const boost::python::object& right)
2372 {
2373 if (isProtected()) {
2374 throw DataException("Error - attempt to update protected Data object.");
2375 }
2376 Data tmp(right,getFunctionSpace(),false);
2377 (*this)+=tmp;
2378 return *this;
2379 }
2380
2381 // Hmmm, operator= makes a deep copy but the copy constructor does not?
2382 Data&
2383 Data::operator=(const Data& other)
2384 {
2385 m_protected=false; // since any changes should be caught by exclusiveWrite();
2386 // m_data=other.m_data;
2387 set_m_data(other.m_data);
2388 return (*this);
2389 }
2390
2391 Data&
2392 Data::operator-=(const Data& right)
2393 {
2394 if (isProtected()) {
2395 throw DataException("Error - attempt to update protected Data object.");
2396 }
2397 MAKELAZYBINSELF(right,SUB)
2398 exclusiveWrite();
2399 binaryOp(right,minus<double>());
2400 return (*this);
2401 }
2402
2403 Data&
2404 Data::operator-=(const boost::python::object& right)
2405 {
2406 if (isProtected()) {
2407 throw DataException("Error - attempt to update protected Data object.");
2408 }
2409 Data tmp(right,getFunctionSpace(),false);
2410 (*this)-=tmp;
2411 return (*this);
2412 }
2413
2414 Data&
2415 Data::operator*=(const Data& right)
2416 {
2417 if (isProtected()) {
2418 throw DataException("Error - attempt to update protected Data object.");
2419 }
2420 MAKELAZYBINSELF(right,MUL)
2421 exclusiveWrite();
2422 binaryOp(right,multiplies<double>());
2423 return (*this);
2424 }
2425
2426 Data&
2427 Data::operator*=(const boost::python::object& right)
2428 {
2429 if (isProtected()) {
2430 throw DataException("Error - attempt to update protected Data object.");
2431 }
2432 Data tmp(right,getFunctionSpace(),false);
2433 (*this)*=tmp;
2434 return (*this);
2435 }
2436
2437 Data&
2438 Data::operator/=(const Data& right)
2439 {
2440 if (isProtected()) {
2441 throw DataException("Error - attempt to update protected Data object.");
2442 }
2443 MAKELAZYBINSELF(right,DIV)
2444 exclusiveWrite();
2445 binaryOp(right,divides<double>());
2446 return (*this);
2447 }
2448
2449 Data&
2450 Data::operator/=(const boost::python::object& right)
2451 {
2452 if (isProtected()) {
2453 throw DataException("Error - attempt to update protected Data object.");
2454 }
2455 Data tmp(right,getFunctionSpace(),false);
2456 (*this)/=tmp;
2457 return (*this);
2458 }
2459
2460 /* Be careful trying to make this operation lazy.
2461 At time of writing, resolve() and resolveSample() do not throw.
2462 Changing this would mean that any resolve call would need to use MPI (to check for global errors)
2463 */
2464 Data
2465 Data::matrixInverse() const
2466 {
2467 if (isLazy())
2468 {
2469 Data d(*this);
2470 d.resolve();
2471 return d.matrixInverse();
2472 }
2473
2474 Data out(0.,getDataPointShape(),getFunctionSpace());
2475 out.typeMatchRight(*this);
2476 int errcode=m_data->matrixInverse(out.getReadyPtr().get());
2477 #ifdef ESYS_MPI
2478 int globalval=0;
2479 MPI_Allreduce( &errcode, &globalval, 1, MPI_INT, MPI_MAX, get_MPIComm() );
2480 errcode=globalval;
2481 #endif
2482 if (errcode)
2483 {
2484 DataMaths::matrixInverseError(errcode); // throws exceptions
2485 }
2486 return out;
2487 }
2488
2489
2490 Data
2491 Data::rpowO(const boost::python::object& left) const
2492 {
2493 Data left_d(left,*this);
2494 return left_d.powD(*this);
2495 }
2496
2497 Data
2498 Data::powO(const boost::python::object& right) const
2499 {
2500 Data tmp(right,getFunctionSpace(),false);
2501 return powD(tmp);
2502 }
2503
2504 Data
2505 Data::powD(const Data& right) const
2506 {
2507 MAKELAZYBIN(right,POW)
2508 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2509 }
2510
2511 //
2512 // NOTE: It is essential to specify the namespace this operator belongs to
2513 Data
2514 escript::operator+(const Data& left, const Data& right)
2515 {
2516 MAKELAZYBIN2(left,right,ADD)
2517 return C_TensorBinaryOperation(left, right, plus<double>());
2518 }
2519
2520 //
2521 // NOTE: It is essential to specify the namespace this operator belongs to
2522 Data
2523 escript::operator-(const Data& left, const Data& right)
2524 {
2525 MAKELAZYBIN2(left,right,SUB)
2526 return C_TensorBinaryOperation(left, right, minus<double>());
2527 }
2528
2529 //
2530 // NOTE: It is essential to specify the namespace this operator belongs to
2531 Data
2532 escript::operator*(const Data& left, const Data& right)
2533 {
2534 MAKELAZYBIN2(left,right,MUL)
2535 return C_TensorBinaryOperation(left, right, multiplies<double>());
2536 }
2537
2538 //
2539 // NOTE: It is essential to specify the namespace this operator belongs to
2540 Data
2541 escript::operator/(const Data& left, const Data& right)
2542 {
2543 MAKELAZYBIN2(left,right,DIV)
2544 return C_TensorBinaryOperation(left, right, divides<double>());
2545 }
2546
2547 //
2548 // NOTE: It is essential to specify the namespace this operator belongs to
2549 Data
2550 escript::operator+(const Data& left, const boost::python::object& right)
2551 {
2552 Data tmp(right,left.getFunctionSpace(),false);
2553 MAKELAZYBIN2(left,tmp,ADD)
2554 return left+tmp;
2555 }
2556
2557 //
2558 // NOTE: It is essential to specify the namespace this operator belongs to
2559 Data
2560 escript::operator-(const Data& left, const boost::python::object& right)
2561 {
2562 Data tmp(right,left.getFunctionSpace(),false);
2563 MAKELAZYBIN2(left,tmp,SUB)
2564 return left-tmp;
2565 }
2566
2567 //
2568 // NOTE: It is essential to specify the namespace this operator belongs to
2569 Data
2570 escript::operator*(const Data& left, const boost::python::object& right)
2571 {
2572 Data tmp(right,left.getFunctionSpace(),false);
2573 MAKELAZYBIN2(left,tmp,MUL)
2574 return left*tmp;
2575 }
2576
2577 //
2578 // NOTE: It is essential to specify the namespace this operator belongs to
2579 Data
2580 escript::operator/(const Data& left, const boost::python::object& right)
2581 {
2582 Data tmp(right,left.getFunctionSpace(),false);
2583 MAKELAZYBIN2(left,tmp,DIV)
2584 return left/tmp;
2585 }
2586
2587 //
2588 // NOTE: It is essential to specify the namespace this operator belongs to
2589 Data
2590 escript::operator+(const boost::python::object& left, const Data& right)
2591 {
2592 Data tmp(left,right.getFunctionSpace(),false);
2593 MAKELAZYBIN2(tmp,right,ADD)
2594 return tmp+right;
2595 }
2596
2597 //
2598 // NOTE: It is essential to specify the namespace this operator belongs to
2599 Data
2600 escript::operator-(const boost::python::object& left, const Data& right)
2601 {
2602 Data tmp(left,right.getFunctionSpace(),false);
2603 MAKELAZYBIN2(tmp,right,SUB)
2604 return tmp-right;
2605 }
2606
2607 //
2608 // NOTE: It is essential to specify the namespace this operator belongs to
2609 Data
2610 escript::operator*(const boost::python::object& left, const Data& right)
2611 {
2612 Data tmp(left,right.getFunctionSpace(),false);
2613 MAKELAZYBIN2(tmp,right,MUL)
2614 return tmp*right;
2615 }
2616
2617 //
2618 // NOTE: It is essential to specify the namespace this operator belongs to
2619 Data
2620 escript::operator/(const boost::python::object& left, const Data& right)
2621 {
2622 Data tmp(left,right.getFunctionSpace(),false);
2623 MAKELAZYBIN2(tmp,right,DIV)
2624 return tmp/right;
2625 }
2626
2627
2628 /* TODO */
2629 /* global reduction */
2630 Data
2631 Data::getItem(const boost::python::object& key) const
2632 {
2633
2634 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2635
2636 if (slice_region.size()!=getDataPointRank()) {
2637 throw DataException("Error - slice size does not match Data rank.");
2638 }
2639
2640 return getSlice(slice_region);
2641 }
2642
2643 /* TODO */
2644 /* global reduction */
2645 Data
2646 Data::getSlice(const DataTypes::RegionType& region) const
2647 {
2648 return Data(*this,region);
2649 }
2650
2651 /* TODO */
2652 /* global reduction */
2653 void
2654 Data::setItemO(const boost::python::object& key,
2655 const boost::python::object& value)
2656 {
2657 Data tempData(value,getFunctionSpace());
2658 setItemD(key,tempData);
2659 }
2660
2661 void
2662 Data::setItemD(const boost::python::object& key,
2663 const Data& value)
2664 {
2665 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2666 if (slice_region.size()!=getDataPointRank()) {
2667 throw DataException("Error - slice size does not match Data rank.");
2668 }
2669 exclusiveWrite();
2670 if (getFunctionSpace()!=value.getFunctionSpace()) {
2671 setSlice(Data(value,getFunctionSpace()),slice_region);
2672 } else {
2673 setSlice(value,slice_region);
2674 }
2675 }
2676
2677 void
2678 Data::setSlice(const Data& value,
2679 const DataTypes::RegionType& region)
2680 {
2681 if (isProtected()) {
2682 throw DataException("Error - attempt to update protected Data object.");
2683 }
2684 forceResolve();
2685 exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2686 Data tempValue(value);
2687 typeMatchLeft(tempValue);
2688 typeMatchRight(tempValue);
2689 getReady()->setSlice(tempValue.m_data.get(),region);
2690 }
2691
2692 void
2693 Data::typeMatchLeft(Data& right) const
2694 {
2695 if (right.isLazy() && !isLazy())
2696 {
2697 right.resolve();
2698 }
2699 if (isExpanded()){
2700 right.expand();
2701 } else if (isTagged()) {
2702 if (right.isConstant()) {
2703 right.tag();
2704 }
2705 }
2706 }
2707
2708 void
2709 Data::typeMatchRight(const Data& right)
2710 {
2711 if (isLazy() && !right.isLazy())
2712 {
2713 resolve();
2714 }
2715 if (isTagged()) {
2716 if (right.isExpanded()) {
2717 expand();
2718 }
2719 } else if (isConstant()) {
2720 if (right.isExpanded()) {
2721 expand();
2722 } else if (right.isTagged()) {
2723 tag();
2724 }
2725 }
2726 }
2727
2728 // The normal TaggedValue adds the tag if it is not already present
2729 // This form does not. It throws instead.
2730 // This is because the names are maintained by the domain and cannot be added
2731 // without knowing the tag number to map it to.
2732 void
2733 Data::setTaggedValueByName(std::string name,
2734 const boost::python::object& value)
2735 {
2736 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2737 forceResolve();
2738 exclusiveWrite();
2739 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2740 setTaggedValue(tagKey,value);
2741 }
2742 else
2743 { // The
2744 throw DataException("Error - unknown tag in setTaggedValueByName.");
2745 }
2746 }
2747
2748 void
2749 Data::setTaggedValue(int tagKey,
2750 const boost::python::object& value)
2751 {
2752 if (isProtected()) {
2753 throw DataException("Error - attempt to update protected Data object.");
2754 }
2755 //
2756 // Ensure underlying data object is of type DataTagged
2757 forceResolve();
2758 exclusiveWrite();
2759 if (isConstant()) tag();
2760 WrappedArray w(value);
2761
2762 DataVector temp_data2;
2763 temp_data2.copyFromArray(w,1);
2764
2765 m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2766 }
2767
2768
2769 void
2770 Data::setTaggedValueFromCPP(int tagKey,
2771 const DataTypes::ShapeType& pointshape,
2772 const DataTypes::ValueType& value,
2773 int dataOffset)
2774 {
2775 if (isProtected()) {
2776 throw DataException("Error - attempt to update protected Data object.");
2777 }
2778 //
2779 // Ensure underlying data object is of type DataTagged
2780 forceResolve();
2781 if (isConstant()) tag();
2782 exclusiveWrite();
2783 //
2784 // Call DataAbstract::setTaggedValue
2785 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2786 }
2787
2788 int
2789 Data::getTagNumber(int dpno)
2790 {
2791 if (isEmpty())
2792 {
2793 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2794 }
2795 return getFunctionSpace().getTagFromDataPointNo(dpno);
2796 }
2797
2798
2799 ostream& escript::operator<<(ostream& o, const Data& data)
2800 {
2801 o << data.toString();
2802 return o;
2803 }
2804
2805 Data
2806 escript::C_GeneralTensorProduct(Data& arg_0,
2807 Data& arg_1,
2808 int axis_offset,
2809 int transpose)
2810 {
2811 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2812 // SM is the product of the last axis_offset entries in arg_0.getShape().
2813
2814 // deal with any lazy data
2815 // if (arg_0.isLazy()) {arg_0.resolve();}
2816 // if (arg_1.isLazy()) {arg_1.resolve();}
2817 if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2818 {
2819 DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2820 return Data(c);
2821 }
2822
2823 // Interpolate if necessary and find an appropriate function space
2824 Data arg_0_Z, arg_1_Z;
2825 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2826 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2827 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2828 arg_1_Z = Data(arg_1);
2829 }
2830 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2831 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2832 arg_0_Z =Data(arg_0);
2833 }
2834 else {
2835 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2836 }
2837 } else {
2838 arg_0_Z = Data(arg_0);
2839 arg_1_Z = Data(arg_1);
2840 }
2841 // Get rank and shape of inputs
2842 int rank0 = arg_0_Z.getDataPointRank();
2843 int rank1 = arg_1_Z.getDataPointRank();
2844 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2845 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2846
2847 // Prepare for the loops of the product and verify compatibility of shapes
2848 int start0=0, start1=0;
2849 if (transpose == 0) {}
2850 else if (transpose == 1) { start0 = axis_offset; }
2851 else if (transpose == 2) { start1 = rank1-axis_offset; }
2852 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2853
2854
2855 // Adjust the shapes for transpose
2856 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2857 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2858 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2859 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2860
2861 #if 0
2862 // For debugging: show shape after transpose
2863 char tmp[100];
2864 std::string shapeStr;
2865 shapeStr = "(";
2866 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2867 shapeStr += ")";
2868 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2869 shapeStr = "(";
2870 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2871 shapeStr += ")";
2872 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2873 #endif
2874
2875 // Prepare for the loops of the product
2876 int SL=1, SM=1, SR=1;
2877 for (int i=0; i<rank0-axis_offset; i++) {
2878 SL *= tmpShape0[i];
2879 }
2880 for (int i=rank0-axis_offset; i<rank0; i++) {
2881 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2882 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2883 }
2884 SM *= tmpShape0[i];
2885 }
2886 for (int i=axis_offset; i<rank1; i++) {
2887 SR *= tmpShape1[i];
2888 }
2889
2890 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2891 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2892 { // block to limit the scope of out_index
2893 int out_index=0;
2894 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2895 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2896 }
2897
2898 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2899 {
2900 ostringstream os;
2901 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2902 throw DataException(os.str());
2903 }
2904
2905 // Declare output Data object
2906 Data res;
2907
2908 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2909 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2910 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2911 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2912 double *ptr_2 = &(res.getDataAtOffsetRW(0));
2913 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2914 }
2915 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2916
2917 // Prepare the DataConstant input
2918 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2919 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2920
2921 // Borrow DataTagged input from Data object
2922 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2923 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2924
2925 // Prepare a DataTagged output 2
2926 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2927 res.tag();
2928 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2929 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2930
2931 // Prepare offset into DataConstant
2932 int offset_0 = tmp_0->getPointOffset(0,0);
2933 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2934
2935 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2936 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2937
2938 // Compute an MVP for the default
2939 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2940 // Compute an MVP for each tag
2941 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2942 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2943 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2944 tmp_2->addTag(i->first);
2945
2946 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2947 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2948
2949 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2950 }
2951
2952 }
2953 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2954
2955 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2956 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2957 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2958 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2959 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2960 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2961 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2962 int sampleNo_1,dataPointNo_1;
2963 int numSamples_1 = arg_1_Z.getNumSamples();
2964 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2965 int offset_0 = tmp_0->getPointOffset(0,0);
2966 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2967 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2968 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2969 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2970 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2971 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2972 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2973 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2974 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2975 }
2976 }
2977
2978 }
2979 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2980
2981 // Borrow DataTagged input from Data object
2982 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2983 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2984
2985 // Prepare the DataConstant input
2986 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2987 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2988
2989 // Prepare a DataTagged output 2
2990 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2991 res.tag();
2992 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2993 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2994
2995 // Prepare offset into DataConstant
2996 int offset_1 = tmp_1->getPointOffset(0,0);
2997 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2998 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2999 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3000
3001 // Compute an MVP for the default
3002 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3003 // Compute an MVP for each tag
3004 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3005 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3006 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3007
3008 tmp_2->addTag(i->first);
3009 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3010 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
3011 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3012 }
3013
3014 }
3015 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
3016
3017 // Borrow DataTagged input from Data object
3018 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3019 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3020
3021 // Borrow DataTagged input from Data object
3022 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3023 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3024
3025 // Prepare a DataTagged output 2
3026 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
3027 res.tag(); // DataTagged output
3028 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
3029 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3030
3031 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
3032 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
3033 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
3034
3035 // Compute an MVP for the default
3036 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3037 // Merge the tags
3038 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
3039 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
3040 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
3041 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
3042 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
3043 }
3044 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
3045 tmp_2->addTag(i->first);
3046 }
3047 // Compute an MVP for each tag
3048 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
3049 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
3050 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
3051 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
3052 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
3053
3054 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3055 }
3056
3057 }
3058 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
3059
3060 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3061 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3062 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
3063 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3064 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3065 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3066 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3067 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3068 int sampleNo_0,dataPointNo_0;
3069 int numSamples_0 = arg_0_Z.getNumSamples();
3070 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3071 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3072 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3073 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
3074 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3075 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3076 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3077 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3078 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3079 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3080 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3081 }
3082 }
3083
3084 }
3085 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
3086
3087 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3088 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3089 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
3090 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3091 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3092 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
3093 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3094 int sampleNo_0,dataPointNo_0;
3095 int numSamples_0 = arg_0_Z.getNumSamples();
3096 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3097 int offset_1 = tmp_1->getPointOffset(0,0);
3098 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3099 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3100 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3101 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3102 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3103 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3104 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3105 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3106 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3107 }
3108 }
3109
3110
3111 }
3112 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
3113
3114 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3115 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3116 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3117 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
3118 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3119 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3120 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
3121 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3122 int sampleNo_0,dataPointNo_0;
3123 int numSamples_0 = arg_0_Z.getNumSamples();
3124 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3125 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3126 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3127 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
3128 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3129 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3130 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3131 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3132 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3133 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3134 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3135 }
3136 }
3137
3138 }
3139 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
3140
3141 // After finding a common function space above the two inputs have the same numSamples and num DPPS
3142 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
3143 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
3144 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
3145 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
3146 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3147 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3148 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
3149 int sampleNo_0,dataPointNo_0;
3150 int numSamples_0 = arg_0_Z.getNumSamples();
3151 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
3152 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
3153 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3154 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3155 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3156 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
3157 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3158 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3159 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
3160 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3161 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
3162 }
3163 }
3164
3165 }
3166 else {
3167 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
3168 }
3169
3170 return res;
3171 }
3172
3173 DataAbstract*
3174 Data::borrowData() const
3175 {
3176 return m_data.get();
3177 }
3178
3179 // Not all that happy about returning a non-const from a const
3180 DataAbstract_ptr
3181 Data::borrowDataPtr() const
3182 {
3183 return m_data;
3184 }
3185
3186 // Not all that happy about returning a non-const from a const
3187 DataReady_ptr
3188 Data::borrowReadyPtr() const
3189 {
3190 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
3191 EsysAssert((dr!=0), "Error - casting to DataReady.");
3192 return dr;
3193 }
3194
3195 std::string
3196 Data::toString() const
3197 {
3198
3199 int localNeedSummary=0;
3200 #ifdef ESYS_MPI
3201 int globalNeedSummary=0;
3202 #endif
3203 if (!m_data->isEmpty() &&
3204 !m_data->isLazy() &&
3205 getLength()>escriptParams.getInt("TOO_MANY_LINES"))
3206 {
3207 localNeedSummary=1;
3208 }
3209
3210 #ifdef ESYS_MPI
3211 MPI_Allreduce( &localNeedSummary, &globalNeedSummary, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3212 localNeedSummary=globalNeedSummary;
3213 #endif
3214
3215 if (localNeedSummary){
3216 stringstream temp;
3217 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
3218 return temp.str();
3219 }
3220 return m_data->toString();
3221 }
3222
3223
3224 // This method is not thread-safe
3225 DataTypes::ValueType::reference
3226 Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
3227 {
3228 checkExclusiveWrite();
3229 return getReady()->getDataAtOffsetRW(i);
3230 }
3231
3232 // This method is not thread-safe
3233 DataTypes::ValueType::const_reference
3234 Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
3235 {
3236 forceResolve();
3237 return getReady()->getDataAtOffsetRO(i);
3238 }
3239
3240
3241 // DataTypes::ValueType::const_reference
3242 // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
3243 // {
3244 // if (isLazy())
3245 // {
3246 // throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
3247 // }
3248 // return getReady()->getDataAtOffsetRO(i);
3249 // }
3250
3251
3252 DataTypes::ValueType::const_reference
3253 Data::getDataPointRO(int sampleNo, int dataPointNo)
3254 {
3255 forceResolve();
3256 if (!isReady())
3257 {
3258 throw DataException("Programmer error -getDataPointRO() not permitted on Lazy Data.");
3259 }
3260 else
3261 {
3262 const DataReady* dr=getReady();
3263 return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
3264 }
3265 }
3266
3267
3268
3269
3270 DataTypes::ValueType::reference
3271 Data::getDataPointRW(int sampleNo, int dataPointNo)
3272 {
3273 checkExclusiveWrite();
3274 DataReady* dr=getReady();
3275 return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
3276 }
3277
3278 Data
3279 Data::interpolateFromTable3DP(boost::python::object table, double Amin, double Astep,
3280 Data& B, double Bmin, double Bstep, Data& C, double Cmin, double Cstep, double undef, bool check_boundaries)
3281 {
3282 WrappedArray t(table);
3283 return interpolateFromTable3D(t, Amin, Astep, undef, B, Bmin, Bstep, C, Cmin, Cstep, check_boundaries);
3284 }
3285
3286 Data
3287 Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3288 Data& B, double Bmin, double Bstep, double undef, bool check_boundaries)
3289 {
3290 WrappedArray t(table);
3291 return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3292 }
3293
3294 Data
3295 Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3296 double undef,bool check_boundaries)
3297 {
3298 WrappedArray t(table);
3299 return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3300 }
3301
3302
3303 Data
3304 Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3305 double undef, bool check_boundaries)
3306 {
3307 table.convertArray(); // critical! Calling getElt on an unconverted array is not thread safe
3308 int error=0;
3309 if ((getDataPointRank()!=0))
3310 {
3311 throw DataException("Input to 1D interpolation must be scalar");
3312 }
3313 if (table.getRank()!=1)
3314 {
3315 throw DataException("Table for 1D interpolation must be 1D");