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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6469 - (show annotations)
Tue Jan 17 07:45:36 2017 UTC (2 years, 3 months ago) by jfenwick
File size: 207132 byte(s)
Add zeroed copy

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