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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1855 - (show annotations)
Tue Oct 7 05:25:30 2008 UTC (11 years ago) by jfenwick
File size: 75842 byte(s)
Modified Data::copyWithMask to have a less cryptic implementation.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26