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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2199 - (hide annotations)
Thu Jan 8 06:10:52 2009 UTC (11 years, 2 months ago) by jfenwick
File size: 81345 byte(s)
Misc fixes:
Added some svn:ignore properties for output files that were cluttering things up.

Lazy fixes:
Fixed shape calculations for TRACE and TRANSPOSE for rank>2.
Adjusted unit test accordingly.

As a Temporary change to DataC.cpp to test for lazy data in DataC's expanded check.
This is wrong but would only affect people using lazy data.
The proper fix will come when the numarray removal code moves over from the branch.

Made tensor product AUTOLAZY capable.
Fixed some bugs resolving tensor products (incorrect offsets in buffers).
Macro'd some stray couts.

- It appears that AUTOLAZY now passes all unit tests.
- It will not be _really_ safe for general use until I can add COW. 
- (Everything's better with COW)
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 jfenwick 2199 if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2297 jfenwick 2066 {
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);