/[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 1889 - (show annotations)
Thu Oct 16 05:57:09 2008 UTC (11 years, 7 months ago) by jfenwick
File size: 78264 byte(s)
Branch commit
Rewrote resolve to take into account Tagged and Constant Data.
Mixing expanded and Tagged does not work yet.

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
1047 boost::python::numeric::array
1048 Data::integrate() const
1049 {
1050 int index;
1051 int rank = getDataPointRank();
1052 DataTypes::ShapeType shape = getDataPointShape();
1053 int dataPointSize = getDataPointSize();
1054
1055 //
1056 // calculate the integral values
1057 vector<double> integrals(dataPointSize);
1058 vector<double> integrals_local(dataPointSize);
1059 #ifdef PASO_MPI
1060 AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
1061 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1062 double *tmp = new double[dataPointSize];
1063 double *tmp_local = new double[dataPointSize];
1064 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1065 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1066 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1067 delete[] tmp;
1068 delete[] tmp_local;
1069 #else
1070 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1071 #endif
1072
1073 //
1074 // create the numeric array to be returned
1075 // and load the array with the integral values
1076 boost::python::numeric::array bp_array(1.0);
1077 if (rank==0) {
1078 bp_array.resize(1);
1079 index = 0;
1080 bp_array[0] = integrals[index];
1081 }
1082 if (rank==1) {
1083 bp_array.resize(shape[0]);
1084 for (int i=0; i<shape[0]; i++) {
1085 index = i;
1086 bp_array[i] = integrals[index];
1087 }
1088 }
1089 if (rank==2) {
1090 bp_array.resize(shape[0],shape[1]);
1091 for (int i=0; i<shape[0]; i++) {
1092 for (int j=0; j<shape[1]; j++) {
1093 index = i + shape[0] * j;
1094 bp_array[make_tuple(i,j)] = integrals[index];
1095 }
1096 }
1097 }
1098 if (rank==3) {
1099 bp_array.resize(shape[0],shape[1],shape[2]);
1100 for (int i=0; i<shape[0]; i++) {
1101 for (int j=0; j<shape[1]; j++) {
1102 for (int k=0; k<shape[2]; k++) {
1103 index = i + shape[0] * ( j + shape[1] * k );
1104 bp_array[make_tuple(i,j,k)] = integrals[index];
1105 }
1106 }
1107 }
1108 }
1109 if (rank==4) {
1110 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1111 for (int i=0; i<shape[0]; i++) {
1112 for (int j=0; j<shape[1]; j++) {
1113 for (int k=0; k<shape[2]; k++) {
1114 for (int l=0; l<shape[3]; l++) {
1115 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1116 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1117 }
1118 }
1119 }
1120 }
1121 }
1122
1123 //
1124 // return the loaded array
1125 return bp_array;
1126 }
1127
1128 Data
1129 Data::sin() const
1130 {
1131 if (isLazy())
1132 {
1133 DataLazy* c=new DataLazy(borrowDataPtr(),SIN);
1134 return Data(c);
1135 }
1136 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1137 }
1138
1139 Data
1140 Data::cos() const
1141 {
1142 if (isLazy())
1143 {
1144 DataLazy* c=new DataLazy(borrowDataPtr(),COS);
1145 return Data(c);
1146 }
1147 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1148 }
1149
1150 Data
1151 Data::tan() const
1152 {
1153 if (isLazy())
1154 {
1155 DataLazy* c=new DataLazy(borrowDataPtr(),TAN);
1156 return Data(c);
1157 }
1158 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1159 }
1160
1161 Data
1162 Data::asin() const
1163 {
1164 if (isLazy())
1165 {
1166 DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);
1167 return Data(c);
1168 }
1169 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1170 }
1171
1172 Data
1173 Data::acos() const
1174 {
1175 if (isLazy())
1176 {
1177 DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);
1178 return Data(c);
1179 }
1180 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1181 }
1182
1183
1184 Data
1185 Data::atan() const
1186 {
1187 if (isLazy())
1188 {
1189 DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);
1190 return Data(c);
1191 }
1192 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1193 }
1194
1195 Data
1196 Data::sinh() const
1197 {
1198 if (isLazy())
1199 {
1200 DataLazy* c=new DataLazy(borrowDataPtr(),SINH);
1201 return Data(c);
1202 }
1203 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1204 }
1205
1206 Data
1207 Data::cosh() const
1208 {
1209 if (isLazy())
1210 {
1211 DataLazy* c=new DataLazy(borrowDataPtr(),COSH);
1212 return Data(c);
1213 }
1214 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1215 }
1216
1217 Data
1218 Data::tanh() const
1219 {
1220 if (isLazy())
1221 {
1222 DataLazy* c=new DataLazy(borrowDataPtr(),TANH);
1223 return Data(c);
1224 }
1225 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1226 }
1227
1228
1229 Data
1230 Data::erf() const
1231 {
1232 #ifdef _WIN32
1233 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1234 #else
1235 if (isLazy())
1236 {
1237 DataLazy* c=new DataLazy(borrowDataPtr(),ERF);
1238 return Data(c);
1239 }
1240 return C_TensorUnaryOperation(*this, ::erf);
1241 #endif
1242 }
1243
1244 Data
1245 Data::asinh() const
1246 {
1247 if (isLazy())
1248 {
1249 DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1250 return Data(c);
1251 }
1252 #ifdef _WIN32
1253 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1254 #else
1255 return C_TensorUnaryOperation(*this, ::asinh);
1256 #endif
1257 }
1258
1259 Data
1260 Data::acosh() const
1261 {
1262 if (isLazy())
1263 {
1264 DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1265 return Data(c);
1266 }
1267 #ifdef _WIN32
1268 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1269 #else
1270 return C_TensorUnaryOperation(*this, ::acosh);
1271 #endif
1272 }
1273
1274 Data
1275 Data::atanh() const
1276 {
1277 if (isLazy())
1278 {
1279 DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1280 return Data(c);
1281 }
1282 #ifdef _WIN32
1283 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1284 #else
1285 return C_TensorUnaryOperation(*this, ::atanh);
1286 #endif
1287 }
1288
1289 Data
1290 Data::log10() const
1291 { if (isLazy())
1292 {
1293 DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);
1294 return Data(c);
1295 }
1296 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1297 }
1298
1299 Data
1300 Data::log() const
1301 {
1302 if (isLazy())
1303 {
1304 DataLazy* c=new DataLazy(borrowDataPtr(),LOG);
1305 return Data(c);
1306 }
1307 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1308 }
1309
1310 Data
1311 Data::sign() const
1312 {
1313 if (isLazy())
1314 {
1315 DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);
1316 return Data(c);
1317 }
1318 return C_TensorUnaryOperation(*this, escript::fsign);
1319 }
1320
1321 Data
1322 Data::abs() const
1323 {
1324 if (isLazy())
1325 {
1326 DataLazy* c=new DataLazy(borrowDataPtr(),ABS);
1327 return Data(c);
1328 }
1329 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1330 }
1331
1332 Data
1333 Data::neg() const
1334 {
1335 if (isLazy())
1336 {
1337 DataLazy* c=new DataLazy(borrowDataPtr(),NEG);
1338 return Data(c);
1339 }
1340 return C_TensorUnaryOperation(*this, negate<double>());
1341 }
1342
1343 Data
1344 Data::pos() const
1345 {
1346 // not doing lazy check here is deliberate.
1347 // since a deep copy of lazy data should be cheap, I'll just let it happen now
1348 Data result;
1349 // perform a deep copy
1350 result.copy(*this);
1351 return result;
1352 }
1353
1354 Data
1355 Data::exp() const
1356 {
1357 if (isLazy())
1358 {
1359 DataLazy* c=new DataLazy(borrowDataPtr(),EXP);
1360 return Data(c);
1361 }
1362 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1363 }
1364
1365 Data
1366 Data::sqrt() const
1367 {
1368 if (isLazy())
1369 {
1370 DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);
1371 return Data(c);
1372 }
1373 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1374 }
1375
1376 double
1377 Data::Lsup_const() const
1378 {
1379 if (isLazy())
1380 {
1381 throw DataException("Error - cannot compute Lsup for constant lazy data.");
1382 }
1383 return LsupWorker();
1384 }
1385
1386 double
1387 Data::Lsup()
1388 {
1389 if (isLazy())
1390 {
1391 expand();
1392 }
1393 return LsupWorker();
1394 }
1395
1396 double
1397 Data::sup_const() const
1398 {
1399 if (isLazy())
1400 {
1401 throw DataException("Error - cannot compute sup for constant lazy data.");
1402 }
1403 return supWorker();
1404 }
1405
1406 double
1407 Data::sup()
1408 {
1409 if (isLazy())
1410 {
1411 expand();
1412 }
1413 return supWorker();
1414 }
1415
1416 double
1417 Data::inf_const() const
1418 {
1419 if (isLazy())
1420 {
1421 throw DataException("Error - cannot compute inf for constant lazy data.");
1422 }
1423 return infWorker();
1424 }
1425
1426 double
1427 Data::inf()
1428 {
1429 if (isLazy())
1430 {
1431 expand();
1432 }
1433 return infWorker();
1434 }
1435
1436 double
1437 Data::LsupWorker() const
1438 {
1439 double localValue;
1440 //
1441 // set the initial absolute maximum value to zero
1442
1443 AbsMax abs_max_func;
1444 localValue = algorithm(abs_max_func,0);
1445 #ifdef PASO_MPI
1446 double globalValue;
1447 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1448 return globalValue;
1449 #else
1450 return localValue;
1451 #endif
1452 }
1453
1454 double
1455 Data::supWorker() const
1456 {
1457 double localValue;
1458 //
1459 // set the initial maximum value to min possible double
1460 FMax fmax_func;
1461 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1462 #ifdef PASO_MPI
1463 double globalValue;
1464 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1465 return globalValue;
1466 #else
1467 return localValue;
1468 #endif
1469 }
1470
1471 double
1472 Data::infWorker() const
1473 {
1474 double localValue;
1475 //
1476 // set the initial minimum value to max possible double
1477 FMin fmin_func;
1478 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1479 #ifdef PASO_MPI
1480 double globalValue;
1481 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1482 return globalValue;
1483 #else
1484 return localValue;
1485 #endif
1486 }
1487
1488 /* TODO */
1489 /* global reduction */
1490 Data
1491 Data::maxval() const
1492 {
1493 //
1494 // set the initial maximum value to min possible double
1495 FMax fmax_func;
1496 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1497 }
1498
1499 Data
1500 Data::minval() const
1501 {
1502 //
1503 // set the initial minimum value to max possible double
1504 FMin fmin_func;
1505 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1506 }
1507
1508 Data
1509 Data::swapaxes(const int axis0, const int axis1) const
1510 {
1511 int axis0_tmp,axis1_tmp;
1512 DataTypes::ShapeType s=getDataPointShape();
1513 DataTypes::ShapeType ev_shape;
1514 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1515 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1516 int rank=getDataPointRank();
1517 if (rank<2) {
1518 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1519 }
1520 if (axis0<0 || axis0>rank-1) {
1521 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1522 }
1523 if (axis1<0 || axis1>rank-1) {
1524 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1525 }
1526 if (axis0 == axis1) {
1527 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1528 }
1529 if (axis0 > axis1) {
1530 axis0_tmp=axis1;
1531 axis1_tmp=axis0;
1532 } else {
1533 axis0_tmp=axis0;
1534 axis1_tmp=axis1;
1535 }
1536 for (int i=0; i<rank; i++) {
1537 if (i == axis0_tmp) {
1538 ev_shape.push_back(s[axis1_tmp]);
1539 } else if (i == axis1_tmp) {
1540 ev_shape.push_back(s[axis0_tmp]);
1541 } else {
1542 ev_shape.push_back(s[i]);
1543 }
1544 }
1545 Data ev(0.,ev_shape,getFunctionSpace());
1546 ev.typeMatchRight(*this);
1547 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1548 return ev;
1549
1550 }
1551
1552 Data
1553 Data::symmetric() const
1554 {
1555 // check input
1556 DataTypes::ShapeType s=getDataPointShape();
1557 if (getDataPointRank()==2) {
1558 if(s[0] != s[1])
1559 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1560 }
1561 else if (getDataPointRank()==4) {
1562 if(!(s[0] == s[2] && s[1] == s[3]))
1563 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1564 }
1565 else {
1566 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1567 }
1568 Data ev(0.,getDataPointShape(),getFunctionSpace());
1569 ev.typeMatchRight(*this);
1570 m_data->symmetric(ev.m_data.get());
1571 return ev;
1572 }
1573
1574 Data
1575 Data::nonsymmetric() const
1576 {
1577 // check input
1578 DataTypes::ShapeType s=getDataPointShape();
1579 if (getDataPointRank()==2) {
1580 if(s[0] != s[1])
1581 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1582 DataTypes::ShapeType ev_shape;
1583 ev_shape.push_back(s[0]);
1584 ev_shape.push_back(s[1]);
1585 Data ev(0.,ev_shape,getFunctionSpace());
1586 ev.typeMatchRight(*this);
1587 m_data->nonsymmetric(ev.m_data.get());
1588 return ev;
1589 }
1590 else if (getDataPointRank()==4) {
1591 if(!(s[0] == s[2] && s[1] == s[3]))
1592 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1593 DataTypes::ShapeType ev_shape;
1594 ev_shape.push_back(s[0]);
1595 ev_shape.push_back(s[1]);
1596 ev_shape.push_back(s[2]);
1597 ev_shape.push_back(s[3]);
1598 Data ev(0.,ev_shape,getFunctionSpace());
1599 ev.typeMatchRight(*this);
1600 m_data->nonsymmetric(ev.m_data.get());
1601 return ev;
1602 }
1603 else {
1604 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1605 }
1606 }
1607
1608 Data
1609 Data::trace(int axis_offset) const
1610 {
1611 DataTypes::ShapeType s=getDataPointShape();
1612 if (getDataPointRank()==2) {
1613 DataTypes::ShapeType ev_shape;
1614 Data ev(0.,ev_shape,getFunctionSpace());
1615 ev.typeMatchRight(*this);
1616 m_data->trace(ev.m_data.get(), axis_offset);
1617 return ev;
1618 }
1619 if (getDataPointRank()==3) {
1620 DataTypes::ShapeType ev_shape;
1621 if (axis_offset==0) {
1622 int s2=s[2];
1623 ev_shape.push_back(s2);
1624 }
1625 else if (axis_offset==1) {
1626 int s0=s[0];
1627 ev_shape.push_back(s0);
1628 }
1629 Data ev(0.,ev_shape,getFunctionSpace());
1630 ev.typeMatchRight(*this);
1631 m_data->trace(ev.m_data.get(), axis_offset);
1632 return ev;
1633 }
1634 if (getDataPointRank()==4) {
1635 DataTypes::ShapeType ev_shape;
1636 if (axis_offset==0) {
1637 ev_shape.push_back(s[2]);
1638 ev_shape.push_back(s[3]);
1639 }
1640 else if (axis_offset==1) {
1641 ev_shape.push_back(s[0]);
1642 ev_shape.push_back(s[3]);
1643 }
1644 else if (axis_offset==2) {
1645 ev_shape.push_back(s[0]);
1646 ev_shape.push_back(s[1]);
1647 }
1648 Data ev(0.,ev_shape,getFunctionSpace());
1649 ev.typeMatchRight(*this);
1650 m_data->trace(ev.m_data.get(), axis_offset);
1651 return ev;
1652 }
1653 else {
1654 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1655 }
1656 }
1657
1658 Data
1659 Data::transpose(int axis_offset) const
1660 {
1661 DataTypes::ShapeType s=getDataPointShape();
1662 DataTypes::ShapeType ev_shape;
1663 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1664 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1665 int rank=getDataPointRank();
1666 if (axis_offset<0 || axis_offset>rank) {
1667 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1668 }
1669 for (int i=0; i<rank; i++) {
1670 int index = (axis_offset+i)%rank;
1671 ev_shape.push_back(s[index]); // Append to new shape
1672 }
1673 Data ev(0.,ev_shape,getFunctionSpace());
1674 ev.typeMatchRight(*this);
1675 m_data->transpose(ev.m_data.get(), axis_offset);
1676 return ev;
1677 }
1678
1679 Data
1680 Data::eigenvalues() const
1681 {
1682 // check input
1683 DataTypes::ShapeType s=getDataPointShape();
1684 if (getDataPointRank()!=2)
1685 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1686 if(s[0] != s[1])
1687 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1688 // create return
1689 DataTypes::ShapeType ev_shape(1,s[0]);
1690 Data ev(0.,ev_shape,getFunctionSpace());
1691 ev.typeMatchRight(*this);
1692 m_data->eigenvalues(ev.m_data.get());
1693 return ev;
1694 }
1695
1696 const boost::python::tuple
1697 Data::eigenvalues_and_eigenvectors(const double tol) const
1698 {
1699 DataTypes::ShapeType s=getDataPointShape();
1700 if (getDataPointRank()!=2)
1701 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1702 if(s[0] != s[1])
1703 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1704 // create return
1705 DataTypes::ShapeType ev_shape(1,s[0]);
1706 Data ev(0.,ev_shape,getFunctionSpace());
1707 ev.typeMatchRight(*this);
1708 DataTypes::ShapeType V_shape(2,s[0]);
1709 Data V(0.,V_shape,getFunctionSpace());
1710 V.typeMatchRight(*this);
1711 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1712 return make_tuple(boost::python::object(ev),boost::python::object(V));
1713 }
1714
1715 const boost::python::tuple
1716 Data::minGlobalDataPoint() const
1717 {
1718 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1719 // abort (for unknown reasons) if there are openmp directives with it in the
1720 // surrounding function
1721
1722 int DataPointNo;
1723 int ProcNo;
1724 calc_minGlobalDataPoint(ProcNo,DataPointNo);
1725 return make_tuple(ProcNo,DataPointNo);
1726 }
1727
1728 void
1729 Data::calc_minGlobalDataPoint(int& ProcNo,
1730 int& DataPointNo) const
1731 {
1732 int i,j;
1733 int lowi=0,lowj=0;
1734 double min=numeric_limits<double>::max();
1735
1736 Data temp=minval();
1737
1738 int numSamples=temp.getNumSamples();
1739 int numDPPSample=temp.getNumDataPointsPerSample();
1740
1741 double next,local_min;
1742 int local_lowi,local_lowj;
1743
1744 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1745 {
1746 local_min=min;
1747 #pragma omp for private(i,j) schedule(static)
1748 for (i=0; i<numSamples; i++) {
1749 for (j=0; j<numDPPSample; j++) {
1750 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1751 if (next<local_min) {
1752 local_min=next;
1753 local_lowi=i;
1754 local_lowj=j;
1755 }
1756 }
1757 }
1758 #pragma omp critical
1759 if (local_min<min) {
1760 min=local_min;
1761 lowi=local_lowi;
1762 lowj=local_lowj;
1763 }
1764 }
1765
1766 #ifdef PASO_MPI
1767 // determine the processor on which the minimum occurs
1768 next = temp.getDataPoint(lowi,lowj);
1769 int lowProc = 0;
1770 double *globalMins = new double[get_MPISize()+1];
1771 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1772
1773 if( get_MPIRank()==0 ){
1774 next = globalMins[lowProc];
1775 for( i=1; i<get_MPISize(); i++ )
1776 if( next>globalMins[i] ){
1777 lowProc = i;
1778 next = globalMins[i];
1779 }
1780 }
1781 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1782
1783 delete [] globalMins;
1784 ProcNo = lowProc;
1785 #else
1786 ProcNo = 0;
1787 #endif
1788 DataPointNo = lowj + lowi * numDPPSample;
1789 }
1790
1791 void
1792 Data::saveDX(std::string fileName) const
1793 {
1794 if (isEmpty())
1795 {
1796 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1797 }
1798 boost::python::dict args;
1799 args["data"]=boost::python::object(this);
1800 getDomain()->saveDX(fileName,args);
1801 return;
1802 }
1803
1804 void
1805 Data::saveVTK(std::string fileName) const
1806 {
1807 if (isEmpty())
1808 {
1809 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1810 }
1811 boost::python::dict args;
1812 args["data"]=boost::python::object(this);
1813 getDomain()->saveVTK(fileName,args);
1814 return;
1815 }
1816
1817 Data&
1818 Data::operator+=(const Data& right)
1819 {
1820 if (isProtected()) {
1821 throw DataException("Error - attempt to update protected Data object.");
1822 }
1823 binaryOp(right,plus<double>());
1824 return (*this);
1825 }
1826
1827 Data&
1828 Data::operator+=(const boost::python::object& right)
1829 {
1830 Data tmp(right,getFunctionSpace(),false);
1831 binaryOp(tmp,plus<double>());
1832 return (*this);
1833 }
1834 Data&
1835 Data::operator=(const Data& other)
1836 {
1837 copy(other);
1838 return (*this);
1839 }
1840
1841 Data&
1842 Data::operator-=(const Data& right)
1843 {
1844 if (isProtected()) {
1845 throw DataException("Error - attempt to update protected Data object.");
1846 }
1847 binaryOp(right,minus<double>());
1848 return (*this);
1849 }
1850
1851 Data&
1852 Data::operator-=(const boost::python::object& right)
1853 {
1854 Data tmp(right,getFunctionSpace(),false);
1855 binaryOp(tmp,minus<double>());
1856 return (*this);
1857 }
1858
1859 Data&
1860 Data::operator*=(const Data& right)
1861 {
1862 if (isProtected()) {
1863 throw DataException("Error - attempt to update protected Data object.");
1864 }
1865 binaryOp(right,multiplies<double>());
1866 return (*this);
1867 }
1868
1869 Data&
1870 Data::operator*=(const boost::python::object& right)
1871 {
1872 Data tmp(right,getFunctionSpace(),false);
1873 binaryOp(tmp,multiplies<double>());
1874 return (*this);
1875 }
1876
1877 Data&
1878 Data::operator/=(const Data& right)
1879 {
1880 if (isProtected()) {
1881 throw DataException("Error - attempt to update protected Data object.");
1882 }
1883 binaryOp(right,divides<double>());
1884 return (*this);
1885 }
1886
1887 Data&
1888 Data::operator/=(const boost::python::object& right)
1889 {
1890 Data tmp(right,getFunctionSpace(),false);
1891 binaryOp(tmp,divides<double>());
1892 return (*this);
1893 }
1894
1895 Data
1896 Data::rpowO(const boost::python::object& left) const
1897 {
1898 Data left_d(left,*this);
1899 return left_d.powD(*this);
1900 }
1901
1902 Data
1903 Data::powO(const boost::python::object& right) const
1904 {
1905 Data tmp(right,getFunctionSpace(),false);
1906 return powD(tmp);
1907 }
1908
1909 Data
1910 Data::powD(const Data& right) const
1911 {
1912 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1913 }
1914
1915 //
1916 // NOTE: It is essential to specify the namespace this operator belongs to
1917 Data
1918 escript::operator+(const Data& left, const Data& right)
1919 {
1920 if (left.isLazy() || right.isLazy())
1921 {
1922 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
1923 return Data(c);
1924 }
1925 return C_TensorBinaryOperation(left, right, plus<double>());
1926 }
1927
1928 //
1929 // NOTE: It is essential to specify the namespace this operator belongs to
1930 Data
1931 escript::operator-(const Data& left, const Data& right)
1932 {
1933 if (left.isLazy() || right.isLazy())
1934 {
1935 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
1936 return Data(c);
1937 }
1938 return C_TensorBinaryOperation(left, right, minus<double>());
1939 }
1940
1941 //
1942 // NOTE: It is essential to specify the namespace this operator belongs to
1943 Data
1944 escript::operator*(const Data& left, const Data& right)
1945 {
1946 if (left.isLazy() || right.isLazy())
1947 {
1948 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
1949 return Data(c);
1950 }
1951 return C_TensorBinaryOperation(left, right, multiplies<double>());
1952 }
1953
1954 //
1955 // NOTE: It is essential to specify the namespace this operator belongs to
1956 Data
1957 escript::operator/(const Data& left, const Data& right)
1958 {
1959 if (left.isLazy() || right.isLazy())
1960 {
1961 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),DIV);
1962 return Data(c);
1963 }
1964 return C_TensorBinaryOperation(left, right, divides<double>());
1965 }
1966
1967 //
1968 // NOTE: It is essential to specify the namespace this operator belongs to
1969 Data
1970 escript::operator+(const Data& left, const boost::python::object& right)
1971 {
1972 return left+Data(right,left.getFunctionSpace(),false);
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 boost::python::object& right)
1979 {
1980 return left-Data(right,left.getFunctionSpace(),false);
1981 }
1982
1983 //
1984 // NOTE: It is essential to specify the namespace this operator belongs to
1985 Data
1986 escript::operator*(const Data& left, const boost::python::object& right)
1987 {
1988 return left*Data(right,left.getFunctionSpace(),false);
1989 }
1990
1991 //
1992 // NOTE: It is essential to specify the namespace this operator belongs to
1993 Data
1994 escript::operator/(const Data& left, const boost::python::object& right)
1995 {
1996 return left/Data(right,left.getFunctionSpace(),false);
1997 }
1998
1999 //
2000 // NOTE: It is essential to specify the namespace this operator belongs to
2001 Data
2002 escript::operator+(const boost::python::object& left, const Data& right)
2003 {
2004 return Data(left,right.getFunctionSpace(),false)+right;
2005 }
2006
2007 //
2008 // NOTE: It is essential to specify the namespace this operator belongs to
2009 Data
2010 escript::operator-(const boost::python::object& left, const Data& right)
2011 {
2012 return Data(left,right.getFunctionSpace(),false)-right;
2013 }
2014
2015 //
2016 // NOTE: It is essential to specify the namespace this operator belongs to
2017 Data
2018 escript::operator*(const boost::python::object& left, const Data& right)
2019 {
2020 return Data(left,right.getFunctionSpace(),false)*right;
2021 }
2022
2023 //
2024 // NOTE: It is essential to specify the namespace this operator belongs to
2025 Data
2026 escript::operator/(const boost::python::object& left, const Data& right)
2027 {
2028 return Data(left,right.getFunctionSpace(),false)/right;
2029 }
2030
2031
2032 /* TODO */
2033 /* global reduction */
2034 Data
2035 Data::getItem(const boost::python::object& key) const
2036 {
2037 // const DataArrayView& view=getPointDataView();
2038
2039 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2040
2041 if (slice_region.size()!=getDataPointRank()) {
2042 throw DataException("Error - slice size does not match Data rank.");
2043 }
2044
2045 return getSlice(slice_region);
2046 }
2047
2048 /* TODO */
2049 /* global reduction */
2050 Data
2051 Data::getSlice(const DataTypes::RegionType& region) const
2052 {
2053 return Data(*this,region);
2054 }
2055
2056 /* TODO */
2057 /* global reduction */
2058 void
2059 Data::setItemO(const boost::python::object& key,
2060 const boost::python::object& value)
2061 {
2062 Data tempData(value,getFunctionSpace());
2063 setItemD(key,tempData);
2064 }
2065
2066 void
2067 Data::setItemD(const boost::python::object& key,
2068 const Data& value)
2069 {
2070 // const DataArrayView& view=getPointDataView();
2071
2072 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2073 if (slice_region.size()!=getDataPointRank()) {
2074 throw DataException("Error - slice size does not match Data rank.");
2075 }
2076 if (getFunctionSpace()!=value.getFunctionSpace()) {
2077 setSlice(Data(value,getFunctionSpace()),slice_region);
2078 } else {
2079 setSlice(value,slice_region);
2080 }
2081 }
2082
2083 void
2084 Data::setSlice(const Data& value,
2085 const DataTypes::RegionType& region)
2086 {
2087 if (isProtected()) {
2088 throw DataException("Error - attempt to update protected Data object.");
2089 }
2090 if (isLazy())
2091 {
2092 throw DataException("Error - setSlice not permitted on lazy data.");
2093 }
2094 Data tempValue(value);
2095 typeMatchLeft(tempValue);
2096 typeMatchRight(tempValue);
2097 getReady()->setSlice(tempValue.m_data.get(),region);
2098 }
2099
2100 void
2101 Data::typeMatchLeft(Data& right) const
2102 {
2103 if (isExpanded()){
2104 right.expand();
2105 } else if (isTagged()) {
2106 if (right.isConstant()) {
2107 right.tag();
2108 }
2109 }
2110 }
2111
2112 void
2113 Data::typeMatchRight(const Data& right)
2114 {
2115 if (isTagged()) {
2116 if (right.isExpanded()) {
2117 expand();
2118 }
2119 } else if (isConstant()) {
2120 if (right.isExpanded()) {
2121 expand();
2122 } else if (right.isTagged()) {
2123 tag();
2124 }
2125 }
2126 }
2127
2128 void
2129 Data::setTaggedValueByName(std::string name,
2130 const boost::python::object& value)
2131 {
2132 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2133 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2134 setTaggedValue(tagKey,value);
2135 }
2136 }
2137 void
2138 Data::setTaggedValue(int tagKey,
2139 const boost::python::object& value)
2140 {
2141 if (isProtected()) {
2142 throw DataException("Error - attempt to update protected Data object.");
2143 }
2144 //
2145 // Ensure underlying data object is of type DataTagged
2146 if (isConstant()) tag();
2147
2148 numeric::array asNumArray(value);
2149
2150
2151 // extract the shape of the numarray
2152 DataTypes::ShapeType tempShape;
2153 for (int i=0; i < asNumArray.getrank(); i++) {
2154 tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2155 }
2156
2157 // get the space for the data vector
2158 // int len = DataTypes::noValues(tempShape);
2159 // DataVector temp_data(len, 0.0, len);
2160 // DataArrayView temp_dataView(temp_data, tempShape);
2161 // temp_dataView.copy(asNumArray);
2162
2163 DataVector temp_data2;
2164 temp_data2.copyFromNumArray(asNumArray);
2165
2166 //
2167 // Call DataAbstract::setTaggedValue
2168 //m_data->setTaggedValue(tagKey,temp_dataView);
2169
2170 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2171 }
2172
2173
2174 void
2175 Data::setTaggedValueFromCPP(int tagKey,
2176 const DataTypes::ShapeType& pointshape,
2177 const DataTypes::ValueType& value,
2178 int dataOffset)
2179 {
2180 if (isProtected()) {
2181 throw DataException("Error - attempt to update protected Data object.");
2182 }
2183 //
2184 // Ensure underlying data object is of type DataTagged
2185 if (isConstant()) tag();
2186
2187 //
2188 // Call DataAbstract::setTaggedValue
2189 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2190 }
2191
2192 int
2193 Data::getTagNumber(int dpno)
2194 {
2195 if (isEmpty())
2196 {
2197 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2198 }
2199 return getFunctionSpace().getTagFromDataPointNo(dpno);
2200 }
2201
2202
2203 ostream& escript::operator<<(ostream& o, const Data& data)
2204 {
2205 o << data.toString();
2206 return o;
2207 }
2208
2209 Data
2210 escript::C_GeneralTensorProduct(Data& arg_0,
2211 Data& arg_1,
2212 int axis_offset,
2213 int transpose)
2214 {
2215 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2216 // SM is the product of the last axis_offset entries in arg_0.getShape().
2217
2218 // Interpolate if necessary and find an appropriate function space
2219 Data arg_0_Z, arg_1_Z;
2220 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2221 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2222 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2223 arg_1_Z = Data(arg_1);
2224 }
2225 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2226 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2227 arg_0_Z =Data(arg_0);
2228 }
2229 else {
2230 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2231 }
2232 } else {
2233 arg_0_Z = Data(arg_0);
2234 arg_1_Z = Data(arg_1);
2235 }
2236 // Get rank and shape of inputs
2237 int rank0 = arg_0_Z.getDataPointRank();
2238 int rank1 = arg_1_Z.getDataPointRank();
2239 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2240 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2241
2242 // Prepare for the loops of the product and verify compatibility of shapes
2243 int start0=0, start1=0;
2244 if (transpose == 0) {}
2245 else if (transpose == 1) { start0 = axis_offset; }
2246 else if (transpose == 2) { start1 = rank1-axis_offset; }
2247 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2248
2249
2250 // Adjust the shapes for transpose
2251 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2252 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2253 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2254 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2255
2256 #if 0
2257 // For debugging: show shape after transpose
2258 char tmp[100];
2259 std::string shapeStr;
2260 shapeStr = "(";
2261 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2262 shapeStr += ")";
2263 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2264 shapeStr = "(";
2265 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2266 shapeStr += ")";
2267 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2268 #endif
2269
2270 // Prepare for the loops of the product
2271 int SL=1, SM=1, SR=1;
2272 for (int i=0; i<rank0-axis_offset; i++) {
2273 SL *= tmpShape0[i];
2274 }
2275 for (int i=rank0-axis_offset; i<rank0; i++) {
2276 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2277 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2278 }
2279 SM *= tmpShape0[i];
2280 }
2281 for (int i=axis_offset; i<rank1; i++) {
2282 SR *= tmpShape1[i];
2283 }
2284
2285 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2286 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2287 { // block to limit the scope of out_index
2288 int out_index=0;
2289 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2290 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2291 }
2292
2293 // Declare output Data object
2294 Data res;
2295
2296 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2297 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2298 double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2299 double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2300 double *ptr_2 = &(res.getDataAtOffset(0));
2301 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2302 }
2303 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2304
2305 // Prepare the DataConstant input
2306 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2307 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2308
2309 // Borrow DataTagged input from Data object
2310 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2311 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2312
2313 // Prepare a DataTagged output 2
2314 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2315 res.tag();
2316 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2317 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2318
2319 // Prepare offset into DataConstant
2320 int offset_0 = tmp_0->getPointOffset(0,0);
2321 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2322 // Get the views
2323 // DataArrayView view_1 = tmp_1->getDefaultValue();
2324 // DataArrayView view_2 = tmp_2->getDefaultValue();
2325 // // Get the pointers to the actual data
2326 // double *ptr_1 = &((view_1.getData())[0]);
2327 // double *ptr_2 = &((view_2.getData())[0]);
2328
2329 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2330 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2331
2332
2333 // Compute an MVP for the default
2334 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2335 // Compute an MVP for each tag
2336 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2337 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2338 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2339 tmp_2->addTag(i->first);
2340 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2341 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2342 // double *ptr_1 = &view_1.getData(0);
2343 // double *ptr_2 = &view_2.getData(0);
2344
2345 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2346 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2347
2348 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2349 }
2350
2351 }
2352 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2353
2354 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2355 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2356 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2357 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2358 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2359 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2360 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2361 int sampleNo_1,dataPointNo_1;
2362 int numSamples_1 = arg_1_Z.getNumSamples();
2363 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2364 int offset_0 = tmp_0->getPointOffset(0,0);
2365 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2366 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2367 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2368 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2369 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2370 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2371 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2372 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2373 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2374 }
2375 }
2376
2377 }
2378 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2379
2380 // Borrow DataTagged input from Data object
2381 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2382 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2383
2384 // Prepare the DataConstant input
2385 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2386 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2387
2388 // Prepare a DataTagged output 2
2389 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2390 res.tag();
2391 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2392 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2393
2394 // Prepare offset into DataConstant
2395 int offset_1 = tmp_1->getPointOffset(0,0);
2396 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2397 // Get the views
2398 // DataArrayView view_0 = tmp_0->getDefaultValue();
2399 // DataArrayView view_2 = tmp_2->getDefaultValue();
2400 // // Get the pointers to the actual data
2401 // double *ptr_0 = &((view_0.getData())[0]);
2402 // double *ptr_2 = &((view_2.getData())[0]);
2403
2404 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2405 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2406
2407 // Compute an MVP for the default
2408 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2409 // Compute an MVP for each tag
2410 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2411 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2412 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2413 // tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2414 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2415 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2416 // double *ptr_0 = &view_0.getData(0);
2417 // double *ptr_2 = &view_2.getData(0);
2418
2419 tmp_2->addTag(i->first);
2420 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2421 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2422 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2423 }
2424
2425 }
2426 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2427
2428 // Borrow DataTagged input from Data object
2429 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2430 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2431
2432 // Borrow DataTagged input from Data object
2433 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2434 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2435
2436 // Prepare a DataTagged output 2
2437 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2438 res.tag(); // DataTagged output
2439 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2440 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2441
2442 // // Get the views
2443 // DataArrayView view_0 = tmp_0->getDefaultValue();
2444 // DataArrayView view_1 = tmp_1->getDefaultValue();
2445 // DataArrayView view_2 = tmp_2->getDefaultValue();
2446 // // Get the pointers to the actual data
2447 // double *ptr_0 = &((view_0.getData())[0]);
2448 // double *ptr_1 = &((view_1.getData())[0]);
2449 // double *ptr_2 = &((view_2.getData())[0]);
2450
2451 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2452 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2453 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2454
2455
2456 // Compute an MVP for the default
2457 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2458 // Merge the tags
2459 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2460 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2461 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2462 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2463 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2464 }
2465 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2466 tmp_2->addTag(i->first);
2467 }
2468 // Compute an MVP for each tag
2469 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2470 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2471 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2472 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2473 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2474 // double *ptr_0 = &view_0.getData(0);
2475 // double *ptr_1 = &view_1.getData(0);
2476 // double *ptr_2 = &view_2.getData(0);
2477
2478 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2479 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2480 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2481
2482 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2483 }
2484
2485 }
2486 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2487
2488 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2489 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2490 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2491 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2492 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2493 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2494 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2495 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2496 int sampleNo_0,dataPointNo_0;
2497 int numSamples_0 = arg_0_Z.getNumSamples();
2498 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2499 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2500 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2501 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2502 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2503 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2504 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2505 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2506 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2507 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2508 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2509 }
2510 }
2511
2512 }
2513 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2514
2515 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2516 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2517 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2518 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2519 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2520 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2521 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2522 int sampleNo_0,dataPointNo_0;
2523 int numSamples_0 = arg_0_Z.getNumSamples();
2524 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2525 int offset_1 = tmp_1->getPointOffset(0,0);
2526 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2527 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2528 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2529 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2530 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2531 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2532 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2533 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2534 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2535 }
2536 }
2537
2538
2539 }
2540 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2541
2542 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2543 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2544 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2545 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2546 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2547 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2548 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2549 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2550 int sampleNo_0,dataPointNo_0;
2551 int numSamples_0 = arg_0_Z.getNumSamples();
2552 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2553 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2554 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2555 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2556 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2557 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2558 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2559 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2560 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2561 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2562 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2563 }
2564 }
2565
2566 }
2567 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2568
2569 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2570 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2571 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2572 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2573 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2574 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2575 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2576 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2577 int sampleNo_0,dataPointNo_0;
2578 int numSamples_0 = arg_0_Z.getNumSamples();
2579 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2580 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2581 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2582 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2583 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2584 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2585 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2586 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2587 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2588 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2589 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2590 }
2591 }
2592
2593 }
2594 else {
2595 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2596 }
2597
2598 return res;
2599 }
2600
2601 DataAbstract*
2602 Data::borrowData() const
2603 {
2604 return m_data.get();
2605 }
2606
2607 // Not all that happy about returning a non-const from a const
2608 DataAbstract_ptr
2609 Data::borrowDataPtr() const
2610 {
2611 return m_data;
2612 }
2613
2614 // Not all that happy about returning a non-const from a const
2615 DataReady_ptr
2616 Data::borrowReadyPtr() const
2617 {
2618 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
2619 EsysAssert((dr!=0), "Error - casting to DataReady.");
2620 return dr;
2621 }
2622
2623 std::string
2624 Data::toString() const
2625 {
2626 static const DataTypes::ValueType::size_type TOO_MANY_POINTS=80;
2627 if (isLazy())
2628 { // This needs to change back to printing out something useful once the summary ops
2629 return m_data->toString(); // are defined
2630 }
2631 if (getNumDataPoints()*getDataPointSize()>TOO_MANY_POINTS)
2632 {
2633 stringstream temp;
2634 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
2635 return temp.str();
2636 }
2637 return m_data->toString();
2638 }
2639
2640
2641
2642 DataTypes::ValueType::const_reference
2643 Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
2644 {
2645 if (isLazy())
2646 {
2647 throw DataException("getDataAtOffset not permitted on lazy data.");
2648 }
2649 return getReady()->getDataAtOffset(i);
2650 }
2651
2652
2653 DataTypes::ValueType::reference
2654 Data::getDataAtOffset(DataTypes::ValueType::size_type i)
2655 {
2656 if (isLazy())
2657 {
2658 throw DataException("getDataAtOffset not permitted on lazy data.");
2659 }
2660 return getReady()->getDataAtOffset(i);
2661 }
2662
2663 DataTypes::ValueType::const_reference
2664 Data::getDataPoint(int sampleNo, int dataPointNo) const
2665 {
2666 if (!isReady())
2667 {
2668 throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");
2669 }
2670 else
2671 {
2672 const DataReady* dr=getReady();
2673 return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
2674 }
2675 }
2676
2677
2678 DataTypes::ValueType::reference
2679 Data::getDataPoint(int sampleNo, int dataPointNo)
2680 {
2681 if (!isReady())
2682 {
2683 throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");
2684 }
2685 else
2686 {
2687 DataReady* dr=getReady();
2688 return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
2689 }
2690 }
2691
2692
2693 /* Member functions specific to the MPI implementation */
2694
2695 void
2696 Data::print()
2697 {
2698 int i,j;
2699
2700 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2701 for( i=0; i<getNumSamples(); i++ )
2702 {
2703 printf( "[%6d]", i );
2704 for( j=0; j<getNumDataPointsPerSample(); j++ )
2705 printf( "\t%10.7g", (getSampleData(i))[j] );
2706 printf( "\n" );
2707 }
2708 }
2709 void
2710 Data::dump(const std::string fileName) const
2711 {
2712 try
2713 {
2714 return m_data->dump(fileName);
2715 }
2716 catch (exception& e)
2717 {
2718 cout << e.what() << endl;
2719 }
2720 }
2721
2722 int
2723 Data::get_MPISize() const
2724 {
2725 int size;
2726 #ifdef PASO_MPI
2727 int error;
2728 error = MPI_Comm_size( get_MPIComm(), &size );
2729 #else
2730 size = 1;
2731 #endif
2732 return size;
2733 }
2734
2735 int
2736 Data::get_MPIRank() const
2737 {
2738 int rank;
2739 #ifdef PASO_MPI
2740 int error;
2741 error = MPI_Comm_rank( get_MPIComm(), &rank );
2742 #else
2743 rank = 0;
2744 #endif
2745 return rank;
2746 }
2747
2748 MPI_Comm
2749 Data::get_MPIComm() const
2750 {
2751 #ifdef PASO_MPI
2752 return MPI_COMM_WORLD;
2753 #else
2754 return -1;
2755 #endif
2756 }
2757
2758

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26