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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2721 - (show annotations)
Fri Oct 16 05:40:12 2009 UTC (10 years, 8 months ago) by jfenwick
File size: 96572 byte(s)
minval and maxval are now lazy operations (they weren't before).
Whether or not Lsup, sup and inf resolve their arguments before computing answers is controlled by the escriptParam 'RESOLVE_COLLECTIVE'.
Note: integrate() still forces a resolve.

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