/[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 1899 - (show annotations)
Mon Oct 20 05:13:24 2008 UTC (10 years, 6 months ago) by jfenwick
File size: 78579 byte(s)
Branch commit.
Added some doco to DataLazy.
Made Data::integrate aware of DataLazy.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26