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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2723 - (show annotations)
Sun Oct 18 23:44:37 2009 UTC (10 years, 8 months ago) by jfenwick
File size: 96637 byte(s)
Use correct types for MPI op parameters

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,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, 0);
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(), 0);
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 /* TODO */
1735 /* global reduction */
1736 Data
1737 Data::maxval() const
1738 {
1739 MAKELAZYOP(MAXVAL)
1740 //
1741 // set the initial maximum value to min possible double
1742 FMax fmax_func;
1743 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1744 }
1745
1746 Data
1747 Data::minval() const
1748 {
1749 MAKELAZYOP(MINVAL)
1750 //
1751 // set the initial minimum value to max possible double
1752 FMin fmin_func;
1753 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1754 }
1755
1756 Data
1757 Data::swapaxes(const int axis0, const int axis1) const
1758 {
1759 int axis0_tmp,axis1_tmp;
1760 DataTypes::ShapeType s=getDataPointShape();
1761 DataTypes::ShapeType ev_shape;
1762 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1763 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1764 int rank=getDataPointRank();
1765 if (rank<2) {
1766 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1767 }
1768 if (axis0<0 || axis0>rank-1) {
1769 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1770 }
1771 if (axis1<0 || axis1>rank-1) {
1772 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1773 }
1774 if (axis0 == axis1) {
1775 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1776 }
1777 MAKELAZYOP2(SWAP,axis0,axis1)
1778 if (axis0 > axis1)
1779 {
1780 axis0_tmp=axis1;
1781 axis1_tmp=axis0;
1782 }
1783 else
1784 {
1785 axis0_tmp=axis0;
1786 axis1_tmp=axis1;
1787 }
1788 for (int i=0; i<rank; i++)
1789 {
1790 if (i == axis0_tmp)
1791 {
1792 ev_shape.push_back(s[axis1_tmp]);
1793 }
1794 else if (i == axis1_tmp)
1795 {
1796 ev_shape.push_back(s[axis0_tmp]);
1797 }
1798 else
1799 {
1800 ev_shape.push_back(s[i]);
1801 }
1802 }
1803 Data ev(0.,ev_shape,getFunctionSpace());
1804 ev.typeMatchRight(*this);
1805 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1806 return ev;
1807 }
1808
1809 Data
1810 Data::symmetric() const
1811 {
1812 // check input
1813 DataTypes::ShapeType s=getDataPointShape();
1814 if (getDataPointRank()==2) {
1815 if(s[0] != s[1])
1816 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1817 }
1818 else if (getDataPointRank()==4) {
1819 if(!(s[0] == s[2] && s[1] == s[3]))
1820 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1821 }
1822 else {
1823 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1824 }
1825 MAKELAZYOP(SYM)
1826 Data ev(0.,getDataPointShape(),getFunctionSpace());
1827 ev.typeMatchRight(*this);
1828 m_data->symmetric(ev.m_data.get());
1829 return ev;
1830 }
1831
1832 Data
1833 Data::nonsymmetric() const
1834 {
1835 MAKELAZYOP(NSYM)
1836 // check input
1837 DataTypes::ShapeType s=getDataPointShape();
1838 if (getDataPointRank()==2) {
1839 if(s[0] != s[1])
1840 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1841 DataTypes::ShapeType ev_shape;
1842 ev_shape.push_back(s[0]);
1843 ev_shape.push_back(s[1]);
1844 Data ev(0.,ev_shape,getFunctionSpace());
1845 ev.typeMatchRight(*this);
1846 m_data->nonsymmetric(ev.m_data.get());
1847 return ev;
1848 }
1849 else if (getDataPointRank()==4) {
1850 if(!(s[0] == s[2] && s[1] == s[3]))
1851 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1852 DataTypes::ShapeType ev_shape;
1853 ev_shape.push_back(s[0]);
1854 ev_shape.push_back(s[1]);
1855 ev_shape.push_back(s[2]);
1856 ev_shape.push_back(s[3]);
1857 Data ev(0.,ev_shape,getFunctionSpace());
1858 ev.typeMatchRight(*this);
1859 m_data->nonsymmetric(ev.m_data.get());
1860 return ev;
1861 }
1862 else {
1863 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1864 }
1865 }
1866
1867 Data
1868 Data::trace(int axis_offset) const
1869 {
1870 MAKELAZYOPOFF(TRACE,axis_offset)
1871 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1872 {
1873 throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
1874 }
1875 DataTypes::ShapeType s=getDataPointShape();
1876 if (getDataPointRank()==2) {
1877 DataTypes::ShapeType ev_shape;
1878 Data ev(0.,ev_shape,getFunctionSpace());
1879 ev.typeMatchRight(*this);
1880 m_data->trace(ev.m_data.get(), axis_offset);
1881 return ev;
1882 }
1883 if (getDataPointRank()==3) {
1884 DataTypes::ShapeType ev_shape;
1885 if (axis_offset==0) {
1886 int s2=s[2];
1887 ev_shape.push_back(s2);
1888 }
1889 else if (axis_offset==1) {
1890 int s0=s[0];
1891 ev_shape.push_back(s0);
1892 }
1893 Data ev(0.,ev_shape,getFunctionSpace());
1894 ev.typeMatchRight(*this);
1895 m_data->trace(ev.m_data.get(), axis_offset);
1896 return ev;
1897 }
1898 if (getDataPointRank()==4) {
1899 DataTypes::ShapeType ev_shape;
1900 if (axis_offset==0) {
1901 ev_shape.push_back(s[2]);
1902 ev_shape.push_back(s[3]);
1903 }
1904 else if (axis_offset==1) {
1905 ev_shape.push_back(s[0]);
1906 ev_shape.push_back(s[3]);
1907 }
1908 else if (axis_offset==2) {
1909 ev_shape.push_back(s[0]);
1910 ev_shape.push_back(s[1]);
1911 }
1912 Data ev(0.,ev_shape,getFunctionSpace());
1913 ev.typeMatchRight(*this);
1914 m_data->trace(ev.m_data.get(), axis_offset);
1915 return ev;
1916 }
1917 else {
1918 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1919 }
1920 }
1921
1922 Data
1923 Data::transpose(int axis_offset) const
1924 {
1925 MAKELAZYOPOFF(TRANS,axis_offset)
1926 DataTypes::ShapeType s=getDataPointShape();
1927 DataTypes::ShapeType ev_shape;
1928 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1929 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1930 int rank=getDataPointRank();
1931 if (axis_offset<0 || axis_offset>rank) {
1932 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1933 }
1934 for (int i=0; i<rank; i++) {
1935
1936 int index = (axis_offset+i)%rank;
1937 ev_shape.push_back(s[index]); // Append to new shape
1938 }
1939 Data ev(0.,ev_shape,getFunctionSpace());
1940 ev.typeMatchRight(*this);
1941 m_data->transpose(ev.m_data.get(), axis_offset);
1942 return ev;
1943 }
1944
1945 Data
1946 Data::eigenvalues() const
1947 {
1948 if (isLazy())
1949 {
1950 Data temp(*this); // to get around the fact that you can't resolve a const Data
1951 temp.resolve();
1952 return temp.eigenvalues();
1953 }
1954 // check input
1955 DataTypes::ShapeType s=getDataPointShape();
1956 if (getDataPointRank()!=2)
1957 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1958 if(s[0] != s[1])
1959 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1960 // create return
1961 DataTypes::ShapeType ev_shape(1,s[0]);
1962 Data ev(0.,ev_shape,getFunctionSpace());
1963 ev.typeMatchRight(*this);
1964 m_data->eigenvalues(ev.m_data.get());
1965 return ev;
1966 }
1967
1968 const boost::python::tuple
1969 Data::eigenvalues_and_eigenvectors(const double tol) const
1970 {
1971 if (isLazy())
1972 {
1973 Data temp(*this); // to get around the fact that you can't resolve a const Data
1974 temp.resolve();
1975 return temp.eigenvalues_and_eigenvectors(tol);
1976 }
1977 DataTypes::ShapeType s=getDataPointShape();
1978 if (getDataPointRank()!=2)
1979 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1980 if(s[0] != s[1])
1981 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1982 // create return
1983 DataTypes::ShapeType ev_shape(1,s[0]);
1984 Data ev(0.,ev_shape,getFunctionSpace());
1985 ev.typeMatchRight(*this);
1986 DataTypes::ShapeType V_shape(2,s[0]);
1987 Data V(0.,V_shape,getFunctionSpace());
1988 V.typeMatchRight(*this);
1989 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1990 return make_tuple(boost::python::object(ev),boost::python::object(V));
1991 }
1992
1993 const boost::python::tuple
1994 Data::minGlobalDataPoint() const
1995 {
1996 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1997 // abort (for unknown reasons) if there are openmp directives with it in the
1998 // surrounding function
1999
2000 int DataPointNo;
2001 int ProcNo;
2002 calc_minGlobalDataPoint(ProcNo,DataPointNo);
2003 return make_tuple(ProcNo,DataPointNo);
2004 }
2005
2006 void
2007 Data::calc_minGlobalDataPoint(int& ProcNo,
2008 int& DataPointNo) const
2009 {
2010 if (isLazy())
2011 {
2012 Data temp(*this); // to get around the fact that you can't resolve a const Data
2013 temp.resolve();
2014 return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2015 }
2016 int i,j;
2017 int lowi=0,lowj=0;
2018 double min=numeric_limits<double>::max();
2019
2020 Data temp=minval();
2021
2022 int numSamples=temp.getNumSamples();
2023 int numDPPSample=temp.getNumDataPointsPerSample();
2024
2025 double next,local_min;
2026 int local_lowi=0,local_lowj=0;
2027
2028 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
2029 {
2030 local_min=min;
2031 #pragma omp for private(i,j) schedule(static)
2032 for (i=0; i<numSamples; i++) {
2033 for (j=0; j<numDPPSample; j++) {
2034 next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2035 if (next<local_min) {
2036 local_min=next;
2037 local_lowi=i;
2038 local_lowj=j;
2039 }
2040 }
2041 }
2042 #pragma omp critical
2043 if (local_min<min) { // If we found a smaller value than our sentinel
2044 min=local_min;
2045 lowi=local_lowi;
2046 lowj=local_lowj;
2047 }
2048 }
2049
2050 #ifdef PASO_MPI
2051 // determine the processor on which the minimum occurs
2052 next = temp.getDataPointRO(lowi,lowj);
2053 int lowProc = 0;
2054 double *globalMins = new double[get_MPISize()+1];
2055 int error;
2056 error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
2057
2058 if( get_MPIRank()==0 ){
2059 next = globalMins[lowProc];
2060 for( i=1; i<get_MPISize(); i++ )
2061 if( next>globalMins[i] ){
2062 lowProc = i;
2063 next = globalMins[i];
2064 }
2065 }
2066 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2067 DataPointNo = lowj + lowi * numDPPSample;
2068 MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2069 delete [] globalMins;
2070 ProcNo = lowProc;
2071 #else
2072 ProcNo = 0;
2073 DataPointNo = lowj + lowi * numDPPSample;
2074 #endif
2075 }
2076
2077
2078 const boost::python::tuple
2079 Data::maxGlobalDataPoint() const
2080 {
2081 int DataPointNo;
2082 int ProcNo;
2083 calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2084 return make_tuple(ProcNo,DataPointNo);
2085 }
2086
2087 void
2088 Data::calc_maxGlobalDataPoint(int& ProcNo,
2089 int& DataPointNo) const
2090 {
2091 if (isLazy())
2092 {
2093 Data temp(*this); // to get around the fact that you can't resolve a const Data
2094 temp.resolve();
2095 return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2096 }
2097 int i,j;
2098 int highi=0,highj=0;
2099 //-------------
2100 double max= -numeric_limits<double>::max();
2101
2102 Data temp=maxval();
2103
2104 int numSamples=temp.getNumSamples();
2105 int numDPPSample=temp.getNumDataPointsPerSample();
2106
2107 double next,local_max;
2108 int local_highi=0,local_highj=0;
2109
2110 #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2111 {
2112 local_max=max;
2113 #pragma omp for private(i,j) schedule(static)
2114 for (i=0; i<numSamples; i++) {
2115 for (j=0; j<numDPPSample; j++) {
2116 next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2117 if (next>local_max) {
2118 local_max=next;
2119 local_highi=i;
2120 local_highj=j;
2121 }
2122 }
2123 }
2124 #pragma omp critical
2125 if (local_max>max) { // If we found a larger value than our sentinel
2126 max=local_max;
2127 highi=local_highi;
2128 highj=local_highj;
2129 }
2130 }
2131 #ifdef PASO_MPI
2132 // determine the processor on which the maximum occurs
2133 next = temp.getDataPointRO(highi,highj);
2134 int highProc = 0;
2135 double *globalMaxs = new double[get_MPISize()+1];
2136 int error;
2137 error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2138 if( get_MPIRank()==0 ){
2139 next = globalMaxs[highProc];
2140 for( i=1; i<get_MPISize(); i++ )
2141 {
2142 if( next<globalMaxs[i] )
2143 {
2144 highProc = i;
2145 next = globalMaxs[i];
2146 }
2147 }
2148 }
2149 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2150 DataPointNo = highj + highi * numDPPSample;
2151 MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2152
2153 delete [] globalMaxs;
2154 ProcNo = highProc;
2155 #else
2156 ProcNo = 0;
2157 DataPointNo = highj + highi * numDPPSample;
2158 #endif
2159 }
2160
2161 void
2162 Data::saveDX(std::string fileName) const
2163 {
2164 if (isEmpty())
2165 {
2166 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2167 }
2168 if (isLazy())
2169 {
2170 Data temp(*this); // to get around the fact that you can't resolve a const Data
2171 temp.resolve();
2172 temp.saveDX(fileName);
2173 return;
2174 }
2175 boost::python::dict args;
2176 args["data"]=boost::python::object(this);
2177 getDomain()->saveDX(fileName,args);
2178 return;
2179 }
2180
2181 void
2182 Data::saveVTK(std::string fileName) const
2183 {
2184 if (isEmpty())
2185 {
2186 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2187 }
2188 if (isLazy())
2189 {
2190 Data temp(*this); // to get around the fact that you can't resolve a const Data
2191 temp.resolve();
2192 temp.saveVTK(fileName);
2193 return;
2194 }
2195 boost::python::dict args;
2196 args["data"]=boost::python::object(this);
2197 getDomain()->saveVTK(fileName,args,"","");
2198 return;
2199 }
2200
2201
2202
2203 Data&
2204 Data::operator+=(const Data& right)
2205 {
2206 if (isProtected()) {
2207 throw DataException("Error - attempt to update protected Data object.");
2208 }
2209 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
2210 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2211 binaryOp(right,plus<double>());
2212 return (*this);
2213 }
2214
2215 Data&
2216 Data::operator+=(const boost::python::object& right)
2217 {
2218 if (isProtected()) {
2219 throw DataException("Error - attempt to update protected Data object.");
2220 }
2221 Data tmp(right,getFunctionSpace(),false);
2222 (*this)+=tmp;
2223 return *this;
2224 }
2225
2226 // Hmmm, operator= makes a deep copy but the copy constructor does not?
2227 Data&
2228 Data::operator=(const Data& other)
2229 {
2230 m_protected=false; // since any changes should be caught by exclusiveWrite();
2231 // m_data=other.m_data;
2232 set_m_data(other.m_data);
2233 return (*this);
2234 }
2235
2236 Data&
2237 Data::operator-=(const Data& right)
2238 {
2239 if (isProtected()) {
2240 throw DataException("Error - attempt to update protected Data object.");
2241 }
2242 MAKELAZYBINSELF(right,SUB)
2243 exclusiveWrite();
2244 binaryOp(right,minus<double>());
2245 return (*this);
2246 }
2247
2248 Data&
2249 Data::operator-=(const boost::python::object& right)
2250 {
2251 if (isProtected()) {
2252 throw DataException("Error - attempt to update protected Data object.");
2253 }
2254 Data tmp(right,getFunctionSpace(),false);
2255 (*this)-=tmp;
2256 return (*this);
2257 }
2258
2259 Data&
2260 Data::operator*=(const Data& right)
2261 {
2262 if (isProtected()) {
2263 throw DataException("Error - attempt to update protected Data object.");
2264 }
2265 MAKELAZYBINSELF(right,MUL)
2266 exclusiveWrite();
2267 binaryOp(right,multiplies<double>());
2268 return (*this);
2269 }
2270
2271 Data&
2272 Data::operator*=(const boost::python::object& right)
2273 {
2274 if (isProtected()) {
2275 throw DataException("Error - attempt to update protected Data object.");
2276 }
2277 Data tmp(right,getFunctionSpace(),false);
2278 (*this)*=tmp;
2279 return (*this);
2280 }
2281
2282 Data&
2283 Data::operator/=(const Data& right)
2284 {
2285 if (isProtected()) {
2286 throw DataException("Error - attempt to update protected Data object.");
2287 }
2288 MAKELAZYBINSELF(right,DIV)
2289 exclusiveWrite();
2290 binaryOp(right,divides<double>());
2291 return (*this);
2292 }
2293
2294 Data&
2295 Data::operator/=(const boost::python::object& right)
2296 {
2297 if (isProtected()) {
2298 throw DataException("Error - attempt to update protected Data object.");
2299 }
2300 Data tmp(right,getFunctionSpace(),false);
2301 (*this)/=tmp;
2302 return (*this);
2303 }
2304
2305 Data
2306 Data::rpowO(const boost::python::object& left) const
2307 {
2308 Data left_d(left,*this);
2309 return left_d.powD(*this);
2310 }
2311
2312 Data
2313 Data::powO(const boost::python::object& right) const
2314 {
2315 Data tmp(right,getFunctionSpace(),false);
2316 return powD(tmp);
2317 }
2318
2319 Data
2320 Data::powD(const Data& right) const
2321 {
2322 MAKELAZYBIN(right,POW)
2323 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2324 }
2325
2326 //
2327 // NOTE: It is essential to specify the namespace this operator belongs to
2328 Data
2329 escript::operator+(const Data& left, const Data& right)
2330 {
2331 MAKELAZYBIN2(left,right,ADD)
2332 return C_TensorBinaryOperation(left, right, plus<double>());
2333 }
2334
2335 //
2336 // NOTE: It is essential to specify the namespace this operator belongs to
2337 Data
2338 escript::operator-(const Data& left, const Data& right)
2339 {
2340 MAKELAZYBIN2(left,right,SUB)
2341 return C_TensorBinaryOperation(left, right, minus<double>());
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,MUL)
2350 return C_TensorBinaryOperation(left, right, multiplies<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,DIV)
2359 return C_TensorBinaryOperation(left, right, divides<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 boost::python::object& right)
2366 {
2367 Data tmp(right,left.getFunctionSpace(),false);
2368 MAKELAZYBIN2(left,tmp,ADD)
2369 return left+tmp;
2370 }
2371
2372 //
2373 // NOTE: It is essential to specify the namespace this operator belongs to
2374 Data
2375 escript::operator-(const Data& left, const boost::python::object& right)
2376 {
2377 Data tmp(right,left.getFunctionSpace(),false);
2378 MAKELAZYBIN2(left,tmp,SUB)
2379 return left-tmp;
2380 }
2381
2382 //
2383 // NOTE: It is essential to specify the namespace this operator belongs to
2384 Data
2385 escript::operator*(const Data& left, const boost::python::object& right)
2386 {
2387 Data tmp(right,left.getFunctionSpace(),false);
2388 MAKELAZYBIN2(left,tmp,MUL)
2389 return left*tmp;
2390 }
2391
2392 //
2393 // NOTE: It is essential to specify the namespace this operator belongs to
2394 Data
2395 escript::operator/(const Data& left, const boost::python::object& right)
2396 {
2397 Data tmp(right,left.getFunctionSpace(),false);
2398 MAKELAZYBIN2(left,tmp,DIV)
2399 return left/tmp;
2400 }
2401
2402 //
2403 // NOTE: It is essential to specify the namespace this operator belongs to
2404 Data
2405 escript::operator+(const boost::python::object& left, const Data& right)
2406 {
2407 Data tmp(left,right.getFunctionSpace(),false);
2408 MAKELAZYBIN2(tmp,right,ADD)
2409 return tmp+right;
2410 }
2411
2412 //
2413 // NOTE: It is essential to specify the namespace this operator belongs to
2414 Data
2415 escript::operator-(const boost::python::object& left, const Data& right)
2416 {
2417 Data tmp(left,right.getFunctionSpace(),false);
2418 MAKELAZYBIN2(tmp,right,SUB)
2419 return tmp-right;
2420 }
2421
2422 //
2423 // NOTE: It is essential to specify the namespace this operator belongs to
2424 Data
2425 escript::operator*(const boost::python::object& left, const Data& right)
2426 {
2427 Data tmp(left,right.getFunctionSpace(),false);
2428 MAKELAZYBIN2(tmp,right,MUL)
2429 return tmp*right;
2430 }
2431
2432 //
2433 // NOTE: It is essential to specify the namespace this operator belongs to
2434 Data
2435 escript::operator/(const boost::python::object& left, const Data& right)
2436 {
2437 Data tmp(left,right.getFunctionSpace(),false);
2438 MAKELAZYBIN2(tmp,right,DIV)
2439 return tmp/right;
2440 }
2441
2442
2443 /* TODO */
2444 /* global reduction */
2445 Data
2446 Data::getItem(const boost::python::object& key) const
2447 {
2448
2449 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2450
2451 if (slice_region.size()!=getDataPointRank()) {
2452 throw DataException("Error - slice size does not match Data rank.");
2453 }
2454
2455 return getSlice(slice_region);
2456 }
2457
2458 /* TODO */
2459 /* global reduction */
2460 Data
2461 Data::getSlice(const DataTypes::RegionType& region) const
2462 {
2463 return Data(*this,region);
2464 }
2465
2466 /* TODO */
2467 /* global reduction */
2468 void
2469 Data::setItemO(const boost::python::object& key,
2470 const boost::python::object& value)
2471 {
2472 Data tempData(value,getFunctionSpace());
2473 setItemD(key,tempData);
2474 }
2475
2476 void
2477 Data::setItemD(const boost::python::object& key,
2478 const Data& value)
2479 {
2480 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2481 if (slice_region.size()!=getDataPointRank()) {
2482 throw DataException("Error - slice size does not match Data rank.");
2483 }
2484 exclusiveWrite();
2485 if (getFunctionSpace()!=value.getFunctionSpace()) {
2486 setSlice(Data(value,getFunctionSpace()),slice_region);
2487 } else {
2488 setSlice(value,slice_region);
2489 }
2490 }
2491
2492 void
2493 Data::setSlice(const Data& value,
2494 const DataTypes::RegionType& region)
2495 {
2496 if (isProtected()) {
2497 throw DataException("Error - attempt to update protected Data object.");
2498 }
2499 forceResolve();
2500 exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2501 Data tempValue(value);
2502 typeMatchLeft(tempValue);
2503 typeMatchRight(tempValue);
2504 getReady()->setSlice(tempValue.m_data.get(),region);
2505 }
2506
2507 void
2508 Data::typeMatchLeft(Data& right) const
2509 {
2510 if (right.isLazy() && !isLazy())
2511 {
2512 right.resolve();
2513 }
2514 if (isExpanded()){
2515 right.expand();
2516 } else if (isTagged()) {
2517 if (right.isConstant()) {
2518 right.tag();
2519 }
2520 }
2521 }
2522
2523 void
2524 Data::typeMatchRight(const Data& right)
2525 {
2526 if (isLazy() && !right.isLazy())
2527 {
2528 resolve();
2529 }
2530 if (isTagged()) {
2531 if (right.isExpanded()) {
2532 expand();
2533 }
2534 } else if (isConstant()) {
2535 if (right.isExpanded()) {
2536 expand();
2537 } else if (right.isTagged()) {
2538 tag();
2539 }
2540 }
2541 }
2542
2543 // The normal TaggedValue adds the tag if it is not already present
2544 // This form does not. It throws instead.
2545 // This is because the names are maintained by the domain and cannot be added
2546 // without knowing the tag number to map it to.
2547 void
2548 Data::setTaggedValueByName(std::string name,
2549 const boost::python::object& value)
2550 {
2551 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2552 forceResolve();
2553 exclusiveWrite();
2554 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2555 setTaggedValue(tagKey,value);
2556 }
2557 else
2558 { // The
2559 throw DataException("Error - unknown tag in setTaggedValueByName.");
2560 }
2561 }
2562
2563 void
2564 Data::setTaggedValue(int tagKey,
2565 const boost::python::object& value)
2566 {
2567 if (isProtected()) {
2568 throw DataException("Error - attempt to update protected Data object.");
2569 }
2570 //
2571 // Ensure underlying data object is of type DataTagged
2572 forceResolve();
2573 exclusiveWrite();
2574 if (isConstant()) tag();
2575 WrappedArray w(value);
2576
2577 DataVector temp_data2;
2578 temp_data2.copyFromArray(w,1);
2579
2580 m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2581 }
2582
2583
2584 void
2585 Data::setTaggedValueFromCPP(int tagKey,
2586 const DataTypes::ShapeType& pointshape,
2587 const DataTypes::ValueType& value,
2588 int dataOffset)
2589 {
2590 if (isProtected()) {
2591 throw DataException("Error - attempt to update protected Data object.");
2592 }
2593 //
2594 // Ensure underlying data object is of type DataTagged
2595 forceResolve();
2596 if (isConstant()) tag();
2597 exclusiveWrite();
2598 //
2599 // Call DataAbstract::setTaggedValue
2600 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2601 }
2602
2603 int
2604 Data::getTagNumber(int dpno)
2605 {
2606 if (isEmpty())
2607 {
2608 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2609 }
2610 return getFunctionSpace().getTagFromDataPointNo(dpno);
2611 }
2612
2613
2614 ostream& escript::operator<<(ostream& o, const Data& data)
2615 {
2616 o << data.toString();
2617 return o;
2618 }
2619
2620 Data
2621 escript::C_GeneralTensorProduct(Data& arg_0,
2622 Data& arg_1,
2623 int axis_offset,
2624 int transpose)
2625 {
2626 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2627 // SM is the product of the last axis_offset entries in arg_0.getShape().
2628
2629 // deal with any lazy data
2630 // if (arg_0.isLazy()) {arg_0.resolve();}
2631 // if (arg_1.isLazy()) {arg_1.resolve();}
2632 if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2633 {
2634 DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2635 return Data(c);
2636 }
2637
2638 // Interpolate if necessary and find an appropriate function space
2639 Data arg_0_Z, arg_1_Z;
2640 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2641 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2642 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2643 arg_1_Z = Data(arg_1);
2644 }
2645 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2646 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2647 arg_0_Z =Data(arg_0);
2648 }
2649 else {
2650 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2651 }
2652 } else {
2653 arg_0_Z = Data(arg_0);
2654 arg_1_Z = Data(arg_1);
2655 }
2656 // Get rank and shape of inputs
2657 int rank0 = arg_0_Z.getDataPointRank();
2658 int rank1 = arg_1_Z.getDataPointRank();
2659 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2660 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2661
2662 // Prepare for the loops of the product and verify compatibility of shapes
2663 int start0=0, start1=0;
2664 if (transpose == 0) {}
2665 else if (transpose == 1) { start0 = axis_offset; }
2666 else if (transpose == 2) { start1 = rank1-axis_offset; }
2667 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2668
2669
2670 // Adjust the shapes for transpose
2671 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2672 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2673 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2674 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2675
2676 #if 0
2677 // For debugging: show shape after transpose
2678 char tmp[100];
2679 std::string shapeStr;
2680 shapeStr = "(";
2681 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2682 shapeStr += ")";
2683 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2684 shapeStr = "(";
2685 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2686 shapeStr += ")";
2687 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2688 #endif
2689
2690 // Prepare for the loops of the product
2691 int SL=1, SM=1, SR=1;
2692 for (int i=0; i<rank0-axis_offset; i++) {
2693 SL *= tmpShape0[i];
2694 }
2695 for (int i=rank0-axis_offset; i<rank0; i++) {
2696 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2697 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2698 }
2699 SM *= tmpShape0[i];
2700 }
2701 for (int i=axis_offset; i<rank1; i++) {
2702 SR *= tmpShape1[i];
2703 }
2704
2705 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2706 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2707 { // block to limit the scope of out_index
2708 int out_index=0;
2709 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2710 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2711 }
2712
2713 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2714 {
2715 ostringstream os;
2716 os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2717 throw DataException(os.str());
2718 }
2719
2720 // Declare output Data object
2721 Data res;
2722
2723 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2724 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2725 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2726 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2727 double *ptr_2 = &(res.getDataAtOffsetRW(0));
2728 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2729 }
2730 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2731
2732 // Prepare the DataConstant input
2733 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2734 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2735
2736 // Borrow DataTagged input from Data object
2737 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2738 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2739
2740 // Prepare a DataTagged output 2
2741 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2742 res.tag();
2743 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2744 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2745
2746 // Prepare offset into DataConstant
2747 int offset_0 = tmp_0->getPointOffset(0,0);
2748 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2749
2750 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2751 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2752
2753 // Compute an MVP for the default
2754 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2755 // Compute an MVP for each tag
2756 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2757 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2758 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2759 tmp_2->addTag(i->first);
2760
2761 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2762 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2763
2764 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2765 }
2766
2767 }
2768 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2769
2770 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2771 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2772 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2773 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2774 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2775 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2776 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2777 int sampleNo_1,dataPointNo_1;
2778 int numSamples_1 = arg_1_Z.getNumSamples();
2779 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2780 int offset_0 = tmp_0->getPointOffset(0,0);
2781 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2782 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2783 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2784 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2785 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2786 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2787 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2788 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2789 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2790 }
2791 }
2792
2793 }
2794 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2795
2796 // Borrow DataTagged input from Data object
2797 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2798 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2799
2800 // Prepare the DataConstant input
2801 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2802 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2803
2804 // Prepare a DataTagged output 2
2805 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2806 res.tag();
2807 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2808 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2809
2810 // Prepare offset into DataConstant
2811 int offset_1 = tmp_1->getPointOffset(0,0);
2812 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2813 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2814 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2815
2816 // Compute an MVP for the default
2817 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2818 // Compute an MVP for each tag
2819 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2820 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2821 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2822
2823 tmp_2->addTag(i->first);
2824 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2825 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2826 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2827 }
2828
2829 }
2830 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2831
2832 // Borrow DataTagged input from Data object
2833 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2834 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2835
2836 // Borrow DataTagged input from Data object
2837 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2838 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2839
2840 // Prepare a DataTagged output 2
2841 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2842 res.tag(); // DataTagged output
2843 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2844 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2845
2846 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2847 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2848 double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2849
2850 // Compute an MVP for the default
2851 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2852 // Merge the tags
2853 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2854 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2855 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2856 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2857 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2858 }
2859 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2860 tmp_2->addTag(i->first);
2861 }
2862 // Compute an MVP for each tag
2863 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2864 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2865 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2866 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2867 double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2868
2869 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2870 }
2871
2872 }
2873 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2874
2875 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2876 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2877 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2878 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2879 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2880 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2881 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2882 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2883 int sampleNo_0,dataPointNo_0;
2884 int numSamples_0 = arg_0_Z.getNumSamples();
2885 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2886 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2887 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2888 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2889 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2890 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2891 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2892 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2893 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2894 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2895 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2896 }
2897 }
2898
2899 }
2900 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2901
2902 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2903 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2904 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2905 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2906 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2907 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2908 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2909 int sampleNo_0,dataPointNo_0;
2910 int numSamples_0 = arg_0_Z.getNumSamples();
2911 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2912 int offset_1 = tmp_1->getPointOffset(0,0);
2913 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2914 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2915 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2916 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2917 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2918 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2919 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2920 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2921 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2922 }
2923 }
2924
2925
2926 }
2927 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2928
2929 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2930 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2931 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2932 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2933 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2934 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2935 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2936 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2937 int sampleNo_0,dataPointNo_0;
2938 int numSamples_0 = arg_0_Z.getNumSamples();
2939 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2940 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2941 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2942 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2943 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2944 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2945 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2946 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2947 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2948 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2949 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2950 }
2951 }
2952
2953 }
2954 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2955
2956 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2957 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2958 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2959 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2960 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2961 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2962 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2963 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2964 int sampleNo_0,dataPointNo_0;
2965 int numSamples_0 = arg_0_Z.getNumSamples();
2966 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2967 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2968 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2969 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2970 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2971 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2972 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2973 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2974 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2975 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2976 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2977 }
2978 }
2979
2980 }
2981 else {
2982 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2983 }
2984
2985 return res;
2986 }
2987
2988 DataAbstract*
2989 Data::borrowData() const
2990 {
2991 return m_data.get();
2992 }
2993
2994 // Not all that happy about returning a non-const from a const
2995 DataAbstract_ptr
2996 Data::borrowDataPtr() const
2997 {
2998 return m_data;
2999 }
3000
3001 // Not all that happy about returning a non-const from a const
3002 DataReady_ptr
3003 Data::borrowReadyPtr() const
3004 {
3005 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
3006 EsysAssert((dr!=0), "Error - casting to DataReady.");
3007 return dr;
3008 }
3009
3010 std::string
3011 Data::toString() const
3012 {
3013 if (!m_data->isEmpty() &&
3014 !m_data->isLazy() &&
3015 getLength()>escriptParams.getInt("TOO_MANY_LINES"))
3016 {
3017 stringstream temp;
3018 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
3019 return temp.str();
3020 }
3021 return m_data->toString();
3022 }
3023
3024
3025 // This method is not thread-safe
3026 DataTypes::ValueType::reference
3027 Data::getDataAtOffsetRW(DataTypes::ValueType::size_type i)
3028 {
3029 checkExclusiveWrite();
3030 return getReady()->getDataAtOffsetRW(i);
3031 }
3032
3033 // This method is not thread-safe
3034 DataTypes::ValueType::const_reference
3035 Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i)
3036 {
3037 forceResolve();
3038 return getReady()->getDataAtOffsetRO(i);
3039 }
3040
3041
3042 // DataTypes::ValueType::const_reference
3043 // Data::getDataAtOffsetRO(DataTypes::ValueType::size_type i) const
3044 // {
3045 // if (isLazy())
3046 // {
3047 // throw DataException("Programmer error - getDataAtOffsetRO() not permitted on Lazy Data (object is const which prevents resolving).");
3048 // }
3049 // return getReady()->getDataAtOffsetRO(i);
3050 // }
3051
3052
3053 DataTypes::ValueType::const_reference
3054 Data::getDataPointRO(int sampleNo, int dataPointNo)
3055 {
3056 forceResolve();
3057 if (!isReady())
3058 {
3059 throw DataException("Programmer error -getDataPointRO() not permitted on Lazy Data.");
3060 }
3061 else
3062 {
3063 const DataReady* dr=getReady();
3064 return dr->getDataAtOffsetRO(dr->getPointOffset(sampleNo, dataPointNo));
3065 }
3066 }
3067
3068
3069
3070
3071 DataTypes::ValueType::reference
3072 Data::getDataPointRW(int sampleNo, int dataPointNo)
3073 {
3074 checkExclusiveWrite();
3075 DataReady* dr=getReady();
3076 return dr->getDataAtOffsetRW(dr->getPointOffset(sampleNo, dataPointNo));
3077 }
3078
3079 BufferGroup*
3080 Data::allocSampleBuffer() const
3081 {
3082 if (isLazy())
3083 {
3084 #ifdef _OPENMP
3085 int tnum=omp_get_max_threads();
3086 #else
3087 int tnum=1;
3088 #endif
3089 return new BufferGroup(getSampleBufferSize(),tnum);
3090 }
3091 else
3092 {
3093 return NULL;
3094 }
3095 }
3096
3097 void
3098 Data::freeSampleBuffer(BufferGroup* bufferg)
3099 {
3100 if (bufferg!=0)
3101 {
3102 delete bufferg;
3103 }
3104 }
3105
3106
3107 Data
3108 Data::interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
3109 Data& B, double Bmin, double Bstep, double undef, bool check_boundaries)
3110 {
3111 WrappedArray t(table);
3112 return interpolateFromTable2D(t, Amin, Astep, undef, B, Bmin, Bstep,check_boundaries);
3113 }
3114
3115 Data
3116 Data::interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
3117 double undef,bool check_boundaries)
3118 {
3119 WrappedArray t(table);
3120 return interpolateFromTable1D(t, Amin, Astep, undef, check_boundaries);
3121 }
3122
3123
3124 Data
3125 Data::interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
3126 double undef, bool check_boundaries)
3127 {
3128 table.convertArray(); // critical! Calling getElt on an unconverted array is not thread safe
3129 int error=0;
3130 if ((getDataPointRank()!=0))
3131 {
3132 throw DataException("Input to 1D interpolation must be scalar");
3133 }
3134 if (table.getRank()!=1)
3135 {
3136 throw DataException("Table for 1D interpolation must be 1D");
3137 }
3138 if (Astep<=0)
3139 {
3140 throw DataException("Astep must be positive");
3141 }
3142 if (!isExpanded())
3143 {
3144 expand();
3145 }
3146 Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3147 try
3148 {
3149 int numpts=getNumDataPoints();
3150 const DataVector& adat=getReady()->getVectorRO();
3151 DataVector& rdat=res.getReady()->getVectorRW();
3152 int twidth=table.getShape()[0]-1;
3153 bool haserror=false;
3154 int l=0;
3155 #pragma omp parallel for private(l) schedule(static)
3156 for (l=0;l<numpts; ++l)
3157 {
3158 #pragma omp flush(haserror) // In case haserror was in register
3159 if (!haserror)
3160 {
3161 int lerror=0;
3162 try
3163 {
3164 do // so we can use break
3165 {
3166 double a=adat[l];
3167 int x=static_cast<int>(((a-Amin)/Astep));
3168 if (check_boundaries) {
3169 if ((a<Amin) || (x<0))
3170 {
3171 lerror=1;
3172 break;
3173 }
3174 if (a>Amin+Astep*twidth)
3175 {
3176 lerror=4;
3177 break;
3178 }
3179 }
3180 if (x<0) x=0;
3181 if (x>twidth) x=twidth;
3182
3183 if (x==twidth) // value is on the far end of the table
3184 {
3185 double e=table.getElt(x);
3186 if (e>undef)
3187 {
3188 lerror=2;
3189 break;
3190 }
3191 rdat[l]=e;
3192 }
3193 else // x and y are in bounds
3194 {
3195 double e=table.getElt(x);
3196 double w=table.getElt(x+1);
3197 if ((e>undef) || (w>undef))
3198 {
3199 lerror=2;
3200 break;
3201 }
3202 // map x*Astep <= a << (x+1)*Astep to [-1,1]
3203 double la = 2.0*(a-Amin-(x*Astep))/Astep-1;
3204 rdat[l]=((1-la)*e + (1+la)*w)/2;
3205 }
3206 } while (false);
3207 } catch (DataException d)
3208 {
3209 lerror=3;
3210 }
3211 if (lerror!=0)
3212 {
3213 #pragma omp critical // Doco says there is a flush associated with critical
3214 {
3215 haserror=true; // We only care that one error is recorded. We don't care which
3216 error=lerror; // one
3217 }
3218 }
3219 } // if (!error)
3220 } // parallelised for
3221 } catch (DataException d)
3222 {
3223 error=3; // this is outside the parallel region so assign directly
3224 }
3225 #ifdef PASO_MPI
3226 int rerror=0;
3227 MPI_Allreduce( &error, &rerror, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3228 error=rerror;
3229 #endif
3230 if (error)
3231 {
3232 switch (error)
3233 {
3234 case 1: throw DataException("Value below lower table range.");
3235 case 2: throw DataException("Interpolated value too large");
3236 case 4: throw DataException("Value greater than upper table range.");
3237 default:
3238 throw DataException("Unknown error in interpolation");
3239 }
3240 }
3241 return res;
3242 }
3243
3244
3245 Data
3246 Data::interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
3247 double undef, Data& B, double Bmin, double Bstep, bool check_boundaries)
3248 {
3249 table.convertArray(); // critical! Calling getElt on an unconverted array is not thread safe
3250 int error=0;
3251 if ((getDataPointRank()!=0) || (B.getDataPointRank()!=0))
3252 {
3253 throw DataException("Inputs to 2D interpolation must be scalar");
3254 }
3255 if (table.getRank()!=2)
3256 {
3257 throw DataException("Table for 2D interpolation must be 2D");
3258 }
3259 if ((Astep<=0) || (Bstep<=0))
3260 {
3261 throw DataException("Astep and Bstep must be postive");
3262 }
3263 if (getFunctionSpace()!=B.getFunctionSpace())
3264 {
3265 Data n=B.interpolate(getFunctionSpace());
3266 return interpolateFromTable2D(table, Amin, Astep, undef,
3267 n , Bmin, Bstep, check_boundaries);
3268 }
3269 if (!isExpanded())
3270 {
3271 expand();
3272 }
3273 if (!B.isExpanded())
3274 {
3275 B.expand();
3276 }
3277 Data res(0, DataTypes::scalarShape, getFunctionSpace(), true);
3278 try
3279 {
3280 int numpts=getNumDataPoints();
3281 const DataVector& adat=getReady()->getVectorRO();
3282 const DataVector& bdat=B.getReady()->getVectorRO();
3283 DataVector& rdat=res.getReady()->getVectorRW();
3284 const DataTypes::ShapeType& ts=table.getShape();
3285 int twx=ts[0]-1; // table width x
3286 int twy=ts[1]-1; // table width y
3287 bool haserror=false;
3288 int l=0;
3289 #pragma omp parallel for private(l) schedule(static)
3290 for (l=0; l<numpts; ++l)
3291 {
3292 #pragma omp flush(haserror) // In case haserror was in register
3293 if (!haserror)
3294 {
3295 int lerror=0;
3296 try
3297 {
3298 do
3299 {
3300