/[escript]/branches/more_shared_ptrs_from_1812/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/more_shared_ptrs_from_1812/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1828 - (show annotations)
Thu Oct 2 04:52:11 2008 UTC (11 years, 4 months ago) by jfenwick
File size: 73046 byte(s)
Branch commit.
Added getPtr to DataAbstract.
Passes all unit tests.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26