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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1799 - (hide annotations)
Wed Sep 17 06:33:18 2008 UTC (11 years ago) by jfenwick
File size: 72960 byte(s)
Added Data::copySelf() [Note: this is exposed as copy() in python].
This method returns a pointer to a deep copy of the target.
There are c++ tests but no python tests for this yet.

All DataAbstracts now have a deepCopy() which simplifies the 
implementation of the compy methods.


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