/[escript]/trunk/escript/src/Data.cpp
ViewVC logotype

Contents of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1921 - (show annotations)
Thu Oct 23 11:32:24 2008 UTC (11 years, 10 months ago) by phornby
File size: 75284 byte(s)
Jump through hoops to meet the OPENMP 2.5 restrictions on loop variables.
It's all ready for OPENMP 3.0 improvments, and should compile cleanly on other platforms.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26