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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


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