/[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 1864 - (show annotations)
Thu Oct 9 03:09:30 2008 UTC (11 years, 5 months ago) by jfenwick
File size: 74130 byte(s)
Branch commit
It compiles but doesn't do much.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26