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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26