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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1796 - (show annotations)
Wed Sep 17 01:45:46 2008 UTC (10 years, 7 months ago) by jfenwick
File size: 73653 byte(s)
Merged noarrayview branch onto trunk.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26