/[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 1712 - (show annotations)
Wed Aug 20 05:04:28 2008 UTC (10 years, 9 months ago) by jfenwick
File size: 77255 byte(s)
Branch commit.

Finished first pass of Data.h - There is still a constructor which takes 
a DataArrayView as a parameter. Apart from that, there are no direct 
references to DataArrayView.

DataTagged has a new constructor for copying just the tags from an 
existing object.
DataTypes:: now has a scalarShape constant (to avoid creating one 
everytime you create a scalar).



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26