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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 8 months ago) by ksteube
File size: 75618 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26