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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26