/[escript]/branches/more_shared_ptrs_from_1812/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/more_shared_ptrs_from_1812/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1828 - (hide annotations)
Thu Oct 2 04:52:11 2008 UTC (11 years, 11 months ago) by jfenwick
File size: 73046 byte(s)
Branch commit.
Added getPtr to DataAbstract.
Passes all unit tests.


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