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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1803 - (hide annotations)
Wed Sep 24 06:20:29 2008 UTC (11 years, 2 months ago) by jfenwick
File size: 73760 byte(s)
All about making DataEmpty instances throw.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Exposed getDim from AbstractDomain to python to fix bug.

Added isEmpty member to DataAbstract to allow it to throw is queries are 
made about a DataEmpty instance.


Added exceptions to DataAbstract, DataEmpty and Data to prevent calls 
being made against DataEmpty objects.
The following still work as expected on DataEmpty instances

copy, getDomain, getFunctionSpace, isEmpty, isExpanded, isProtected, 
isTagged, setprotection.

You can also call interpolate, however it should throw if you try to 
change FunctionSpaces.


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