/[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 1868 - (show annotations)
Thu Oct 9 06:30:49 2008 UTC (11 years, 10 months ago) by jfenwick
File size: 75555 byte(s)
Branch commit.
Bulk resolve for + - * /

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26