/[escript]/trunk-mpi-branch/escript/src/Data.cpp
ViewVC logotype

Contents of /trunk-mpi-branch/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1311 - (show annotations)
Mon Sep 24 06:13:58 2007 UTC (14 years, 10 months ago) by matt
File size: 76030 byte(s)
Replaced DataArray in several constructors.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26