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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26