/[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 1943 - (show annotations)
Wed Oct 29 04:05:14 2008 UTC (10 years, 6 months ago) by jfenwick
File size: 83824 byte(s)
Branch commit
Added interpolation.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26