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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 6 months ago) by jfenwick
File size: 104402 byte(s)
Don't panic.
Updating copyright stamps

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