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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1897 - (show annotations)
Mon Oct 20 00:32:30 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 75172 byte(s)
Modified Data::toString() so it doesn't throw on DataEmpty.
Added setEscriptParamInt and getEscriptParamInt as free functions.
At the moment all they do is allow you to set the param TOO_MANY_LINES.
This is used to determine when printing a Data object will show you the 
points and when it will print a summary.

I've set the default value back to 80 lines.
If you need to see more lines use (in python):

setEscriptParamInt("TOO_MANY_LINES",80000)


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 long i;
545 #pragma omp parallel for private(i) schedule(static)
546 for (i=0;i<num_points;++i)
547 {
548 if (mvec[i]>0)
549 {
550 self[i]=ovec[i];
551 }
552 }
553 }
554
555
556
557 bool
558 Data::isExpanded() const
559 {
560 DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
561 return (temp!=0);
562 }
563
564 bool
565 Data::isTagged() const
566 {
567 DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
568 return (temp!=0);
569 }
570
571 bool
572 Data::isEmpty() const
573 {
574 DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
575 return (temp!=0);
576 }
577
578 bool
579 Data::isConstant() const
580 {
581 DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
582 return (temp!=0);
583 }
584
585 void
586 Data::setProtection()
587 {
588 m_protected=true;
589 }
590
591 bool
592 Data::isProtected() const
593 {
594 return m_protected;
595 }
596
597
598
599 void
600 Data::expand()
601 {
602 if (isConstant()) {
603 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
604 DataAbstract* temp=new DataExpanded(*tempDataConst);
605 // shared_ptr<DataAbstract> temp_data(temp);
606 // m_data=temp_data;
607 m_data=temp->getPtr();
608 } else if (isTagged()) {
609 DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
610 DataAbstract* temp=new DataExpanded(*tempDataTag);
611 // shared_ptr<DataAbstract> temp_data(temp);
612 // m_data=temp_data;
613 m_data=temp->getPtr();
614 } else if (isExpanded()) {
615 //
616 // do nothing
617 } else if (isEmpty()) {
618 throw DataException("Error - Expansion of DataEmpty not possible.");
619 } else {
620 throw DataException("Error - Expansion not implemented for this Data type.");
621 }
622 }
623
624 void
625 Data::tag()
626 {
627 if (isConstant()) {
628 DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
629 DataAbstract* temp=new DataTagged(*tempDataConst);
630 // shared_ptr<DataAbstract> temp_data(temp);
631 // m_data=temp_data;
632 m_data=temp->getPtr();
633 } else if (isTagged()) {
634 // do nothing
635 } else if (isExpanded()) {
636 throw DataException("Error - Creating tag data from DataExpanded not possible.");
637 } else if (isEmpty()) {
638 throw DataException("Error - Creating tag data from DataEmpty not possible.");
639 } else {
640 throw DataException("Error - Tagging not implemented for this Data type.");
641 }
642 }
643
644 Data
645 Data::oneOver() const
646 {
647 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
648 }
649
650 Data
651 Data::wherePositive() const
652 {
653 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
654 }
655
656 Data
657 Data::whereNegative() const
658 {
659 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
660 }
661
662 Data
663 Data::whereNonNegative() const
664 {
665 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
666 }
667
668 Data
669 Data::whereNonPositive() const
670 {
671 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
672 }
673
674 Data
675 Data::whereZero(double tol) const
676 {
677 Data dataAbs=abs();
678 return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
679 }
680
681 Data
682 Data::whereNonZero(double tol) const
683 {
684 Data dataAbs=abs();
685 return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
686 }
687
688 Data
689 Data::interpolate(const FunctionSpace& functionspace) const
690 {
691 return Data(*this,functionspace);
692 }
693
694 bool
695 Data::probeInterpolation(const FunctionSpace& functionspace) const
696 {
697 if (getFunctionSpace()==functionspace) {
698 return true;
699 } else {
700 const_Domain_ptr domain=getDomain();
701 if (*domain==*functionspace.getDomain()) {
702 return domain->probeInterpolationOnDomain(getFunctionSpace().getTypeCode(),functionspace.getTypeCode());
703 } else {
704 return domain->probeInterpolationACross(getFunctionSpace().getTypeCode(),*(functionspace.getDomain()),functionspace.getTypeCode());
705 }
706 }
707 }
708
709 Data
710 Data::gradOn(const FunctionSpace& functionspace) const
711 {
712 if (isEmpty())
713 {
714 throw DataException("Error - operation not permitted on instances of DataEmpty.");
715 }
716 double blocktimer_start = blocktimer_time();
717 if (functionspace.getDomain()!=getDomain())
718 throw DataException("Error - gradient cannot be calculated on different domains.");
719 DataTypes::ShapeType grad_shape=getDataPointShape();
720 grad_shape.push_back(functionspace.getDim());
721 Data out(0.0,grad_shape,functionspace,true);
722 getDomain()->setToGradient(out,*this);
723 blocktimer_increment("grad()", blocktimer_start);
724 return out;
725 }
726
727 Data
728 Data::grad() const
729 {
730 if (isEmpty())
731 {
732 throw DataException("Error - operation not permitted on instances of DataEmpty.");
733 }
734 return gradOn(escript::function(*getDomain()));
735 }
736
737 int
738 Data::getDataPointSize() const
739 {
740 return m_data->getNoValues();
741 }
742
743 DataTypes::ValueType::size_type
744 Data::getLength() const
745 {
746 return m_data->getLength();
747 }
748
749 const
750 boost::python::numeric::array
751 Data:: getValueOfDataPoint(int dataPointNo)
752 {
753 size_t length=0;
754 int i, j, k, l;
755 //
756 // determine the rank and shape of each data point
757 int dataPointRank = getDataPointRank();
758 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
759
760 //
761 // create the numeric array to be returned
762 boost::python::numeric::array numArray(0.0);
763
764 //
765 // the shape of the returned numeric array will be the same
766 // as that of the data point
767 int arrayRank = dataPointRank;
768 const DataTypes::ShapeType& arrayShape = dataPointShape;
769
770 //
771 // resize the numeric array to the shape just calculated
772 if (arrayRank==0) {
773 numArray.resize(1);
774 }
775 if (arrayRank==1) {
776 numArray.resize(arrayShape[0]);
777 }
778 if (arrayRank==2) {
779 numArray.resize(arrayShape[0],arrayShape[1]);
780 }
781 if (arrayRank==3) {
782 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
783 }
784 if (arrayRank==4) {
785 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
786 }
787
788 if (getNumDataPointsPerSample()>0) {
789 int sampleNo = dataPointNo/getNumDataPointsPerSample();
790 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
791 //
792 // Check a valid sample number has been supplied
793 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
794 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
795 }
796
797 //
798 // Check a valid data point number has been supplied
799 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
800 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
801 }
802 // TODO: global error handling
803 // create a view of the data if it is stored locally
804 // DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
805 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
806
807
808 switch( dataPointRank ){
809 case 0 :
810 numArray[0] = getDataAtOffset(offset);
811 break;
812 case 1 :
813 for( i=0; i<dataPointShape[0]; i++ )
814 numArray[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
815 break;
816 case 2 :
817 for( i=0; i<dataPointShape[0]; i++ )
818 for( j=0; j<dataPointShape[1]; j++)
819 numArray[make_tuple(i,j)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
820 break;
821 case 3 :
822 for( i=0; i<dataPointShape[0]; i++ )
823 for( j=0; j<dataPointShape[1]; j++ )
824 for( k=0; k<dataPointShape[2]; k++)
825 numArray[make_tuple(i,j,k)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
826 break;
827 case 4 :
828 for( i=0; i<dataPointShape[0]; i++ )
829 for( j=0; j<dataPointShape[1]; j++ )
830 for( k=0; k<dataPointShape[2]; k++ )
831 for( l=0; l<dataPointShape[3]; l++)
832 numArray[make_tuple(i,j,k,l)]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
833 break;
834 }
835 }
836 //
837 // return the array
838 return numArray;
839
840 }
841
842 void
843 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
844 {
845 // this will throw if the value cannot be represented
846 boost::python::numeric::array num_array(py_object);
847 setValueOfDataPointToArray(dataPointNo,num_array);
848 }
849
850 void
851 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
852 {
853 if (isProtected()) {
854 throw DataException("Error - attempt to update protected Data object.");
855 }
856 //
857 // check rank
858 if (num_array.getrank()<getDataPointRank())
859 throw DataException("Rank of numarray does not match Data object rank");
860
861 //
862 // check shape of num_array
863 for (int i=0; i<getDataPointRank(); i++) {
864 if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
865 throw DataException("Shape of numarray does not match Data object rank");
866 }
867 //
868 // make sure data is expanded:
869 //
870 if (!isExpanded()) {
871 expand();
872 }
873 if (getNumDataPointsPerSample()>0) {
874 int sampleNo = dataPointNo/getNumDataPointsPerSample();
875 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
876 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
877 } else {
878 m_data->copyToDataPoint(-1, 0,num_array);
879 }
880 }
881
882 void
883 Data::setValueOfDataPoint(int dataPointNo, const double value)
884 {
885 if (isProtected()) {
886 throw DataException("Error - attempt to update protected Data object.");
887 }
888 //
889 // make sure data is expanded:
890 if (!isExpanded()) {
891 expand();
892 }
893 if (getNumDataPointsPerSample()>0) {
894 int sampleNo = dataPointNo/getNumDataPointsPerSample();
895 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
896 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
897 } else {
898 m_data->copyToDataPoint(-1, 0,value);
899 }
900 }
901
902 const
903 boost::python::numeric::array
904 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
905 {
906 size_t length=0;
907 int i, j, k, l, pos;
908 //
909 // determine the rank and shape of each data point
910 int dataPointRank = getDataPointRank();
911 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
912
913 //
914 // create the numeric array to be returned
915 boost::python::numeric::array numArray(0.0);
916
917 //
918 // the shape of the returned numeric array will be the same
919 // as that of the data point
920 int arrayRank = dataPointRank;
921 const DataTypes::ShapeType& arrayShape = dataPointShape;
922
923 //
924 // resize the numeric array to the shape just calculated
925 if (arrayRank==0) {
926 numArray.resize(1);
927 }
928 if (arrayRank==1) {
929 numArray.resize(arrayShape[0]);
930 }
931 if (arrayRank==2) {
932 numArray.resize(arrayShape[0],arrayShape[1]);
933 }
934 if (arrayRank==3) {
935 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
936 }
937 if (arrayRank==4) {
938 numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
939 }
940
941 // added for the MPI communication
942 length=1;
943 for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
944 double *tmpData = new double[length];
945
946 //
947 // load the values for the data point into the numeric array.
948
949 // updated for the MPI case
950 if( get_MPIRank()==procNo ){
951 if (getNumDataPointsPerSample()>0) {
952 int sampleNo = dataPointNo/getNumDataPointsPerSample();
953 int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
954 //
955 // Check a valid sample number has been supplied
956 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
957 throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
958 }
959
960 //
961 // Check a valid data point number has been supplied
962 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
963 throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
964 }
965 // TODO: global error handling
966 // create a view of the data if it is stored locally
967 //DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
968 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
969
970 // pack the data from the view into tmpData for MPI communication
971 pos=0;
972 switch( dataPointRank ){
973 case 0 :
974 tmpData[0] = getDataAtOffset(offset);
975 break;
976 case 1 :
977 for( i=0; i<dataPointShape[0]; i++ )
978 tmpData[i]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i));
979 break;
980 case 2 :
981 for( i=0; i<dataPointShape[0]; i++ )
982 for( j=0; j<dataPointShape[1]; j++, pos++ )
983 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j));
984 break;
985 case 3 :
986 for( i=0; i<dataPointShape[0]; i++ )
987 for( j=0; j<dataPointShape[1]; j++ )
988 for( k=0; k<dataPointShape[2]; k++, pos++ )
989 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k));
990 break;
991 case 4 :
992 for( i=0; i<dataPointShape[0]; i++ )
993 for( j=0; j<dataPointShape[1]; j++ )
994 for( k=0; k<dataPointShape[2]; k++ )
995 for( l=0; l<dataPointShape[3]; l++, pos++ )
996 tmpData[pos]=getDataAtOffset(offset+DataTypes::getRelIndex(dataPointShape, i,j,k,l));
997 break;
998 }
999 }
1000 }
1001 #ifdef PASO_MPI
1002 // broadcast the data to all other processes
1003 MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1004 #endif
1005
1006 // unpack the data
1007 switch( dataPointRank ){
1008 case 0 :
1009 numArray[0]=tmpData[0];
1010 break;
1011 case 1 :
1012 for( i=0; i<dataPointShape[0]; i++ )
1013 numArray[i]=tmpData[i];
1014 break;
1015 case 2 :
1016 for( i=0; i<dataPointShape[0]; i++ )
1017 for( j=0; j<dataPointShape[1]; j++ )
1018 numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
1019 break;
1020 case 3 :
1021 for( i=0; i<dataPointShape[0]; i++ )
1022 for( j=0; j<dataPointShape[1]; j++ )
1023 for( k=0; k<dataPointShape[2]; k++ )
1024 numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
1025 break;
1026 case 4 :
1027 for( i=0; i<dataPointShape[0]; i++ )
1028 for( j=0; j<dataPointShape[1]; j++ )
1029 for( k=0; k<dataPointShape[2]; k++ )
1030 for( l=0; l<dataPointShape[3]; l++ )
1031 numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
1032 break;
1033 }
1034
1035 delete [] tmpData;
1036 //
1037 // return the loaded array
1038 return numArray;
1039 }
1040
1041
1042
1043 boost::python::numeric::array
1044 Data::integrate() const
1045 {
1046 int index;
1047 int rank = getDataPointRank();
1048 DataTypes::ShapeType shape = getDataPointShape();
1049 int dataPointSize = getDataPointSize();
1050
1051 //
1052 // calculate the integral values
1053 vector<double> integrals(dataPointSize);
1054 vector<double> integrals_local(dataPointSize);
1055 #ifdef PASO_MPI
1056 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);
1057 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1058 double *tmp = new double[dataPointSize];
1059 double *tmp_local = new double[dataPointSize];
1060 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1061 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1062 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1063 delete[] tmp;
1064 delete[] tmp_local;
1065 #else
1066 AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);
1067 #endif
1068
1069 //
1070 // create the numeric array to be returned
1071 // and load the array with the integral values
1072 boost::python::numeric::array bp_array(1.0);
1073 if (rank==0) {
1074 bp_array.resize(1);
1075 index = 0;
1076 bp_array[0] = integrals[index];
1077 }
1078 if (rank==1) {
1079 bp_array.resize(shape[0]);
1080 for (int i=0; i<shape[0]; i++) {
1081 index = i;
1082 bp_array[i] = integrals[index];
1083 }
1084 }
1085 if (rank==2) {
1086 bp_array.resize(shape[0],shape[1]);
1087 for (int i=0; i<shape[0]; i++) {
1088 for (int j=0; j<shape[1]; j++) {
1089 index = i + shape[0] * j;
1090 bp_array[make_tuple(i,j)] = integrals[index];
1091 }
1092 }
1093 }
1094 if (rank==3) {
1095 bp_array.resize(shape[0],shape[1],shape[2]);
1096 for (int i=0; i<shape[0]; i++) {
1097 for (int j=0; j<shape[1]; j++) {
1098 for (int k=0; k<shape[2]; k++) {
1099 index = i + shape[0] * ( j + shape[1] * k );
1100 bp_array[make_tuple(i,j,k)] = integrals[index];
1101 }
1102 }
1103 }
1104 }
1105 if (rank==4) {
1106 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
1107 for (int i=0; i<shape[0]; i++) {
1108 for (int j=0; j<shape[1]; j++) {
1109 for (int k=0; k<shape[2]; k++) {
1110 for (int l=0; l<shape[3]; l++) {
1111 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1112 bp_array[make_tuple(i,j,k,l)] = integrals[index];
1113 }
1114 }
1115 }
1116 }
1117 }
1118
1119 //
1120 // return the loaded array
1121 return bp_array;
1122 }
1123
1124 Data
1125 Data::sin() const
1126 {
1127 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1128 }
1129
1130 Data
1131 Data::cos() const
1132 {
1133 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1134 }
1135
1136 Data
1137 Data::tan() const
1138 {
1139 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1140 }
1141
1142 Data
1143 Data::asin() const
1144 {
1145 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1146 }
1147
1148 Data
1149 Data::acos() const
1150 {
1151 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1152 }
1153
1154
1155 Data
1156 Data::atan() const
1157 {
1158 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1159 }
1160
1161 Data
1162 Data::sinh() const
1163 {
1164 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1165
1166 }
1167
1168 Data
1169 Data::cosh() const
1170 {
1171 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1172 }
1173
1174 Data
1175 Data::tanh() const
1176 {
1177 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1178 }
1179
1180
1181 Data
1182 Data::erf() const
1183 {
1184 #ifdef _WIN32
1185 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1186 #else
1187 return C_TensorUnaryOperation(*this, ::erf);
1188 #endif
1189 }
1190
1191 Data
1192 Data::asinh() const
1193 {
1194 #ifdef _WIN32
1195 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1196 #else
1197 return C_TensorUnaryOperation(*this, ::asinh);
1198 #endif
1199 }
1200
1201 Data
1202 Data::acosh() const
1203 {
1204 #ifdef _WIN32
1205 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1206 #else
1207 return C_TensorUnaryOperation(*this, ::acosh);
1208 #endif
1209 }
1210
1211 Data
1212 Data::atanh() const
1213 {
1214 #ifdef _WIN32
1215 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1216 #else
1217 return C_TensorUnaryOperation(*this, ::atanh);
1218 #endif
1219 }
1220
1221 Data
1222 Data::log10() const
1223 {
1224 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1225 }
1226
1227 Data
1228 Data::log() const
1229 {
1230 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1231 }
1232
1233 Data
1234 Data::sign() const
1235 {
1236 return C_TensorUnaryOperation(*this, escript::fsign);
1237 }
1238
1239 Data
1240 Data::abs() const
1241 {
1242 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1243 }
1244
1245 Data
1246 Data::neg() const
1247 {
1248 return C_TensorUnaryOperation(*this, negate<double>());
1249 }
1250
1251 Data
1252 Data::pos() const
1253 {
1254 Data result;
1255 // perform a deep copy
1256 result.copy(*this);
1257 return result;
1258 }
1259
1260 Data
1261 Data::exp() const
1262 {
1263 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1264 }
1265
1266 Data
1267 Data::sqrt() const
1268 {
1269 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1270 }
1271
1272 double
1273 Data::Lsup() const
1274 {
1275 double localValue;
1276 //
1277 // set the initial absolute maximum value to zero
1278
1279 AbsMax abs_max_func;
1280 localValue = algorithm(abs_max_func,0);
1281 #ifdef PASO_MPI
1282 double globalValue;
1283 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1284 return globalValue;
1285 #else
1286 return localValue;
1287 #endif
1288 }
1289
1290 double
1291 Data::sup() const
1292 {
1293 double localValue;
1294 //
1295 // set the initial maximum value to min possible double
1296 FMax fmax_func;
1297 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1298 #ifdef PASO_MPI
1299 double globalValue;
1300 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1301 return globalValue;
1302 #else
1303 return localValue;
1304 #endif
1305 }
1306
1307 double
1308 Data::inf() const
1309 {
1310 double localValue;
1311 //
1312 // set the initial minimum value to max possible double
1313 FMin fmin_func;
1314 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1315 #ifdef PASO_MPI
1316 double globalValue;
1317 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1318 return globalValue;
1319 #else
1320 return localValue;
1321 #endif
1322 }
1323
1324 /* TODO */
1325 /* global reduction */
1326 Data
1327 Data::maxval() const
1328 {
1329 //
1330 // set the initial maximum value to min possible double
1331 FMax fmax_func;
1332 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1333 }
1334
1335 Data
1336 Data::minval() const
1337 {
1338 //
1339 // set the initial minimum value to max possible double
1340 FMin fmin_func;
1341 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1342 }
1343
1344 Data
1345 Data::swapaxes(const int axis0, const int axis1) const
1346 {
1347 int axis0_tmp,axis1_tmp;
1348 DataTypes::ShapeType s=getDataPointShape();
1349 DataTypes::ShapeType ev_shape;
1350 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1351 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1352 int rank=getDataPointRank();
1353 if (rank<2) {
1354 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1355 }
1356 if (axis0<0 || axis0>rank-1) {
1357 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1358 }
1359 if (axis1<0 || axis1>rank-1) {
1360 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1361 }
1362 if (axis0 == axis1) {
1363 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1364 }
1365 if (axis0 > axis1) {
1366 axis0_tmp=axis1;
1367 axis1_tmp=axis0;
1368 } else {
1369 axis0_tmp=axis0;
1370 axis1_tmp=axis1;
1371 }
1372 for (int i=0; i<rank; i++) {
1373 if (i == axis0_tmp) {
1374 ev_shape.push_back(s[axis1_tmp]);
1375 } else if (i == axis1_tmp) {
1376 ev_shape.push_back(s[axis0_tmp]);
1377 } else {
1378 ev_shape.push_back(s[i]);
1379 }
1380 }
1381 Data ev(0.,ev_shape,getFunctionSpace());
1382 ev.typeMatchRight(*this);
1383 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1384 return ev;
1385
1386 }
1387
1388 Data
1389 Data::symmetric() const
1390 {
1391 // check input
1392 DataTypes::ShapeType s=getDataPointShape();
1393 if (getDataPointRank()==2) {
1394 if(s[0] != s[1])
1395 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1396 }
1397 else if (getDataPointRank()==4) {
1398 if(!(s[0] == s[2] && s[1] == s[3]))
1399 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1400 }
1401 else {
1402 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1403 }
1404 Data ev(0.,getDataPointShape(),getFunctionSpace());
1405 ev.typeMatchRight(*this);
1406 m_data->symmetric(ev.m_data.get());
1407 return ev;
1408 }
1409
1410 Data
1411 Data::nonsymmetric() const
1412 {
1413 // check input
1414 DataTypes::ShapeType s=getDataPointShape();
1415 if (getDataPointRank()==2) {
1416 if(s[0] != s[1])
1417 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1418 DataTypes::ShapeType ev_shape;
1419 ev_shape.push_back(s[0]);
1420 ev_shape.push_back(s[1]);
1421 Data ev(0.,ev_shape,getFunctionSpace());
1422 ev.typeMatchRight(*this);
1423 m_data->nonsymmetric(ev.m_data.get());
1424 return ev;
1425 }
1426 else if (getDataPointRank()==4) {
1427 if(!(s[0] == s[2] && s[1] == s[3]))
1428 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1429 DataTypes::ShapeType ev_shape;
1430 ev_shape.push_back(s[0]);
1431 ev_shape.push_back(s[1]);
1432 ev_shape.push_back(s[2]);
1433 ev_shape.push_back(s[3]);
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 {
1440 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1441 }
1442 }
1443
1444 Data
1445 Data::trace(int axis_offset) const
1446 {
1447 DataTypes::ShapeType s=getDataPointShape();
1448 if (getDataPointRank()==2) {
1449 DataTypes::ShapeType ev_shape;
1450 Data ev(0.,ev_shape,getFunctionSpace());
1451 ev.typeMatchRight(*this);
1452 m_data->trace(ev.m_data.get(), axis_offset);
1453 return ev;
1454 }
1455 if (getDataPointRank()==3) {
1456 DataTypes::ShapeType ev_shape;
1457 if (axis_offset==0) {
1458 int s2=s[2];
1459 ev_shape.push_back(s2);
1460 }
1461 else if (axis_offset==1) {
1462 int s0=s[0];
1463 ev_shape.push_back(s0);
1464 }
1465 Data ev(0.,ev_shape,getFunctionSpace());
1466 ev.typeMatchRight(*this);
1467 m_data->trace(ev.m_data.get(), axis_offset);
1468 return ev;
1469 }
1470 if (getDataPointRank()==4) {
1471 DataTypes::ShapeType ev_shape;
1472 if (axis_offset==0) {
1473 ev_shape.push_back(s[2]);
1474 ev_shape.push_back(s[3]);
1475 }
1476 else if (axis_offset==1) {
1477 ev_shape.push_back(s[0]);
1478 ev_shape.push_back(s[3]);
1479 }
1480 else if (axis_offset==2) {
1481 ev_shape.push_back(s[0]);
1482 ev_shape.push_back(s[1]);
1483 }
1484 Data ev(0.,ev_shape,getFunctionSpace());
1485 ev.typeMatchRight(*this);
1486 m_data->trace(ev.m_data.get(), axis_offset);
1487 return ev;
1488 }
1489 else {
1490 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1491 }
1492 }
1493
1494 Data
1495 Data::transpose(int axis_offset) const
1496 {
1497 DataTypes::ShapeType s=getDataPointShape();
1498 DataTypes::ShapeType ev_shape;
1499 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1500 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1501 int rank=getDataPointRank();
1502 if (axis_offset<0 || axis_offset>rank) {
1503 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1504 }
1505 for (int i=0; i<rank; i++) {
1506 int index = (axis_offset+i)%rank;
1507 ev_shape.push_back(s[index]); // Append to new shape
1508 }
1509 Data ev(0.,ev_shape,getFunctionSpace());
1510 ev.typeMatchRight(*this);
1511 m_data->transpose(ev.m_data.get(), axis_offset);
1512 return ev;
1513 }
1514
1515 Data
1516 Data::eigenvalues() const
1517 {
1518 // check input
1519 DataTypes::ShapeType s=getDataPointShape();
1520 if (getDataPointRank()!=2)
1521 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1522 if(s[0] != s[1])
1523 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1524 // create return
1525 DataTypes::ShapeType ev_shape(1,s[0]);
1526 Data ev(0.,ev_shape,getFunctionSpace());
1527 ev.typeMatchRight(*this);
1528 m_data->eigenvalues(ev.m_data.get());
1529 return ev;
1530 }
1531
1532 const boost::python::tuple
1533 Data::eigenvalues_and_eigenvectors(const double tol) const
1534 {
1535 DataTypes::ShapeType s=getDataPointShape();
1536 if (getDataPointRank()!=2)
1537 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1538 if(s[0] != s[1])
1539 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1540 // create return
1541 DataTypes::ShapeType ev_shape(1,s[0]);
1542 Data ev(0.,ev_shape,getFunctionSpace());
1543 ev.typeMatchRight(*this);
1544 DataTypes::ShapeType V_shape(2,s[0]);
1545 Data V(0.,V_shape,getFunctionSpace());
1546 V.typeMatchRight(*this);
1547 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1548 return make_tuple(boost::python::object(ev),boost::python::object(V));
1549 }
1550
1551 const boost::python::tuple
1552 Data::minGlobalDataPoint() const
1553 {
1554 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1555 // abort (for unknown reasons) if there are openmp directives with it in the
1556 // surrounding function
1557
1558 int DataPointNo;
1559 int ProcNo;
1560 calc_minGlobalDataPoint(ProcNo,DataPointNo);
1561 return make_tuple(ProcNo,DataPointNo);
1562 }
1563
1564 void
1565 Data::calc_minGlobalDataPoint(int& ProcNo,
1566 int& DataPointNo) const
1567 {
1568 int i,j;
1569 int lowi=0,lowj=0;
1570 double min=numeric_limits<double>::max();
1571
1572 Data temp=minval();
1573
1574 int numSamples=temp.getNumSamples();
1575 int numDPPSample=temp.getNumDataPointsPerSample();
1576
1577 double next,local_min;
1578 int local_lowi,local_lowj;
1579
1580 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1581 {
1582 local_min=min;
1583 #pragma omp for private(i,j) schedule(static)
1584 for (i=0; i<numSamples; i++) {
1585 for (j=0; j<numDPPSample; j++) {
1586 next=temp.getDataAtOffset(temp.getDataOffset(i,j));
1587 if (next<local_min) {
1588 local_min=next;
1589 local_lowi=i;
1590 local_lowj=j;
1591 }
1592 }
1593 }
1594 #pragma omp critical
1595 if (local_min<min) {
1596 min=local_min;
1597 lowi=local_lowi;
1598 lowj=local_lowj;
1599 }
1600 }
1601
1602 #ifdef PASO_MPI
1603 // determine the processor on which the minimum occurs
1604 next = temp.getDataPoint(lowi,lowj);
1605 int lowProc = 0;
1606 double *globalMins = new double[get_MPISize()+1];
1607 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1608
1609 if( get_MPIRank()==0 ){
1610 next = globalMins[lowProc];
1611 for( i=1; i<get_MPISize(); i++ )
1612 if( next>globalMins[i] ){
1613 lowProc = i;
1614 next = globalMins[i];
1615 }
1616 }
1617 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1618
1619 delete [] globalMins;
1620 ProcNo = lowProc;
1621 #else
1622 ProcNo = 0;
1623 #endif
1624 DataPointNo = lowj + lowi * numDPPSample;
1625 }
1626
1627 void
1628 Data::saveDX(std::string fileName) const
1629 {
1630 if (isEmpty())
1631 {
1632 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1633 }
1634 boost::python::dict args;
1635 args["data"]=boost::python::object(this);
1636 getDomain()->saveDX(fileName,args);
1637 return;
1638 }
1639
1640 void
1641 Data::saveVTK(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()->saveVTK(fileName,args);
1650 return;
1651 }
1652
1653 Data&
1654 Data::operator+=(const Data& right)
1655 {
1656 if (isProtected()) {
1657 throw DataException("Error - attempt to update protected Data object.");
1658 }
1659 binaryOp(right,plus<double>());
1660 return (*this);
1661 }
1662
1663 Data&
1664 Data::operator+=(const boost::python::object& right)
1665 {
1666 Data tmp(right,getFunctionSpace(),false);
1667 binaryOp(tmp,plus<double>());
1668 return (*this);
1669 }
1670 Data&
1671 Data::operator=(const Data& other)
1672 {
1673 copy(other);
1674 return (*this);
1675 }
1676
1677 Data&
1678 Data::operator-=(const Data& right)
1679 {
1680 if (isProtected()) {
1681 throw DataException("Error - attempt to update protected Data object.");
1682 }
1683 binaryOp(right,minus<double>());
1684 return (*this);
1685 }
1686
1687 Data&
1688 Data::operator-=(const boost::python::object& right)
1689 {
1690 Data tmp(right,getFunctionSpace(),false);
1691 binaryOp(tmp,minus<double>());
1692 return (*this);
1693 }
1694
1695 Data&
1696 Data::operator*=(const Data& right)
1697 {
1698 if (isProtected()) {
1699 throw DataException("Error - attempt to update protected Data object.");
1700 }
1701 binaryOp(right,multiplies<double>());
1702 return (*this);
1703 }
1704
1705 Data&
1706 Data::operator*=(const boost::python::object& right)
1707 {
1708 Data tmp(right,getFunctionSpace(),false);
1709 binaryOp(tmp,multiplies<double>());
1710 return (*this);
1711 }
1712
1713 Data&
1714 Data::operator/=(const Data& right)
1715 {
1716 if (isProtected()) {
1717 throw DataException("Error - attempt to update protected Data object.");
1718 }
1719 binaryOp(right,divides<double>());
1720 return (*this);
1721 }
1722
1723 Data&
1724 Data::operator/=(const boost::python::object& right)
1725 {
1726 Data tmp(right,getFunctionSpace(),false);
1727 binaryOp(tmp,divides<double>());
1728 return (*this);
1729 }
1730
1731 Data
1732 Data::rpowO(const boost::python::object& left) const
1733 {
1734 Data left_d(left,*this);
1735 return left_d.powD(*this);
1736 }
1737
1738 Data
1739 Data::powO(const boost::python::object& right) const
1740 {
1741 Data tmp(right,getFunctionSpace(),false);
1742 return powD(tmp);
1743 }
1744
1745 Data
1746 Data::powD(const Data& right) const
1747 {
1748 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
1749 }
1750
1751 //
1752 // NOTE: It is essential to specify the namespace this operator belongs to
1753 Data
1754 escript::operator+(const Data& left, const Data& right)
1755 {
1756 return C_TensorBinaryOperation(left, right, plus<double>());
1757 }
1758
1759 //
1760 // NOTE: It is essential to specify the namespace this operator belongs to
1761 Data
1762 escript::operator-(const Data& left, const Data& right)
1763 {
1764 return C_TensorBinaryOperation(left, right, minus<double>());
1765 }
1766
1767 //
1768 // NOTE: It is essential to specify the namespace this operator belongs to
1769 Data
1770 escript::operator*(const Data& left, const Data& right)
1771 {
1772 return C_TensorBinaryOperation(left, right, multiplies<double>());
1773 }
1774
1775 //
1776 // NOTE: It is essential to specify the namespace this operator belongs to
1777 Data
1778 escript::operator/(const Data& left, const Data& right)
1779 {
1780 return C_TensorBinaryOperation(left, right, divides<double>());
1781 }
1782
1783 //
1784 // NOTE: It is essential to specify the namespace this operator belongs to
1785 Data
1786 escript::operator+(const Data& left, const boost::python::object& right)
1787 {
1788 return left+Data(right,left.getFunctionSpace(),false);
1789 }
1790
1791 //
1792 // NOTE: It is essential to specify the namespace this operator belongs to
1793 Data
1794 escript::operator-(const Data& left, const boost::python::object& right)
1795 {
1796 return left-Data(right,left.getFunctionSpace(),false);
1797 }
1798
1799 //
1800 // NOTE: It is essential to specify the namespace this operator belongs to
1801 Data
1802 escript::operator*(const Data& left, const boost::python::object& right)
1803 {
1804 return left*Data(right,left.getFunctionSpace(),false);
1805 }
1806
1807 //
1808 // NOTE: It is essential to specify the namespace this operator belongs to
1809 Data
1810 escript::operator/(const Data& left, const boost::python::object& right)
1811 {
1812 return left/Data(right,left.getFunctionSpace(),false);
1813 }
1814
1815 //
1816 // NOTE: It is essential to specify the namespace this operator belongs to
1817 Data
1818 escript::operator+(const boost::python::object& left, const Data& right)
1819 {
1820 return Data(left,right.getFunctionSpace(),false)+right;
1821 }
1822
1823 //
1824 // NOTE: It is essential to specify the namespace this operator belongs to
1825 Data
1826 escript::operator-(const boost::python::object& left, const Data& right)
1827 {
1828 return Data(left,right.getFunctionSpace(),false)-right;
1829 }
1830
1831 //
1832 // NOTE: It is essential to specify the namespace this operator belongs to
1833 Data
1834 escript::operator*(const boost::python::object& left, const Data& right)
1835 {
1836 return Data(left,right.getFunctionSpace(),false)*right;
1837 }
1838
1839 //
1840 // NOTE: It is essential to specify the namespace this operator belongs to
1841 Data
1842 escript::operator/(const boost::python::object& left, const Data& right)
1843 {
1844 return Data(left,right.getFunctionSpace(),false)/right;
1845 }
1846
1847 //
1848 //bool escript::operator==(const Data& left, const Data& right)
1849 //{
1850 // /*
1851 // NB: this operator does very little at this point, and isn't to
1852 // be relied on. Requires further implementation.
1853 // */
1854 //
1855 // bool ret;
1856 //
1857 // if (left.isEmpty()) {
1858 // if(!right.isEmpty()) {
1859 // ret = false;
1860 // } else {
1861 // ret = true;
1862 // }
1863 // }
1864 //
1865 // if (left.isConstant()) {
1866 // if(!right.isConstant()) {
1867 // ret = false;
1868 // } else {
1869 // ret = true;
1870 // }
1871 // }
1872 //
1873 // if (left.isTagged()) {
1874 // if(!right.isTagged()) {
1875 // ret = false;
1876 // } else {
1877 // ret = true;
1878 // }
1879 // }
1880 //
1881 // if (left.isExpanded()) {
1882 // if(!right.isExpanded()) {
1883 // ret = false;
1884 // } else {
1885 // ret = true;
1886 // }
1887 // }
1888 //
1889 // return ret;
1890 //}
1891
1892 /* TODO */
1893 /* global reduction */
1894 Data
1895 Data::getItem(const boost::python::object& key) const
1896 {
1897 // const DataArrayView& view=getPointDataView();
1898
1899 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1900
1901 if (slice_region.size()!=getDataPointRank()) {
1902 throw DataException("Error - slice size does not match Data rank.");
1903 }
1904
1905 return getSlice(slice_region);
1906 }
1907
1908 /* TODO */
1909 /* global reduction */
1910 Data
1911 Data::getSlice(const DataTypes::RegionType& region) const
1912 {
1913 return Data(*this,region);
1914 }
1915
1916 /* TODO */
1917 /* global reduction */
1918 void
1919 Data::setItemO(const boost::python::object& key,
1920 const boost::python::object& value)
1921 {
1922 Data tempData(value,getFunctionSpace());
1923 setItemD(key,tempData);
1924 }
1925
1926 void
1927 Data::setItemD(const boost::python::object& key,
1928 const Data& value)
1929 {
1930 // const DataArrayView& view=getPointDataView();
1931
1932 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
1933 if (slice_region.size()!=getDataPointRank()) {
1934 throw DataException("Error - slice size does not match Data rank.");
1935 }
1936 if (getFunctionSpace()!=value.getFunctionSpace()) {
1937 setSlice(Data(value,getFunctionSpace()),slice_region);
1938 } else {
1939 setSlice(value,slice_region);
1940 }
1941 }
1942
1943 void
1944 Data::setSlice(const Data& value,
1945 const DataTypes::RegionType& region)
1946 {
1947 if (isProtected()) {
1948 throw DataException("Error - attempt to update protected Data object.");
1949 }
1950 Data tempValue(value);
1951 typeMatchLeft(tempValue);
1952 typeMatchRight(tempValue);
1953 m_data->setSlice(tempValue.m_data.get(),region);
1954 }
1955
1956 void
1957 Data::typeMatchLeft(Data& right) const
1958 {
1959 if (isExpanded()){
1960 right.expand();
1961 } else if (isTagged()) {
1962 if (right.isConstant()) {
1963 right.tag();
1964 }
1965 }
1966 }
1967
1968 void
1969 Data::typeMatchRight(const Data& right)
1970 {
1971 if (isTagged()) {
1972 if (right.isExpanded()) {
1973 expand();
1974 }
1975 } else if (isConstant()) {
1976 if (right.isExpanded()) {
1977 expand();
1978 } else if (right.isTagged()) {
1979 tag();
1980 }
1981 }
1982 }
1983
1984 void
1985 Data::setTaggedValueByName(std::string name,
1986 const boost::python::object& value)
1987 {
1988 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
1989 int tagKey=getFunctionSpace().getDomain()->getTag(name);
1990 setTaggedValue(tagKey,value);
1991 }
1992 }
1993 void
1994 Data::setTaggedValue(int tagKey,
1995 const boost::python::object& value)
1996 {
1997 if (isProtected()) {
1998 throw DataException("Error - attempt to update protected Data object.");
1999 }
2000 //
2001 // Ensure underlying data object is of type DataTagged
2002 if (isConstant()) tag();
2003
2004 numeric::array asNumArray(value);
2005
2006
2007 // extract the shape of the numarray
2008 DataTypes::ShapeType tempShape;
2009 for (int i=0; i < asNumArray.getrank(); i++) {
2010 tempShape.push_back(extract<int>(asNumArray.getshape()[i]));
2011 }
2012
2013 // get the space for the data vector
2014 // int len = DataTypes::noValues(tempShape);
2015 // DataVector temp_data(len, 0.0, len);
2016 // DataArrayView temp_dataView(temp_data, tempShape);
2017 // temp_dataView.copy(asNumArray);
2018
2019 DataVector temp_data2;
2020 temp_data2.copyFromNumArray(asNumArray);
2021
2022 //
2023 // Call DataAbstract::setTaggedValue
2024 //m_data->setTaggedValue(tagKey,temp_dataView);
2025
2026 m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2027 }
2028
2029
2030 void
2031 Data::setTaggedValueFromCPP(int tagKey,
2032 const DataTypes::ShapeType& pointshape,
2033 const DataTypes::ValueType& value,
2034 int dataOffset)
2035 {
2036 if (isProtected()) {
2037 throw DataException("Error - attempt to update protected Data object.");
2038 }
2039 //
2040 // Ensure underlying data object is of type DataTagged
2041 if (isConstant()) tag();
2042
2043 //
2044 // Call DataAbstract::setTaggedValue
2045 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2046 }
2047
2048 int
2049 Data::getTagNumber(int dpno)
2050 {
2051 if (isEmpty())
2052 {
2053 throw DataException("Error - operation not permitted on instances of DataEmpty.");
2054 }
2055 return getFunctionSpace().getTagFromDataPointNo(dpno);
2056 }
2057
2058
2059 ostream& escript::operator<<(ostream& o, const Data& data)
2060 {
2061 o << data.toString();
2062 return o;
2063 }
2064
2065 Data
2066 escript::C_GeneralTensorProduct(Data& arg_0,
2067 Data& arg_1,
2068 int axis_offset,
2069 int transpose)
2070 {
2071 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2072 // SM is the product of the last axis_offset entries in arg_0.getShape().
2073
2074 // Interpolate if necessary and find an appropriate function space
2075 Data arg_0_Z, arg_1_Z;
2076 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2077 if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2078 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2079 arg_1_Z = Data(arg_1);
2080 }
2081 else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2082 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2083 arg_0_Z =Data(arg_0);
2084 }
2085 else {
2086 throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2087 }
2088 } else {
2089 arg_0_Z = Data(arg_0);
2090 arg_1_Z = Data(arg_1);
2091 }
2092 // Get rank and shape of inputs
2093 int rank0 = arg_0_Z.getDataPointRank();
2094 int rank1 = arg_1_Z.getDataPointRank();
2095 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2096 const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2097
2098 // Prepare for the loops of the product and verify compatibility of shapes
2099 int start0=0, start1=0;
2100 if (transpose == 0) {}
2101 else if (transpose == 1) { start0 = axis_offset; }
2102 else if (transpose == 2) { start1 = rank1-axis_offset; }
2103 else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2104
2105
2106 // Adjust the shapes for transpose
2107 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2108 DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2109 for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2110 for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2111
2112 #if 0
2113 // For debugging: show shape after transpose
2114 char tmp[100];
2115 std::string shapeStr;
2116 shapeStr = "(";
2117 for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2118 shapeStr += ")";
2119 cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2120 shapeStr = "(";
2121 for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2122 shapeStr += ")";
2123 cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2124 #endif
2125
2126 // Prepare for the loops of the product
2127 int SL=1, SM=1, SR=1;
2128 for (int i=0; i<rank0-axis_offset; i++) {
2129 SL *= tmpShape0[i];
2130 }
2131 for (int i=rank0-axis_offset; i<rank0; i++) {
2132 if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2133 throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2134 }
2135 SM *= tmpShape0[i];
2136 }
2137 for (int i=axis_offset; i<rank1; i++) {
2138 SR *= tmpShape1[i];
2139 }
2140
2141 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2142 DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2143 { // block to limit the scope of out_index
2144 int out_index=0;
2145 for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2146 for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2147 }
2148
2149 // Declare output Data object
2150 Data res;
2151
2152 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2153 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2154 double *ptr_0 = &(arg_0_Z.getDataAtOffset(0));
2155 double *ptr_1 = &(arg_1_Z.getDataAtOffset(0));
2156 double *ptr_2 = &(res.getDataAtOffset(0));
2157 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2158 }
2159 else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2160
2161 // Prepare the DataConstant input
2162 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2163 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2164
2165 // Borrow DataTagged input from Data object
2166 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2167 if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2168
2169 // Prepare a DataTagged output 2
2170 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataTagged output
2171 res.tag();
2172 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2173 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2174
2175 // Prepare offset into DataConstant
2176 int offset_0 = tmp_0->getPointOffset(0,0);
2177 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2178 // Get the views
2179 // DataArrayView view_1 = tmp_1->getDefaultValue();
2180 // DataArrayView view_2 = tmp_2->getDefaultValue();
2181 // // Get the pointers to the actual data
2182 // double *ptr_1 = &((view_1.getData())[0]);
2183 // double *ptr_2 = &((view_2.getData())[0]);
2184
2185 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2186 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2187
2188
2189 // Compute an MVP for the default
2190 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2191 // Compute an MVP for each tag
2192 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2193 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2194 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2195 tmp_2->addTag(i->first);
2196 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2197 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2198 // double *ptr_1 = &view_1.getData(0);
2199 // double *ptr_2 = &view_2.getData(0);
2200
2201 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2202 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2203
2204 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2205 }
2206
2207 }
2208 else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2209
2210 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2211 DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2212 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2213 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2214 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2215 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2216 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2217 int sampleNo_1,dataPointNo_1;
2218 int numSamples_1 = arg_1_Z.getNumSamples();
2219 int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2220 int offset_0 = tmp_0->getPointOffset(0,0);
2221 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2222 for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2223 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2224 int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2225 int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2226 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2227 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2228 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2229 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2230 }
2231 }
2232
2233 }
2234 else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2235
2236 // Borrow DataTagged input from Data object
2237 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2238 if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2239
2240 // Prepare the DataConstant input
2241 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2242 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2243
2244 // Prepare a DataTagged output 2
2245 res = Data(0.0, shape2, arg_0_Z.getFunctionSpace()); // DataTagged output
2246 res.tag();
2247 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2248 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2249
2250 // Prepare offset into DataConstant
2251 int offset_1 = tmp_1->getPointOffset(0,0);
2252 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2253 // Get the views
2254 // DataArrayView view_0 = tmp_0->getDefaultValue();
2255 // DataArrayView view_2 = tmp_2->getDefaultValue();
2256 // // Get the pointers to the actual data
2257 // double *ptr_0 = &((view_0.getData())[0]);
2258 // double *ptr_2 = &((view_2.getData())[0]);
2259
2260 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2261 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2262
2263 // Compute an MVP for the default
2264 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2265 // Compute an MVP for each tag
2266 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2267 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2268 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2269 // tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2270 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2271 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2272 // double *ptr_0 = &view_0.getData(0);
2273 // double *ptr_2 = &view_2.getData(0);
2274
2275 tmp_2->addTag(i->first);
2276 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2277 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2278 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2279 }
2280
2281 }
2282 else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2283
2284 // Borrow DataTagged input from Data object
2285 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2286 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2287
2288 // Borrow DataTagged input from Data object
2289 DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2290 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2291
2292 // Prepare a DataTagged output 2
2293 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2294 res.tag(); // DataTagged output
2295 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2296 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2297
2298 // // Get the views
2299 // DataArrayView view_0 = tmp_0->getDefaultValue();
2300 // DataArrayView view_1 = tmp_1->getDefaultValue();
2301 // DataArrayView view_2 = tmp_2->getDefaultValue();
2302 // // Get the pointers to the actual data
2303 // double *ptr_0 = &((view_0.getData())[0]);
2304 // double *ptr_1 = &((view_1.getData())[0]);
2305 // double *ptr_2 = &((view_2.getData())[0]);
2306
2307 double *ptr_0 = &(tmp_0->getDefaultValue(0));
2308 double *ptr_1 = &(tmp_1->getDefaultValue(0));
2309 double *ptr_2 = &(tmp_2->getDefaultValue(0));
2310
2311
2312 // Compute an MVP for the default
2313 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2314 // Merge the tags
2315 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2316 const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2317 const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2318 for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2319 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2320 }
2321 for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2322 tmp_2->addTag(i->first);
2323 }
2324 // Compute an MVP for each tag
2325 const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2326 for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2327 // DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2328 // DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2329 // DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2330 // double *ptr_0 = &view_0.getData(0);
2331 // double *ptr_1 = &view_1.getData(0);
2332 // double *ptr_2 = &view_2.getData(0);
2333
2334 double *ptr_0 = &(tmp_0->getDataByTag(i->first,0));
2335 double *ptr_1 = &(tmp_1->getDataByTag(i->first,0));
2336 double *ptr_2 = &(tmp_2->getDataByTag(i->first,0));
2337
2338 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2339 }
2340
2341 }
2342 else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2343
2344 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2345 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2346 DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2347 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2348 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2349 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2350 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2351 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2352 int sampleNo_0,dataPointNo_0;
2353 int numSamples_0 = arg_0_Z.getNumSamples();
2354 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2355 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2356 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2357 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2358 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2359 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2360 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2361 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2362 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2363 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2364 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2365 }
2366 }
2367
2368 }
2369 else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2370
2371 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2372 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2373 DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2374 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2375 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2376 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2377 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2378 int sampleNo_0,dataPointNo_0;
2379 int numSamples_0 = arg_0_Z.getNumSamples();
2380 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2381 int offset_1 = tmp_1->getPointOffset(0,0);
2382 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2383 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2384 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2385 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2386 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2387 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2388 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2389 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2390 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2391 }
2392 }
2393
2394
2395 }
2396 else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2397
2398 // After finding a common function space above the two inputs have the same numSamples and num DPPS
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 DataTagged* tmp_1=dynamic_cast<DataTagged*>(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 DataTagged."); }
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 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2410 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2411 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2412 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2413 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2414 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2415 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2416 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
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 else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2424
2425 // After finding a common function space above the two inputs have the same numSamples and num DPPS
2426 res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2427 DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2428 DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2429 DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2430 if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2431 if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2432 if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2433 int sampleNo_0,dataPointNo_0;
2434 int numSamples_0 = arg_0_Z.getNumSamples();
2435 int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2436 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2437 for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2438 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2439 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2440 int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2441 int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2442 double *ptr_0 = &(arg_0_Z.getDataAtOffset(offset_0));
2443 double *ptr_1 = &(arg_1_Z.getDataAtOffset(offset_1));
2444 double *ptr_2 = &(res.getDataAtOffset(offset_2));
2445 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2446 }
2447 }
2448
2449 }
2450 else {
2451 throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2452 }
2453
2454 return res;
2455 }
2456
2457 DataAbstract*
2458 Data::borrowData() const
2459 {
2460 return m_data.get();
2461 }
2462
2463
2464 std::string
2465 Data::toString() const
2466 {
2467 if (!m_data->isEmpty() &&
2468 getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))
2469 {
2470 stringstream temp;
2471 temp << "Summary: inf="<< inf() << " sup=" << sup() << " data points=" << getNumDataPoints();
2472 return temp.str();
2473 }
2474 return m_data->toString();
2475 }
2476
2477
2478
2479 DataTypes::ValueType::const_reference
2480 Data::getDataAtOffset(DataTypes::ValueType::size_type i) const
2481 {
2482 return m_data->getDataAtOffset(i);
2483 }
2484
2485
2486 DataTypes::ValueType::reference
2487 Data::getDataAtOffset(DataTypes::ValueType::size_type i)
2488 {
2489 return m_data->getDataAtOffset(i);
2490 }
2491
2492 DataTypes::ValueType::const_reference
2493 Data::getDataPoint(int sampleNo, int dataPointNo) const
2494 {
2495 return m_data->getDataAtOffset(m_data->getPointOffset(sampleNo, dataPointNo));
2496 }
2497
2498
2499 DataTypes::ValueType::reference
2500 Data::getDataPoint(int sampleNo, int dataPointNo)
2501 {
2502 return m_data->getDataAtOffset(m_data->getPointOffset(sampleNo, dataPointNo));
2503 }
2504
2505
2506 /* Member functions specific to the MPI implementation */
2507
2508 void
2509 Data::print()
2510 {
2511 int i,j;
2512
2513 printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2514 for( i=0; i<getNumSamples(); i++ )
2515 {
2516 printf( "[%6d]", i );
2517 for( j=0; j<getNumDataPointsPerSample(); j++ )
2518 printf( "\t%10.7g", (getSampleData(i))[j] );
2519 printf( "\n" );
2520 }
2521 }
2522 void
2523 Data::dump(const std::string fileName) const
2524 {
2525 try
2526 {
2527 return m_data->dump(fileName);
2528 }
2529 catch (exception& e)
2530 {
2531 cout << e.what() << endl;
2532 }
2533 }
2534
2535 int
2536 Data::get_MPISize() const
2537 {
2538 int size;
2539 #ifdef PASO_MPI
2540 int error;
2541 error = MPI_Comm_size( get_MPIComm(), &size );
2542 #else
2543 size = 1;
2544 #endif
2545 return size;
2546 }
2547
2548 int
2549 Data::get_MPIRank() const
2550 {
2551 int rank;
2552 #ifdef PASO_MPI
2553 int error;
2554 error = MPI_Comm_rank( get_MPIComm(), &rank );
2555 #else
2556 rank = 0;
2557 #endif
2558 return rank;
2559 }
2560
2561 MPI_Comm
2562 Data::get_MPIComm() const
2563 {
2564 #ifdef PASO_MPI
2565 return MPI_COMM_WORLD;
2566 #else
2567 return -1;
2568 #endif
2569 }
2570
2571

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26