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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2787 - (show annotations)
Fri Nov 27 05:03:09 2009 UTC (10 years, 7 months ago) by jfenwick
File size: 99064 byte(s)
Worked around icc11 problem in table interpolate.
Best guess is that a stackframe was getting mauled somehow.


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