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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (hide annotations)
Wed Sep 17 01:45:46 2008 UTC (11 years, 6 months ago) by jfenwick
File size: 73653 byte(s)
Merged noarrayview branch onto trunk.


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