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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3506 - (show annotations)
Wed May 11 01:59:45 2011 UTC (8 years, 2 months ago) by jfenwick
File size: 123306 byte(s)
Use VSL if MKLRANDOM is defined.
Use different seeds for different ranks/threads.
Use different seeds for successive calls with no seed given.

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