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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26