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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 6 months ago) by ksteube
File size: 73725 byte(s)
Copyright updated in all files

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