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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2653 - (hide annotations)
Tue Sep 8 04:26:30 2009 UTC (10 years, 6 months ago) by jfenwick
File size: 90989 byte(s)
Fix bug in maxGlobalDataPoint and minGlobalDataPoint.
They now give the correct answers and the datapoint ids returned are globally
correct.

Removed some #defines from before COW
Removed hasNoSamples() - I don't trust myself to use that properly let alone anybody else.


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