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

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

Parent Directory Parent Directory | Revision Log Revision Log


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