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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1855 - (hide annotations)
Tue Oct 7 05:25:30 2008 UTC (11 years, 5 months ago) by jfenwick
File size: 75842 byte(s)
Modified Data::copyWithMask to have a less cryptic implementation.

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