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

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

Parent Directory Parent Directory | Revision Log Revision Log


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