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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1865 - (show annotations)
Thu Oct 9 03:53:57 2008 UTC (12 years, 8 months ago) by jfenwick
File size: 74374 byte(s)
Branch commit
Added some missing files.

In python can ask a data object if it-  isReady(), isConstant(), 
isLazy().



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26