/[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 1935 - (show annotations)
Mon Oct 27 06:06:39 2008 UTC (11 years, 8 months ago) by jfenwick
File size: 83903 byte(s)
Branch commit
More cleanup of DataTestCase - still don't have all the LazyTests in 
there yet.
Added tests to make sure the resolve() operation does what it is 
supposed to.
Added non-constant versions of getPointOffset to DataAbstract 
classes.
Fixed a bug in deepCopy on DataLazy.
Changed setToZero to not check the type of the data.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26