/[escript]/branches/arrayview_from_1695_trunk/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1714 - (show annotations)
Thu Aug 21 00:01:55 2008 UTC (12 years, 1 month ago) by jfenwick
File size: 77941 byte(s)
Branch commit

Moved getSliceRegion() and getSliceRange() into DataTypes
Data.cpp - modified not to rely on operator() from DataArrayView
         - Used more const& to avoid copies


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26