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

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

Parent Directory Parent Directory | Revision Log Revision Log


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