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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1319 - (show annotations)
Thu Sep 27 00:27:51 2007 UTC (12 years, 6 months ago) by matt
File size: 76648 byte(s)
DataArray is no longer needed. However, the unit tests still require it.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26