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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2078 - (hide annotations)
Thu Nov 20 16:10:10 2008 UTC (11 years, 4 months ago) by phornby
File size: 86256 byte(s)
Two changes.

1. Move blocktimer from escript to esysUtils.
2. Make it possible to link to paso as a DLL or .so.

Should have no effect on 'nix's

In respect of 1., blocktimer had begun to spring up everywhere, so
for the moment I thought it best to move it to the only other library that
pops up all over the place.

In respect of 2., paso needed to be a DLL in order to use the windows intelc /fast
option, which does aggressive multi-file optimisations. Even in its current form, it either
vectorises or parallelises  hundreds more loops in the esys system than appear in the pragmas.

In achieving 2. I have not been too delicate in adding

PASO_DLL_API

declarations to the .h files in paso/src. Only toward the end of the process of
the conversion, when the number of linker errors dropped below 20, say, did I choosy about what
functions in a header I declared PASO_DLL_API. As a result, there are likely to be many routines
declared as external functions symbols that are in fact internal to the paso DLL. 
Why is this an issue?  It prevents the intelc compiler from getting aggressive on the paso module.
With pain there is sometimes gain. At least all the DLL rules in windows give good
(non-microsoft) compiler writers a chance to really shine.

So, if you should see a PASO_DLL_API on a function in a paso header file,
and think to yourself, "that function is only called in paso, why export it?", then feel free to
delete the PASO_DLL_API export declaration.

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