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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2492 - (hide annotations)
Wed Jun 24 01:28:57 2009 UTC (10 years, 9 months ago) by gross
File size: 86474 byte(s)
MPI_DOUBLE changed to MPI_INT as lowProc is int
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 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
156 jfenwick 2459 }
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 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
191 jfenwick 2459 }
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 2481
1072 jfenwick 2459 // 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 jfenwick 2476 if (local_min<min) { // If we found a smaller value than our sentinel
1960 jgs 148 min=local_min;
1961     lowi=local_lowi;
1962     lowj=local_lowj;
1963     }
1964 jgs 121 }
1965    
1966 bcumming 782 #ifdef PASO_MPI
1967 jfenwick 2476 // determine the processor on which the minimum occurs
1968     next = temp.getDataPointRO(lowi,lowj);
1969     int lowProc = 0;
1970     double *globalMins = new double[get_MPISize()+1];
1971     int error;
1972     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1973 ksteube 1312
1974 jfenwick 2476 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 gross 2492 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
1983 bcumming 782
1984 jfenwick 2476 delete [] globalMins;
1985     ProcNo = lowProc;
1986 bcumming 790 #else
1987 jfenwick 2476 ProcNo = 0;
1988 bcumming 782 #endif
1989 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1990 jgs 121 }
1991    
1992 jfenwick 2476
1993     const boost::python::tuple
1994     Data::maxGlobalDataPoint() const
1995     {
1996     int DataPointNo;
1997     int ProcNo;
1998     calc_maxGlobalDataPoint(ProcNo,DataPointNo);
1999     return make_tuple(ProcNo,DataPointNo);
2000     }
2001    
2002 jgs 104 void
2003 jfenwick 2476 Data::calc_maxGlobalDataPoint(int& ProcNo,
2004     int& DataPointNo) const
2005     {
2006     if (isLazy())
2007     {
2008     Data temp(*this); // to get around the fact that you can't resolve a const Data
2009     temp.resolve();
2010     return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2011     }
2012     int i,j;
2013     int highi=0,highj=0;
2014     //-------------
2015     double max=numeric_limits<double>::min();
2016    
2017     Data temp=maxval();
2018    
2019     int numSamples=temp.getNumSamples();
2020     int numDPPSample=temp.getNumDataPointsPerSample();
2021    
2022     double next,local_max;
2023     int local_highi=0,local_highj=0;
2024    
2025     #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2026     {
2027     local_max=max;
2028     #pragma omp for private(i,j) schedule(static)
2029     for (i=0; i<numSamples; i++) {
2030     for (j=0; j<numDPPSample; j++) {
2031     next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2032     if (next>local_max) {
2033     local_max=next;
2034     local_highi=i;
2035     local_highj=j;
2036     }
2037     }
2038     }
2039     #pragma omp critical
2040     if (local_max>max) { // If we found a larger value than our sentinel
2041     max=local_max;
2042     highi=local_highi;
2043     highj=local_highj;
2044     }
2045     }
2046    
2047     #ifdef PASO_MPI
2048     // determine the processor on which the maximum occurs
2049     next = temp.getDataPointRO(highi,highj);
2050     int highProc = 0;
2051     double *globalMaxs = new double[get_MPISize()+1];
2052     int error;
2053     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2054    
2055     if( get_MPIRank()==0 ){
2056     next = globalMaxs[highProc];
2057     for( i=1; i<get_MPISize(); i++ )
2058     if( next>globalMaxs[i] ){
2059     highProc = i;
2060     next = globalMaxs[i];
2061     }
2062     }
2063 gross 2492 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2064 jfenwick 2476 delete [] globalMaxs;
2065     ProcNo = highProc;
2066     #else
2067     ProcNo = 0;
2068     #endif
2069     DataPointNo = highj + highi * numDPPSample;
2070     }
2071    
2072     void
2073 jgs 104 Data::saveDX(std::string fileName) const
2074     {
2075 jfenwick 1803 if (isEmpty())
2076     {
2077     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2078     }
2079 jfenwick 2005 if (isLazy())
2080     {
2081     Data temp(*this); // to get around the fact that you can't resolve a const Data
2082     temp.resolve();
2083     temp.saveDX(fileName);
2084     return;
2085     }
2086 jgs 153 boost::python::dict args;
2087     args["data"]=boost::python::object(this);
2088 jfenwick 1872 getDomain()->saveDX(fileName,args);
2089 jgs 104 return;
2090     }
2091    
2092 jgs 110 void
2093     Data::saveVTK(std::string fileName) const
2094     {
2095 jfenwick 1803 if (isEmpty())
2096     {
2097     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2098     }
2099 jfenwick 2005 if (isLazy())
2100     {
2101     Data temp(*this); // to get around the fact that you can't resolve a const Data
2102     temp.resolve();
2103     temp.saveVTK(fileName);
2104     return;
2105     }
2106 jgs 153 boost::python::dict args;
2107     args["data"]=boost::python::object(this);
2108 gross 2421 getDomain()->saveVTK(fileName,args,"","");
2109 jgs 110 return;
2110     }
2111    
2112 jgs 102 Data&
2113     Data::operator+=(const Data& right)
2114     {
2115 gross 783 if (isProtected()) {
2116     throw DataException("Error - attempt to update protected Data object.");
2117     }
2118 jfenwick 2146 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
2119 jfenwick 2271 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2120 jfenwick 2146 binaryOp(right,plus<double>());
2121     return (*this);
2122 jgs 94 }
2123    
2124 jgs 102 Data&
2125     Data::operator+=(const boost::python::object& right)
2126 jgs 94 {
2127 jfenwick 2005 if (isProtected()) {
2128     throw DataException("Error - attempt to update protected Data object.");
2129     }
2130 gross 854 Data tmp(right,getFunctionSpace(),false);
2131 jfenwick 2271 (*this)+=tmp;
2132     return *this;
2133 jgs 94 }
2134 jfenwick 2005
2135     // Hmmm, operator= makes a deep copy but the copy constructor does not?
2136 ksteube 1312 Data&
2137     Data::operator=(const Data& other)
2138     {
2139 jfenwick 2271 #if defined ASSIGNMENT_MEANS_DEEPCOPY
2140     // This should not be used.
2141     // Just leaving this here until I have completed transition
2142 ksteube 1312 copy(other);
2143 jfenwick 2271 #else
2144     m_protected=false; // since any changes should be caught by exclusiveWrite();
2145     // m_data=other.m_data;
2146     set_m_data(other.m_data);
2147     #endif
2148 ksteube 1312 return (*this);
2149     }
2150 jgs 94
2151 jgs 102 Data&
2152     Data::operator-=(const Data& right)
2153 jgs 94 {
2154 gross 783 if (isProtected()) {
2155     throw DataException("Error - attempt to update protected Data object.");
2156     }
2157 jfenwick 2146 MAKELAZYBINSELF(right,SUB)
2158 jfenwick 2271 exclusiveWrite();
2159 jfenwick 2146 binaryOp(right,minus<double>());
2160     return (*this);
2161 jgs 94 }
2162    
2163 jgs 102 Data&
2164     Data::operator-=(const boost::python::object& right)
2165 jgs 94 {
2166 jfenwick 2005 if (isProtected()) {
2167     throw DataException("Error - attempt to update protected Data object.");
2168     }
2169 gross 854 Data tmp(right,getFunctionSpace(),false);
2170 jfenwick 2271 (*this)-=tmp;
2171 jfenwick 2146 return (*this);
2172 jgs 94 }
2173    
2174 jgs 102 Data&
2175     Data::operator*=(const Data& right)
2176 jgs 94 {
2177 gross 783 if (isProtected()) {
2178     throw DataException("Error - attempt to update protected Data object.");
2179     }
2180 jfenwick 2146 MAKELAZYBINSELF(right,MUL)
2181 jfenwick 2271 exclusiveWrite();
2182 jfenwick 2146 binaryOp(right,multiplies<double>());
2183     return (*this);
2184 jgs 94 }
2185    
2186 jgs 102 Data&
2187     Data::operator*=(const boost::python::object& right)
2188 jfenwick 2005 {
2189     if (isProtected()) {
2190     throw DataException("Error - attempt to update protected Data object.");
2191     }
2192 gross 854 Data tmp(right,getFunctionSpace(),false);
2193 jfenwick 2271 (*this)*=tmp;
2194 jfenwick 2146 return (*this);
2195 jgs 94 }
2196    
2197 jgs 102 Data&
2198     Data::operator/=(const Data& right)
2199 jgs 94 {
2200 gross 783 if (isProtected()) {
2201     throw DataException("Error - attempt to update protected Data object.");
2202     }
2203 jfenwick 2146 MAKELAZYBINSELF(right,DIV)
2204 jfenwick 2271 exclusiveWrite();
2205 jfenwick 2146 binaryOp(right,divides<double>());
2206     return (*this);
2207 jgs 94 }
2208    
2209 jgs 102 Data&
2210     Data::operator/=(const boost::python::object& right)
2211 jgs 94 {
2212 jfenwick 2005 if (isProtected()) {
2213     throw DataException("Error - attempt to update protected Data object.");
2214     }
2215 gross 854 Data tmp(right,getFunctionSpace(),false);
2216 jfenwick 2271 (*this)/=tmp;
2217 jfenwick 2146 return (*this);
2218 jgs 94 }
2219    
2220 jgs 102 Data
2221 gross 699 Data::rpowO(const boost::python::object& left) const
2222     {
2223     Data left_d(left,*this);
2224     return left_d.powD(*this);
2225     }
2226    
2227     Data
2228 jgs 102 Data::powO(const boost::python::object& right) const
2229 jgs 94 {
2230 gross 854 Data tmp(right,getFunctionSpace(),false);
2231     return powD(tmp);
2232 jgs 94 }
2233    
2234 jgs 102 Data
2235     Data::powD(const Data& right) const
2236 jgs 94 {
2237 jfenwick 2146 MAKELAZYBIN(right,POW)
2238 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2239 jgs 94 }
2240    
2241     //
2242 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2243 jgs 102 Data
2244     escript::operator+(const Data& left, const Data& right)
2245 jgs 94 {
2246 jfenwick 2146 MAKELAZYBIN2(left,right,ADD)
2247 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
2248 jgs 94 }
2249    
2250     //
2251 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2252 jgs 102 Data
2253     escript::operator-(const Data& left, const Data& right)
2254 jgs 94 {
2255 jfenwick 2146 MAKELAZYBIN2(left,right,SUB)
2256 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2257 jgs 94 }
2258    
2259     //
2260 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2261 jgs 102 Data
2262     escript::operator*(const Data& left, const Data& right)
2263 jgs 94 {
2264 jfenwick 2146 MAKELAZYBIN2(left,right,MUL)
2265 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2266 jgs 94 }
2267    
2268     //
2269 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2270 jgs 102 Data
2271     escript::operator/(const Data& left, const Data& right)
2272 jgs 94 {
2273 jfenwick 2146 MAKELAZYBIN2(left,right,DIV)
2274 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2275 jgs 94 }
2276    
2277     //
2278 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2279 jgs 102 Data
2280     escript::operator+(const Data& left, const boost::python::object& right)
2281 jgs 94 {
2282 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2283     MAKELAZYBIN2(left,tmp,ADD)
2284     return left+tmp;
2285 jgs 94 }
2286    
2287     //
2288 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2289 jgs 102 Data
2290     escript::operator-(const Data& left, const boost::python::object& right)
2291 jgs 94 {
2292 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2293     MAKELAZYBIN2(left,tmp,SUB)
2294     return left-tmp;
2295 jgs 94 }
2296    
2297     //
2298 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2299 jgs 102 Data
2300     escript::operator*(const Data& left, const boost::python::object& right)
2301 jgs 94 {
2302 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2303     MAKELAZYBIN2(left,tmp,MUL)
2304     return left*tmp;
2305 jgs 94 }
2306    
2307     //
2308 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2309 jgs 102 Data
2310     escript::operator/(const Data& left, const boost::python::object& right)
2311 jgs 94 {
2312 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2313     MAKELAZYBIN2(left,tmp,DIV)
2314     return left/tmp;
2315 jgs 94 }
2316    
2317     //
2318 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2319 jgs 102 Data
2320     escript::operator+(const boost::python::object& left, const Data& right)
2321 jgs 94 {
2322 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2323     MAKELAZYBIN2(tmp,right,ADD)
2324     return tmp+right;
2325 jgs 94 }
2326    
2327     //
2328 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2329 jgs 102 Data
2330     escript::operator-(const boost::python::object& left, const Data& right)
2331 jgs 94 {
2332 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2333     MAKELAZYBIN2(tmp,right,SUB)
2334     return tmp-right;
2335 jgs 94 }
2336    
2337     //
2338 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2339 jgs 102 Data
2340     escript::operator*(const boost::python::object& left, const Data& right)
2341 jgs 94 {
2342 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2343     MAKELAZYBIN2(tmp,right,MUL)
2344     return tmp*right;
2345 jgs 94 }
2346    
2347     //
2348 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2349 jgs 102 Data
2350     escript::operator/(const boost::python::object& left, const Data& right)
2351 jgs 94 {
2352 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2353     MAKELAZYBIN2(tmp,right,DIV)
2354     return tmp/right;
2355 jgs 94 }
2356    
2357 jgs 102
2358 bcumming 751 /* TODO */
2359     /* global reduction */
2360 jgs 102 Data
2361 ksteube 1312 Data::getItem(const boost::python::object& key) const
2362 jgs 94 {
2363    
2364 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2365 jgs 94
2366 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2367 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2368 jgs 94 }
2369    
2370 jgs 102 return getSlice(slice_region);
2371 jgs 94 }
2372    
2373 bcumming 751 /* TODO */
2374     /* global reduction */
2375 jgs 94 Data
2376 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2377 jgs 94 {
2378 jgs 102 return Data(*this,region);
2379 jgs 94 }
2380    
2381 bcumming 751 /* TODO */
2382     /* global reduction */
2383 jgs 94 void
2384 jgs 102 Data::setItemO(const boost::python::object& key,
2385     const boost::python::object& value)
2386 jgs 94 {
2387 jgs 102 Data tempData(value,getFunctionSpace());
2388     setItemD(key,tempData);
2389     }
2390    
2391     void
2392     Data::setItemD(const boost::python::object& key,
2393     const Data& value)
2394     {
2395 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2396     if (slice_region.size()!=getDataPointRank()) {
2397 jgs 94 throw DataException("Error - slice size does not match Data rank.");
2398     }
2399 jfenwick 2271 exclusiveWrite();
2400 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2401     setSlice(Data(value,getFunctionSpace()),slice_region);
2402     } else {
2403     setSlice(value,slice_region);
2404     }
2405 jgs 94 }
2406    
2407     void
2408 jgs 102 Data::setSlice(const Data& value,
2409 jfenwick 1796 const DataTypes::RegionType& region)
2410 jgs 94 {
2411 gross 783 if (isProtected()) {
2412     throw DataException("Error - attempt to update protected Data object.");
2413     }
2414 jfenwick 2271 forceResolve();
2415     exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2416 jgs 102 Data tempValue(value);
2417     typeMatchLeft(tempValue);
2418     typeMatchRight(tempValue);
2419 jfenwick 2005 getReady()->setSlice(tempValue.m_data.get(),region);
2420 jgs 102 }
2421    
2422     void
2423     Data::typeMatchLeft(Data& right) const
2424     {
2425 jfenwick 2005 if (right.isLazy() && !isLazy())
2426     {
2427     right.resolve();
2428     }
2429 jgs 102 if (isExpanded()){
2430     right.expand();
2431     } else if (isTagged()) {
2432     if (right.isConstant()) {
2433     right.tag();
2434     }
2435     }
2436     }
2437    
2438     void
2439     Data::typeMatchRight(const Data& right)
2440     {
2441 jfenwick 2005 if (isLazy() && !right.isLazy())
2442     {
2443     resolve();
2444     }
2445 jgs 94 if (isTagged()) {
2446     if (right.isExpanded()) {
2447     expand();
2448     }
2449     } else if (isConstant()) {
2450     if (right.isExpanded()) {
2451     expand();
2452     } else if (right.isTagged()) {
2453     tag();
2454     }
2455     }
2456     }
2457    
2458 jfenwick 2282 // The normal TaggedValue adds the tag if it is not already present
2459     // This form does not. It throws instead.
2460     // This is because the names are maintained by the domain and cannot be added
2461     // without knowing the tag number to map it to.
2462 jgs 94 void
2463 gross 1044 Data::setTaggedValueByName(std::string name,
2464 ksteube 1312 const boost::python::object& value)
2465 gross 1044 {
2466 jfenwick 1872 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2467 jfenwick 2271 forceResolve();
2468     exclusiveWrite();
2469 jfenwick 1872 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2470 gross 1044 setTaggedValue(tagKey,value);
2471     }
2472 jfenwick 2282 else
2473     { // The