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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 1 month ago) by ksteube
File size: 73725 byte(s)
Copyright updated in all files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26