/[escript]/branches/complex/escriptcore/src/Data.cpp
ViewVC logotype

Contents of /branches/complex/escriptcore/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5881 - (show annotations)
Thu Jan 21 04:21:25 2016 UTC (2 years, 5 months ago) by jfenwick
File size: 144331 byte(s)
TensorOps functions are now defined in terms of template parameters rather
than using "double" directly.
This uses the typeinfo and c++11 feature decltype.
This commit does not necessarily fix the functions which TensorOps 
call but it passes unit tests for now.

Also replaced all <double (*)(double) > instantiations with function objects.

This commit does not force c++11 compilation in options files so
if your options are not configured properly you may have compile issues.


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