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

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

Parent Directory Parent Directory | Revision Log Revision Log


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