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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2089 - (show annotations)
Mon Nov 24 06:07:29 2008 UTC (10 years, 9 months ago) by jfenwick
File size: 85302 byte(s)
Resolve mantis issue 216.

Removed the asAbstractContinuousDomain(X) member which did not take cast 
failure into account.
Rewrote the few places where it was used.



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26