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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1953 - (hide annotations)
Thu Oct 30 07:01:11 2008 UTC (11 years, 5 months ago) by phornby
File size: 75229 byte(s)
removed one signed/unsigned value comparison. The remaining one at line 863 is far more difficult to deal with.
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 1921
545     // OPENMP 3.0 allows unsigned loop vars.
546     #if defined(_OPENMP) && (_OPENMP < 200805)
547     long i;
548     #else
549 phornby 1918 size_t i;
550 phornby 1921 #endif
551 jfenwick 1855 #pragma omp parallel for private(i) schedule(static)
552     for (i=0;i<num_points;++i)
553     {
554     if (mvec[i]>0)
555     {
556     self[i]=ovec[i];
557     }
558     }
559     }
560 jgs 94
561    
562    
563     bool
564     Data::isExpanded() const
565     {
566     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
567     return (temp!=0);
568     }
569    
570     bool
571     Data::isTagged() const
572     {
573     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
574     return (temp!=0);
575     }
576    
577     bool
578     Data::isEmpty() const
579     {
580     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
581     return (temp!=0);
582     }
583    
584     bool
585     Data::isConstant() const
586     {
587     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
588     return (temp!=0);
589     }
590    
591     void
592 ksteube 1312 Data::setProtection()
593     {
594 gross 783 m_protected=true;
595     }
596    
597     bool
598 ksteube 1312 Data::isProtected() const
599     {
600 gross 783 return m_protected;
601     }
602    
603    
604    
605     void
606 jgs 94 Data::expand()
607     {
608     if (isConstant()) {
609     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
610     DataAbstract* temp=new DataExpanded(*tempDataConst);
611 jfenwick 1872 // shared_ptr<DataAbstract> temp_data(temp);
612     // m_data=temp_data;
613     m_data=temp->getPtr();
614 jgs 94 } else if (isTagged()) {
615     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
616     DataAbstract* temp=new DataExpanded(*tempDataTag);
617 jfenwick 1872 // shared_ptr<DataAbstract> temp_data(temp);
618     // m_data=temp_data;
619     m_data=temp->getPtr();
620 jgs 94 } else if (isExpanded()) {
621     //
622     // do nothing
623     } else if (isEmpty()) {
624     throw DataException("Error - Expansion of DataEmpty not possible.");
625     } else {
626     throw DataException("Error - Expansion not implemented for this Data type.");
627     }
628     }
629    
630     void
631     Data::tag()
632     {
633     if (isConstant()) {
634     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
635     DataAbstract* temp=new DataTagged(*tempDataConst);
636 jfenwick 1872 // shared_ptr<DataAbstract> temp_data(temp);
637     // m_data=temp_data;
638     m_data=temp->getPtr();
639 jgs 94 } else if (isTagged()) {
640     // do nothing
641     } else if (isExpanded()) {
642     throw DataException("Error - Creating tag data from DataExpanded not possible.");
643     } else if (isEmpty()) {
644     throw DataException("Error - Creating tag data from DataEmpty not possible.");
645     } else {
646     throw DataException("Error - Tagging not implemented for this Data type.");
647     }
648     }
649    
650 gross 854 Data
651     Data::oneOver() const
652 jgs 102 {
653 matt 1334 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
654 jgs 102 }
655    
656 jgs 94 Data
657 gross 698 Data::wherePositive() const
658 jgs 94 {
659 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
660 jgs 94 }
661    
662     Data
663 gross 698 Data::whereNegative() const
664 jgs 102 {
665 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
666 jgs 102 }
667    
668     Data
669 gross 698 Data::whereNonNegative() const
670 jgs 94 {
671 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
672 jgs 94 }
673    
674     Data
675 gross 698 Data::whereNonPositive() const
676 jgs 94 {
677 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
678 jgs 94 }
679    
680     Data
681 jgs 571 Data::whereZero(double tol) const
682 jgs 94 {
683 jgs 571 Data dataAbs=abs();
684 matt 1334 return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
685 jgs 94 }
686    
687     Data
688 jgs 571 Data::whereNonZero(double tol) const
689 jgs 102 {
690 jgs 571 Data dataAbs=abs();
691 matt 1334 return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
692 jgs 102 }
693    
694     Data
695 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
696     {
697     return Data(*this,functionspace);
698     }
699    
700     bool
701     Data::probeInterpolation(const FunctionSpace& functionspace) const
702     {
703     if (getFunctionSpace()==functionspace) {
704     return true;
705     } else {
706 jfenwick 1872 const_Domain_ptr domain=getDomain();
707     if (*domain==*functionspace.getDomain()) {
708     return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
709 jgs 94 } else {
710 jfenwick 1872 return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
711 jgs 94 }
712     }
713     }
714    
715     Data
716     Data::gradOn(const FunctionSpace& functionspace) const
717     {
718 jfenwick 1803 if (isEmpty())
719     {
720     throw DataException("Error - operation not permitted on instances of DataEmpty.");
721     }
722 matt 1353 double blocktimer_start = blocktimer_time();
723 jgs 94 if (functionspace.getDomain()!=getDomain())
724     throw DataException("Error - gradient cannot be calculated on different domains.");
725 jfenwick 1796 DataTypes::ShapeType grad_shape=getDataPointShape();
726 jgs 94 grad_shape.push_back(functionspace.getDim());
727     Data out(0.0,grad_shape,functionspace,true);
728 jfenwick 1872 getDomain()->setToGradient(out,*this);
729 matt 1353 blocktimer_increment("grad()", blocktimer_start);
730 jgs 94 return out;
731     }
732    
733     Data
734     Data::grad() const
735     {
736 jfenwick 1803 if (isEmpty())
737     {
738     throw DataException("Error - operation not permitted on instances of DataEmpty.");
739     }
740 jfenwick 1872 return gradOn(escript::function(*getDomain()));
741 jgs 94 }
742    
743     int
744     Data::getDataPointSize() const
745     {
746 jfenwick 1796 return m_data->getNoValues();
747 jgs 94 }
748    
749 jfenwick 1796 DataTypes::ValueType::size_type
750 jgs 94 Data::getLength() const
751     {
752     return m_data->getLength();
753     }
754    
755 ksteube 1312 const
756 jgs 121 boost::python::numeric::array
757 ksteube 1312 Data:: getValueOfDataPoint(int dataPointNo)
758 jgs 121 {
759 gross 921 int i, j, k, l;
760 jgs 121 //
761     // determine the rank and shape of each data point
762     int dataPointRank = getDataPointRank();
763 jfenwick 1796 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
764 jgs 121
765     //
766     // create the numeric array to be returned
767     boost::python::numeric::array numArray(0.0);
768    
769     //
770 gross 921 // the shape of the returned numeric array will be the same
771     // as that of the data point
772     int arrayRank = dataPointRank;
773 jfenwick 1796 const DataTypes::ShapeType& arrayShape = dataPointShape;
774 jgs 121
775     //
776     // resize the numeric array to the shape just calculated
777 gross 921 if (arrayRank==0) {
778     numArray.resize(1);
779     }
780 jgs 121 if (arrayRank==1) {
781     numArray.resize(arrayShape[0]);
782     }
783     if (arrayRank==2) {
784     numArray.resize(arrayShape[0],arrayShape[1]);
785     }
786     if (arrayRank==3) {
787     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
788     }
789     if (arrayRank==4) {
790     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
791     }
792    
793 gross 921 if (getNumDataPointsPerSample()>0) {
794     int sampleNo = dataPointNo/getNumDataPointsPerSample();
795     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
796     //
797     // Check a valid sample number has been supplied
798 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
799 gross 921 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
800     }
801 ksteube 1312
802 gross 921 //
803     // Check a valid data point number has been supplied
804 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
805 gross 921 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
806     }
807     // TODO: global error handling
808     // create a view of the data if it is stored locally
809 jfenwick 1796 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
810     DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
811 ksteube 1312
812 jfenwick 1796
813 gross 921 switch( dataPointRank ){
814     case 0 :
815 jfenwick 1796 numArray[0] = getDataAtOffset(offset);
816 gross 921 break;
817 ksteube 1312 case 1 :
818 gross 921 for( i=0; i<dataPointShape[0]; i++ )
819 jfenwick 1796 numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
820 gross 921 break;
821 ksteube 1312 case 2 :
822 gross 921 for( i=0; i<dataPointShape[0]; i++ )
823     for( j=0; j<dataPointShape[1]; j++)
824 jfenwick 1796 numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
825 gross 921 break;
826 ksteube 1312 case 3 :
827 gross 921 for( i=0; i<dataPointShape[0]; i++ )
828     for( j=0; j<dataPointShape[1]; j++ )
829     for( k=0; k<dataPointShape[2]; k++)
830 jfenwick 1796 numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
831 gross 921 break;
832     case 4 :
833     for( i=0; i<dataPointShape[0]; i++ )
834     for( j=0; j<dataPointShape[1]; j++ )
835     for( k=0; k<dataPointShape[2]; k++ )
836     for( l=0; l<dataPointShape[3]; l++)
837 jfenwick 1796 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
838 gross 921 break;
839     }
840 jgs 117 }
841     //
842 gross 921 // return the array
843 jgs 117 return numArray;
844 gross 921
845 jgs 117 }
846 jfenwick 1796
847 gross 921 void
848 ksteube 1312 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
849 jgs 121 {
850 gross 1034 // this will throw if the value cannot be represented
851     boost::python::numeric::array num_array(py_object);
852     setValueOfDataPointToArray(dataPointNo,num_array);
853     }
854    
855     void
856     Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
857     {
858 gross 921 if (isProtected()) {
859     throw DataException("Error - attempt to update protected Data object.");
860     }
861     //
862     // check rank
863 ksteube 1312 if (num_array.getrank()<getDataPointRank())
864 gross 921 throw DataException("Rank of numarray does not match Data object rank");
865 bcumming 790
866 jgs 121 //
867 gross 921 // check shape of num_array
868 phornby 1953 for (unsigned int i=0; i<getDataPointRank(); i++) {
869 gross 921 if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
870     throw DataException("Shape of numarray does not match Data object rank");
871 jgs 121 }
872     //
873 gross 921 // make sure data is expanded:
874 gross 1859 //
875 gross 921 if (!isExpanded()) {
876     expand();
877 jgs 121 }
878 gross 921 if (getNumDataPointsPerSample()>0) {
879     int sampleNo = dataPointNo/getNumDataPointsPerSample();
880     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
881     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
882     } else {
883     m_data->copyToDataPoint(-1, 0,num_array);
884     }
885     }
886 jgs 121
887 gross 922 void
888     Data::setValueOfDataPoint(int dataPointNo, const double value)
889     {
890     if (isProtected()) {
891     throw DataException("Error - attempt to update protected Data object.");
892     }
893     //
894     // make sure data is expanded:
895     if (!isExpanded()) {
896     expand();
897     }
898     if (getNumDataPointsPerSample()>0) {
899     int sampleNo = dataPointNo/getNumDataPointsPerSample();
900     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
901     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
902     } else {
903     m_data->copyToDataPoint(-1, 0,value);
904     }
905     }
906    
907 ksteube 1312 const
908 gross 921 boost::python::numeric::array
909 ksteube 1312 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
910 gross 921 {
911     size_t length=0;
912     int i, j, k, l, pos;
913 jgs 121 //
914     // determine the rank and shape of each data point
915     int dataPointRank = getDataPointRank();
916 jfenwick 1796 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
917 jgs 121
918     //
919     // create the numeric array to be returned
920     boost::python::numeric::array numArray(0.0);
921    
922     //
923     // the shape of the returned numeric array will be the same
924     // as that of the data point
925     int arrayRank = dataPointRank;
926 jfenwick 1796 const DataTypes::ShapeType& arrayShape = dataPointShape;
927 jgs 121
928     //
929     // resize the numeric array to the shape just calculated
930     if (arrayRank==0) {
931     numArray.resize(1);
932     }
933     if (arrayRank==1) {
934     numArray.resize(arrayShape[0]);
935     }
936     if (arrayRank==2) {
937     numArray.resize(arrayShape[0],arrayShape[1]);
938     }
939     if (arrayRank==3) {
940     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
941     }
942     if (arrayRank==4) {
943     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
944     }
945    
946 gross 921 // added for the MPI communication
947     length=1;
948     for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
949     double *tmpData = new double[length];
950 bcumming 790
951 jgs 121 //
952     // load the values for the data point into the numeric array.
953 bcumming 790
954     // updated for the MPI case
955     if( get_MPIRank()==procNo ){
956 gross 921 if (getNumDataPointsPerSample()>0) {
957     int sampleNo = dataPointNo/getNumDataPointsPerSample();
958     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
959     //
960     // Check a valid sample number has been supplied
961 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
962 gross 921 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
963     }
964 ksteube 1312
965 gross 921 //
966     // Check a valid data point number has been supplied
967 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
968 gross 921 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
969     }
970     // TODO: global error handling
971 bcumming 790 // create a view of the data if it is stored locally
972 jfenwick 1796 //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
973     DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
974 ksteube 1312
975 bcumming 790 // pack the data from the view into tmpData for MPI communication
976     pos=0;
977     switch( dataPointRank ){
978     case 0 :
979 jfenwick 1796 tmpData[0] = getDataAtOffset(offset);
980 bcumming 790 break;
981 ksteube 1312 case 1 :
982 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
983 jfenwick 1796 tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
984 bcumming 790 break;
985 ksteube 1312 case 2 :
986 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
987     for( j=0; j<dataPointShape[1]; j++, pos++ )
988 jfenwick 1796 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
989 bcumming 790 break;
990 ksteube 1312 case 3 :
991 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
992     for( j=0; j<dataPointShape[1]; j++ )
993     for( k=0; k<dataPointShape[2]; k++, pos++ )
994 jfenwick 1796 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
995 bcumming 790 break;
996     case 4 :
997     for( i=0; i<dataPointShape[0]; i++ )
998     for( j=0; j<dataPointShape[1]; j++ )
999     for( k=0; k<dataPointShape[2]; k++ )
1000     for( l=0; l<dataPointShape[3]; l++, pos++ )
1001 jfenwick 1796 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
1002 bcumming 790 break;
1003     }
1004 gross 921 }
1005 bcumming 790 }
1006 ksteube 1312 #ifdef PASO_MPI
1007 gross 921 // broadcast the data to all other processes
1008     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1009     #endif
1010 bcumming 790
1011     // unpack the data
1012     switch( dataPointRank ){
1013     case 0 :
1014 gross 921 numArray[0]=tmpData[0];
1015 bcumming 790 break;
1016 ksteube 1312 case 1 :
1017 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
1018     numArray[i]=tmpData[i];
1019     break;
1020 ksteube 1312 case 2 :
1021 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
1022     for( j=0; j<dataPointShape[1]; j++ )
1023 gross 921 numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1024 bcumming 790 break;
1025 ksteube 1312 case 3 :
1026 bcumming 790 for( i=0; i<dataPointShape[0]; i++ )
1027     for( j=0; j<dataPointShape[1]; j++ )
1028     for( k=0; k<dataPointShape[2]; k++ )
1029 gross 921 numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1030 bcumming 790 break;
1031     case 4 :
1032     for( i=0; i<dataPointShape[0]; i++ )
1033     for( j=0; j<dataPointShape[1]; j++ )
1034     for( k=0; k<dataPointShape[2]; k++ )
1035     for( l=0; l<dataPointShape[3]; l++ )
1036 gross 921 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1037 bcumming 790 break;
1038     }
1039    
1040 ksteube 1312 delete [] tmpData;
1041 jgs 121 //
1042     // return the loaded array
1043     return numArray;
1044     }
1045    
1046 gross 921
1047    
1048 jgs 121 boost::python::numeric::array
1049 jgs 94 Data::integrate() const
1050     {
1051     int index;
1052     int rank = getDataPointRank();
1053 jfenwick 1796 DataTypes::ShapeType shape = getDataPointShape();
1054 ksteube 1312 int dataPointSize = getDataPointSize();
1055 jgs 94
1056     //
1057     // calculate the integral values
1058 ksteube 1312 vector<double> integrals(dataPointSize);
1059     vector<double> integrals_local(dataPointSize);
1060     #ifdef PASO_MPI
1061 jfenwick 1873 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);
1062 ksteube 1312 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1063     double *tmp = new double[dataPointSize];
1064     double *tmp_local = new double[dataPointSize];
1065     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1066     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1067     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1068     delete[] tmp;
1069     delete[] tmp_local;
1070     #else
1071 jfenwick 1872 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1072 ksteube 1312 #endif
1073 jgs 94
1074     //
1075     // create the numeric array to be returned
1076     // and load the array with the integral values
1077     boost::python::numeric::array bp_array(1.0);
1078     if (rank==0) {
1079 jgs 108 bp_array.resize(1);
1080 jgs 94 index = 0;
1081     bp_array[0] = integrals[index];
1082     }
1083     if (rank==1) {
1084     bp_array.resize(shape[0]);
1085     for (int i=0; i<shape[0]; i++) {
1086     index = i;
1087     bp_array[i] = integrals[index];
1088     }
1089     }
1090     if (rank==2) {
1091 gross 436 bp_array.resize(shape[0],shape[1]);
1092     for (int i=0; i<shape[0]; i++) {
1093     for (int j=0; j<shape[1]; j++) {
1094     index = i + shape[0] * j;
1095     bp_array[make_tuple(i,j)] = integrals[index];
1096     }
1097     }
1098 jgs 94 }
1099     if (rank==3) {
1100     bp_array.resize(shape[0],shape[1],shape[2]);
1101     for (int i=0; i<shape[0]; i++) {
1102     for (int j=0; j<shape[1]; j++) {
1103     for (int k=0; k<shape[2]; k++) {
1104     index = i + shape[0] * ( j + shape[1] * k );
1105 gross 436 bp_array[make_tuple(i,j,k)] = integrals[index];
1106 jgs 94 }
1107     }
1108     }
1109     }
1110     if (rank==4) {
1111     bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1112     for (int i=0; i<shape[0]; i++) {
1113     for (int j=0; j<shape[1]; j++) {
1114     for (int k=0; k<shape[2]; k++) {
1115     for (int l=0; l<shape[3]; l++) {
1116     index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1117 gross 436 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1118 jgs 94 }
1119     }
1120     }
1121     }
1122     }
1123    
1124     //
1125     // return the loaded array
1126     return bp_array;
1127     }
1128    
1129     Data
1130     Data::sin() const
1131     {
1132 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1133 jgs 94 }
1134    
1135     Data
1136     Data::cos() const
1137     {
1138 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1139 jgs 94 }
1140    
1141     Data
1142     Data::tan() const
1143     {
1144 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1145 jgs 94 }
1146    
1147     Data
1148 jgs 150 Data::asin() const
1149     {
1150 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1151 jgs 150 }
1152    
1153     Data
1154     Data::acos() const
1155     {
1156 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1157 jgs 150 }
1158    
1159 phornby 1026
1160 jgs 150 Data
1161     Data::atan() const
1162     {
1163 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1164 jgs 150 }
1165    
1166     Data
1167     Data::sinh() const
1168     {
1169 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1170 matt 1334
1171 jgs 150 }
1172    
1173     Data
1174     Data::cosh() const
1175     {
1176 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1177 jgs 150 }
1178    
1179     Data
1180     Data::tanh() const
1181     {
1182 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1183 jgs 150 }
1184    
1185 phornby 1026
1186 jgs 150 Data
1187 phornby 1026 Data::erf() const
1188     {
1189 gross 1028 #ifdef _WIN32
1190     throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1191     #else
1192 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
1193 phornby 1026 #endif
1194     }
1195    
1196     Data
1197 jgs 150 Data::asinh() const
1198     {
1199 phornby 1032 #ifdef _WIN32
1200 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1201 phornby 1032 #else
1202 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
1203 phornby 1032 #endif
1204 jgs 150 }
1205    
1206     Data
1207     Data::acosh() const
1208     {
1209 phornby 1032 #ifdef _WIN32
1210 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1211 phornby 1032 #else
1212 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
1213 phornby 1032 #endif
1214 jgs 150 }
1215    
1216     Data
1217     Data::atanh() const
1218     {
1219 phornby 1032 #ifdef _WIN32
1220 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1221 phornby 1032 #else
1222 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
1223 phornby 1032 #endif
1224 jgs 150 }
1225    
1226     Data
1227 gross 286 Data::log10() const
1228 jgs 94 {
1229 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1230 jgs 94 }
1231    
1232     Data
1233 gross 286 Data::log() const
1234 jgs 94 {
1235 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1236 jgs 94 }
1237    
1238 jgs 106 Data
1239     Data::sign() const
1240 jgs 94 {
1241 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1242 jgs 94 }
1243    
1244 jgs 106 Data
1245     Data::abs() const
1246 jgs 94 {
1247 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1248 jgs 94 }
1249    
1250 jgs 106 Data
1251     Data::neg() const
1252 jgs 94 {
1253 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1254 jgs 94 }
1255    
1256 jgs 102 Data
1257 jgs 106 Data::pos() const
1258 jgs 94 {
1259 jgs 148 Data result;
1260     // perform a deep copy
1261     result.copy(*this);
1262     return result;
1263 jgs 102 }
1264    
1265     Data
1266 jgs 106 Data::exp() const
1267 jgs 102 {
1268 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1269 jgs 102 }
1270    
1271     Data
1272 jgs 106 Data::sqrt() const
1273 jgs 102 {
1274 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1275 jgs 102 }
1276    
1277 jgs 106 double
1278     Data::Lsup() const
1279 jgs 102 {
1280 phornby 1455 double localValue;
1281 jgs 106 //
1282     // set the initial absolute maximum value to zero
1283 bcumming 751
1284 jgs 147 AbsMax abs_max_func;
1285 bcumming 751 localValue = algorithm(abs_max_func,0);
1286     #ifdef PASO_MPI
1287 phornby 1455 double globalValue;
1288 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1289     return globalValue;
1290     #else
1291     return localValue;
1292     #endif
1293 jgs 117 }
1294    
1295     double
1296 jgs 106 Data::sup() const
1297 jgs 102 {
1298 phornby 1455 double localValue;
1299 jgs 106 //
1300     // set the initial maximum value to min possible double
1301 jgs 147 FMax fmax_func;
1302 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1303     #ifdef PASO_MPI
1304 phornby 1455 double globalValue;
1305 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1306     return globalValue;
1307     #else
1308     return localValue;
1309     #endif
1310 jgs 102 }
1311    
1312 jgs 106 double
1313     Data::inf() const
1314 jgs 102 {
1315 phornby 1455 double localValue;
1316 jgs 106 //
1317     // set the initial minimum value to max possible double
1318 jgs 147 FMin fmin_func;
1319 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1320     #ifdef PASO_MPI
1321 phornby 1455 double globalValue;
1322 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1323     return globalValue;
1324     #else
1325     return localValue;
1326     #endif
1327 jgs 102 }
1328    
1329 bcumming 751 /* TODO */
1330     /* global reduction */
1331 jgs 102 Data
1332 jgs 106 Data::maxval() const
1333 jgs 102 {
1334 jgs 113 //
1335     // set the initial maximum value to min possible double
1336 jgs 147 FMax fmax_func;
1337     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1338 jgs 102 }
1339    
1340     Data
1341 jgs 106 Data::minval() const
1342 jgs 102 {
1343 jgs 113 //
1344     // set the initial minimum value to max possible double
1345 jgs 147 FMin fmin_func;
1346     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1347 jgs 102 }
1348    
1349 jgs 123 Data
1350 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1351 jgs 123 {
1352 gross 804 int axis0_tmp,axis1_tmp;
1353 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1354     DataTypes::ShapeType ev_shape;
1355 gross 800 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1356     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1357     int rank=getDataPointRank();
1358 gross 804 if (rank<2) {
1359     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1360 gross 800 }
1361 gross 804 if (axis0<0 || axis0>rank-1) {
1362     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1363     }
1364     if (axis1<0 || axis1>rank-1) {
1365     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1366     }
1367     if (axis0 == axis1) {
1368     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1369     }
1370     if (axis0 > axis1) {
1371     axis0_tmp=axis1;
1372     axis1_tmp=axis0;
1373     } else {
1374     axis0_tmp=axis0;
1375     axis1_tmp=axis1;
1376     }
1377 gross 800 for (int i=0; i<rank; i++) {
1378 gross 804 if (i == axis0_tmp) {
1379 ksteube 1312 ev_shape.push_back(s[axis1_tmp]);
1380 gross 804 } else if (i == axis1_tmp) {
1381 ksteube 1312 ev_shape.push_back(s[axis0_tmp]);
1382 gross 800 } else {
1383 ksteube 1312 ev_shape.push_back(s[i]);
1384 gross 800 }
1385     }
1386     Data ev(0.,ev_shape,getFunctionSpace());
1387     ev.typeMatchRight(*this);
1388 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1389 gross 800 return ev;
1390    
1391 jgs 123 }
1392    
1393     Data
1394 ksteube 775 Data::symmetric() const
1395 jgs 123 {
1396 ksteube 775 // check input
1397 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1398 ksteube 775 if (getDataPointRank()==2) {
1399 ksteube 1312 if(s[0] != s[1])
1400 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1401     }
1402     else if (getDataPointRank()==4) {
1403     if(!(s[0] == s[2] && s[1] == s[3]))
1404     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1405     }
1406     else {
1407     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1408     }
1409     Data ev(0.,getDataPointShape(),getFunctionSpace());
1410     ev.typeMatchRight(*this);
1411     m_data->symmetric(ev.m_data.get());
1412     return ev;
1413     }
1414    
1415     Data
1416     Data::nonsymmetric() const
1417     {
1418     // check input
1419 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1420 ksteube 775 if (getDataPointRank()==2) {
1421 ksteube 1312 if(s[0] != s[1])
1422 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1423 jfenwick 1796 DataTypes::ShapeType ev_shape;
1424 ksteube 775 ev_shape.push_back(s[0]);
1425     ev_shape.push_back(s[1]);
1426     Data ev(0.,ev_shape,getFunctionSpace());
1427     ev.typeMatchRight(*this);
1428     m_data->nonsymmetric(ev.m_data.get());
1429     return ev;
1430     }
1431     else if (getDataPointRank()==4) {
1432     if(!(s[0] == s[2] && s[1] == s[3]))
1433     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1434 jfenwick 1796 DataTypes::ShapeType ev_shape;
1435 ksteube 775 ev_shape.push_back(s[0]);
1436     ev_shape.push_back(s[1]);
1437     ev_shape.push_back(s[2]);
1438     ev_shape.push_back(s[3]);
1439     Data ev(0.,ev_shape,getFunctionSpace());
1440     ev.typeMatchRight(*this);
1441     m_data->nonsymmetric(ev.m_data.get());
1442     return ev;
1443     }
1444     else {
1445     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1446     }
1447     }
1448    
1449     Data
1450 gross 800 Data::trace(int axis_offset) const
1451 ksteube 775 {
1452 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1453 ksteube 775 if (getDataPointRank()==2) {
1454 jfenwick 1796 DataTypes::ShapeType ev_shape;
1455 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1456     ev.typeMatchRight(*this);
1457 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1458 ksteube 775 return ev;
1459     }
1460     if (getDataPointRank()==3) {
1461 jfenwick 1796 DataTypes::ShapeType ev_shape;
1462 ksteube 775 if (axis_offset==0) {
1463     int s2=s[2];
1464     ev_shape.push_back(s2);
1465     }
1466     else if (axis_offset==1) {
1467     int s0=s[0];
1468     ev_shape.push_back(s0);
1469     }
1470     Data ev(0.,ev_shape,getFunctionSpace());
1471     ev.typeMatchRight(*this);
1472 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1473 ksteube 775 return ev;
1474     }
1475     if (getDataPointRank()==4) {
1476 jfenwick 1796 DataTypes::ShapeType ev_shape;
1477 ksteube 775 if (axis_offset==0) {
1478     ev_shape.push_back(s[2]);
1479     ev_shape.push_back(s[3]);
1480     }
1481     else if (axis_offset==1) {
1482     ev_shape.push_back(s[0]);
1483     ev_shape.push_back(s[3]);
1484     }
1485     else if (axis_offset==2) {
1486     ev_shape.push_back(s[0]);
1487     ev_shape.push_back(s[1]);
1488     }
1489     Data ev(0.,ev_shape,getFunctionSpace());
1490     ev.typeMatchRight(*this);
1491 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1492 ksteube 775 return ev;
1493     }
1494     else {
1495 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1496 ksteube 775 }
1497     }
1498    
1499     Data
1500     Data::transpose(int axis_offset) const
1501     {
1502 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1503     DataTypes::ShapeType ev_shape;
1504 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1505     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1506     int rank=getDataPointRank();
1507     if (axis_offset<0 || axis_offset>rank) {
1508     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1509     }
1510     for (int i=0; i<rank; i++) {
1511     int index = (axis_offset+i)%rank;
1512     ev_shape.push_back(s[index]); // Append to new shape
1513     }
1514     Data ev(0.,ev_shape,getFunctionSpace());
1515     ev.typeMatchRight(*this);
1516     m_data->transpose(ev.m_data.get(), axis_offset);
1517     return ev;
1518 jgs 123 }
1519    
1520 gross 576 Data
1521     Data::eigenvalues() const
1522     {
1523     // check input
1524 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1525 ksteube 1312 if (getDataPointRank()!=2)
1526 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1527 ksteube 1312 if(s[0] != s[1])
1528 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1529     // create return
1530 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1531 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1532     ev.typeMatchRight(*this);
1533     m_data->eigenvalues(ev.m_data.get());
1534     return ev;
1535     }
1536    
1537 jgs 121 const boost::python::tuple
1538 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1539     {
1540 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1541 ksteube 1312 if (getDataPointRank()!=2)
1542 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1543 ksteube 1312 if(s[0] != s[1])
1544 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1545     // create return
1546 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1547 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1548     ev.typeMatchRight(*this);
1549 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1550 gross 576 Data V(0.,V_shape,getFunctionSpace());
1551     V.typeMatchRight(*this);
1552     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1553     return make_tuple(boost::python::object(ev),boost::python::object(V));
1554     }
1555    
1556     const boost::python::tuple
1557 gross 921 Data::minGlobalDataPoint() const
1558 jgs 121 {
1559 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1560 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1561     // surrounding function
1562    
1563     int DataPointNo;
1564 gross 921 int ProcNo;
1565     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1566     return make_tuple(ProcNo,DataPointNo);
1567 jgs 148 }
1568    
1569     void
1570 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1571     int& DataPointNo) const
1572 jgs 148 {
1573     int i,j;
1574     int lowi=0,lowj=0;
1575     double min=numeric_limits<double>::max();
1576    
1577 jgs 121 Data temp=minval();
1578    
1579     int numSamples=temp.getNumSamples();
1580     int numDPPSample=temp.getNumDataPointsPerSample();
1581    
1582 jgs 148 double next,local_min;
1583 jfenwick 1946 int local_lowi=0,local_lowj=0;
1584 jgs 121
1585 jgs 148 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1586     {
1587     local_min=min;
1588     #pragma omp for private(i,j) schedule(static)
1589     for (i=0; i<numSamples; i++) {
1590     for (j=0; j<numDPPSample; j++) {
1591 jfenwick 1796 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1592 jgs 148 if (next<local_min) {
1593     local_min=next;
1594     local_lowi=i;
1595     local_lowj=j;
1596     }
1597 jgs 121 }
1598     }
1599 jgs 148 #pragma omp critical
1600     if (local_min<min) {
1601     min=local_min;
1602     lowi=local_lowi;
1603     lowj=local_lowj;
1604     }
1605 jgs 121 }
1606    
1607 bcumming 782 #ifdef PASO_MPI
1608     // determine the processor on which the minimum occurs
1609 jfenwick 1796 next = temp.getDataPoint(lowi,lowj);
1610 bcumming 782 int lowProc = 0;
1611     double *globalMins = new double[get_MPISize()+1];
1612     int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1613 ksteube 1312
1614 bcumming 782 if( get_MPIRank()==0 ){
1615     next = globalMins[lowProc];
1616     for( i=1; i<get_MPISize(); i++ )
1617     if( next>globalMins[i] ){
1618     lowProc = i;
1619     next = globalMins[i];
1620     }
1621     }
1622     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1623    
1624     delete [] globalMins;
1625     ProcNo = lowProc;
1626 bcumming 790 #else
1627     ProcNo = 0;
1628 bcumming 782 #endif
1629 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1630 jgs 121 }
1631    
1632 jgs 104 void
1633     Data::saveDX(std::string fileName) const
1634     {
1635 jfenwick 1803 if (isEmpty())
1636     {
1637     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1638     }
1639 jgs 153 boost::python::dict args;
1640     args["data"]=boost::python::object(this);
1641 jfenwick 1872 getDomain()->saveDX(fileName,args);
1642 jgs 104 return;
1643     }
1644    
1645 jgs 110 void
1646     Data::saveVTK(std::string fileName) const
1647     {
1648 jfenwick 1803 if (isEmpty())
1649     {
1650     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1651     }
1652 jgs 153 boost::python::dict args;
1653     args["data"]=boost::python::object(this);
1654 jfenwick 1872 getDomain()->saveVTK(fileName,args);
1655 jgs 110 return;
1656     }
1657    
1658 jgs 102 Data&
1659     Data::operator+=(const Data& right)
1660     {
1661 gross 783 if (isProtected()) {
1662     throw DataException("Error - attempt to update protected Data object.");
1663     }
1664 jgs 94 binaryOp(right,plus<double>());
1665     return (*this);
1666     }
1667    
1668 jgs 102 Data&
1669     Data::operator+=(const boost::python::object& right)
1670 jgs 94 {
1671 gross 854 Data tmp(right,getFunctionSpace(),false);
1672     binaryOp(tmp,plus<double>());
1673 jgs 94 return (*this);
1674     }
1675 ksteube 1312 Data&
1676     Data::operator=(const Data& other)
1677     {
1678     copy(other);
1679     return (*this);
1680     }
1681 jgs 94
1682 jgs 102 Data&
1683     Data::operator-=(const Data& right)
1684 jgs 94 {
1685 gross 783 if (isProtected()) {
1686     throw DataException("Error - attempt to update protected Data object.");
1687     }
1688 jgs 94 binaryOp(right,minus<double>());
1689     return (*this);
1690     }
1691    
1692 jgs 102 Data&
1693     Data::operator-=(const boost::python::object& right)
1694 jgs 94 {
1695 gross 854 Data tmp(right,getFunctionSpace(),false);
1696     binaryOp(tmp,minus<double>());
1697 jgs 94 return (*this);
1698     }
1699    
1700 jgs 102 Data&
1701     Data::operator*=(const Data& right)
1702 jgs 94 {
1703 gross 783 if (isProtected()) {
1704     throw DataException("Error - attempt to update protected Data object.");
1705     }
1706 jgs 94 binaryOp(right,multiplies<double>());
1707     return (*this);
1708     }
1709    
1710 jgs 102 Data&
1711     Data::operator*=(const boost::python::object& right)
1712 jgs 94 {
1713 gross 854 Data tmp(right,getFunctionSpace(),false);
1714     binaryOp(tmp,multiplies<double>());
1715 jgs 94 return (*this);
1716     }
1717    
1718 jgs 102 Data&
1719     Data::operator/=(const Data& right)
1720 jgs 94 {
1721 gross 783 if (isProtected()) {
1722     throw DataException("Error - attempt to update protected Data object.");
1723     }
1724 jgs 94 binaryOp(right,divides<double>());
1725     return (*this);
1726     }
1727    
1728 jgs 102 Data&
1729     Data::operator/=(const boost::python::object& right)
1730 jgs 94 {
1731 gross 854 Data tmp(right,getFunctionSpace(),false);
1732     binaryOp(tmp,divides<double>());
1733 jgs 94 return (*this);
1734     }
1735    
1736 jgs 102 Data
1737 gross 699 Data::rpowO(const boost::python::object& left) const
1738     {
1739     Data left_d(left,*this);
1740     return left_d.powD(*this);
1741     }
1742    
1743     Data
1744 jgs 102 Data::powO(const boost::python::object& right) const
1745 jgs 94 {
1746 gross 854 Data tmp(right,getFunctionSpace(),false);
1747     return powD(tmp);
1748 jgs 94 }
1749    
1750 jgs 102 Data
1751     Data::powD(const Data& right) const
1752 jgs 94 {
1753 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1754 jgs 94 }
1755    
1756     //
1757 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
1758 jgs 102 Data
1759     escript::operator+(const Data& left, const Data& right)
1760 jgs 94 {
1761 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
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, minus<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, multiplies<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, divides<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 boost::python::object& right)
1792 jgs 94 {
1793 gross 854 return left+Data(right,left.getFunctionSpace(),false);
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 boost::python::object& left, const Data& right)
1824 jgs 94 {
1825 gross 854 return Data(left,right.getFunctionSpace(),false)+right;
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 102 //bool escript::operator==(const Data& left, const Data& right)
1854     //{
1855     // /*
1856     // NB: this operator does very little at this point, and isn't to
1857     // be relied on. Requires further implementation.
1858     // */
1859     //
1860     // bool ret;
1861     //
1862     // if (left.isEmpty()) {
1863     // if(!right.isEmpty()) {
1864     // ret = false;
1865     // } else {
1866     // ret = true;
1867     // }
1868     // }
1869     //
1870     // if (left.isConstant()) {
1871     // if(!right.isConstant()) {
1872     // ret = false;
1873     // } else {
1874     // ret = true;
1875     // }
1876     // }
1877     //
1878     // if (left.isTagged()) {
1879     // if(!right.isTagged()) {
1880     // ret = false;
1881     // } else {
1882     // ret = true;
1883     // }
1884     // }
1885     //
1886     // if (left.isExpanded()) {
1887     // if(!right.isExpanded()) {
1888     // ret = false;
1889     // } else {
1890     // ret = true;
1891     // }
1892     // }
1893     //
1894     // return ret;
1895     //}
1896    
1897 bcumming 751 /* TODO */
1898     /* global reduction */
1899 jgs 102 Data
1900 ksteube 1312 Data::getItem(const boost::python::object& key) const
1901 jgs 94 {
1902    
1903 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1904 jgs 94
1905 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
1906 jgs 102 throw DataException("Error - slice size does not match Data rank.");
1907 jgs 94 }
1908    
1909 jgs 102 return getSlice(slice_region);
1910 jgs 94 }
1911    
1912 bcumming 751 /* TODO */
1913     /* global reduction */
1914 jgs 94 Data
1915 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
1916 jgs 94 {
1917 jgs 102 return Data(*this,region);
1918 jgs 94 }
1919    
1920 bcumming 751 /* TODO */
1921     /* global reduction */
1922 jgs 94 void
1923 jgs 102 Data::setItemO(const boost::python::object& key,
1924     const boost::python::object& value)
1925 jgs 94 {
1926 jgs 102 Data tempData(value,getFunctionSpace());
1927     setItemD(key,tempData);
1928     }
1929    
1930     void
1931     Data::setItemD(const boost::python::object& key,
1932     const Data& value)
1933     {
1934 jfenwick 1796 // const DataArrayView& view=getPointDataView();
1935 jgs 123
1936 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1937     if (slice_region.size()!=getDataPointRank()) {
1938 jgs 94 throw DataException("Error - slice size does not match Data rank.");
1939     }
1940 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
1941     setSlice(Data(value,getFunctionSpace()),slice_region);
1942     } else {
1943     setSlice(value,slice_region);
1944     }
1945 jgs 94 }
1946    
1947     void
1948 jgs 102 Data::setSlice(const Data& value,
1949 jfenwick 1796 const DataTypes::RegionType& region)
1950 jgs 94 {
1951 gross 783 if (isProtected()) {
1952     throw DataException("Error - attempt to update protected Data object.");
1953     }
1954 jgs 102 Data tempValue(value);
1955     typeMatchLeft(tempValue);
1956     typeMatchRight(tempValue);
1957     m_data->setSlice(tempValue.m_data.get(),region);
1958     }
1959    
1960     void
1961     Data::typeMatchLeft(Data& right) const
1962     {
1963     if (isExpanded()){
1964     right.expand();
1965     } else if (isTagged()) {
1966     if (right.isConstant()) {
1967     right.tag();
1968     }
1969     }
1970     }
1971    
1972     void
1973     Data::typeMatchRight(const Data& right)
1974     {
1975 jgs 94 if (isTagged()) {
1976     if (right.isExpanded()) {
1977     expand();
1978     }
1979     } else if (isConstant()) {
1980     if (right.isExpanded()) {
1981     expand();
1982     } else if (right.isTagged()) {
1983     tag();
1984     }
1985     }
1986     }
1987    
1988     void
1989 gross 1044 Data::setTaggedValueByName(std::string name,
1990 ksteube 1312 const boost::python::object& value)
1991 gross 1044 {
1992 jfenwick 1872 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
1993     int tagKey=getFunctionSpace().getDomain()->getTag(name);
1994 gross 1044 setTaggedValue(tagKey,value);
1995     }
1996     }
1997     void
1998 jgs 94 Data::setTaggedValue(int tagKey,
1999     const boost::python::object& value)
2000     {
2001 gross 783 if (isProtected()) {
2002     throw DataException("Error - attempt to update protected Data object.");
2003     }
2004 jgs 94 //
2005     // Ensure underlying data object is of type DataTagged
2006 gross 1358 if (isConstant()) tag();
2007 jgs 94
2008 matt 1319 numeric::array asNumArray(value);
2009 jgs 94
2010 matt 1319
2011     // extract the shape of the numarray
2012 jfenwick 1796 DataTypes::ShapeType tempShape;
2013 matt 1319 for (int i=0; i < asNumArray.getrank(); i++) {
2014     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2015     }
2016    
2017     // get the space for the data vector
2018 jfenwick 1796 // int len = DataTypes::noValues(tempShape);
2019     // DataVector temp_data(len, 0.0, len);
2020     // DataArrayView temp_dataView(temp_data, tempShape);
2021     // temp_dataView.copy(asNumArray);
2022 matt 1319
2023 jfenwick 1796 DataVector temp_data2;
2024     temp_data2.copyFromNumArray(asNumArray);
2025    
2026 jgs 94 //
2027     // Call DataAbstract::setTaggedValue
2028 jfenwick 1796 //m_data->setTaggedValue(tagKey,temp_dataView);
2029    
2030     m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2031 jgs 94 }
2032    
2033 jfenwick 1796
2034 jgs 110 void
2035 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2036 jfenwick 1796 const DataTypes::ShapeType& pointshape,
2037     const DataTypes::ValueType& value,
2038     int dataOffset)
2039 jgs 121 {
2040 gross 783 if (isProtected()) {
2041     throw DataException("Error - attempt to update protected Data object.");
2042     }
2043 jgs 121 //
2044     // Ensure underlying data object is of type DataTagged
2045 gross 1358 if (isConstant()) tag();
2046 jgs 121
2047     //
2048     // Call DataAbstract::setTaggedValue
2049 jfenwick 1796 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2050 jgs 121 }
2051    
2052 jgs 149 int
2053     Data::getTagNumber(int dpno)
2054     {
2055 jfenwick 1803 if (isEmpty())
2056     {
2057     throw DataException("Error - operation not permitted on instances of DataEmpty.");
2058     }
2059 gross 1409 return getFunctionSpace().getTagFromDataPointNo(dpno);
2060 jgs 149 }
2061    
2062 jgs 119
2063 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2064     {
2065     o << data.toString();
2066     return o;
2067     }
2068 bcumming 782
2069 ksteube 813 Data
2070     escript::C_GeneralTensorProduct(Data& arg_0,
2071     Data& arg_1,
2072     int axis_offset,
2073     int transpose)
2074     {
2075 gross 826 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2076 ksteube 813 // SM is the product of the last axis_offset entries in arg_0.getShape().
2077    
2078     // Interpolate if necessary and find an appropriate function space
2079 gross 826 Data arg_0_Z, arg_1_Z;
2080 ksteube 813 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2081     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2082 gross 826 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2083     arg_1_Z = Data(arg_1);
2084 ksteube 813 }
2085     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2086 gross 826 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2087     arg_0_Z =Data(arg_0);
2088 ksteube 813 }
2089     else {
2090     throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2091     }
2092 gross 826 } else {
2093     arg_0_Z = Data(arg_0);
2094     arg_1_Z = Data(arg_1);
2095 ksteube 813 }
2096     // Get rank and shape of inputs
2097 gross 826 int rank0 = arg_0_Z.getDataPointRank();
2098     int rank1 = arg_1_Z.getDataPointRank();
2099 jfenwick 1796 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2100     const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2101 ksteube 813
2102     // Prepare for the loops of the product and verify compatibility of shapes
2103     int start0=0, start1=0;
2104     if (transpose == 0) {}
2105     else if (transpose == 1) { start0 = axis_offset; }
2106     else if (transpose == 2) { start1 = rank1-axis_offset; }
2107     else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2108    
2109 jfenwick 1796
2110 ksteube 813 // Adjust the shapes for transpose
2111 jfenwick 1796 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2112     DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2113     for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2114     for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2115 ksteube 813
2116     #if 0
2117     // For debugging: show shape after transpose
2118     char tmp[100];
2119     std::string shapeStr;
2120     shapeStr = "(";
2121     for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2122     shapeStr += ")";
2123     cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2124     shapeStr = "(";
2125     for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2126     shapeStr += ")";
2127     cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2128     #endif
2129    
2130     // Prepare for the loops of the product
2131     int SL=1, SM=1, SR=1;
2132     for (int i=0; i<rank0-axis_offset; i++) {
2133     SL *= tmpShape0[i];
2134     }
2135     for (int i=rank0-axis_offset; i<rank0; i++) {
2136     if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2137     throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2138     }
2139     SM *= tmpShape0[i];
2140     }
2141     for (int i=axis_offset; i<rank1; i++) {
2142     SR *= tmpShape1[i];
2143     }
2144    
2145 jfenwick 1796 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2146     DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2147     { // block to limit the scope of out_index
2148     int out_index=0;
2149     for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2150     for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2151     }
2152 ksteube 813
2153     // Declare output Data object
2154 gross 826 Data res;
2155 ksteube 813
2156 gross 826 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2157     res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2158 jfenwick 1796 double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2159     double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2160     double *ptr_2 = &(res.getDataAtOffset(0));
2161 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2162     }
2163 gross 826 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2164 ksteube 813
2165     // Prepare the DataConstant input
2166 gross 826 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2167 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2168    
2169     // Borrow DataTagged input from Data object
2170 gross 826 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2171 ksteube 813 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2172    
2173     // Prepare a DataTagged output 2
2174 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2175     res.tag();
2176     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2177 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2178    
2179     // Prepare offset into DataConstant
2180     int offset_0 = tmp_0->getPointOffset(0,0);
2181 jfenwick 1796 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2182 ksteube 813 // Get the views
2183 jfenwick 1796 // DataArrayView view_1 = tmp_1->getDefaultValue();
2184     // DataArrayView view_2 = tmp_2->getDefaultValue();
2185     // // Get the pointers to the actual data
2186     // double *ptr_1 = &((view_1.getData())[0]);
2187     // double *ptr_2 = &((view_2.getData())[0]);
2188    
2189     double *ptr_1 = &(tmp_1->getDefaultValue(0));
2190     double *ptr_2 = &(tmp_2->getDefaultValue(0));
2191    
2192    
2193 ksteube 813 // Compute an MVP for the default
2194     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2195     // Compute an MVP for each tag
2196     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2197     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2198     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2199 jfenwick 1796 tmp_2->addTag(i->first);
2200     // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2201     // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2202     // double *ptr_1 = &view_1.getData(0);
2203     // double *ptr_2 = &view_2.getData(0);
2204    
2205     double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2206     double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2207    
2208 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2209     }
2210    
2211     }
2212 gross 826 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2213 ksteube 813
2214 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2215     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2216     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2217     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2218 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2219     if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2220     if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2221     int sampleNo_1,dataPointNo_1;
2222 gross 826 int numSamples_1 = arg_1_Z.getNumSamples();
2223     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2224 ksteube 813 int offset_0 = tmp_0->getPointOffset(0,0);
2225     #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2226     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2227     for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2228     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2229     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2230 jfenwick 1796 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2231     double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2232     double *ptr_2 = &(res.getDataAtOffset(offset_2));
2233 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2234     }
2235     }
2236    
2237     }
2238 gross 826 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2239 ksteube 813
2240     // Borrow DataTagged input from Data object
2241 gross 826 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2242 ksteube 813 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2243    
2244     // Prepare the DataConstant input
2245 gross 826 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2246 ksteube 813 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2247    
2248     // Prepare a DataTagged output 2
2249 gross 826 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2250     res.tag();
2251     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2252 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2253    
2254     // Prepare offset into DataConstant
2255     int offset_1 = tmp_1->getPointOffset(0,0);
2256 jfenwick 1796 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2257 ksteube 813 // Get the views
2258 jfenwick 1796 // DataArrayView view_0 = tmp_0->getDefaultValue();
2259     // DataArrayView view_2 = tmp_2->getDefaultValue();
2260     // // Get the pointers to the actual data
2261     // double *ptr_0 = &((view_0.getData())[0]);
2262     // double *ptr_2 = &((view_2.getData())[0]);
2263    
2264     double *ptr_0 = &(tmp_0->getDefaultValue(0));
2265     double *ptr_2 = &(tmp_2->getDefaultValue(0));
2266    
2267 ksteube 813 // Compute an MVP for the default
2268     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2269     // Compute an MVP for each tag
2270     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2271     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2272     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2273 jfenwick 1796 // tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2274     // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2275     // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2276     // double *ptr_0 = &view_0.getData(0);
2277     // double *ptr_2 = &view_2.getData(0);
2278    
2279     tmp_2->addTag(i->first);
2280     double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2281     double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2282 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2283     }
2284    
2285     }
2286 gross 826 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2287 ksteube 813
2288     // Borrow DataTagged input from Data object
2289 gross 826 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2290 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2291    
2292     // Borrow DataTagged input from Data object
2293 gross 826 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2294 ksteube 813 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2295    
2296     // Prepare a DataTagged output 2
2297 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2298     res.tag(); // DataTagged output
2299     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2300 ksteube 813 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2301    
2302 jfenwick 1796 // // Get the views
2303     // DataArrayView view_0 = tmp_0->getDefaultValue();
2304     // DataArrayView view_1 = tmp_1->getDefaultValue();
2305     // DataArrayView view_2 = tmp_2->getDefaultValue();
2306     // // Get the pointers to the actual data
2307     // double *ptr_0 = &((view_0.getData())[0]);
2308     // double *ptr_1 = &((view_1.getData())[0]);
2309     // double *ptr_2 = &((view_2.getData())[0]);
2310    
2311     double *ptr_0 = &(tmp_0->getDefaultValue(0));
2312     double *ptr_1 = &(tmp_1->getDefaultValue(0));
2313     double *ptr_2 = &(tmp_2->getDefaultValue(0));
2314    
2315    
2316 ksteube 813 // Compute an MVP for the default
2317     matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2318     // Merge the tags
2319     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2320     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2321     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2322     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2323 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2324 ksteube 813 }
2325     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2326 jfenwick 1796 tmp_2->addTag(i->first);
2327 ksteube 813 }
2328     // Compute an MVP for each tag
2329     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2330     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2331 jfenwick 1796 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2332     // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2333     // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2334     // double *ptr_0 = &view_0.getData(0);
2335     // double *ptr_1 = &view_1.getData(0);
2336     // double *ptr_2 = &view_2.getData(0);
2337    
2338     double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2339     double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2340     double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2341    
2342 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2343     }
2344    
2345     }
2346 gross 826 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2347 ksteube 813
2348     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2349 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2350     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2351     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2352     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2353 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2354     if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2355     if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2356     int sampleNo_0,dataPointNo_0;
2357 gross 826 int numSamples_0 = arg_0_Z.getNumSamples();
2358     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2359 ksteube 813 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2360     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2361     int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2362 jfenwick 1796 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2363 ksteube 813 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2364     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2365     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2366 jfenwick 1796 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2367     double *ptr_2 = &(res.getDataAtOffset(offset_2));
2368 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2369     }
2370     }
2371    
2372     }
2373 gross 826 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2374 ksteube 813
2375 gross 826 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2376     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2377     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2378     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2379 ksteube 813 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2380     if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2381     if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2382     int sampleNo_0,dataPointNo_0;
2383 gross 826 int numSamples_0 = arg_0_Z.getNumSamples();
2384     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2385 ksteube 813 int offset_1 = tmp_1->getPointOffset(0,0);
2386     #pragma omp parallel for private(sampleNo_0,data