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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2644 - (hide annotations)
Wed Sep 2 04:14:03 2009 UTC (10 years, 7 months ago) by jfenwick
File size: 89214 byte(s)
Add unit tests for saveDataCSV which should be ready for use now.
Keyword args are now output in sorted order.

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