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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2084 - (hide annotations)
Fri Nov 21 05:20:42 2008 UTC (11 years, 4 months ago) by jfenwick
File size: 85826 byte(s)
Fixed a warning in cpp unit tests under dodebug
Pointed the url for python doco at shake200 rather than iservo.
Added support for trace and transpose to LazyData.
Fixed bug in trace to initialise running totals.
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     Data
1726 gross 800 Data::trace(int axis_offset) const
1727 jfenwick 2084 {
1728 jfenwick 2005 if (isLazy())
1729     {
1730 jfenwick 2084 DataLazy* c=new DataLazy(borrowDataPtr(),TRACE,axis_offset);
1731     return Data(c);
1732 jfenwick 2005 }
1733 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1734 ksteube 775 if (getDataPointRank()==2) {
1735 jfenwick 1796 DataTypes::ShapeType ev_shape;
1736 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1737     ev.typeMatchRight(*this);
1738 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1739 ksteube 775 return ev;
1740     }
1741     if (getDataPointRank()==3) {
1742 jfenwick 1796 DataTypes::ShapeType ev_shape;
1743 ksteube 775 if (axis_offset==0) {
1744     int s2=s[2];
1745     ev_shape.push_back(s2);
1746     }
1747     else if (axis_offset==1) {
1748     int s0=s[0];
1749     ev_shape.push_back(s0);
1750     }
1751     Data ev(0.,ev_shape,getFunctionSpace());
1752     ev.typeMatchRight(*this);
1753 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1754 ksteube 775 return ev;
1755     }
1756     if (getDataPointRank()==4) {
1757 jfenwick 1796 DataTypes::ShapeType ev_shape;
1758 ksteube 775 if (axis_offset==0) {
1759     ev_shape.push_back(s[2]);
1760     ev_shape.push_back(s[3]);
1761     }
1762     else if (axis_offset==1) {
1763     ev_shape.push_back(s[0]);
1764     ev_shape.push_back(s[3]);
1765     }
1766     else if (axis_offset==2) {
1767     ev_shape.push_back(s[0]);
1768     ev_shape.push_back(s[1]);
1769     }
1770     Data ev(0.,ev_shape,getFunctionSpace());
1771     ev.typeMatchRight(*this);
1772 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1773 ksteube 775 return ev;
1774     }
1775     else {
1776 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1777 ksteube 775 }
1778     }
1779    
1780     Data
1781     Data::transpose(int axis_offset) const
1782 jfenwick 2005 {
1783     if (isLazy())
1784     {
1785 jfenwick 2084 DataLazy* c=new DataLazy(borrowDataPtr(),TRANS,axis_offset);
1786     return Data(c);
1787 jfenwick 2005 }
1788 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1789     DataTypes::ShapeType ev_shape;
1790 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1791     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1792     int rank=getDataPointRank();
1793     if (axis_offset<0 || axis_offset>rank) {
1794     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1795     }
1796     for (int i=0; i<rank; i++) {
1797     int index = (axis_offset+i)%rank;
1798     ev_shape.push_back(s[index]); // Append to new shape
1799     }
1800     Data ev(0.,ev_shape,getFunctionSpace());
1801     ev.typeMatchRight(*this);
1802     m_data->transpose(ev.m_data.get(), axis_offset);
1803     return ev;
1804 jgs 123 }
1805    
1806 gross 576 Data
1807     Data::eigenvalues() const
1808     {
1809 jfenwick 2005 if (isLazy())
1810     {
1811     Data temp(*this); // to get around the fact that you can't resolve a const Data
1812     temp.resolve();
1813     return temp.eigenvalues();
1814     }
1815 gross 576 // check input
1816 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1817 ksteube 1312 if (getDataPointRank()!=2)
1818 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1819 ksteube 1312 if(s[0] != s[1])
1820 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1821     // create return
1822 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1823 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1824     ev.typeMatchRight(*this);
1825     m_data->eigenvalues(ev.m_data.get());
1826     return ev;
1827     }
1828    
1829 jgs 121 const boost::python::tuple
1830 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1831     {
1832 jfenwick 2005 if (isLazy())
1833     {
1834     Data temp(*this); // to get around the fact that you can't resolve a const Data
1835     temp.resolve();
1836     return temp.eigenvalues_and_eigenvectors(tol);
1837     }
1838 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1839 ksteube 1312 if (getDataPointRank()!=2)
1840 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1841 ksteube 1312 if(s[0] != s[1])
1842 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1843     // create return
1844 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1845 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1846     ev.typeMatchRight(*this);
1847 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1848 gross 576 Data V(0.,V_shape,getFunctionSpace());
1849     V.typeMatchRight(*this);
1850     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1851     return make_tuple(boost::python::object(ev),boost::python::object(V));
1852     }
1853    
1854     const boost::python::tuple
1855 gross 921 Data::minGlobalDataPoint() const
1856 jgs 121 {
1857 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1858 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1859     // surrounding function
1860    
1861     int DataPointNo;
1862 gross 921 int ProcNo;
1863     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1864     return make_tuple(ProcNo,DataPointNo);
1865 jgs 148 }
1866    
1867     void
1868 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1869     int& DataPointNo) const
1870 jgs 148 {
1871 jfenwick 2005 if (isLazy())
1872     {
1873     Data temp(*this); // to get around the fact that you can't resolve a const Data
1874     temp.resolve();
1875     return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
1876     }
1877 jgs 148 int i,j;
1878     int lowi=0,lowj=0;
1879     double min=numeric_limits<double>::max();
1880    
1881 jgs 121 Data temp=minval();
1882    
1883     int numSamples=temp.getNumSamples();
1884     int numDPPSample=temp.getNumDataPointsPerSample();
1885    
1886 jgs 148 double next,local_min;
1887 jfenwick 1946 int local_lowi=0,local_lowj=0;
1888 jgs 121
1889 caltinay 2081 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1890 jgs 148 {
1891     local_min=min;
1892     #pragma omp for private(i,j) schedule(static)
1893     for (i=0; i<numSamples; i++) {
1894     for (j=0; j<numDPPSample; j++) {
1895 jfenwick 1796 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1896 jgs 148 if (next<local_min) {
1897     local_min=next;
1898     local_lowi=i;
1899     local_lowj=j;
1900     }
1901 jgs 121 }
1902     }
1903 jgs 148 #pragma omp critical
1904     if (local_min<min) {
1905     min=local_min;
1906     lowi=local_lowi;
1907     lowj=local_lowj;
1908     }
1909 jgs 121 }
1910    
1911 bcumming 782 #ifdef PASO_MPI
1912     // determine the processor on which the minimum occurs
1913 jfenwick 1796 next = temp.getDataPoint(lowi,lowj);
1914 bcumming 782 int lowProc = 0;
1915     double *globalMins = new double[get_MPISize()+1];
1916 caltinay 2081 int error;
1917     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1918 ksteube 1312
1919 bcumming 782 if( get_MPIRank()==0 ){
1920     next = globalMins[lowProc];
1921     for( i=1; i<get_MPISize(); i++ )
1922     if( next>globalMins[i] ){
1923     lowProc = i;
1924     next = globalMins[i];
1925     }
1926     }
1927     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1928    
1929     delete [] globalMins;
1930     ProcNo = lowProc;
1931 bcumming 790 #else
1932     ProcNo = 0;
1933 bcumming 782 #endif
1934 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1935 jgs 121 }
1936    
1937 jgs 104 void
1938     Data::saveDX(std::string fileName) const
1939     {
1940 jfenwick 1803 if (isEmpty())
1941     {
1942     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1943     }
1944 jfenwick 2005 if (isLazy())
1945     {
1946     Data temp(*this); // to get around the fact that you can't resolve a const Data
1947     temp.resolve();
1948     temp.saveDX(fileName);
1949     return;
1950     }
1951 jgs 153 boost::python::dict args;
1952     args["data"]=boost::python::object(this);
1953 jfenwick 1872 getDomain()->saveDX(fileName,args);
1954 jgs 104 return;
1955     }
1956    
1957 jgs 110 void
1958     Data::saveVTK(std::string fileName) const
1959     {
1960 jfenwick 1803 if (isEmpty())
1961     {
1962     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1963     }
1964 jfenwick 2005 if (isLazy())
1965     {
1966     Data temp(*this); // to get around the fact that you can't resolve a const Data
1967     temp.resolve();
1968     temp.saveVTK(fileName);
1969     return;
1970     }
1971 jgs 153 boost::python::dict args;
1972     args["data"]=boost::python::object(this);
1973 jfenwick 1872 getDomain()->saveVTK(fileName,args);
1974 jgs 110 return;
1975     }
1976    
1977 jgs 102 Data&
1978     Data::operator+=(const Data& right)
1979     {
1980 gross 783 if (isProtected()) {
1981     throw DataException("Error - attempt to update protected Data object.");
1982     }
1983 jfenwick 2005 if (isLazy() || right.isLazy())
1984     {
1985     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
1986     m_data=c->getPtr();
1987     return (*this);
1988     }
1989     else
1990     {
1991     binaryOp(right,plus<double>());
1992     return (*this);
1993     }
1994 jgs 94 }
1995    
1996 jgs 102 Data&
1997     Data::operator+=(const boost::python::object& right)
1998 jgs 94 {
1999 jfenwick 2005 if (isProtected()) {
2000     throw DataException("Error - attempt to update protected Data object.");
2001     }
2002 gross 854 Data tmp(right,getFunctionSpace(),false);
2003 jfenwick 2005 if (isLazy())
2004     {
2005     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),ADD); // for lazy + is equivalent to +=
2006     m_data=c->getPtr();
2007     return (*this);
2008     }
2009     else
2010     {
2011     binaryOp(tmp,plus<double>());
2012     return (*this);
2013     }
2014 jgs 94 }
2015 jfenwick 2005
2016     // Hmmm, operator= makes a deep copy but the copy constructor does not?
2017 ksteube 1312 Data&
2018     Data::operator=(const Data& other)
2019     {
2020     copy(other);
2021     return (*this);
2022     }
2023 jgs 94
2024 jgs 102 Data&
2025     Data::operator-=(const Data& right)
2026 jgs 94 {
2027 gross 783 if (isProtected()) {
2028     throw DataException("Error - attempt to update protected Data object.");
2029     }
2030 jfenwick 2005 if (isLazy() || right.isLazy())
2031     {
2032     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2033     m_data=c->getPtr();
2034     return (*this);
2035     }
2036     else
2037     {
2038     binaryOp(right,minus<double>());
2039     return (*this);
2040     }
2041 jgs 94 }
2042    
2043 jgs 102 Data&
2044     Data::operator-=(const boost::python::object& right)
2045 jgs 94 {
2046 jfenwick 2005 if (isProtected()) {
2047     throw DataException("Error - attempt to update protected Data object.");
2048     }
2049 gross 854 Data tmp(right,getFunctionSpace(),false);
2050 jfenwick 2005 if (isLazy())
2051     {
2052     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),SUB); // for lazy - is equivalent to -=
2053     m_data=c->getPtr();
2054     return (*this);
2055     }
2056     else
2057     {
2058     binaryOp(tmp,minus<double>());
2059     return (*this);
2060     }
2061 jgs 94 }
2062    
2063 jgs 102 Data&
2064     Data::operator*=(const Data& right)
2065 jgs 94 {
2066 gross 783 if (isProtected()) {
2067     throw DataException("Error - attempt to update protected Data object.");
2068     }
2069 jfenwick 2005 if (isLazy() || right.isLazy())
2070     {
2071     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2072     m_data=c->getPtr();
2073     return (*this);
2074     }
2075     else
2076     {
2077     binaryOp(right,multiplies<double>());
2078     return (*this);
2079     }
2080 jgs 94 }
2081    
2082 jgs 102 Data&
2083     Data::operator*=(const boost::python::object& right)
2084 jfenwick 2005 {
2085     if (isProtected()) {
2086     throw DataException("Error - attempt to update protected Data object.");
2087     }
2088 gross 854 Data tmp(right,getFunctionSpace(),false);
2089 jfenwick 2005 if (isLazy())
2090     {
2091     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),MUL); // for lazy * is equivalent to *=
2092     m_data=c->getPtr();
2093     return (*this);
2094     }
2095     else
2096     {
2097     binaryOp(tmp,multiplies<double>());
2098     return (*this);
2099     }
2100 jgs 94 }
2101    
2102 jgs 102 Data&
2103     Data::operator/=(const Data& right)
2104 jgs 94 {
2105 gross 783 if (isProtected()) {
2106     throw DataException("Error - attempt to update protected Data object.");
2107     }
2108 jfenwick 2005 if (isLazy() || right.isLazy())
2109     {
2110     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2111     m_data=c->getPtr();
2112     return (*this);
2113     }
2114     else
2115     {
2116     binaryOp(right,divides<double>());
2117     return (*this);
2118     }
2119 jgs 94 }
2120    
2121 jgs 102 Data&
2122     Data::operator/=(const boost::python::object& right)
2123 jgs 94 {
2124 jfenwick 2005 if (isProtected()) {
2125     throw DataException("Error - attempt to update protected Data object.");
2126     }
2127 gross 854 Data tmp(right,getFunctionSpace(),false);
2128 jfenwick 2005 if (isLazy())
2129     {
2130     DataLazy* c=new DataLazy(m_data,tmp.borrowDataPtr(),DIV); // for lazy / is equivalent to /=
2131     m_data=c->getPtr();
2132     return (*this);
2133     }
2134     else
2135     {
2136     binaryOp(tmp,divides<double>());
2137     return (*this);
2138     }
2139 jgs 94 }
2140    
2141 jgs 102 Data
2142 gross 699 Data::rpowO(const boost::python::object& left) const
2143     {
2144     Data left_d(left,*this);
2145     return left_d.powD(*this);
2146     }
2147    
2148     Data
2149 jgs 102 Data::powO(const boost::python::object& right) const
2150 jgs 94 {
2151 gross 854 Data tmp(right,getFunctionSpace(),false);
2152     return powD(tmp);
2153 jgs 94 }
2154    
2155 jgs 102 Data
2156     Data::powD(const Data& right) const
2157 jgs 94 {
2158 jfenwick 2005 if (isLazy() || right.isLazy())
2159     {
2160     DataLazy* c=new DataLazy(m_data,right.borrowDataPtr(),POW);
2161     return Data(c);
2162     }
2163 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2164 jgs 94 }
2165    
2166     //
2167 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2168 jgs 102 Data
2169     escript::operator+(const Data& left, const Data& right)
2170 jgs 94 {
2171 jfenwick 2005 if (left.isLazy() || right.isLazy())
2172     {
2173     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
2174     return Data(c);
2175     }
2176 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
2177 jgs 94 }
2178    
2179     //
2180 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2181 jgs 102 Data
2182     escript::operator-(const Data& left, const Data& right)
2183 jgs 94 {
2184 jfenwick 2005 if (left.isLazy() || right.isLazy())
2185     {
2186     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
2187     return Data(c);
2188     }
2189 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2190 jgs 94 }
2191    
2192     //
2193 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2194 jgs 102 Data
2195     escript::operator*(const Data& left, const Data& right)
2196 jgs 94 {
2197 jfenwick 2005 if (left.isLazy() || right.isLazy())
2198     {
2199     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
2200     return Data(c);
2201     }
2202 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2203 jgs 94 }
2204    
2205     //
2206 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2207 jgs 102 Data
2208     escript::operator/(const Data& left, const Data& right)
2209 jgs 94 {
2210 jfenwick 2005 if (left.isLazy() || right.isLazy())
2211     {
2212     DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);
2213     return Data(c);
2214     }
2215 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2216 jgs 94 }
2217    
2218     //
2219 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2220 jgs 102 Data
2221     escript::operator+(const Data& left, const boost::python::object& right)
2222 jgs 94 {
2223 jfenwick 2005 if (left.isLazy())
2224     {
2225     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),ADD);
2226     return Data(c);
2227     }
2228 gross 854 return left+Data(right,left.getFunctionSpace(),false);
2229 jgs 94 }
2230    
2231     //
2232 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2233 jgs 102 Data
2234     escript::operator-(const Data& left, const boost::python::object& right)
2235 jgs 94 {
2236 jfenwick 2005 if (left.isLazy())
2237     {
2238     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),SUB);
2239     return Data(c);
2240     }
2241 gross 854 return left-Data(right,left.getFunctionSpace(),false);
2242 jgs 94 }
2243    
2244     //
2245 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2246 jgs 102 Data
2247     escript::operator*(const Data& left, const boost::python::object& right)
2248 jgs 94 {
2249 jfenwick 2005 if (left.isLazy())
2250     {
2251     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),MUL);
2252     return Data(c);
2253     }
2254 gross 854 return left*Data(right,left.getFunctionSpace(),false);
2255 jgs 94 }
2256    
2257     //
2258 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2259 jgs 102 Data
2260     escript::operator/(const Data& left, const boost::python::object& right)
2261 jgs 94 {
2262 jfenwick 2005 if (left.isLazy())
2263     {
2264     DataLazy* c=new DataLazy(left.borrowDataPtr(),Data(right,left.getFunctionSpace(),false).borrowDataPtr(),DIV);
2265     return Data(c);
2266     }
2267 gross 854 return left/Data(right,left.getFunctionSpace(),false);
2268 jgs 94 }
2269    
2270     //
2271 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2272 jgs 102 Data
2273     escript::operator+(const boost::python::object& left, const Data& right)
2274 jgs 94 {
2275 jfenwick 2005 if (right.isLazy())
2276     {
2277     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),ADD);
2278     return Data(c);
2279     }
2280 gross 854 return Data(left,right.getFunctionSpace(),false)+right;
2281 jgs 94 }
2282    
2283     //
2284 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2285 jgs 102 Data
2286     escript::operator-(const boost::python::object& left, const Data& right)
2287 jgs 94 {
2288 jfenwick 2005 if (right.isLazy())
2289     {
2290     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),SUB);
2291     return Data(c);
2292     }
2293 gross 854 return Data(left,right.getFunctionSpace(),false)-right;
2294 jgs 94 }
2295    
2296     //
2297 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2298 jgs 102 Data
2299     escript::operator*(const boost::python::object& left, const Data& right)
2300 jgs 94 {
2301 jfenwick 2005 if (right.isLazy())
2302     {
2303     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),MUL);
2304     return Data(c);
2305     }
2306 gross 854 return Data(left,right.getFunctionSpace(),false)*right;
2307 jgs 94 }
2308    
2309     //
2310 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2311 jgs 102 Data
2312     escript::operator/(const boost::python::object& left, const Data& right)
2313 jgs 94 {
2314 jfenwick 2005 if (right.isLazy())
2315     {
2316     DataLazy* c=new DataLazy(Data(left,right.getFunctionSpace(),false).borrowDataPtr(),right.borrowDataPtr(),DIV);
2317     return Data(c);
2318     }
2319 gross 854 return Data(left,right.getFunctionSpace(),false)/right;
2320 jgs 94 }
2321    
2322 jgs 102
2323 bcumming 751 /* TODO */
2324     /* global reduction */
2325 jgs 102 Data
2326 ksteube 1312 Data::getItem(const boost::python::object& key) const
2327 jgs 94 {
2328    
2329 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2330 jgs 94
2331 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2332 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2333 jgs 94 }
2334    
2335 jgs 102 return getSlice(slice_region);
2336 jgs 94 }
2337    
2338 bcumming 751 /* TODO */
2339     /* global reduction */
2340 jgs 94 Data
2341 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2342 jgs 94 {
2343 jgs 102 return Data(*this,region);
2344 jgs 94 }
2345    
2346 bcumming 751 /* TODO */
2347     /* global reduction */
2348 jgs 94 void
2349 jgs 102 Data::setItemO(const boost::python::object& key,
2350     const boost::python::object& value)
2351 jgs 94 {
2352 jgs 102 Data tempData(value,getFunctionSpace());
2353     setItemD(key,tempData);
2354     }
2355    
2356     void
2357     Data::setItemD(const boost::python::object& key,
2358     const Data& value)
2359     {
2360 jfenwick 1796 // const DataArrayView& view=getPointDataView();
2361 jgs 123
2362 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2363     if (slice_region.size()!=getDataPointRank()) {
2364 jgs 94 throw DataException("Error - slice size does not match Data rank.");
2365     }
2366 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2367     setSlice(Data(value,getFunctionSpace()),slice_region);
2368     } else {
2369     setSlice(value,slice_region);
2370     }
2371 jgs 94 }
2372    
2373     void
2374 jgs 102 Data::setSlice(const Data& value,
2375 jfenwick 1796 const DataTypes::RegionType& region)
2376 jgs 94 {
2377 gross 783 if (isProtected()) {
2378     throw DataException("Error - attempt to update protected Data object.");
2379     }
2380 jfenwick 2005 FORCERESOLVE;
2381     /* if (isLazy())
2382     {
2383     throw DataException("Error - setSlice not permitted on lazy data.");
2384     }*/
2385 jgs 102 Data tempValue(value);
2386     typeMatchLeft(tempValue);
2387     typeMatchRight(tempValue);
2388 jfenwick 2005 getReady()->setSlice(tempValue.m_data.get(),region);
2389 jgs 102 }
2390    
2391     void
2392     Data::typeMatchLeft(Data& right) const
2393     {
2394 jfenwick 2005 if (right.isLazy() && !isLazy())
2395     {
2396     right.resolve();
2397     }
2398 jgs 102 if (isExpanded()){
2399     right.expand();
2400     } else if (isTagged()) {
2401     if (right.isConstant()) {
2402     right.tag();
2403     }
2404     }
2405     }
2406    
2407     void
2408     Data::typeMatchRight(const Data& right)
2409     {
2410 jfenwick 2005 if (isLazy() && !right.isLazy())
2411     {
2412     resolve();
2413     }
2414 jgs 94 if (isTagged()) {
2415     if (right.isExpanded()) {
2416     expand();
2417     }
2418     } else if (isConstant()) {
2419     if (right.isExpanded()) {
2420     expand();
2421     } else if (right.isTagged()) {
2422     tag();
2423     }
2424     }
2425     }
2426    
2427     void
2428 gross 1044 Data::setTaggedValueByName(std::string name,
2429 ksteube 1312 const boost::python::object& value)
2430 gross 1044 {
2431 jfenwick 1872 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2432 jfenwick 2005 FORCERESOLVE;
2433 jfenwick 1872 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2434 gross 1044 setTaggedValue(tagKey,value);
2435     }
2436     }
2437     void
2438 jgs 94 Data::setTaggedValue(int tagKey,
2439     const boost::python::object& value)
2440     {
2441 gross 783 if (isProtected()) {
2442     throw DataException("Error - attempt to update protected Data object.");
2443     }
2444 jgs 94 //
2445     // Ensure underlying data object is of type DataTagged
2446 jfenwick 2005 FORCERESOLVE;
2447 gross 1358 if (isConstant()) tag();
2448 matt 1319 numeric::array asNumArray(value);
2449 jgs 94
2450 matt 1319 // extract the shape of the numarray
2451 jfenwick 1796 DataTypes::ShapeType tempShape;
2452 matt 1319 for (int i=0; i < asNumArray.getrank(); i++) {
2453     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2454     }
2455    
2456 jfenwick 1796 DataVector temp_data2;
2457     temp_data2.copyFromNumArray(asNumArray);
2458    
2459 jfenwick 2005 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2460 jgs 94 }
2461    
2462 jfenwick 1796
2463 jgs 110 void
2464 jgs 121