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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1799 - (show annotations)
Wed Sep 17 06:33:18 2008 UTC (10 years, 9 months ago) by jfenwick
File size: 72960 byte(s)
Added Data::copySelf() [Note: this is exposed as copy() in python].
This method returns a pointer to a deep copy of the target.
There are c++ tests but no python tests for this yet.

All DataAbstracts now have a deepCopy() which simplifies the 
implementation of the compy methods.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26