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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6517 - (show annotations)
Mon Mar 6 06:35:31 2017 UTC (2 years, 1 month ago) by jfenwick
File size: 207527 byte(s)
Work towards complex lazy.

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 if (isLazy())
1700 {
1701 throw DataException("Error - cannot integrate for constant lazy data.");
1702 }
1703 if (isComplex()) {
1704 return integrateWorker<cplx_t>();
1705 } else {
1706 return integrateWorker<real_t>();
1707 }
1708 }
1709
1710 bp::object
1711 Data::integrateToTuple()
1712 {
1713 if (isLazy())
1714 {
1715 expand(); // Can't do a non-resolving version of this without changing the domain object
1716 } // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1717 if (isComplex()) {
1718 return integrateWorker<cplx_t>();
1719 } else {
1720 return integrateWorker<real_t>();
1721 }
1722 }
1723
1724 template<typename Scalar>
1725 bp::object
1726 Data::integrateWorker() const
1727 {
1728 DataTypes::ShapeType shape = getDataPointShape();
1729 int dataPointSize = getDataPointSize();
1730
1731 //
1732 // calculate the integral values
1733 vector<Scalar> integrals(dataPointSize);
1734 vector<Scalar> integrals_local(dataPointSize);
1735 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1736 if (dom == 0)
1737 {
1738 throw DataException("Can not integrate over non-continuous domains.");
1739 }
1740 #ifdef ESYS_MPI
1741 dom->setToIntegrals(integrals_local, *this);
1742 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1743 Scalar* tmp = new Scalar[dataPointSize];
1744 Scalar* tmp_local = new Scalar[dataPointSize];
1745 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1746 if (sizeof(Scalar) == sizeof(double)) {
1747 MPI_Allreduce(&tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, getDomain()->getMPIComm());
1748 } else {
1749 MPI_Allreduce(&tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE_COMPLEX, MPI_SUM, getDomain()->getMPIComm());
1750 }
1751 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1752 bp::tuple result = pointToTuple(shape, tmp);
1753 delete[] tmp;
1754 delete[] tmp_local;
1755 #else
1756 dom->setToIntegrals(integrals, *this);
1757 bp::tuple result = pointToTuple(shape, integrals);
1758 #endif
1759
1760 return result;
1761 }
1762
1763 Data
1764 Data::besselFirstKind(int order)
1765 {
1766 THROWONCOMPLEX
1767 return bessel(order,boost::math::cyl_bessel_j);
1768 }
1769
1770 Data
1771 Data::besselSecondKind(int order)
1772 {
1773 THROWONCOMPLEX
1774 return bessel(order,boost::math::cyl_neumann);
1775 }
1776
1777 Data
1778 Data::bessel(int order, real_t (*besselfunc) (int,real_t) )
1779 {
1780 THROWONCOMPLEX
1781 DataTypes::real_t wantreal=0;
1782 if (isEmpty()) // do this before we attempt to interpolate
1783 {
1784 throw DataException("Error - Operations (bessel) not permitted on instances of DataEmpty.");
1785 }
1786 if (isLazy())
1787 {
1788 resolve();
1789 }
1790 // Interpolate if necessary and find an appropriate function space
1791 Data arg_0_Z = Data(*this);
1792
1793 // Get rank and shape of inputs
1794 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
1795 int size0 = arg_0_Z.getDataPointSize();
1796
1797 // Declare output Data object
1798 Data res;
1799
1800 if (arg_0_Z.isConstant()) {
1801 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),false); // DataConstant output
1802 const real_t *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0, wantreal));
1803 real_t *ptr_2 = &(res.getDataAtOffsetRW(0, wantreal));
1804 for (int i = 0; i < size0; ++i) {
1805 ptr_2[i] = besselfunc(order, ptr_0[i]);
1806 }
1807 }
1808 else if (arg_0_Z.isTagged()) {
1809
1810 // Borrow DataTagged input from Data object
1811 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
1812
1813 // Prepare a DataTagged output 2
1814 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(), false); // DataTagged output
1815 res.tag();
1816 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
1817
1818 // Get the pointers to the actual data
1819 const real_t *ptr_0 = &(tmp_0->getDefaultValueRO(0, wantreal));
1820 real_t *ptr_2 = &(tmp_2->getDefaultValueRW(0, wantreal));
1821 // Compute a result for the default
1822 for (int i = 0; i < size0; ++i) {
1823 ptr_2[i] = besselfunc(order, ptr_0[i]);
1824 }
1825 // Compute a result for each tag
1826 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
1827 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
1828 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
1829 tmp_2->addTag(i->first);
1830 const real_t *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
1831 real_t *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
1832 for (int i = 0; i < size0; ++i) {
1833 ptr_2[i] = besselfunc(order, ptr_0[i]);
1834 }
1835
1836 }
1837
1838 }
1839 else if (arg_0_Z.isExpanded()) {
1840
1841 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
1842 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
1843 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
1844
1845 int sampleNo_0,dataPointNo_0;
1846 int numSamples_0 = arg_0_Z.getNumSamples();
1847 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
1848 DataTypes::real_t wantreal=0;
1849 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
1850 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
1851 dataPointNo_0=0;
1852 // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
1853 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
1854 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
1855 const real_t *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0, wantreal));
1856 real_t *ptr_2 = &(res.getDataAtOffsetRW(offset_2, wantreal));
1857 for (int i = 0; i < size0*numDataPointsPerSample_0; ++i) {
1858 ptr_2[i] = besselfunc(order, ptr_0[i]);
1859 }
1860 // }
1861 }
1862 }
1863 else {
1864 throw DataException("Error - Bessel function: unknown combination of inputs");
1865 }
1866
1867 return res;
1868 }
1869
1870 Data
1871 Data::conjugate() const
1872 {
1873 if (isLazy())
1874 {
1875 Data temp(*this);
1876 temp.resolve();
1877 return temp.conjugate();
1878 }
1879 if (isComplex())
1880 {
1881 return C_TensorUnaryOperation(*this, escript::ES_optype::CONJ);
1882 }
1883 else
1884 {
1885 return copySelf();
1886 }
1887 }
1888
1889
1890 Data
1891 Data::real() const
1892 {
1893 if (isLazy())
1894 {
1895 Data temp(*this);
1896 temp.resolve();
1897 return temp.real();
1898 }
1899 if (isComplex())
1900 {
1901 return C_TensorUnaryOperation(*this, escript::ES_optype::REAL);
1902 }
1903 else
1904 {
1905 return copySelf();
1906 }
1907 }
1908
1909 Data
1910 Data::imag() const
1911 {
1912 if (isLazy())
1913 {
1914 Data temp(*this);
1915 temp.resolve();
1916 return temp.imag();
1917 }
1918 if (isComplex())
1919 {
1920 return C_TensorUnaryOperation(*this, escript::ES_optype::IMAG);
1921 }
1922 else
1923 {
1924 return copySelf()*Data(0, m_data->getShape(), getFunctionSpace(),false); // return an object with same tags etc but all values 0
1925 // This is not efficient, but why are you taking imag of R anyway?
1926 }
1927 }
1928
1929 Data
1930 Data::phase() const
1931 {
1932 if (isLazy())
1933 {
1934 Data temp(*this);
1935 temp.resolve();
1936 return temp.phase();
1937 }
1938 if (isComplex())
1939 {
1940 return C_TensorUnaryOperation(*this, escript::ES_optype::PHS);
1941 }
1942 else
1943 {
1944 return copySelf()*Data(0, m_data->getShape(), getFunctionSpace(),false); // return an object with same tags etc but all values 0
1945 // This is not efficient, but why are you taking imag of R anyway?
1946 }
1947 }
1948
1949
1950
1951
1952
1953 Data
1954 Data::sin() const
1955 {
1956 MAKELAZYOP(SIN);
1957 return C_TensorUnaryOperation(*this, escript::ES_optype::SIN);
1958 }
1959
1960 Data
1961 Data::cos() const
1962 {
1963 MAKELAZYOP(COS);
1964 return C_TensorUnaryOperation(*this, escript::ES_optype::COS);
1965 }
1966
1967 Data
1968 Data::tan() const
1969 {
1970 MAKELAZYOP(TAN);
1971 return C_TensorUnaryOperation(*this, escript::ES_optype::TAN);
1972 }
1973
1974 Data
1975 Data::asin() const
1976 {
1977 MAKELAZYOP(ASIN);
1978 return C_TensorUnaryOperation(*this, escript::ES_optype::ASIN);
1979 }
1980
1981 Data
1982 Data::acos() const
1983 {
1984 MAKELAZYOP(ACOS);
1985 return C_TensorUnaryOperation(*this, escript::ES_optype::ACOS);
1986 }
1987
1988
1989 Data
1990 Data::atan() const
1991 {
1992 MAKELAZYOP(ATAN);
1993 return C_TensorUnaryOperation(*this, escript::ES_optype::ATAN);
1994 }
1995
1996 Data
1997 Data::sinh() const
1998 {
1999 MAKELAZYOP(SINH);
2000 return C_TensorUnaryOperation(*this, escript::ES_optype::SINH);
2001 }
2002
2003 Data
2004 Data::cosh() const
2005 {
2006 MAKELAZYOP(COSH);
2007 return C_TensorUnaryOperation(*this, escript::ES_optype::COSH);
2008 }
2009
2010 Data
2011 Data::tanh() const
2012 {
2013 MAKELAZYOP(TANH);
2014 return C_TensorUnaryOperation(*this, escript::ES_optype::TANH);
2015 }
2016
2017
2018 Data
2019 Data::erf() const
2020 {
2021 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
2022 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
2023 #else
2024 MAKELAZYOP(ERF);
2025 return C_TensorUnaryOperation(*this, escript::ES_optype::ERF);
2026 #endif
2027 }
2028
2029 Data
2030 Data::asinh() const
2031 {
2032 MAKELAZYOP(ASINH);
2033 return C_TensorUnaryOperation(*this, escript::ES_optype::ASINH);
2034 }
2035
2036 Data
2037 Data::acosh() const
2038 {
2039 MAKELAZYOP(ACOSH);
2040 return C_TensorUnaryOperation(*this, escript::ES_optype::ACOSH);
2041 }
2042
2043 Data
2044 Data::atanh() const
2045 {
2046 MAKELAZYOP(ATANH);
2047 return C_TensorUnaryOperation(*this, escript::ES_optype::ATANH);
2048 }
2049
2050 Data
2051 Data::log10() const
2052 {
2053 MAKELAZYOP(LOG10);
2054 return C_TensorUnaryOperation(*this, escript::ES_optype::LOG10);
2055 }
2056
2057 Data
2058 Data::log() const
2059 {
2060 MAKELAZYOP(LOG);
2061 return C_TensorUnaryOperation(*this, escript::ES_optype::LOG);
2062 }
2063
2064 Data
2065 Data::sign() const
2066 {
2067 MAKELAZYOP(SIGN);
2068 return C_TensorUnaryOperation(*this, escript::ES_optype::SIGN);
2069 }
2070
2071 Data
2072 Data::abs() const
2073 {
2074 MAKELAZYOP(ABS);
2075 return C_TensorUnaryOperation(*this, escript::ES_optype::ABS);
2076 }
2077
2078 Data
2079 Data::neg() const
2080 {
2081 MAKELAZYOP(NEG);
2082 return C_TensorUnaryOperation(*this, escript::ES_optype::NEG);
2083 }
2084
2085 Data
2086 Data::pos() const
2087 {
2088 THROWONCOMPLEX
2089 // not doing lazy check here is deliberate.
2090 // since a deep copy of lazy data should be cheap, I'll just let it happen now
2091 Data result;
2092 // perform a deep copy
2093 result.copy(*this);
2094 return result;
2095 }
2096
2097 Data
2098 Data::exp() const
2099 {
2100 MAKELAZYOP(EXP);
2101 return C_TensorUnaryOperation(*this, escript::ES_optype::EXP);
2102 }
2103
2104 Data
2105 Data::sqrt() const
2106 {
2107 MAKELAZYOP(SQRT);
2108 return C_TensorUnaryOperation(*this, escript::ES_optype::SQRT);
2109 }
2110
2111 real_t
2112 Data::Lsup_const() const
2113 {
2114 if (isLazy())
2115 {
2116 throw DataException("Error - cannot compute Lsup for constant lazy data.");
2117 }
2118 return LsupWorker();
2119 }
2120
2121 real_t
2122 Data::Lsup()
2123 {
2124 if (isLazy())
2125 {
2126 if (!actsExpanded() || CHECK_DO_CRES)
2127 {
2128 resolve();
2129 }
2130 else
2131 {
2132 #ifdef ESYS_MPI
2133 if (isComplex())
2134 {
2135 return lazyAlgWorker<AbsMax<cplx_t> >(0,MPI_MAX);
2136 }
2137 else
2138 {
2139 return lazyAlgWorker<AbsMax<cplx_t> >(0,MPI_MAX);
2140 }
2141 #else
2142 if (isComplex())
2143 {
2144 return lazyAlgWorker<AbsMax<real_t> >(0);
2145 }
2146 else
2147 {
2148 return lazyAlgWorker<AbsMax<real_t> >(0);
2149 }
2150 #endif
2151 }
2152 }
2153 return LsupWorker();
2154 }
2155
2156 real_t
2157 Data::sup_const() const
2158 {
2159 if (isComplex())
2160 {
2161 throw DataException("Error Cannot compute sup() for complex data.");
2162 }
2163 if (isLazy())
2164 {
2165 throw DataException("Error - cannot compute sup for constant lazy data.");
2166 }
2167 return supWorker();
2168 }
2169
2170 real_t
2171 Data::sup()
2172 {
2173 if (isComplex())
2174 {
2175 throw DataException("Error Cannot compute sup() for complex data.");
2176 }
2177 if (isLazy())
2178 {
2179 if (!actsExpanded() || CHECK_DO_CRES)
2180 {
2181 resolve();
2182 }
2183 else
2184 {
2185 #ifdef ESYS_MPI
2186 return lazyAlgWorker<FMax>(numeric_limits<real_t>::max()*-1, MPI_MAX);
2187 #else
2188 return lazyAlgWorker<FMax>(numeric_limits<real_t>::max()*-1);
2189 #endif
2190 }
2191 }
2192 return supWorker();
2193 }
2194
2195 real_t
2196 Data::inf_const() const
2197 {
2198 if (isComplex())
2199 {
2200 throw DataException("Error Cannot compute inf() for complex data.");
2201 }
2202 if (isLazy())
2203 {
2204 throw DataException("Error - cannot compute inf for constant lazy data.");
2205 }
2206 return infWorker();
2207 }
2208
2209 real_t
2210 Data::inf()
2211 {
2212 if (isComplex())
2213 {
2214 throw DataException("Error Cannot compute inf() for complex data.");
2215 }
2216 if (isLazy())
2217 {
2218 if (!actsExpanded() || CHECK_DO_CRES)
2219 {
2220 resolve();
2221 }
2222 else
2223 {
2224 #ifdef ESYS_MPI
2225 return lazyAlgWorker<FMin>(numeric_limits<real_t>::max(), MPI_MIN);
2226 #else
2227 return lazyAlgWorker<FMin>(numeric_limits<real_t>::max());
2228 #endif
2229 }
2230 }
2231 return infWorker();
2232 }
2233
2234 template <class BinaryOp>
2235 real_t
2236 #ifdef ESYS_MPI
2237 Data::lazyAlgWorker(real_t init, MPI_Op mpiop_type)
2238 #else
2239 Data::lazyAlgWorker(real_t init)
2240 #endif
2241 {
2242 if (!isLazy() || !m_data->actsExpanded())
2243 {
2244 throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
2245 }
2246 DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
2247 ESYS_ASSERT(dl!=NULL, "Programming error - casting to DataLazy.");
2248 real_t val=init;
2249 int i=0;
2250 const size_t numsamples=getNumSamples();
2251 const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
2252 BinaryOp operation;
2253 real_t localValue=0, globalValue;
2254 #pragma omp parallel private(i)
2255 {
2256 real_t localtot=init;
2257 #pragma omp for schedule(static) private(i)
2258 for (i=0;i<numsamples;++i)
2259 {
2260 size_t roffset=0;
2261 const DataTypes::RealVectorType* v=dl->resolveSample(i, roffset);
2262 // Now we have the sample, run operation on all points
2263 for (size_t j=0;j<samplesize;++j)
2264 {
2265 localtot=operation(localtot,(*v)[j+roffset]);
2266 }
2267 if (escript::vectorHasNaN(*v,roffset, samplesize))
2268 {
2269 #pragma omp critical
2270 {
2271 localValue=1.0;
2272 }
2273 }
2274 }
2275 #pragma omp critical
2276 val=operation(val,localtot);
2277 }
2278 #ifdef ESYS_MPI
2279 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2280 #else
2281 globalValue=localValue;
2282 #endif
2283 if (globalValue!=0)
2284 {
2285 return makeNaN();
2286 }
2287 #ifdef ESYS_MPI
2288 MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, getDomain()->getMPIComm() );
2289 return globalValue;
2290 #else
2291 return val;
2292 #endif
2293 }
2294
2295 bool
2296 Data::hasNaN()
2297 {
2298 if (isLazy())
2299 {
2300 resolve();
2301 }
2302 bool haveNaN=getReady()->hasNaN();
2303 return haveNaN;
2304 }
2305
2306 void
2307 Data::replaceNaN(real_t value)
2308 {
2309 if (isLazy())
2310 {
2311 resolve();
2312 }
2313 getReady()->replaceNaN(value);
2314 }
2315
2316 void
2317 Data::replaceNaN(cplx_t value)
2318 {
2319 if (isLazy())
2320 {
2321 resolve();
2322 }
2323 getReady()->replaceNaN(value);
2324 }
2325
2326 void
2327 Data::replaceNaNPython(boost::python::object obj)
2328 {
2329 boost::python::extract<DataTypes::real_t> exr(obj);
2330 if (exr.check())
2331 {
2332 replaceNaN(exr());
2333 }
2334 else
2335 {
2336 replaceNaN(boost::python::extract<DataTypes::cplx_t>(obj)());
2337 }
2338 }
2339
2340 // Do not call this on Lazy Data use the proper entry point
2341 real_t
2342 Data::LsupWorker() const
2343 {
2344 bool haveNaN=getReady()->hasNaN();
2345
2346
2347 #ifdef ESYS_MPI
2348 int nanchecker=0;
2349 if (haveNaN)
2350 {
2351 nanchecker=1.0;
2352 }
2353 int globalnan;
2354 MPI_Allreduce( &nanchecker, &globalnan, 1, MPI_INT, MPI_MAX, getDomain()->getMPIComm() );
2355 if (globalnan!=0)
2356 {
2357 return makeNaN();
2358 }
2359 #else
2360 if (haveNaN)
2361 {
2362 return makeNaN();
2363 }
2364 #endif
2365
2366 //
2367 // set the initial absolute maximum value to zero
2368 if (isComplex())
2369 {
2370 AbsMax<cplx_t> abs_max_func;
2371 real_t localValue=0;
2372 localValue = reduction(abs_max_func,0);
2373
2374 #ifdef ESYS_MPI
2375 real_t globalValue=0;
2376 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2377 return globalValue;
2378 #else
2379 return localValue;
2380 #endif
2381 }
2382 else
2383 {
2384 AbsMax<real_t> abs_max_func;
2385 real_t localValue=0;
2386 localValue = reduction(abs_max_func,0);
2387
2388 #ifdef ESYS_MPI
2389 real_t globalValue=0;
2390 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2391 return globalValue;
2392 #else
2393 return localValue;
2394 #endif
2395 }
2396 }
2397
2398 real_t
2399 Data::supWorker() const
2400 {
2401 bool haveNaN=getReady()->hasNaN();
2402 real_t localValue=0;
2403
2404 #ifdef ESYS_MPI
2405 if (haveNaN)
2406 {
2407 localValue=1.0;
2408 }
2409 real_t globalValue;
2410 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2411 if (globalValue!=0)
2412 {
2413 return makeNaN();
2414 }
2415 #else
2416 if (haveNaN)
2417 {
2418 return makeNaN();
2419 }
2420 #endif
2421
2422 //
2423 // set the initial maximum value to min possible real_t
2424 FMax fmax_func;
2425 if (hasNoSamples())
2426 {
2427 localValue = numeric_limits<real_t>::infinity()*-1;
2428 }
2429 else
2430 {
2431 localValue = reduction(fmax_func,numeric_limits<real_t>::infinity()*-1);
2432 }
2433 #ifdef ESYS_MPI
2434 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2435 return globalValue;
2436 #else
2437 return localValue;
2438 #endif
2439 }
2440
2441 real_t
2442 Data::infWorker() const
2443 {
2444 bool haveNaN=getReady()->hasNaN();
2445 real_t localValue=0;
2446
2447 #ifdef ESYS_MPI
2448 if (haveNaN)
2449 {
2450 localValue=1.0;
2451 }
2452 real_t globalValue;
2453 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, getDomain()->getMPIComm() );
2454 if (globalValue!=0)
2455 {
2456 return makeNaN();
2457 }
2458 #else
2459 if (haveNaN)
2460 {
2461 return makeNaN();
2462 }
2463 #endif
2464 //
2465 // set the initial minimum value to max possible real_t
2466 FMin fmin_func;
2467 if (hasNoSamples())
2468 {
2469 localValue = numeric_limits<real_t>::infinity();
2470 }
2471 else
2472 {
2473 localValue = reduction(fmin_func,numeric_limits<real_t>::infinity());
2474 }
2475 #ifdef ESYS_MPI
2476 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, getDomain()->getMPIComm() );
2477 return globalValue;
2478 #else
2479 return localValue;
2480 #endif
2481 }
2482
2483 /* global reduction */
2484
2485
2486 inline Data
2487 Data::minval_nonlazy() const
2488 {
2489 THROWONCOMPLEX
2490 //
2491 // set the initial minimum value to max possible real_t
2492 FMin fmin_func;
2493 return dp_algorithm(fmin_func,numeric_limits<real_t>::max());
2494 }
2495
2496
2497 inline Data
2498 Data::maxval_nonlazy() const
2499 {
2500 THROWONCOMPLEX
2501 //
2502 // set the initial maximum value to min possible real_t
2503 FMax fmax_func;
2504 return dp_algorithm(fmax_func,numeric_limits<real_t>::max()*-1);
2505 }
2506
2507
2508 Data
2509 Data::maxval() const
2510 {
2511 THROWONCOMPLEX
2512 MAKELAZYOP(MAXVAL);
2513 return maxval_nonlazy();
2514 }
2515
2516
2517 Data
2518 Data::minval() const
2519 {
2520 THROWONCOMPLEX
2521 MAKELAZYOP(MINVAL);
2522 return minval_nonlazy();
2523 }
2524
2525
2526 Data
2527 Data::swapaxes(const int axis0, const int axis1) const
2528 {
2529 int axis0_tmp,axis1_tmp;
2530 DataTypes::ShapeType s=getDataPointShape();
2531 DataTypes::ShapeType ev_shape;
2532 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
2533 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
2534 int rank=getDataPointRank();
2535 if (rank<2) {
2536 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
2537 }
2538 if (axis0<0 || axis0>rank-1) {
2539 stringstream e;
2540 e << "Error - Data::swapaxes: axis0 must be between 0 and rank-1=" << (rank-1);
2541 throw DataException(e.str());
2542 }
2543 if (axis1<0 || axis1>rank-1) {
2544 stringstream e;
2545 e << "Error - Data::swapaxes: axis1 must be between 0 and rank-1=" << (rank-1);
2546 throw DataException(e.str());
2547 }
2548 if (axis0 == axis1) {
2549 throw DataException("Error - Data::swapaxes: axis indices must be different.");
2550 }
2551 MAKELAZYOP2(SWAP,axis0,axis1);
2552 if (axis0 > axis1)
2553 {
2554 axis0_tmp=axis1;
2555 axis1_tmp=axis0;
2556 }
2557 else
2558 {
2559 axis0_tmp=axis0;
2560 axis1_tmp=axis1;
2561 }
2562 for (int i=0; i<rank; i++)
2563 {
2564 if (i == axis0_tmp)
2565 {
2566 ev_shape.push_back(s[axis1_tmp]);
2567 }
2568 else if (i == axis1_tmp)
2569 {
2570 ev_shape.push_back(s[axis0_tmp]);
2571 }
2572 else
2573 {
2574 ev_shape.push_back(s[i]);
2575 }
2576 }
2577 Data ev(0.,ev_shape,getFunctionSpace(), false);
2578 ev.typeMatchRight(*this);
2579 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
2580 return ev;
2581 }
2582
2583 Data
2584 Data::symmetric() const
2585 {
2586 // check input
2587 DataTypes::ShapeType s=getDataPointShape();
2588 if (getDataPointRank()==2) {
2589 if(s[0] != s[1])
2590 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
2591 }
2592 else if (getDataPointRank()==4) {
2593 if(!(s[0] == s[2] && s[1] == s[3]))
2594 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2595 }
2596 else {
2597 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
2598 }
2599 MAKELAZYOP(SYM);
2600 Data ev(0.,getDataPointShape(),getFunctionSpace(), false);
2601 ev.typeMatchRight(*this);
2602 m_data->symmetric(ev.m_data.get());
2603 return ev;
2604 }
2605
2606 Data
2607 Data::antisymmetric() const
2608 {
2609 // check input
2610 DataTypes::ShapeType s=getDataPointShape();
2611 if (getDataPointRank()==2) {
2612 if(s[0] != s[1])
2613 throw DataException("Error - Data::antisymmetric can only be calculated for rank 2 object with equal first and second dimension.");
2614 MAKELAZYOP(NSYM);
2615 DataTypes::ShapeType ev_shape;
2616 ev_shape.push_back(s[0]);
2617 ev_shape.push_back(s[1]);
2618 Data ev(0.,ev_shape,getFunctionSpace(),false);
2619 ev.typeMatchRight(*this);
2620 m_data->antisymmetric(ev.m_data.get());
2621 return ev;
2622 }
2623 else if (getDataPointRank()==4) {
2624 if(!(s[0] == s[2] && s[1] == s[3]))
2625 throw DataException("Error - Data::antisymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2626 MAKELAZYOP(NSYM);
2627 DataTypes::ShapeType ev_shape;
2628 ev_shape.push_back(s[0]);
2629 ev_shape.push_back(s[1]);
2630 ev_shape.push_back(s[2]);
2631 ev_shape.push_back(s[3]);
2632 Data ev(0.,ev_shape,getFunctionSpace(),false);
2633 ev.typeMatchRight(*this);
2634 m_data->antisymmetric(ev.m_data.get());
2635 return ev;
2636 }
2637 else {
2638 throw DataException("Error - Data::antisymmetric can only be calculated for rank 2 or 4 object.");
2639 }
2640 }
2641
2642 Data
2643 Data::hermitian() const
2644 {
2645 if (!isComplex())
2646 {
2647 return symmetric();
2648 }
2649 // check input
2650 DataTypes::ShapeType s=getDataPointShape();
2651 if (getDataPointRank()==2) {
2652 if(s[0] != s[1])
2653 throw DataException("Error - Data::hermitian can only be calculated for rank 2 object with equal first and second dimension.");
2654 }
2655 else if (getDataPointRank()==4) {
2656 if(!(s[0] == s[2] && s[1] == s[3]))
2657 throw DataException("Error - Data::hermitian can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2658 }
2659 else {
2660 throw DataException("Error - Data::hermitian can only be calculated for rank 2 or 4 object.");
2661 }
2662 MAKELAZYOP(HER);
2663 Data ev(0.,getDataPointShape(),getFunctionSpace(), false);
2664 ev.typeMatchRight(*this);
2665 m_data->hermitian(ev.m_data.get());
2666 return ev;
2667 }
2668
2669 Data
2670 Data::antihermitian() const
2671 {
2672 if (!isComplex())
2673 {
2674 return antisymmetric();
2675 }
2676 // check input
2677 DataTypes::ShapeType s=getDataPointShape();
2678 if (getDataPointRank()==2) {
2679 if(s[0] != s[1])
2680 throw DataException("Error - Data::antisymmetric can only be calculated for rank 2 object with equal first and second dimension.");
2681 MAKELAZYOP(NHER);
2682 DataTypes::ShapeType ev_shape;
2683 ev_shape.push_back(s[0]);
2684 ev_shape.push_back(s[1]);
2685 Data ev(0.,ev_shape,getFunctionSpace(), false);
2686 ev.typeMatchRight(*this);
2687 m_data->antihermitian(ev.m_data.get());
2688 return ev;
2689 }
2690 else if (getDataPointRank()==4) {
2691 if(!(s[0] == s[2] && s[1] == s[3]))
2692 throw DataException("Error - Data::antisymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
2693 MAKELAZYOP(NHER);
2694 DataTypes::ShapeType ev_shape;
2695 ev_shape.push_back(s[0]);
2696 ev_shape.push_back(s[1]);
2697 ev_shape.push_back(s[2]);
2698 ev_shape.push_back(s[3]);
2699 Data ev(0.,ev_shape,getFunctionSpace(), false);
2700 ev.typeMatchRight(*this);
2701 m_data->antihermitian(ev.m_data.get());
2702 return ev;
2703 }
2704 else {
2705 throw DataException("Error - Data::antisymmetric can only be calculated for rank 2 or 4 object.");
2706 }
2707 }
2708
2709 Data
2710 Data::trace(int axis_offset) const
2711 {
2712 MAKELAZYOPOFF(TRACE,axis_offset);
2713 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
2714 {
2715 throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
2716 }
2717 DataTypes::ShapeType s=getDataPointShape();
2718 if (getDataPointRank()==2) {
2719 DataTypes::ShapeType ev_shape;
2720 Data ev(0.,ev_shape,getFunctionSpace(),false);
2721 ev.typeMatchRight(*this);
2722 m_data->trace(ev.m_data.get(), axis_offset);
2723 return ev;
2724 }
2725 if (getDataPointRank()==3) {
2726 DataTypes::ShapeType ev_shape;
2727 if (axis_offset==0) {
2728 int s2=s[2];
2729 ev_shape.push_back(s2);
2730 }
2731 else if (axis_offset==1) {
2732 int s0=s[0];
2733 ev_shape.push_back(s0);
2734 }
2735 Data ev(0.,ev_shape,getFunctionSpace(),false);
2736 ev.typeMatchRight(*this);
2737 m_data->trace(ev.m_data.get(), axis_offset);
2738 return ev;
2739 }
2740 if (getDataPointRank()==4) {
2741 DataTypes::ShapeType ev_shape;
2742 if (axis_offset==0) {
2743 ev_shape.push_back(s[2]);
2744 ev_shape.push_back(s[3]);
2745 }
2746 else if (axis_offset==1) {
2747 ev_shape.push_back(s[0]);
2748 ev_shape.push_back(s[3]);
2749 }
2750 else if (axis_offset==2) {
2751 ev_shape.push_back(s[0]);
2752 ev_shape.push_back(s[1]);
2753 }
2754 Data ev(0.,ev_shape,getFunctionSpace(),false);
2755 ev.typeMatchRight(*this);
2756 m_data->trace(ev.m_data.get(), axis_offset);
2757 return ev;
2758 }
2759 else {
2760 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
2761 }
2762 }
2763
2764 Data
2765 Data::transpose(int axis_offset) const
2766 {
2767 MAKELAZYOPOFF(TRANS,axis_offset);
2768 DataTypes::ShapeType s=getDataPointShape();
2769 DataTypes::ShapeType ev_shape;
2770 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
2771 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
2772 int rank=getDataPointRank();
2773 if (axis_offset<0 || axis_offset>rank) {
2774 stringstream e;
2775 e << "Error - Data::transpose must have 0 <= axis_offset <= rank=" << rank;
2776 throw DataException(e.str());
2777 }
2778 for (int i=0; i<rank; i++) {
2779 int index = (axis_offset+i)%rank;
2780 ev_shape.push_back(s[index]); // Append to new shape
2781 }
2782 Data ev(0.,ev_shape,getFunctionSpace(), false);
2783 ev.typeMatchRight(*this);
2784 m_data->transpose(ev.m_data.get(), axis_offset);
2785 return ev;
2786 }
2787
2788 Data
2789 Data::eigenvalues() const
2790 {
2791 if (isLazy())
2792 {
2793 Data temp(*this); // to get around the fact that you can't resolve a const Data
2794 temp.resolve();
2795 return temp.eigenvalues();
2796 }
2797 // check input
2798 DataTypes::ShapeType s=getDataPointShape();
2799 if (getDataPointRank()!=2)
2800 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
2801 if(s[0] != s[1])
2802 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
2803 if (isComplex() && (s[0]>2))
2804 {
2805 throw DataException("Error - Data::eigenvalues not supported for complex 3x3.");
2806 }
2807 // create return
2808 DataTypes::ShapeType ev_shape(1,s[0]);
2809 Data ev(0.,ev_shape,getFunctionSpace(),false);
2810 ev.typeMatchRight(*this);
2811 m_data->eigenvalues(ev.m_data.get());
2812 return ev;
2813 }
2814
2815 const bp::tuple
2816 Data::eigenvalues_and_eigenvectors(const real_t tol) const
2817 {
2818 THROWONCOMPLEX
2819 if (isLazy())
2820 {
2821 Data temp(*this); // to get around the fact that you can't resolve a const Data
2822 temp.resolve();
2823 return temp.eigenvalues_and_eigenvectors(tol);
2824 }
2825 DataTypes::ShapeType s=getDataPointShape();
2826 if (getDataPointRank()!=2)
2827 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
2828 if(s[0] != s[1])
2829 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
2830 // create return
2831 DataTypes::ShapeType ev_shape(1,s[0]);
2832 Data ev(0.,ev_shape,getFunctionSpace(), false);
2833 ev.typeMatchRight(*this);
2834 DataTypes::ShapeType V_shape(2,s[0]);
2835 Data V(0.,V_shape,getFunctionSpace(), false);
2836 V.typeMatchRight(*this);
2837 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
2838 return bp::make_tuple(bp::object(ev),bp::object(V));
2839 }
2840
2841 const bp::tuple
2842 Data::minGlobalDataPoint() const
2843 {
2844 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
2845 // abort (for unknown reasons) if there are openmp directives with it in the
2846 // surrounding function
2847
2848 THROWONCOMPLEX
2849 int DataPointNo;
2850 int ProcNo;
2851 calc_minGlobalDataPoint(ProcNo,DataPointNo);
2852 return bp::make_tuple(ProcNo,DataPointNo);
2853 }
2854
2855 void
2856 Data::calc_minGlobalDataPoint(int& ProcNo,
2857 int& DataPointNo) const
2858 {
2859 THROWONCOMPLEX
2860 if (isLazy())
2861 {
2862 Data temp(*this); // to get around the fact that you can't resolve a const Data
2863 temp.resolve();
2864 return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2865 }
2866 int i,j;
2867 int lowi=0,lowj=0;
2868 real_t min=numeric_limits<real_t>::max();
2869
2870 Data temp=minval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2871
2872 int numSamples=temp.getNumSamples();
2873 int numDPPSample=temp.getNumDataPointsPerSample();
2874
2875 real_t local_val, local_min;
2876 #ifdef ESYS_MPI
2877 real_t next[2];
2878 #endif
2879 int local_lowi=0,local_lowj=0;
2880
2881 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(local_val,local_min)
2882 {
2883 local_min=min;
2884 DataTypes::real_t wantreal=0;
2885 #pragma omp for private(i,j) schedule(static)
2886 for (i=0; i<numSamples; i++) {
2887 for (j=0; j<numDPPSample; j++) {
2888 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j), wantreal);
2889 if (local_val<local_min) {
2890 local_min=local_val;
2891 local_lowi=i;
2892 local_lowj=j;
2893 }
2894 }
2895 }
2896 #pragma omp critical
2897 if (local_min<min) { // If we found a smaller value than our sentinel
2898 min=local_min;
2899 lowi=local_lowi;
2900 lowj=local_lowj;
2901 }
2902 } // parallel section
2903
2904 #ifdef ESYS_MPI
2905 // determine the processor on which the minimum occurs
2906 next[0] = min;
2907 next[1] = numSamples;
2908 int lowProc = 0;
2909 real_t *globalMins = new real_t[get_MPISize()*2+1];
2910 /*int error =*/ MPI_Gather (next, 2, MPI_DOUBLE, globalMins, 2, MPI_DOUBLE, 0, get_MPIComm() );
2911
2912 if ( get_MPIRank()==0 ) {
2913 for (lowProc=0; lowProc<get_MPISize(); lowProc++)
2914 if (globalMins[lowProc*2+1] > 0) break;
2915 min = globalMins[lowProc*2];
2916 for( i=lowProc+1; i<get_MPISize(); i++ )
2917 if( globalMins[i*2+1]>0 && min>globalMins[i*2] ) {
2918 lowProc = i;
2919 min = globalMins[i*2];
2920 }
2921 }
2922 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2923 DataPointNo = lowj + lowi * numDPPSample;
2924 MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2925 delete [] globalMins;
2926 ProcNo = lowProc;
2927 #else
2928 ProcNo = 0;
2929 DataPointNo = lowj + lowi * numDPPSample;
2930 #endif
2931 }
2932
2933
2934 const bp::tuple
2935 Data::maxGlobalDataPoint() const
2936 {
2937 THROWONCOMPLEX
2938 int DataPointNo;
2939 int ProcNo;
2940 calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2941 return bp::make_tuple(ProcNo,DataPointNo);
2942 }
2943
2944 void
2945 Data::calc_maxGlobalDataPoint(int& ProcNo,
2946 int& DataPointNo) const
2947 {
2948 if (isLazy())
2949 {
2950 Data temp(*this); // to get around the fact that you can't resolve a const Data
2951 temp.resolve();
2952 return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2953 }
2954 THROWONCOMPLEX
2955 int i,j;
2956 int highi=0,highj=0;
2957 //-------------
2958 real_t max= -numeric_limits<real_t>::max();
2959
2960 Data temp=maxval_nonlazy(); // need to do this to prevent autolazy from reintroducing laziness
2961
2962 int numSamples=temp.getNumSamples();
2963 int numDPPSample=temp.getNumDataPointsPerSample();
2964
2965 real_t local_val, local_max;
2966 #ifdef ESYS_MPI
2967 real_t next[2];
2968 #endif
2969 int local_highi=0,local_highj=0;
2970
2971 #pragma omp parallel firstprivate(local_highi,local_highj) private(local_val,local_max)
2972 {
2973 local_max=max;
2974 DataTypes::real_t wantreal=0;
2975 #pragma omp for private(i,j) schedule(static)
2976 for (i=0; i<numSamples; i++) {
2977 for (j=0; j<numDPPSample; j++) {
2978 local_val=temp.getDataAtOffsetRO(temp.getDataOffset(i,j), wantreal);
2979 if (local_val>local_max) {
2980 local_max=local_val;
2981 local_highi=i;
2982 local_highj=j;
2983 }
2984 }
2985 }
2986 #pragma omp critical
2987 if (local_max>max) { // If we found a larger value than our sentinel
2988 max=local_max;
2989 highi=local_highi;
2990 highj=local_highj;
2991 }
2992 }
2993 #ifdef ESYS_MPI
2994 // determine the processor on which the maximum occurs
2995 next[0] = max;
2996 next[1] = numSamples;
2997 int highProc = 0;
2998 real_t *globalMaxs = new real_t[get_MPISize()*2+1];
2999 /*int error =*/ MPI_Gather ( next, 2, MPI_DOUBLE, globalMaxs, 2, MPI_DOUBLE, 0, get_MPIComm() );
3000 if( get_MPIRank()==0 ){
3001 for (highProc=0; highProc<get_MPISize(); highProc++)
3002 if (globalMaxs[highProc*2+1] > 0) break;
3003 max = globalMaxs[highProc*2];
3004 for( i=highProc+1; i<get_MPISize(); i++ )
3005 {
3006 if( globalMaxs[i*2+1]>0 && max<globalMaxs[i*2] )
3007 {
3008 highProc = i;
3009 max = globalMaxs[i*2];
3010 }
3011 }
3012 }
3013 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
3014 DataPointNo = highj + highi * numDPPSample;
3015 MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
3016
3017 delete [] globalMaxs;
3018 ProcNo = highProc;
3019 #else
3020 ProcNo = 0;
3021 DataPointNo = highj + highi * numDPPSample;
3022 #endif
3023 }
3024
3025 Data&
3026 Data::operator+=(const Data& right)
3027 {
3028 if (isProtected()) {
3029 throw DataException("Error - attempt to update protected Data object.");
3030 }
3031 MAKELAZYBINSELF(right,ADD); // for lazy + is equivalent to +=
3032 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
3033 if (!isComplex() && right.isComplex())
3034 {
3035 complicate();
3036 }
3037 TensorSelfUpdateBinaryOperation(right, escript::ES_optype::ADD);
3038 return (*this);
3039 }
3040
3041 Data&
3042 Data::operator+=(const boost::python::object& right)
3043 {
3044 if (isProtected()) {
3045 throw DataException("Error - attempt to update protected Data object.");
3046 }
3047 Data tmp(right,getFunctionSpace(),false);
3048 (*this)+=tmp;
3049 return *this;
3050 }
3051
3052 // Hmmm, operator= makes a deep copy but the copy constructor does not?
3053 Data&
3054 Data::operator=(const Data& other)
3055 {
3056 m_protected=false; // since any changes should be caught by exclusiveWrite()
3057 set_m_data(other.m_data);
3058 return (*this);
3059 }
3060
3061 Data&
3062 Data::operator-=(const Data& right)
3063 {
3064 if (isProtected()) {
3065 throw DataException("Error - attempt to update protected Data object.");
3066 }
3067 MAKELAZYBINSELF(right,SUB);
3068 exclusiveWrite();
3069 if (!isComplex() && right.isComplex())
3070 {
3071 complicate();
3072 }
3073 TensorSelfUpdateBinaryOperation(right, escript::ES_optype::SUB);
3074 return (*this);
3075 }
3076
3077 Data&
3078 Data::operator-=(const boost::python::object& right)
3079 {
3080 if (isProtected()) {
3081 throw DataException("Error - attempt to update protected Data object.");
3082 }
3083 Data tmp(right,getFunctionSpace(),false);
3084 (*this)-=tmp;
3085 return (*this);
3086 }
3087
3088 Data&
3089 Data::operator*=(const Data& right)
3090 {
3091 if (isProtected()) {
3092 throw DataException("Error - attempt to update protected Data object.");
3093 }
3094 MAKELAZYBINSELF(right,MUL);
3095 exclusiveWrite();
3096 if (!isComplex() && right.isComplex())
3097 {
3098 complicate();
3099 }
3100 TensorSelfUpdateBinaryOperation(right, escript::ES_optype::MUL);
3101 return (*this);
3102 }
3103
3104 Data&
3105 Data::operator*=(const boost::python::object& right)
3106 {
3107 if (isProtected()) {
3108 throw DataException("Error - attempt to update protected Data object.");
3109 }
3110 Data tmp(right,getFunctionSpace(),false);
3111 (*this)*=tmp;
3112 return (*this);
3113 }
3114
3115 Data&
3116 Data::operator/=(const Data& right)
3117 {
3118 if (isProtected()) {
3119 throw DataException("Error - attempt to update protected Data object.");
3120 }
3121 MAKELAZYBINSELF(right,DIV);
3122 exclusiveWrite();
3123 if (!isComplex() && right.isComplex())
3124 {
3125 complicate();
3126 }
3127 TensorSelfUpdateBinaryOperation(right, escript::ES_optype::DIV);
3128 return (*this);
3129 }
3130
3131 Data&
3132 Data::operator/=(const boost::python::object& right)
3133 {
3134 if (isProtected()) {
3135 throw DataException("Error - attempt to update protected Data object.");
3136 }
3137 Data tmp(right,getFunctionSpace(),false);
3138 (*this)/=tmp;
3139 return (*this);
3140 }
3141
3142 /* Be careful trying to make this operation lazy.
3143 At time of writing, resolve() and resolveSample() do not throw.
3144 Changing this would mean that any resolve call would need to use MPI (to check for global errors)
3145 */
3146 Data
3147 Data::matrixInverse() const
3148 {
3149 if (isLazy()) // Cannot use lazy for this because individual inversions could throw.
3150 {
3151 Data d(*this);
3152 d.resolve();
3153 return d.matrixInverse();
3154 }
3155 THROWONCOMPLEX
3156 Data out(0.,getDataPointShape(),getFunctionSpace(),false);
3157 out.typeMatchRight(*this);
3158
3159 DataReady* drp=out.getReadyPtr().get();
3160 int errcode=m_data->matrixInverse(drp);
3161 #ifdef ESYS_MPI
3162 int globalval=0;
3163 MPI_Allreduce( &errcode, &globalval, 1, MPI_INT, MPI_MAX, get_MPIComm() );
3164 errcode=globalval;
3165 #endif
3166 if (errcode)
3167 {
3168 escript::matrixInverseError(errcode); // throws exceptions
3169 }
3170 return out;
3171 }
3172
3173
3174 Data
3175 Data::rpowO(const bp::object& left) const
3176 {
3177 Data left_d(left,*this);
3178 return left_d.powD(*this);
3179 }
3180
3181 Data
3182 Data::powO(const bp::object& right) const
3183 {
3184 Data tmp(right,getFunctionSpace(),false);
3185 return powD(tmp);
3186 }
3187
3188 Data
3189 Data::powD(const Data& right) const
3190 {
3191 MAKELAZYBIN(right,POW);
3192
3193 return C_TensorBinaryOperation(*this, right, ES_optype::POW);
3194
3195 }
3196
3197
3198 //
3199 // NOTE: It is essential to specify the namespace this operator belongs to
3200 Data
3201 escript::operator+(const Data& left, const Data& right)
3202 {
3203 MAKELAZYBIN2(left,right,ADD);
3204
3205 return C_TensorBinaryOperation(left, right, ES_optype::ADD);
3206 }
3207
3208 //
3209 // NOTE: It is essential to specify the namespace this operator belongs to
3210 Data
3211 escript::operator-(const Data& left, const Data& right)
3212 {
3213 MAKELAZYBIN2(left,right,SUB);
3214 return C_TensorBinaryOperation(left, right, ES_optype::SUB);
3215 }
3216
3217 //
3218 // NOTE: It is essential to specify the namespace this operator belongs to
3219 Data
3220 escript::operator*(const Data& left, const Data& right)
3221 {
3222 MAKELAZYBIN2(left,right,MUL);
3223
3224 return C_TensorBinaryOperation(left, right, ES_optype::MUL);
3225 }
3226
3227 //
3228 // NOTE: It is essential to specify the namespace this operator belongs to
3229 Data
3230 escript::operator/(const Data& left, const Data& right)
3231 {
3232 MAKELAZYBIN2(left,right,DIV);
3233 return C_TensorBinaryOperation(left, right, ES_optype::DIV);
3234 }
3235
3236 //
3237 // NOTE: It is essential to specify the namespace this operator belongs to
3238 Data
3239 escript::operator+(const Data& left, const boost::python::object& right)
3240 {
3241 Data tmp(right,left.getFunctionSpace(),false);
3242 MAKELAZYBIN2(left,tmp,ADD);
3243 return left+tmp;
3244 }
3245
3246 //
3247 // NOTE: It is essential to specify the namespace this operator belongs to
3248 Data
3249 escript::operator-(const Data& left, const boost::python::object& right)
3250 {
3251 Data tmp(right,left.getFunctionSpace(),false);
3252 MAKELAZYBIN2(left,tmp,SUB);
3253 return left-tmp;
3254 }
3255
3256 //
3257 // NOTE: It is essential to specify the namespace this operator belongs to
3258 Data
3259 escript::operator*(const Data& left, const boost::python::object& right)
3260 {
3261 Data tmp(right,left.getFunctionSpace(),false);
3262 MAKELAZYBIN2(left,tmp,MUL);
3263 return left*tmp;
3264 }
3265
3266 //
3267 // NOTE: It is essential to specify the namespace this operator belongs to
3268 Data
3269 escript::operator/(const Data& left, const boost::python::object& right)
3270 {
3271 Data tmp(right,left.getFunctionSpace(),false);
3272 MAKELAZYBIN2(left,tmp,DIV);
3273 return left/tmp;
3274 }
3275
3276 //
3277 // NOTE: It is essential to specify the namespace this operator belongs to
3278 Data
3279 escript::operator+(const boost::python::object& left, const Data& right)
3280 {
3281 Data tmp(left,right.getFunctionSpace(),false);
3282 MAKELAZYBIN2(tmp,right,ADD);
3283 return tmp+right;
3284 }
3285
3286 //
3287 // NOTE: It is essential to specify the namespace this operator belongs to
3288 Data
3289 escript::operator-(const boost::python::object& left, const Data& right)
3290 {
3291 Data tmp(left,right.getFunctionSpace(),false);
3292 MAKELAZYBIN2(tmp,right,SUB);
3293 return tmp-right;
3294 }
3295
3296 //
3297 // NOTE: It is essential to specify the namespace this operator belongs to
3298 Data
3299 escript::operator*(const boost::python::object& left, const Data& right)
3300 {
3301 Data tmp(left,right.getFunctionSpace(),false);
3302 MAKELAZYBIN2(tmp,right,MUL);
3303 return tmp*right;
3304 }
3305
3306 //
3307 // NOTE: It is essential to specify the namespace this operator belongs to
3308 Data
3309 escript::operator/(const boost::python::object& left, const Data& right)
3310 {
3311 Data tmp(left,right.getFunctionSpace(),false);
3312 MAKELAZYBIN2(tmp,right,DIV);
3313 return tmp/right;
3314 }
3315
3316
3317 /* TODO */
3318 /* global reduction */
3319 Data
3320 Data::getItem(const bp::object& key) const
3321 {
3322 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
3323
3324 if (slice_region.size()!=getDataPointRank()) {
3325 throw DataException("Error - slice size does not match Data rank.");
3326 }
3327
3328 return getSlice(slice_region);
3329 }
3330
3331 /* TODO */
3332 /* global reduction */
3333 Data
3334 Data::getSlice(const DataTypes::RegionType& region) const
3335 {
3336 return Data(*this,region);
3337 }
3338
3339 /* TODO */
3340 /* global reduction */
3341 void
3342 Data::setItemO(const bp::object& key,
3343 const bp::object& value)
3344 {
3345 Data tempData(value,getFunctionSpace(),false);
3346 setItemD(key,tempData);
3347 }
3348
3349 void
3350 Data::setItemD(const bp::object& key,
3351 const Data& value)
3352 {
3353 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
3354 if (slice_region.size()!=getDataPointRank()) {
3355 throw DataException("Error - slice size does not match Data rank.");
3356 }
3357 exclusiveWrite();
3358 if (getFunctionSpace()!=value.getFunctionSpace()) {
3359 setSlice(Data(value,getFunctionSpace()),slice_region);
3360 } else {
3361 setSlice(value,slice_region);
3362 }
3363 }
3364
3365 void
3366 Data::setSlice(const Data& value,
3367 const DataTypes::RegionType& region)
3368 {
3369 if (isProtected()) {
3370 throw DataException("Error - attempt to update protected Data object.");
3371 }
3372 forceResolve();
3373 exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
3374 Data tempValue(value);
3375 typeMatchLeft(tempValue);
3376 typeMatchRight(tempValue);
3377 getReady()->setSlice(tempValue.m_data.get(),region);
3378 }
3379
3380 void
3381 Data::typeMatchLeft(Data& right) const
3382 {
3383 if (right.isLazy() && !isLazy())
3384 {
3385 right.resolve();
3386 }
3387 if (isComplex())
3388 {
3389 right.complicate();
3390 }
3391 if (isExpanded()) {
3392 right.expand();
3393 } else if (isTagged()) {
3394 if (right.isConstant()) {
3395 right.tag();
3396 }
3397 }
3398 }
3399
3400 void
3401 Data::typeMatchRight(const Data& right)
3402 {
3403 if (isLazy() && !right.isLazy())
3404 {
3405 resolve();
3406 }
3407 if (right.isComplex())
3408 {
3409 complicate();
3410 }
3411 if (isTagged()) {
3412 if (right.isExpanded()) {
3413 expand();
3414 }
3415 } else if (isConstant()) {
3416 if (right.isExpanded()) {
3417 expand();
3418 } else if (right.isTagged()) {
3419 tag();
3420 }
3421 }
3422 }
3423
3424 // The normal TaggedValue adds the tag if it is not already present
3425 // This form does not. It throws instead.
3426 // This is because the names are maintained by the domain and cannot be added
3427 // without knowing the tag number to map it to.
3428 void
3429 Data::setTaggedValueByName(std::string name,
3430 const bp::object& value)
3431 {
3432 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
3433 forceResolve();
3434 exclusiveWrite();
3435 int tagKey=getFunctionSpace().getDomain()->getTag(name);
3436 setTaggedValue(tagKey,value);
3437 }
3438 else
3439 {
3440 std::string msg="Error - unknown tag ("+name+") in setTaggedValueByName.";
3441 throw DataException(msg);
3442 }
3443 }
3444
3445 void
3446 Data::setTaggedValue(int tagKey,
3447