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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1910 - (show annotations)
Thu Oct 23 03:05:28 2008 UTC (11 years, 8 months ago) by jfenwick
File size: 82102 byte(s)
Branch commit.
Support for ** added.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26