/[escript]/branches/arrexp_trunk2098/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2143 - (show annotations)
Tue Dec 9 06:28:10 2008 UTC (12 years, 4 months ago) by jfenwick
File size: 84682 byte(s)
Branch commit.
Added scons option usenumarray (defaults to yes).
This controls the NONUMARRAY define which switches off all explicit 
references to boost::python::numeric::array

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26