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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1868 - (hide annotations)
Thu Oct 9 06:30:49 2008 UTC (13 years, 1 month ago) by jfenwick
File size: 75555 byte(s)
Branch commit.
Bulk resolve for + - * /

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