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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2735 - (show annotations)
Mon Nov 2 02:03:24 2009 UTC (10 years ago) by jfenwick
File size: 96924 byte(s)
Fixed bug where calling minval() reintroduced laziness after it had been removed.

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