/[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 1911 - (show annotations)
Thu Oct 23 04:50:35 2008 UTC (11 years, 9 months ago) by jfenwick
File size: 84109 byte(s)
Branch commit.
Added forceresolve to operations like symmetric and eigenvalues (will go 
back and implement them in Lazy later).


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26