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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1921 - (hide annotations)
Thu Oct 23 11:32:24 2008 UTC (11 years, 5 months ago) by phornby
File size: 75284 byte(s)
Jump through hoops to meet the OPENMP 2.5 restrictions on loop variables.
It's all ready for OPENMP 3.0 improvments, and should compile cleanly on other platforms.


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