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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26