/[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 1888 - (show annotations)
Wed Oct 15 04:00:21 2008 UTC (11 years, 10 months ago) by jfenwick
File size: 77927 byte(s)
Branch commit.
Added more binary ops.

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 DataReady_ptr res=m_data->resolve();
578 m_data=res;
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 {
600 throw DataException("Error - Tagging not implemented for this Data type.");
601 }
602 }
603
604 Data
605 Data::oneOver() const
606 {
607 if (isLazy())
608 {
609 DataLazy* c=new DataLazy(borrowDataPtr(),RECIP);
610 return Data(c);
611 }
612 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
613 }
614
615 Data
616 Data::wherePositive() const
617 {
618 if (isLazy())
619 {
620 DataLazy* c=new DataLazy(borrowDataPtr(),GZ);
621 return Data(c);
622 }
623 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
624 }
625
626 Data
627 Data::whereNegative() const
628 {
629 if (isLazy())
630 {
631 DataLazy* c=new DataLazy(borrowDataPtr(),LZ);
632 return Data(c);
633 }
634 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
635 }
636
637 Data
638 Data::whereNonNegative() const
639 {
640 if (isLazy())
641 {
642 DataLazy* c=new DataLazy(borrowDataPtr(),GEZ);
643 return Data(c);
644 }
645 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
646 }
647
648 Data
649 Data::whereNonPositive() const
650 {
651 if (isLazy())
652 {
653 DataLazy* c=new DataLazy(borrowDataPtr(),LEZ);
654 return Data(c);
655 }
656 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
657 }
658
659 Data
660 Data::whereZero(double tol) const
661 {
662 Data dataAbs=abs();
663 return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
664 }
665
666 Data
667 Data::whereNonZero(double tol) const
668 {
669 Data dataAbs=abs();
670 return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
671 }
672
673 Data
674 Data::interpolate(const FunctionSpace& functionspace) const
675 {
676 return Data(*this,functionspace);
677 }
678
679 bool
680 Data::probeInterpolation(const FunctionSpace& functionspace) const
681 {
682 if (getFunctionSpace()==functionspace) {
683 return true;
684 } else {
685 const_Domain_ptr domain=getDomain();
686 if (*domain==*functionspace.getDomain()) {
687 return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
688 } else {
689 return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
690 }
691 }
692 }
693
694 Data
695 Data::gradOn(const FunctionSpace& functionspace) const
696 {
697 if (isEmpty())
698 {
699 throw DataException("Error - operation not permitted on instances of DataEmpty.");
700 }
701 double blocktimer_start = blocktimer_time();
702 if (functionspace.getDomain()!=getDomain())
703 throw DataException("Error - gradient cannot be calculated on different domains.");
704 DataTypes::ShapeType grad_shape=getDataPointShape();
705 grad_shape.push_back(functionspace.getDim());
706 Data out(0.0,grad_shape,functionspace,true);
707 getDomain()->setToGradient(out,*this);
708 blocktimer_increment("grad()", blocktimer_start);
709 return out;
710 }
711
712 Data
713 Data::grad() const
714 {
715 if (isEmpty())
716 {
717 throw DataException("Error - operation not permitted on instances of DataEmpty.");
718 }
719 return gradOn(escript::function(*getDomain()));
720 }
721
722 int
723 Data::getDataPointSize() const
724 {
725 return m_data->getNoValues();
726 }
727
728 DataTypes::ValueType::size_type
729 Data::getLength() const
730 {
731 return m_data->getLength();
732 }
733
734 const
735 boost::python::numeric::array
736 Data:: getValueOfDataPoint(int dataPointNo)
737 {
738 size_t length=0;
739 int i, j, k, l;
740 //
741 // determine the rank and shape of each data point
742 int dataPointRank = getDataPointRank();
743 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
744
745 //
746 // create the numeric array to be returned
747 boost::python::numeric::array numArray(0.0);
748
749 //
750 // the shape of the returned numeric array will be the same
751 // as that of the data point
752 int arrayRank = dataPointRank;
753 const DataTypes::ShapeType& arrayShape = dataPointShape;
754
755 //
756 // resize the numeric array to the shape just calculated
757 if (arrayRank==0) {
758 numArray.resize(1);
759 }
760 if (arrayRank==1) {
761 numArray.resize(arrayShape[0]);
762 }
763 if (arrayRank==2) {
764 numArray.resize(arrayShape[0],arrayShape[1]);
765 }
766 if (arrayRank==3) {
767 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
768 }
769 if (arrayRank==4) {
770 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
771 }
772
773 if (getNumDataPointsPerSample()>0) {
774 int sampleNo = dataPointNo/getNumDataPointsPerSample();
775 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
776 //
777 // Check a valid sample number has been supplied
778 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
779 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
780 }
781
782 //
783 // Check a valid data point number has been supplied
784 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
785 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
786 }
787 // TODO: global error handling
788 // create a view of the data if it is stored locally
789 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
790 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
791
792
793 switch( dataPointRank ){
794 case 0 :
795 numArray[0] = getDataAtOffset(offset);
796 break;
797 case 1 :
798 for( i=0; i<dataPointShape[0]; i++ )
799 numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
800 break;
801 case 2 :
802 for( i=0; i<dataPointShape[0]; i++ )
803 for( j=0; j<dataPointShape[1]; j++)
804 numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
805 break;
806 case 3 :
807 for( i=0; i<dataPointShape[0]; i++ )
808 for( j=0; j<dataPointShape[1]; j++ )
809 for( k=0; k<dataPointShape[2]; k++)
810 numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
811 break;
812 case 4 :
813 for( i=0; i<dataPointShape[0]; i++ )
814 for( j=0; j<dataPointShape[1]; j++ )
815 for( k=0; k<dataPointShape[2]; k++ )
816 for( l=0; l<dataPointShape[3]; l++)
817 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
818 break;
819 }
820 }
821 //
822 // return the array
823 return numArray;
824
825 }
826
827 void
828 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
829 {
830 // this will throw if the value cannot be represented
831 boost::python::numeric::array num_array(py_object);
832 setValueOfDataPointToArray(dataPointNo,num_array);
833
834
835 }
836
837 void
838 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
839 {
840 if (isProtected()) {
841 throw DataException("Error - attempt to update protected Data object.");
842 }
843 //
844 // check rank
845 if (num_array.getrank()<getDataPointRank())
846 throw DataException("Rank of numarray does not match Data object rank");
847
848 //
849 // check shape of num_array
850 for (int i=0; i<getDataPointRank(); i++) {
851 if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
852 throw DataException("Shape of numarray does not match Data object rank");
853 }
854 //
855 // make sure data is expanded:
856 if (!isExpanded()) {
857 expand();
858 }
859 if (getNumDataPointsPerSample()>0) {
860 int sampleNo = dataPointNo/getNumDataPointsPerSample();
861 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
862 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
863 } else {
864 m_data->copyToDataPoint(-1, 0,num_array);
865 }
866 }
867
868 void
869 Data::setValueOfDataPoint(int dataPointNo, const double value)
870 {
871 if (isProtected()) {
872 throw DataException("Error - attempt to update protected Data object.");
873 }
874 //
875 // make sure data is expanded:
876 if (!isExpanded()) {
877 expand();
878 }
879 if (getNumDataPointsPerSample()>0) {
880 int sampleNo = dataPointNo/getNumDataPointsPerSample();
881 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
882 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
883 } else {
884 m_data->copyToDataPoint(-1, 0,value);
885 }
886 }
887
888 const
889 boost::python::numeric::array
890 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
891 {
892 size_t length=0;
893 int i, j, k, l, pos;
894 //
895 // determine the rank and shape of each data point
896 int dataPointRank = getDataPointRank();
897 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
898
899 //
900 // create the numeric array to be returned
901 boost::python::numeric::array numArray(0.0);
902
903 //
904 // the shape of the returned numeric array will be the same
905 // as that of the data point
906 int arrayRank = dataPointRank;
907 const DataTypes::ShapeType& arrayShape = dataPointShape;
908
909 //
910 // resize the numeric array to the shape just calculated
911 if (arrayRank==0) {
912 numArray.resize(1);
913 }
914 if (arrayRank==1) {
915 numArray.resize(arrayShape[0]);
916 }
917 if (arrayRank==2) {
918 numArray.resize(arrayShape[0],arrayShape[1]);
919 }
920 if (arrayRank==3) {
921 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
922 }
923 if (arrayRank==4) {
924 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
925 }
926
927 // added for the MPI communication
928 length=1;
929 for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
930 double *tmpData = new double[length];
931
932 //
933 // load the values for the data point into the numeric array.
934
935 // updated for the MPI case
936 if( get_MPIRank()==procNo ){
937 if (getNumDataPointsPerSample()>0) {
938 int sampleNo = dataPointNo/getNumDataPointsPerSample();
939 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
940 //
941 // Check a valid sample number has been supplied
942 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
943 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
944 }
945
946 //
947 // Check a valid data point number has been supplied
948 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
949 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
950 }
951 // TODO: global error handling
952 // create a view of the data if it is stored locally
953 //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
954 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
955
956 // pack the data from the view into tmpData for MPI communication
957 pos=0;
958 switch( dataPointRank ){
959 case 0 :
960 tmpData[0] = getDataAtOffset(offset);
961 break;
962 case 1 :
963 for( i=0; i<dataPointShape[0]; i++ )
964 tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
965 break;
966 case 2 :
967 for( i=0; i<dataPointShape[0]; i++ )
968 for( j=0; j<dataPointShape[1]; j++, pos++ )
969 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
970 break;
971 case 3 :
972 for( i=0; i<dataPointShape[0]; i++ )
973 for( j=0; j<dataPointShape[1]; j++ )
974 for( k=0; k<dataPointShape[2]; k++, pos++ )
975 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
976 break;
977 case 4 :
978 for( i=0; i<dataPointShape[0]; i++ )
979 for( j=0; j<dataPointShape[1]; j++ )
980 for( k=0; k<dataPointShape[2]; k++ )
981 for( l=0; l<dataPointShape[3]; l++, pos++ )
982 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
983 break;
984 }
985 }
986 }
987 #ifdef PASO_MPI
988 // broadcast the data to all other processes
989 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
990 #endif
991
992 // unpack the data
993 switch( dataPointRank ){
994 case 0 :
995 numArray[0]=tmpData[0];
996 break;
997 case 1 :
998 for( i=0; i<dataPointShape[0]; i++ )
999 numArray[i]=tmpData[i];
1000 break;
1001 case 2 :
1002 for( i=0; i<dataPointShape[0]; i++ )
1003 for( j=0; j<dataPointShape[1]; j++ )
1004 numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1005 break;
1006 case 3 :
1007 for( i=0; i<dataPointShape[0]; i++ )
1008 for( j=0; j<dataPointShape[1]; j++ )
1009 for( k=0; k<dataPointShape[2]; k++ )
1010 numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1011 break;
1012 case 4 :
1013 for( i=0; i<dataPointShape[0]; i++ )
1014 for( j=0; j<dataPointShape[1]; j++ )
1015 for( k=0; k<dataPointShape[2]; k++ )
1016 for( l=0; l<dataPointShape[3]; l++ )
1017 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1018 break;
1019 }
1020
1021 delete [] tmpData;
1022 //
1023 // return the loaded array
1024 return numArray;
1025 }
1026
1027
1028
1029 boost::python::numeric::array
1030 Data::integrate() const
1031 {
1032 int index;
1033 int rank = getDataPointRank();
1034 DataTypes::ShapeType shape = getDataPointShape();
1035 int dataPointSize = getDataPointSize();
1036
1037 //
1038 // calculate the integral values
1039 vector<double> integrals(dataPointSize);
1040 vector<double> integrals_local(dataPointSize);
1041 #ifdef PASO_MPI
1042 AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals_local,*this);
1043 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1044 double *tmp = new double[dataPointSize];
1045 double *tmp_local = new double[dataPointSize];
1046 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1047 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1048 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1049 delete[] tmp;
1050 delete[] tmp_local;
1051 #else
1052 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1053 #endif
1054
1055 //
1056 // create the numeric array to be returned
1057 // and load the array with the integral values
1058 boost::python::numeric::array bp_array(1.0);
1059 if (rank==0) {
1060 bp_array.resize(1);
1061 index = 0;
1062 bp_array[0] = integrals[index];
1063 }
1064 if (rank==1) {
1065 bp_array.resize(shape[0]);
1066 for (int i=0; i<shape[0]; i++) {
1067 index = i;
1068 bp_array[i] = integrals[index];
1069 }
1070 }
1071 if (rank==2) {
1072 bp_array.resize(shape[0],shape[1]);
1073 for (int i=0; i<shape[0]; i++) {
1074 for (int j=0; j<shape[1]; j++) {
1075 index = i + shape[0] * j;
1076 bp_array[make_tuple(i,j)] = integrals[index];
1077 }
1078 }
1079 }
1080 if (rank==3) {
1081 bp_array.resize(shape[0],shape[1],shape[2]);
1082 for (int i=0; i<shape[0]; i++) {
1083 for (int j=0; j<shape[1]; j++) {
1084 for (int k=0; k<shape[2]; k++) {
1085 index = i + shape[0] * ( j + shape[1] * k );
1086 bp_array[make_tuple(i,j,k)] = integrals[index];
1087 }
1088 }
1089 }
1090 }
1091 if (rank==4) {
1092 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1093 for (int i=0; i<shape[0]; i++) {
1094 for (int j=0; j<shape[1]; j++) {
1095 for (int k=0; k<shape[2]; k++) {
1096 for (int l=0; l<shape[3]; l++) {
1097 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1098 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1099 }
1100 }
1101 }
1102 }
1103 }
1104
1105 //
1106 // return the loaded array
1107 return bp_array;
1108 }
1109
1110 Data
1111 Data::sin() const
1112 {
1113 if (isLazy())
1114 {
1115 DataLazy* c=new DataLazy(borrowDataPtr(),SIN);
1116 return Data(c);
1117 }
1118 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1119 }
1120
1121 Data
1122 Data::cos() const
1123 {
1124 if (isLazy())
1125 {
1126 DataLazy* c=new DataLazy(borrowDataPtr(),COS);
1127 return Data(c);
1128 }
1129 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1130 }
1131
1132 Data
1133 Data::tan() const
1134 {
1135 if (isLazy())
1136 {
1137 DataLazy* c=new DataLazy(borrowDataPtr(),TAN);
1138 return Data(c);
1139 }
1140 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1141 }
1142
1143 Data
1144 Data::asin() const
1145 {
1146 if (isLazy())
1147 {
1148 DataLazy* c=new DataLazy(borrowDataPtr(),ASIN);
1149 return Data(c);
1150 }
1151 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1152 }
1153
1154 Data
1155 Data::acos() const
1156 {
1157 if (isLazy())
1158 {
1159 DataLazy* c=new DataLazy(borrowDataPtr(),ACOS);
1160 return Data(c);
1161 }
1162 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1163 }
1164
1165
1166 Data
1167 Data::atan() const
1168 {
1169 if (isLazy())
1170 {
1171 DataLazy* c=new DataLazy(borrowDataPtr(),ATAN);
1172 return Data(c);
1173 }
1174 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1175 }
1176
1177 Data
1178 Data::sinh() const
1179 {
1180 if (isLazy())
1181 {
1182 DataLazy* c=new DataLazy(borrowDataPtr(),SINH);
1183 return Data(c);
1184 }
1185 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1186 }
1187
1188 Data
1189 Data::cosh() const
1190 {
1191 if (isLazy())
1192 {
1193 DataLazy* c=new DataLazy(borrowDataPtr(),COSH);
1194 return Data(c);
1195 }
1196 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1197 }
1198
1199 Data
1200 Data::tanh() const
1201 {
1202 if (isLazy())
1203 {
1204 DataLazy* c=new DataLazy(borrowDataPtr(),TANH);
1205 return Data(c);
1206 }
1207 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1208 }
1209
1210
1211 Data
1212 Data::erf() const
1213 {
1214 #ifdef _WIN32
1215 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1216 #else
1217 if (isLazy())
1218 {
1219 DataLazy* c=new DataLazy(borrowDataPtr(),ERF);
1220 return Data(c);
1221 }
1222 return C_TensorUnaryOperation(*this, ::erf);
1223 #endif
1224 }
1225
1226 Data
1227 Data::asinh() const
1228 {
1229 if (isLazy())
1230 {
1231 DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1232 return Data(c);
1233 }
1234 #ifdef _WIN32
1235 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1236 #else
1237 return C_TensorUnaryOperation(*this, ::asinh);
1238 #endif
1239 }
1240
1241 Data
1242 Data::acosh() const
1243 {
1244 if (isLazy())
1245 {
1246 DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1247 return Data(c);
1248 }
1249 #ifdef _WIN32
1250 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1251 #else
1252 return C_TensorUnaryOperation(*this, ::acosh);
1253 #endif
1254 }
1255
1256 Data
1257 Data::atanh() const
1258 {
1259 if (isLazy())
1260 {
1261 DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1262 return Data(c);
1263 }
1264 #ifdef _WIN32
1265 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1266 #else
1267 return C_TensorUnaryOperation(*this, ::atanh);
1268 #endif
1269 }
1270
1271 Data
1272 Data::log10() const
1273 { if (isLazy())
1274 {
1275 DataLazy* c=new DataLazy(borrowDataPtr(),LOG10);
1276 return Data(c);
1277 }
1278 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1279 }
1280
1281 Data
1282 Data::log() const
1283 {
1284 if (isLazy())
1285 {
1286 DataLazy* c=new DataLazy(borrowDataPtr(),LOG);
1287 return Data(c);
1288 }
1289 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1290 }
1291
1292 Data
1293 Data::sign() const
1294 {
1295 if (isLazy())
1296 {
1297 DataLazy* c=new DataLazy(borrowDataPtr(),SIGN);
1298 return Data(c);
1299 }
1300 return C_TensorUnaryOperation(*this, escript::fsign);
1301 }
1302
1303 Data
1304 Data::abs() const
1305 {
1306 if (isLazy())
1307 {
1308 DataLazy* c=new DataLazy(borrowDataPtr(),ABS);
1309 return Data(c);
1310 }
1311 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1312 }
1313
1314 Data
1315 Data::neg() const
1316 {
1317 if (isLazy())
1318 {
1319 DataLazy* c=new DataLazy(borrowDataPtr(),NEG);
1320 return Data(c);
1321 }
1322 return C_TensorUnaryOperation(*this, negate<double>());
1323 }
1324
1325 Data
1326 Data::pos() const
1327 {
1328 // not doing lazy check here is deliberate.
1329 // since a deep copy of lazy data should be cheap, I'll just let it happen now
1330 Data result;
1331 // perform a deep copy
1332 result.copy(*this);
1333 return result;
1334 }
1335
1336 Data
1337 Data::exp() const
1338 {
1339 if (isLazy())
1340 {
1341 DataLazy* c=new DataLazy(borrowDataPtr(),EXP);
1342 return Data(c);
1343 }
1344 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1345 }
1346
1347 Data
1348 Data::sqrt() const
1349 {
1350 if (isLazy())
1351 {
1352 DataLazy* c=new DataLazy(borrowDataPtr(),SQRT);
1353 return Data(c);
1354 }
1355 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1356 }
1357
1358 double
1359 Data::Lsup_const() const
1360 {
1361 if (isLazy())
1362 {
1363 throw DataException("Error - cannot compute Lsup for constant lazy data.");
1364 }
1365 return LsupWorker();
1366 }
1367
1368 double
1369 Data::Lsup()
1370 {
1371 if (isLazy())
1372 {
1373 expand();
1374 }
1375 return LsupWorker();
1376 }
1377
1378 double
1379 Data::sup_const() const
1380 {
1381 if (isLazy())
1382 {
1383 throw DataException("Error - cannot compute sup for constant lazy data.");
1384 }
1385 return supWorker();
1386 }
1387
1388 double
1389 Data::sup()
1390 {
1391 if (isLazy())
1392 {
1393 expand();
1394 }
1395 return supWorker();
1396 }
1397
1398 double
1399 Data::inf_const() const
1400 {
1401 if (isLazy())
1402 {
1403 throw DataException("Error - cannot compute inf for constant lazy data.");
1404 }
1405 return infWorker();
1406 }
1407
1408 double
1409 Data::inf()
1410 {
1411 if (isLazy())
1412 {
1413 expand();
1414 }
1415 return infWorker();
1416 }
1417
1418 double
1419 Data::LsupWorker() const
1420 {
1421 double localValue;
1422 //
1423 // set the initial absolute maximum value to zero
1424
1425 AbsMax abs_max_func;
1426 localValue = algorithm(abs_max_func,0);
1427 #ifdef PASO_MPI
1428 double globalValue;
1429 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1430 return globalValue;
1431 #else
1432 return localValue;
1433 #endif
1434 }
1435
1436 double
1437 Data::supWorker() const
1438 {
1439 double localValue;
1440 //
1441 // set the initial maximum value to min possible double
1442 FMax fmax_func;
1443 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1444 #ifdef PASO_MPI
1445 double globalValue;
1446 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1447 return globalValue;
1448 #else
1449 return localValue;
1450 #endif
1451 }
1452
1453 double
1454 Data::infWorker() const
1455 {
1456 double localValue;
1457 //
1458 // set the initial minimum value to max possible double
1459 FMin fmin_func;
1460 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1461 #ifdef PASO_MPI
1462 double globalValue;
1463 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1464 return globalValue;
1465 #else
1466 return localValue;
1467 #endif
1468 }
1469
1470 /* TODO */
1471 /* global reduction */
1472 Data
1473 Data::maxval() const
1474 {
1475 //
1476 // set the initial maximum value to min possible double
1477 FMax fmax_func;
1478 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1479 }
1480
1481 Data
1482 Data::minval() const
1483 {
1484 //
1485 // set the initial minimum value to max possible double
1486 FMin fmin_func;
1487 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1488 }
1489
1490 Data
1491 Data::swapaxes(const int axis0, const int axis1) const
1492 {
1493 int axis0_tmp,axis1_tmp;
1494 DataTypes::ShapeType s=getDataPointShape();
1495 DataTypes::ShapeType ev_shape;
1496 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1497 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1498 int rank=getDataPointRank();
1499 if (rank<2) {
1500 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1501 }
1502 if (axis0<0 || axis0>rank-1) {
1503 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1504 }
1505 if (axis1<0 || axis1>rank-1) {
1506 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1507 }
1508 if (axis0 == axis1) {
1509 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1510 }
1511 if (axis0 > axis1) {
1512 axis0_tmp=axis1;
1513 axis1_tmp=axis0;
1514 } else {
1515 axis0_tmp=axis0;
1516 axis1_tmp=axis1;
1517 }
1518 for (int i=0; i<rank; i++) {
1519 if (i == axis0_tmp) {
1520 ev_shape.push_back(s[axis1_tmp]);
1521 } else if (i == axis1_tmp) {
1522 ev_shape.push_back(s[axis0_tmp]);
1523 } else {
1524 ev_shape.push_back(s[i]);
1525 }
1526 }
1527 Data ev(0.,ev_shape,getFunctionSpace());
1528 ev.typeMatchRight(*this);
1529 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1530 return ev;
1531
1532 }
1533
1534 Data
1535 Data::symmetric() const
1536 {
1537 // check input
1538 DataTypes::ShapeType s=getDataPointShape();
1539 if (getDataPointRank()==2) {
1540 if(s[0] != s[1])
1541 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1542 }
1543 else if (getDataPointRank()==4) {
1544 if(!(s[0] == s[2] && s[1] == s[3]))
1545 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1546 }
1547 else {
1548 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1549 }
1550 Data ev(0.,getDataPointShape(),getFunctionSpace());
1551 ev.typeMatchRight(*this);
1552 m_data->symmetric(ev.m_data.get());
1553 return ev;
1554 }
1555
1556 Data
1557 Data::nonsymmetric() const
1558 {
1559 // check input
1560 DataTypes::ShapeType s=getDataPointShape();
1561 if (getDataPointRank()==2) {
1562 if(s[0] != s[1])
1563 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1564 DataTypes::ShapeType ev_shape;
1565 ev_shape.push_back(s[0]);
1566 ev_shape.push_back(s[1]);
1567 Data ev(0.,ev_shape,getFunctionSpace());
1568 ev.typeMatchRight(*this);
1569 m_data->nonsymmetric(ev.m_data.get());
1570 return ev;
1571 }
1572 else if (getDataPointRank()==4) {
1573 if(!(s[0] == s[2] && s[1] == s[3]))
1574 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1575 DataTypes::ShapeType ev_shape;
1576 ev_shape.push_back(s[0]);
1577 ev_shape.push_back(s[1]);
1578 ev_shape.push_back(s[2]);
1579 ev_shape.push_back(s[3]);
1580 Data ev(0.,ev_shape,getFunctionSpace());
1581 ev.typeMatchRight(*this);
1582 m_data->nonsymmetric(ev.m_data.get());
1583 return ev;
1584 }
1585 else {
1586 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1587 }
1588 }
1589
1590 Data
1591 Data::trace(int axis_offset) const
1592 {
1593 DataTypes::ShapeType s=getDataPointShape();
1594 if (getDataPointRank()==2) {
1595 DataTypes::ShapeType ev_shape;
1596 Data ev(0.,ev_shape,getFunctionSpace());
1597 ev.typeMatchRight(*this);
1598 m_data->trace(ev.m_data.get(), axis_offset);
1599 return ev;
1600 }
1601 if (getDataPointRank()==3) {
1602 DataTypes::ShapeType ev_shape;
1603 if (axis_offset==0) {
1604 int s2=s[2];
1605 ev_shape.push_back(s2);
1606 }
1607 else if (axis_offset==1) {
1608 int s0=s[0];
1609 ev_shape.push_back(s0);
1610 }
1611 Data ev(0.,ev_shape,getFunctionSpace());
1612 ev.typeMatchRight(*this);
1613 m_data->trace(ev.m_data.get(), axis_offset);
1614 return ev;
1615 }
1616 if (getDataPointRank()==4) {
1617 DataTypes::ShapeType ev_shape;
1618 if (axis_offset==0) {
1619 ev_shape.push_back(s[2]);
1620 ev_shape.push_back(s[3]);
1621 }
1622 else if (axis_offset==1) {
1623 ev_shape.push_back(s[0]);
1624 ev_shape.push_back(s[3]);
1625 }
1626 else if (axis_offset==2) {
1627 ev_shape.push_back(s[0]);
1628 ev_shape.push_back(s[1]);
1629 }
1630 Data ev(0.,ev_shape,getFunctionSpace());
1631 ev.typeMatchRight(*this);
1632 m_data->trace(ev.m_data.get(), axis_offset);
1633 return ev;
1634 }
1635 else {
1636 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1637 }
1638 }
1639
1640 Data
1641 Data::transpose(int axis_offset) const
1642 {
1643 DataTypes::ShapeType s=getDataPointShape();
1644 DataTypes::ShapeType ev_shape;
1645 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1646 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1647 int rank=getDataPointRank();
1648 if (axis_offset<0 || axis_offset>rank) {
1649 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1650 }
1651 for (int i=0; i<rank; i++) {
1652 int index = (axis_offset+i)%rank;
1653 ev_shape.push_back(s[index]); // Append to new shape
1654 }
1655 Data ev(0.,ev_shape,getFunctionSpace());
1656 ev.typeMatchRight(*this);
1657 m_data->transpose(ev.m_data.get(), axis_offset);
1658 return ev;
1659 }
1660
1661 Data
1662 Data::eigenvalues() const
1663 {
1664 // check input
1665 DataTypes::ShapeType s=getDataPointShape();
1666 if (getDataPointRank()!=2)
1667 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1668 if(s[0] != s[1])
1669 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1670 // create return
1671 DataTypes::ShapeType ev_shape(1,s[0]);
1672 Data ev(0.,ev_shape,getFunctionSpace());
1673 ev.typeMatchRight(*this);
1674 m_data->eigenvalues(ev.m_data.get());
1675 return ev;
1676 }
1677
1678 const boost::python::tuple
1679 Data::eigenvalues_and_eigenvectors(const double tol) const
1680 {
1681 DataTypes::ShapeType s=getDataPointShape();
1682 if (getDataPointRank()!=2)
1683 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1684 if(s[0] != s[1])
1685 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1686 // create return
1687 DataTypes::ShapeType ev_shape(1,s[0]);
1688 Data ev(0.,ev_shape,getFunctionSpace());
1689 ev.typeMatchRight(*this);
1690 DataTypes::ShapeType V_shape(2,s[0]);
1691 Data V(0.,V_shape,getFunctionSpace());
1692 V.typeMatchRight(*this);
1693 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1694 return make_tuple(boost::python::object(ev),boost::python::object(V));
1695 }
1696
1697 const boost::python::tuple
1698 Data::minGlobalDataPoint() const
1699 {
1700 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1701 // abort (for unknown reasons) if there are openmp directives with it in the
1702 // surrounding function
1703
1704 int DataPointNo;
1705 int ProcNo;
1706 calc_minGlobalDataPoint(ProcNo,DataPointNo);
1707 return make_tuple(ProcNo,DataPointNo);
1708 }
1709
1710 void
1711 Data::calc_minGlobalDataPoint(int& ProcNo,
1712 int& DataPointNo) const
1713 {
1714 int i,j;
1715 int lowi=0,lowj=0;
1716 double min=numeric_limits<double>::max();
1717
1718 Data temp=minval();
1719
1720 int numSamples=temp.getNumSamples();
1721 int numDPPSample=temp.getNumDataPointsPerSample();
1722
1723 double next,local_min;
1724 int local_lowi,local_lowj;
1725
1726 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1727 {
1728 local_min=min;
1729 #pragma omp for private(i,j) schedule(static)
1730 for (i=0; i<numSamples; i++) {
1731 for (j=0; j<numDPPSample; j++) {
1732 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1733 if (next<local_min) {
1734 local_min=next;
1735 local_lowi=i;
1736 local_lowj=j;
1737 }
1738 }
1739 }
1740 #pragma omp critical
1741 if (local_min<min) {
1742 min=local_min;
1743 lowi=local_lowi;
1744 lowj=local_lowj;
1745 }
1746 }
1747
1748 #ifdef PASO_MPI
1749 // determine the processor on which the minimum occurs
1750 next = temp.getDataPoint(lowi,lowj);
1751 int lowProc = 0;
1752 double *globalMins = new double[get_MPISize()+1];
1753 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1754
1755 if( get_MPIRank()==0 ){
1756 next = globalMins[lowProc];
1757 for( i=1; i<get_MPISize(); i++ )
1758 if( next>globalMins[i] ){
1759 lowProc = i;
1760 next = globalMins[i];
1761 }
1762 }
1763 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1764
1765 delete [] globalMins;
1766 ProcNo = lowProc;
1767 #else
1768 ProcNo = 0;
1769 #endif
1770 DataPointNo = lowj + lowi * numDPPSample;
1771 }
1772
1773 void
1774 Data::saveDX(std::string fileName) const
1775 {
1776 if (isEmpty())
1777 {
1778 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1779 }
1780 boost::python::dict args;
1781 args["data"]=boost::python::object(this);
1782 getDomain()->saveDX(fileName,args);
1783 return;
1784 }
1785
1786 void
1787 Data::saveVTK(std::string fileName) const
1788 {
1789 if (isEmpty())
1790 {
1791 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1792 }
1793 boost::python::dict args;
1794 args["data"]=boost::python::object(this);
1795 getDomain()->saveVTK(fileName,args);
1796 return;
1797 }
1798
1799 Data&
1800 Data::operator+=(const Data& right)
1801 {
1802 if (isProtected()) {
1803 throw DataException("Error - attempt to update protected Data object.");
1804 }
1805 binaryOp(right,plus<double>());
1806 return (*this);
1807 }
1808
1809 Data&
1810 Data::operator+=(const boost::python::object& right)
1811 {
1812 Data tmp(right,getFunctionSpace(),false);
1813 binaryOp(tmp,plus<double>());
1814 return (*this);
1815 }
1816 Data&
1817 Data::operator=(const Data& other)
1818 {
1819 copy(other);
1820 return (*this);
1821 }
1822
1823 Data&
1824 Data::operator-=(const Data& right)
1825 {
1826 if (isProtected()) {
1827 throw DataException("Error - attempt to update protected Data object.");
1828 }
1829 binaryOp(right,minus<double>());
1830 return (*this);
1831 }
1832
1833 Data&
1834 Data::operator-=(const boost::python::object& right)
1835 {
1836 Data tmp(right,getFunctionSpace(),false);
1837 binaryOp(tmp,minus<double>());
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,multiplies<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,multiplies<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,divides<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,divides<double>());
1874 return (*this);
1875 }
1876
1877 Data
1878 Data::rpowO(const boost::python::object& left) const
1879 {
1880 Data left_d(left,*this);
1881 return left_d.powD(*this);
1882 }
1883
1884 Data
1885 Data::powO(const boost::python::object& right) const
1886 {
1887 Data tmp(right,getFunctionSpace(),false);
1888 return powD(tmp);
1889 }
1890
1891 Data
1892 Data::powD(const Data& right) const
1893 {
1894 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1895 }
1896
1897 //
1898 // NOTE: It is essential to specify the namespace this operator belongs to
1899 Data
1900 escript::operator+(const Data& left, const Data& right)
1901 {
1902 if (left.isLazy() || right.isLazy())
1903 {
1904 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),ADD);
1905 return Data(c);
1906 }
1907 return C_TensorBinaryOperation(left, right, plus<double>());
1908 }
1909
1910 //
1911 // NOTE: It is essential to specify the namespace this operator belongs to
1912 Data
1913 escript::operator-(const Data& left, const Data& right)
1914 {
1915 if (left.isLazy() || right.isLazy())
1916 {
1917 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),SUB);
1918 return Data(c);
1919 }
1920 return C_TensorBinaryOperation(left, right, minus<double>());
1921 }
1922
1923 //
1924 // NOTE: It is essential to specify the namespace this operator belongs to
1925 Data
1926 escript::operator*(const Data& left, const Data& right)
1927 {
1928 if (left.isLazy() || right.isLazy())
1929 {
1930 DataLazy* c=new DataLazy(left.borrowDataPtr(),right.borrowDataPtr(),MUL);
1931 return Data(c);
1932 }
1933 return C_TensorBinaryOperation(left, right, multiplies<double>());
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(),DIV);
1944 return Data(c);
1945 }
1946 return C_TensorBinaryOperation(left, right, divides<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 boost::python::object& right)
1953 {
1954 return left+Data(right,left.getFunctionSpace(),false);
1955 }
1956
1957 //
1958 // NOTE: It is essential to specify the namespace this operator belongs to
1959 Data
1960 escript::operator-(const Data& left, const boost::python::object& right)
1961 {
1962 return left-Data(right,left.getFunctionSpace(),false);
1963 }
1964
1965 //
1966 // NOTE: It is essential to specify the namespace this operator belongs to
1967 Data
1968 escript::operator*(const Data& left, const boost::python::object& right)
1969 {
1970 return left*Data(right,left.getFunctionSpace(),false);
1971 }
1972
1973 //
1974 // NOTE: It is essential to specify the namespace this operator belongs to
1975 Data
1976 escript::operator/(const Data& left, const boost::python::object& right)
1977 {
1978 return left/Data(right,left.getFunctionSpace(),false);
1979 }
1980
1981 //
1982 // NOTE: It is essential to specify the namespace this operator belongs to
1983 Data
1984 escript::operator+(const boost::python::object& left, const Data& right)
1985 {
1986 return Data(left,right.getFunctionSpace(),false)+right;
1987 }
1988
1989 //
1990 // NOTE: It is essential to specify the namespace this operator belongs to
1991 Data
1992 escript::operator-(const boost::python::object& left, const Data& right)
1993 {
1994 return Data(left,right.getFunctionSpace(),false)-right;
1995 }
1996
1997 //
1998 // NOTE: It is essential to specify the namespace this operator belongs to
1999 Data
2000 escript::operator*(const boost::python::object& left, const Data& right)
2001 {
2002 return Data(left,right.getFunctionSpace(),false)*right;
2003 }
2004
2005 //
2006 // NOTE: It is essential to specify the namespace this operator belongs to
2007 Data
2008 escript::operator/(const boost::python::object& left, const Data& right)
2009 {
2010 return Data(left,right.getFunctionSpace(),false)/right;
2011 }
2012
2013
2014 /* TODO */
2015 /* global reduction */
2016 Data
2017 Data::getItem(const boost::python::object& key) const
2018 {
2019 // const DataArrayView& view=getPointDataView();
2020
2021 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2022
2023 if (slice_region.size()!=getDataPointRank()) {
2024 throw DataException("Error - slice size does not match Data rank.");
2025 }
2026
2027 return getSlice(slice_region);
2028 }
2029
2030 /* TODO */
2031 /* global reduction */
2032 Data
2033 Data::getSlice(const DataTypes::RegionType& region) const
2034 {
2035 return Data(*this,region);
2036 }
2037
2038 /* TODO */
2039 /* global reduction */
2040 void
2041 Data::setItemO(const boost::python::object& key,
2042 const boost::python::object& value)
2043 {
2044 Data tempData(value,getFunctionSpace());
2045 setItemD(key,tempData);
2046 }
2047
2048 void
2049 Data::setItemD(const boost::python::object& key,
2050 const Data& value)
2051 {
2052 // const DataArrayView& view=getPointDataView();
2053
2054 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2055 if (slice_region.size()!=getDataPointRank()) {
2056 throw DataException("Error - slice size does not match Data rank.");
2057 }
2058 if (getFunctionSpace()!=value.getFunctionSpace()) {
2059 setSlice(Data(value,getFunctionSpace()),slice_region);
2060 } else {
2061 setSlice(value,slice_region);
2062 }
2063 }
2064
2065 void
2066 Data::setSlice(const Data& value,
2067 const DataTypes::RegionType& region)
2068 {
2069 if (isProtected()) {
2070 throw DataException("Error - attempt to update protected Data object.");
2071 }
2072 if (isLazy())
2073 {
2074 throw DataException("Error - setSlice not permitted on lazy data.");
2075 }
2076 Data tempValue(value);
2077 typeMatchLeft(tempValue);
2078 typeMatchRight(tempValue);
2079 getReady()->setSlice(tempValue.m_data.get(),region);
2080 }
2081
2082 void
2083 Data::typeMatchLeft(Data& right) const
2084 {
2085 if (isExpanded()){
2086 right.expand();
2087 } else if (isTagged()) {
2088 if (right.isConstant()) {
2089 right.tag();
2090 }
2091 }
2092 }
2093
2094 void
2095 Data::typeMatchRight(const Data& right)
2096 {
2097 if (isTagged()) {
2098 if (right.isExpanded()) {
2099 expand();
2100 }
2101 } else if (isConstant()) {
2102 if (right.isExpanded()) {
2103 expand();
2104 } else if (right.isTagged()) {
2105 tag();
2106 }
2107 }
2108 }
2109
2110 void
2111 Data::setTaggedValueByName(std::string name,
2112 const boost::python::object& value)
2113 {
2114 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2115 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2116 setTaggedValue(tagKey,value);
2117 }
2118 }
2119 void
2120 Data::setTaggedValue(int tagKey,
2121 const boost::python::object& value)
2122 {
2123 if (isProtected()) {
2124 throw DataException("Error - attempt to update protected Data object.");
2125 }
2126 //
2127 // Ensure underlying data object is of type DataTagged
2128 if (isConstant()) tag();
2129
2130 numeric::array asNumArray(value);
2131
2132
2133 // extract the shape of the numarray
2134 DataTypes::ShapeType tempShape;
2135 for (int i=0; i < asNumArray.getrank(); i++) {
2136 tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2137 }
2138
2139 // get the space for the data vector
2140 // int len = DataTypes::noValues(tempShape);
2141 // DataVector temp_data(len, 0.0, len);
2142 // DataArrayView temp_dataView(temp_data, tempShape);
2143 // temp_dataView.copy(asNumArray);
2144
2145 DataVector temp_data2;
2146 temp_data2.copyFromNumArray(asNumArray);
2147
2148 //
2149 // Call DataAbstract::setTaggedValue
2150 //m_data->setTaggedValue(tagKey,temp_dataView);
2151
2152 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2153 }
2154
2155
2156 void
2157 Data::setTaggedValueFromCPP(int tagKey,
2158 const DataTypes::ShapeType& pointshape,
2159 const DataTypes::ValueType& value,
2160 int dataOffset)
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 //
2170 // Call DataAbstract::setTaggedValue
2171 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2172 }
2173
2174 int
2175 Data::getTagNumber(int dpno)
2176 {
2177 if (isEmpty())
2178 {
2179 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2180 }
2181 return getFunctionSpace().getTagFromDataPointNo(dpno);
2182 }
2183
2184
2185 ostream& escript::operator<<(ostream& o, const Data& data)
2186 {
2187 o << data.toString();
2188 return o;
2189 }
2190
2191 Data
2192 escript::C_GeneralTensorProduct(Data& arg_0,
2193 Data& arg_1,
2194 int axis_offset,
2195 int transpose)
2196 {
2197 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2198 // SM is the product of the last axis_offset entries in arg_0.getShape().
2199
2200 // Interpolate if necessary and find an appropriate function space
2201 Data arg_0_Z, arg_1_Z;
2202 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2203 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2204 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2205 arg_1_Z = Data(arg_1);
2206 }
2207 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2208 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2209 arg_0_Z =Data(arg_0);
2210 }
2211 else {
2212 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2213 }
2214 } else {
2215 arg_0_Z = Data(arg_0);
2216 arg_1_Z = Data(arg_1);
2217 }
2218 // Get rank and shape of inputs
2219 int rank0 = arg_0_Z.getDataPointRank();
2220 int rank1 = arg_1_Z.getDataPointRank();
2221 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2222 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2223
2224 // Prepare for the loops of the product and verify compatibility of shapes
2225 int start0=0, start1=0;
2226 if (transpose == 0) {}
2227 else if (transpose == 1) { start0 = axis_offset; }
2228 else if (transpose == 2) { start1 = rank1-axis_offset; }
2229 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2230
2231
2232 // Adjust the shapes for transpose
2233 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2234 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2235 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2236 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2237
2238 #if 0
2239 // For debugging: show shape after transpose
2240 char tmp[100];
2241 std::string shapeStr;
2242 shapeStr = "(";
2243 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2244 shapeStr += ")";
2245 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2246 shapeStr = "(";
2247 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2248 shapeStr += ")";
2249 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2250 #endif
2251
2252 // Prepare for the loops of the product
2253 int SL=1, SM=1, SR=1;
2254 for (int i=0; i<rank0-axis_offset; i++) {
2255 SL *= tmpShape0[i];
2256 }
2257 for (int i=rank0-axis_offset; i<rank0; i++) {
2258 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2259 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2260 }
2261 SM *= tmpShape0[i];
2262 }
2263 for (int i=axis_offset; i<rank1; i++) {
2264 SR *= tmpShape1[i];
2265 }
2266
2267 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2268 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2269 { // block to limit the scope of out_index
2270 int out_index=0;
2271 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2272 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2273 }
2274
2275 // Declare output Data object
2276 Data res;
2277
2278 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2279 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2280 double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2281 double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2282 double *ptr_2 = &(res.getDataAtOffset(0));
2283 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2284 }
2285 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2286
2287 // Prepare the DataConstant input
2288 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2289 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2290
2291 // Borrow DataTagged input from Data object
2292 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2293 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2294
2295 // Prepare a DataTagged output 2
2296 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2297 res.tag();
2298 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2299 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2300
2301 // Prepare offset into DataConstant
2302 int offset_0 = tmp_0->getPointOffset(0,0);
2303 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2304 // Get the views
2305 // DataArrayView view_1 = tmp_1->getDefaultValue();
2306 // DataArrayView view_2 = tmp_2->getDefaultValue();
2307 // // Get the pointers to the actual data
2308 // double *ptr_1 = &((view_1.getData())[0]);
2309 // double *ptr_2 = &((view_2.getData())[0]);
2310
2311 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2312 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2313
2314
2315 // Compute an MVP for the default
2316 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2317 // Compute an MVP for each tag
2318 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2319 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2320 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2321 tmp_2->addTag(i->first);
2322 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2323 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2324 // double *ptr_1 = &view_1.getData(0);
2325 // double *ptr_2 = &view_2.getData(0);
2326
2327 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2328 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2329
2330 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2331 }
2332
2333 }
2334 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2335
2336 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2337 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2338 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2339 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2340 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2341 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2342 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2343 int sampleNo_1,dataPointNo_1;
2344 int numSamples_1 = arg_1_Z.getNumSamples();
2345 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2346 int offset_0 = tmp_0->getPointOffset(0,0);
2347 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2348 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2349 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2350 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2351 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2352 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2353 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2354 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2355 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2356 }
2357 }
2358
2359 }
2360 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2361
2362 // Borrow DataTagged input from Data object
2363 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2364 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2365
2366 // Prepare the DataConstant input
2367 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2368 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2369
2370 // Prepare a DataTagged output 2
2371 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2372 res.tag();
2373 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2374 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2375
2376 // Prepare offset into DataConstant
2377 int offset_1 = tmp_1->getPointOffset(0,0);
2378 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2379 // Get the views
2380 // DataArrayView view_0 = tmp_0->getDefaultValue();
2381 // DataArrayView view_2 = tmp_2->getDefaultValue();
2382 // // Get the pointers to the actual data
2383 // double *ptr_0 = &((view_0.getData())[0]);
2384 // double *ptr_2 = &((view_2.getData())[0]);
2385
2386 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2387 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2388
2389 // Compute an MVP for the default
2390 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2391 // Compute an MVP for each tag
2392 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2393 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2394 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2395 // tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2396 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2397 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2398 // double *ptr_0 = &view_0.getData(0);
2399 // double *ptr_2 = &view_2.getData(0);
2400
2401 tmp_2->addTag(i->first);
2402 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2403 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2404 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2405 }
2406
2407 }
2408 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2409
2410 // Borrow DataTagged input from Data object
2411 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2412 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2413
2414 // Borrow DataTagged input from Data object
2415 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2416 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2417
2418 // Prepare a DataTagged output 2
2419 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2420 res.tag(); // DataTagged output
2421 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2422 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2423
2424 // // Get the views
2425 // DataArrayView view_0 = tmp_0->getDefaultValue();
2426 // DataArrayView view_1 = tmp_1->getDefaultValue();
2427 // DataArrayView view_2 = tmp_2->getDefaultValue();
2428 // // Get the pointers to the actual data
2429 // double *ptr_0 = &((view_0.getData())[0]);
2430 // double *ptr_1 = &((view_1.getData())[0]);
2431 // double *ptr_2 = &((view_2.getData())[0]);
2432
2433 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2434 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2435 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2436
2437
2438 // Compute an MVP for the default
2439 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2440 // Merge the tags
2441 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2442 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2443 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2444 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2445 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2446 }
2447 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2448 tmp_2->addTag(i->first);
2449 }
2450 // Compute an MVP for each tag
2451 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2452 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2453 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2454 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2455 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2456 // double *ptr_0 = &view_0.getData(0);
2457 // double *ptr_1 = &view_1.getData(0);
2458 // double *ptr_2 = &view_2.getData(0);
2459
2460 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2461 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2462 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2463
2464 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2465 }
2466
2467 }
2468 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2469
2470 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2471 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2472 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2473 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2474 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2475 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2476 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2477 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2478 int sampleNo_0,dataPointNo_0;
2479 int numSamples_0 = arg_0_Z.getNumSamples();
2480 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2481 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2482 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2483 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2484 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2485 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2486 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2487 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2488 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2489 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2490 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2491 }
2492 }
2493
2494 }
2495 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2496
2497 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2498 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2499 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2500 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2501 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2502 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2503 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2504 int sampleNo_0,dataPointNo_0;
2505 int numSamples_0 = arg_0_Z.getNumSamples();
2506 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2507 int offset_1 = tmp_1->getPointOffset(0,0);
2508 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2509 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2510 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2511 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2512 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2513 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2514 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2515 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2516 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2517 }
2518 }
2519
2520
2521 }
2522 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2523
2524 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2525 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2526 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2527 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2528 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2529 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2530 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2531 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2532 int sampleNo_0,dataPointNo_0;
2533 int numSamples_0 = arg_0_Z.getNumSamples();
2534 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2535 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2536 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2537 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2538 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2539 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2540 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2541 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2542 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2543 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2544 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2545 }
2546 }
2547
2548 }
2549 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2550
2551 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2552 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2553 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2554 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2555 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2556 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2557 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2558 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2559 int sampleNo_0,dataPointNo_0;
2560 int numSamples_0 = arg_0_Z.getNumSamples();
2561 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2562 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2563 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2564 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2565 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2566 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2567 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2568 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2569 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2570 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2571 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2572 }
2573 }
2574
2575 }
2576 else {
2577 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2578 }
2579
2580 return res;
2581 }
2582
2583 DataAbstract*
2584 Data::borrowData() const
2585 {
2586 return m_data.get();
2587 }
2588
2589 // Not all that happy about returning a non-const from a const
2590 DataAbstract_ptr
2591 Data::borrowDataPtr() const
2592 {
2593 return m_data;
2594 }
2595
2596 // Not all that happy about returning a non-const from a const
2597 DataReady_ptr
2598 Data::borrowReadyPtr() const
2599 {
2600 DataReady_ptr dr=dynamic_pointer_cast<DataReady>(m_data);
2601 EsysAssert((dr!=0), "Error - casting to DataReady.");
2602 return dr;
2603 }
2604
2605 std::string
2606 Data::toString() const
2607 {
2608 static const DataTypes::ValueType::size_type TOO_MANY_POINTS=80;
2609 if (isLazy())
2610 { // This needs to change back to printing out something useful once the summary ops
2611 return m_data->toString(); // are defined
2612 }
2613 if (getNumDataPoints()*getDataPointSize()>TOO_MANY_POINTS)
2614 {
2615 stringstream temp;
2616 temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();
2617 return temp.str();
2618 }
2619 return m_data->toString();
2620 }
2621
2622
2623
2624 DataTypes::ValueType::const_reference
2625 Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
2626 {
2627 if (isLazy())
2628 {
2629 throw DataException("getDataAtOffset not permitted on lazy data.");
2630 }
2631 return getReady()->getDataAtOffset(i);
2632 }
2633
2634
2635 DataTypes::ValueType::reference
2636 Data::getDataAtOffset(DataTypes::ValueType::size_type i)
2637 {
2638 if (isLazy())
2639 {
2640 throw DataException("getDataAtOffset not permitted on lazy data.");
2641 }
2642 return getReady()->getDataAtOffset(i);
2643 }
2644
2645 DataTypes::ValueType::const_reference
2646 Data::getDataPoint(int sampleNo, int dataPointNo) const
2647 {
2648 if (!isReady())
2649 {
2650 throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");
2651 }
2652 else
2653 {
2654 const DataReady* dr=getReady();
2655 return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
2656 }
2657 }
2658
2659
2660 DataTypes::ValueType::reference
2661 Data::getDataPoint(int sampleNo, int dataPointNo)
2662 {
2663 if (!isReady())
2664 {
2665 throw DataException("Programmer error - getDataPoint() not permitted on Lazy Data.");
2666 }
2667 else
2668 {
2669 DataReady* dr=getReady();
2670 return dr->getDataAtOffset(dr->getPointOffset(sampleNo, dataPointNo));
2671 }
2672 }
2673
2674
2675 /* Member functions specific to the MPI implementation */
2676
2677 void
2678 Data::print()
2679 {
2680 int i,j;
2681
2682 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2683 for( i=0; i<getNumSamples(); i++ )
2684 {
2685 printf( "[%6d]", i );
2686 for( j=0; j<getNumDataPointsPerSample(); j++ )
2687 printf( "\t%10.7g", (getSampleData(i))[j] );
2688 printf( "\n" );
2689 }
2690 }
2691 void
2692 Data::dump(const std::string fileName) const
2693 {
2694 try
2695 {
2696 return m_data->dump(fileName);
2697 }
2698 catch (exception& e)
2699 {
2700 cout << e.what() << endl;
2701 }
2702 }
2703
2704 int
2705 Data::get_MPISize() const
2706 {
2707 int size;
2708 #ifdef PASO_MPI
2709 int error;
2710 error = MPI_Comm_size( get_MPIComm(), &size );
2711 #else
2712 size = 1;
2713 #endif
2714 return size;
2715 }
2716
2717 int
2718 Data::get_MPIRank() const
2719 {
2720 int rank;
2721 #ifdef PASO_MPI
2722 int error;
2723 error = MPI_Comm_rank( get_MPIComm(), &rank );
2724 #else
2725 rank = 0;
2726 #endif
2727 return rank;
2728 }
2729
2730 MPI_Comm
2731 Data::get_MPIComm() const
2732 {
2733 #ifdef PASO_MPI
2734 return MPI_COMM_WORLD;
2735 #else
2736 return -1;
2737 #endif
2738 }
2739
2740

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26