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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2086 - (show annotations)
Mon Nov 24 02:38:50 2008 UTC (10 years, 11 months ago) by jfenwick
File size: 85194 byte(s)
Added checks in C_GeneralTensorProduct (Data:: and Delayed forms) as 
well as the DataAbstract Constructor to prevent Objects with Rank>4 
being created.

Moved the relevant #define into systemdep.

Removed some comments.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26