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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1872 - (show annotations)
Mon Oct 13 00:18:55 2008 UTC (11 years ago) by jfenwick
File size: 75166 byte(s)
Closing the moreshared branch

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26