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

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

Parent Directory Parent Directory | Revision Log Revision Log


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