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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1385 - (show annotations)
Fri Jan 11 07:33:30 2008 UTC (13 years, 1 month ago) by trankine
File size: 75706 byte(s)
This is a commit to the trunk of the current windows version undergoing debugging.
This trunk will shortly be moved to the branches area
of the repository, and the saved trunk version currently over in branches restored to the trunk.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26