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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6168 - (show annotations)
Wed Apr 13 03:08:12 2016 UTC (2 years, 6 months ago) by jfenwick
File size: 202553 byte(s)
Making Lazy and the rest of escript use the same operator enumeration


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