/[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 1196 - (show annotations)
Fri Jun 15 03:45:48 2007 UTC (11 years, 10 months ago) by ksteube
File size: 76223 byte(s)
Use of PAPI on solver is now enabled with papi_instrument_solver=1 in scons/ess_options.py.
Can instrument other blocks of code with blockpapi.c.
Added interval timers to grad, integrate and Assemble_PDE.

1 // $Id$
2
3 /*
4 ************************************************************
5 * Copyright 2006 by ACcESS MNRF *
6 * *
7 * http://www.access.edu.au *
8 * Primary Business: Queensland, Australia *
9 * Licensed under the Open Software License version 3.0 *
10 * http://www.opensource.org/licenses/osl-3.0.php *
11 * *
12 ************************************************************
13 */
14 #include "Data.h"
15
16 #include "paso/Common.h"
17 #include "DataExpanded.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20 #include "DataEmpty.h"
21 #include "DataArray.h"
22 #include "DataArrayView.h"
23 #include "FunctionSpaceFactory.h"
24 #include "AbstractContinuousDomain.h"
25 #include "UnaryFuncs.h"
26 extern "C" {
27 #include "escript/blocktimer.h"
28 }
29
30 #include <fstream>
31 #include <algorithm>
32 #include <vector>
33 #include <functional>
34
35 #include <boost/python/dict.hpp>
36 #include <boost/python/extract.hpp>
37 #include <boost/python/long.hpp>
38
39 using namespace std;
40 using namespace boost::python;
41 using namespace boost;
42 using namespace escript;
43
44 Data::Data()
45 {
46 //
47 // Default data is type DataEmpty
48 DataAbstract* temp=new DataEmpty();
49 shared_ptr<DataAbstract> temp_data(temp);
50 m_data=temp_data;
51 m_protected=false;
52 }
53
54 Data::Data(double value,
55 const tuple& shape,
56 const FunctionSpace& what,
57 bool expanded)
58 {
59 DataArrayView::ShapeType dataPointShape;
60 for (int i = 0; i < shape.attr("__len__")(); ++i) {
61 dataPointShape.push_back(extract<const int>(shape[i]));
62 }
63 DataArray temp(dataPointShape,value);
64 initialise(temp.getView(),what,expanded);
65 m_protected=false;
66 }
67
68 Data::Data(double value,
69 const DataArrayView::ShapeType& dataPointShape,
70 const FunctionSpace& what,
71 bool expanded)
72 {
73 DataArray temp(dataPointShape,value);
74 pair<int,int> dataShape=what.getDataShape();
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, *tmp_local;
812 tmp = MEMALLOC(dataPointSize,double);
813 tmp_local = MEMALLOC(dataPointSize,double);
814 for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
815 MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
816 for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
817 MEMFREE(tmp);
818 MEMFREE(tmp_local);
819 #else
820 AbstractContinuousDomain::asAbstractContinuousDomain(getDomain()).setToIntegrals(integrals,*this);
821 #endif
822
823 //
824 // create the numeric array to be returned
825 // and load the array with the integral values
826 boost::python::numeric::array bp_array(1.0);
827 if (rank==0) {
828 bp_array.resize(1);
829 index = 0;
830 bp_array[0] = integrals[index];
831 }
832 if (rank==1) {
833 bp_array.resize(shape[0]);
834 for (int i=0; i<shape[0]; i++) {
835 index = i;
836 bp_array[i] = integrals[index];
837 }
838 }
839 if (rank==2) {
840 bp_array.resize(shape[0],shape[1]);
841 for (int i=0; i<shape[0]; i++) {
842 for (int j=0; j<shape[1]; j++) {
843 index = i + shape[0] * j;
844 bp_array[make_tuple(i,j)] = integrals[index];
845 }
846 }
847 }
848 if (rank==3) {
849 bp_array.resize(shape[0],shape[1],shape[2]);
850 for (int i=0; i<shape[0]; i++) {
851 for (int j=0; j<shape[1]; j++) {
852 for (int k=0; k<shape[2]; k++) {
853 index = i + shape[0] * ( j + shape[1] * k );
854 bp_array[make_tuple(i,j,k)] = integrals[index];
855 }
856 }
857 }
858 }
859 if (rank==4) {
860 bp_array.resize(shape[0],shape[1],shape[2],shape[3]);
861 for (int i=0; i<shape[0]; i++) {
862 for (int j=0; j<shape[1]; j++) {
863 for (int k=0; k<shape[2]; k++) {
864 for (int l=0; l<shape[3]; l++) {
865 index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
866 bp_array[make_tuple(i,j,k,l)] = integrals[index];
867 }
868 }
869 }
870 }
871 }
872
873 //
874 // return the loaded array
875 return bp_array;
876 }
877
878 Data
879 Data::sin() const
880 {
881 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
882 }
883
884 Data
885 Data::cos() const
886 {
887 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
888 }
889
890 Data
891 Data::tan() const
892 {
893 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
894 }
895
896 Data
897 Data::asin() const
898 {
899 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
900 }
901
902 Data
903 Data::acos() const
904 {
905 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
906 }
907
908
909 Data
910 Data::atan() const
911 {
912 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
913 }
914
915 Data
916 Data::sinh() const
917 {
918 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
919 }
920
921 Data
922 Data::cosh() const
923 {
924 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
925 }
926
927 Data
928 Data::tanh() const
929 {
930 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
931 }
932
933
934 Data
935 Data::erf() const
936 {
937 #ifdef _WIN32
938 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
939 #else
940 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
941 #endif
942 }
943
944 Data
945 Data::asinh() const
946 {
947 #ifdef _WIN32
948 return escript::unaryOp(*this,escript::asinh_substitute);
949 #else
950 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
951 #endif
952 }
953
954 Data
955 Data::acosh() const
956 {
957 #ifdef _WIN32
958 return escript::unaryOp(*this,escript::acosh_substitute);
959 #else
960 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
961 #endif
962 }
963
964 Data
965 Data::atanh() const
966 {
967 #ifdef _WIN32
968 return escript::unaryOp(*this,escript::atanh_substitute);
969 #else
970 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
971 #endif
972 }
973
974 Data
975 Data::log10() const
976 {
977 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
978 }
979
980 Data
981 Data::log() const
982 {
983 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
984 }
985
986 Data
987 Data::sign() const
988 {
989 return escript::unaryOp(*this,escript::fsign);
990 }
991
992 Data
993 Data::abs() const
994 {
995 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
996 }
997
998 Data
999 Data::neg() const
1000 {
1001 return escript::unaryOp(*this,negate<double>());
1002 }
1003
1004 Data
1005 Data::pos() const
1006 {
1007 Data result;
1008 // perform a deep copy
1009 result.copy(*this);
1010 return result;
1011 }
1012
1013 Data
1014 Data::exp() const
1015 {
1016 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1017 }
1018
1019 Data
1020 Data::sqrt() const
1021 {
1022 return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1023 }
1024
1025 double
1026 Data::Lsup() const
1027 {
1028 double localValue, globalValue;
1029 //
1030 // set the initial absolute maximum value to zero
1031
1032 AbsMax abs_max_func;
1033 localValue = algorithm(abs_max_func,0);
1034 #ifdef PASO_MPI
1035 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1036 return globalValue;
1037 #else
1038 return localValue;
1039 #endif
1040 }
1041
1042 double
1043 Data::Linf() const
1044 {
1045 double localValue, globalValue;
1046 //
1047 // set the initial absolute minimum value to max double
1048 AbsMin abs_min_func;
1049 localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1050
1051 #ifdef PASO_MPI
1052 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1053 return globalValue;
1054 #else
1055 return localValue;
1056 #endif
1057 }
1058
1059 double
1060 Data::sup() const
1061 {
1062 double localValue, globalValue;
1063 //
1064 // set the initial maximum value to min possible double
1065 FMax fmax_func;
1066 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1067 #ifdef PASO_MPI
1068 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1069 return globalValue;
1070 #else
1071 return localValue;
1072 #endif
1073 }
1074
1075 double
1076 Data::inf() const
1077 {
1078 double localValue, globalValue;
1079 //
1080 // set the initial minimum value to max possible double
1081 FMin fmin_func;
1082 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1083 #ifdef PASO_MPI
1084 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1085 return globalValue;
1086 #else
1087 return localValue;
1088 #endif
1089 }
1090
1091 /* TODO */
1092 /* global reduction */
1093 Data
1094 Data::maxval() const
1095 {
1096 //
1097 // set the initial maximum value to min possible double
1098 FMax fmax_func;
1099 return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1100 }
1101
1102 Data
1103 Data::minval() const
1104 {
1105 //
1106 // set the initial minimum value to max possible double
1107 FMin fmin_func;
1108 return dp_algorithm(fmin_func,numeric_limits<double>::max());
1109 }
1110
1111 Data
1112 Data::swapaxes(const int axis0, const int axis1) const
1113 {
1114 int axis0_tmp,axis1_tmp;
1115 DataArrayView::ShapeType s=getDataPointShape();
1116 DataArrayView::ShapeType ev_shape;
1117 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1118 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1119 int rank=getDataPointRank();
1120 if (rank<2) {
1121 throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1122 }
1123 if (axis0<0 || axis0>rank-1) {
1124 throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1125 }
1126 if (axis1<0 || axis1>rank-1) {
1127 throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1128 }
1129 if (axis0 == axis1) {
1130 throw DataException("Error - Data::swapaxes: axis indices must be different.");
1131 }
1132 if (axis0 > axis1) {
1133 axis0_tmp=axis1;
1134 axis1_tmp=axis0;
1135 } else {
1136 axis0_tmp=axis0;
1137 axis1_tmp=axis1;
1138 }
1139 for (int i=0; i<rank; i++) {
1140 if (i == axis0_tmp) {
1141 ev_shape.push_back(s[axis1_tmp]);
1142 } else if (i == axis1_tmp) {
1143 ev_shape.push_back(s[axis0_tmp]);
1144 } else {
1145 ev_shape.push_back(s[i]);
1146 }
1147 }
1148 Data ev(0.,ev_shape,getFunctionSpace());
1149 ev.typeMatchRight(*this);
1150 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1151 return ev;
1152
1153 }
1154
1155 Data
1156 Data::symmetric() const
1157 {
1158 // check input
1159 DataArrayView::ShapeType s=getDataPointShape();
1160 if (getDataPointRank()==2) {
1161 if(s[0] != s[1])
1162 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1163 }
1164 else if (getDataPointRank()==4) {
1165 if(!(s[0] == s[2] && s[1] == s[3]))
1166 throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1167 }
1168 else {
1169 throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1170 }
1171 Data ev(0.,getDataPointShape(),getFunctionSpace());
1172 ev.typeMatchRight(*this);
1173 m_data->symmetric(ev.m_data.get());
1174 return ev;
1175 }
1176
1177 Data
1178 Data::nonsymmetric() const
1179 {
1180 // check input
1181 DataArrayView::ShapeType s=getDataPointShape();
1182 if (getDataPointRank()==2) {
1183 if(s[0] != s[1])
1184 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1185 DataArrayView::ShapeType ev_shape;
1186 ev_shape.push_back(s[0]);
1187 ev_shape.push_back(s[1]);
1188 Data ev(0.,ev_shape,getFunctionSpace());
1189 ev.typeMatchRight(*this);
1190 m_data->nonsymmetric(ev.m_data.get());
1191 return ev;
1192 }
1193 else if (getDataPointRank()==4) {
1194 if(!(s[0] == s[2] && s[1] == s[3]))
1195 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1196 DataArrayView::ShapeType ev_shape;
1197 ev_shape.push_back(s[0]);
1198 ev_shape.push_back(s[1]);
1199 ev_shape.push_back(s[2]);
1200 ev_shape.push_back(s[3]);
1201 Data ev(0.,ev_shape,getFunctionSpace());
1202 ev.typeMatchRight(*this);
1203 m_data->nonsymmetric(ev.m_data.get());
1204 return ev;
1205 }
1206 else {
1207 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1208 }
1209 }
1210
1211 Data
1212 Data::trace(int axis_offset) const
1213 {
1214 DataArrayView::ShapeType s=getDataPointShape();
1215 if (getDataPointRank()==2) {
1216 DataArrayView::ShapeType ev_shape;
1217 Data ev(0.,ev_shape,getFunctionSpace());
1218 ev.typeMatchRight(*this);
1219 m_data->trace(ev.m_data.get(), axis_offset);
1220 return ev;
1221 }
1222 if (getDataPointRank()==3) {
1223 DataArrayView::ShapeType ev_shape;
1224 if (axis_offset==0) {
1225 int s2=s[2];
1226 ev_shape.push_back(s2);
1227 }
1228 else if (axis_offset==1) {
1229 int s0=s[0];
1230 ev_shape.push_back(s0);
1231 }
1232 Data ev(0.,ev_shape,getFunctionSpace());
1233 ev.typeMatchRight(*this);
1234 m_data->trace(ev.m_data.get(), axis_offset);
1235 return ev;
1236 }
1237 if (getDataPointRank()==4) {
1238 DataArrayView::ShapeType ev_shape;
1239 if (axis_offset==0) {
1240 ev_shape.push_back(s[2]);
1241 ev_shape.push_back(s[3]);
1242 }
1243 else if (axis_offset==1) {
1244 ev_shape.push_back(s[0]);
1245 ev_shape.push_back(s[3]);
1246 }
1247 else if (axis_offset==2) {
1248 ev_shape.push_back(s[0]);
1249 ev_shape.push_back(s[1]);
1250 }
1251 Data ev(0.,ev_shape,getFunctionSpace());
1252 ev.typeMatchRight(*this);
1253 m_data->trace(ev.m_data.get(), axis_offset);
1254 return ev;
1255 }
1256 else {
1257 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1258 }
1259 }
1260
1261 Data
1262 Data::transpose(int axis_offset) const
1263 {
1264 DataArrayView::ShapeType s=getDataPointShape();
1265 DataArrayView::ShapeType ev_shape;
1266 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1267 // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1268 int rank=getDataPointRank();
1269 if (axis_offset<0 || axis_offset>rank) {
1270 throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1271 }
1272 for (int i=0; i<rank; i++) {
1273 int index = (axis_offset+i)%rank;
1274 ev_shape.push_back(s[index]); // Append to new shape
1275 }
1276 Data ev(0.,ev_shape,getFunctionSpace());
1277 ev.typeMatchRight(*this);
1278 m_data->transpose(ev.m_data.get(), axis_offset);
1279 return ev;
1280 }
1281
1282 Data
1283 Data::eigenvalues() const
1284 {
1285 // check input
1286 DataArrayView::ShapeType s=getDataPointShape();
1287 if (getDataPointRank()!=2)
1288 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1289 if(s[0] != s[1])
1290 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1291 // create return
1292 DataArrayView::ShapeType ev_shape(1,s[0]);
1293 Data ev(0.,ev_shape,getFunctionSpace());
1294 ev.typeMatchRight(*this);
1295 m_data->eigenvalues(ev.m_data.get());
1296 return ev;
1297 }
1298
1299 const boost::python::tuple
1300 Data::eigenvalues_and_eigenvectors(const double tol) const
1301 {
1302 DataArrayView::ShapeType s=getDataPointShape();
1303 if (getDataPointRank()!=2)
1304 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1305 if(s[0] != s[1])
1306 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1307 // create return
1308 DataArrayView::ShapeType ev_shape(1,s[0]);
1309 Data ev(0.,ev_shape,getFunctionSpace());
1310 ev.typeMatchRight(*this);
1311 DataArrayView::ShapeType V_shape(2,s[0]);
1312 Data V(0.,V_shape,getFunctionSpace());
1313 V.typeMatchRight(*this);
1314 m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1315 return make_tuple(boost::python::object(ev),boost::python::object(V));
1316 }
1317
1318 const boost::python::tuple
1319 Data::minGlobalDataPoint() const
1320 {
1321 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1322 // abort (for unknown reasons) if there are openmp directives with it in the
1323 // surrounding function
1324
1325 int DataPointNo;
1326 int ProcNo;
1327 calc_minGlobalDataPoint(ProcNo,DataPointNo);
1328 return make_tuple(ProcNo,DataPointNo);
1329 }
1330
1331 void
1332 Data::calc_minGlobalDataPoint(int& ProcNo,
1333 int& DataPointNo) const
1334 {
1335 int i,j;
1336 int lowi=0,lowj=0;
1337 double min=numeric_limits<double>::max();
1338
1339 Data temp=minval();
1340
1341 int numSamples=temp.getNumSamples();
1342 int numDPPSample=temp.getNumDataPointsPerSample();
1343
1344 double next,local_min;
1345 int local_lowi,local_lowj;
1346
1347 #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1348 {
1349 local_min=min;
1350 #pragma omp for private(i,j) schedule(static)
1351 for (i=0; i<numSamples; i++) {
1352 for (j=0; j<numDPPSample; j++) {
1353 next=temp.getDataPoint(i,j)();
1354 if (next<local_min) {
1355 local_min=next;
1356 local_lowi=i;
1357 local_lowj=j;
1358 }
1359 }
1360 }
1361 #pragma omp critical
1362 if (local_min<min) {
1363 min=local_min;
1364 lowi=local_lowi;
1365 lowj=local_lowj;
1366 }
1367 }
1368
1369 #ifdef PASO_MPI
1370 // determine the processor on which the minimum occurs
1371 next = temp.getDataPoint(lowi,lowj)();
1372 int lowProc = 0;
1373 double *globalMins = new double[get_MPISize()+1];
1374 int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1375
1376 if( get_MPIRank()==0 ){
1377 next = globalMins[lowProc];
1378 for( i=1; i<get_MPISize(); i++ )
1379 if( next>globalMins[i] ){
1380 lowProc = i;
1381 next = globalMins[i];
1382 }
1383 }
1384 MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1385
1386 delete [] globalMins;
1387 ProcNo = lowProc;
1388 #else
1389 ProcNo = 0;
1390 #endif
1391 DataPointNo = lowj + lowi * numDPPSample;
1392 }
1393
1394 void
1395 Data::saveDX(std::string fileName) const
1396 {
1397 boost::python::dict args;
1398 args["data"]=boost::python::object(this);
1399 getDomain().saveDX(fileName,args);
1400 return;
1401 }
1402
1403 void
1404 Data::saveVTK(std::string fileName) const
1405 {
1406 boost::python::dict args;
1407 args["data"]=boost::python::object(this);
1408 getDomain().saveVTK(fileName,args);
1409 return;
1410 }
1411
1412 Data&
1413 Data::operator+=(const Data& right)
1414 {
1415 if (isProtected()) {
1416 throw DataException("Error - attempt to update protected Data object.");
1417 }
1418 binaryOp(right,plus<double>());
1419 return (*this);
1420 }
1421
1422 Data&
1423 Data::operator+=(const boost::python::object& right)
1424 {
1425 Data tmp(right,getFunctionSpace(),false);
1426 binaryOp(tmp,plus<double>());
1427 return (*this);
1428 }
1429
1430 Data&
1431 Data::operator-=(const Data& right)
1432 {
1433 if (isProtected()) {
1434 throw DataException("Error - attempt to update protected Data object.");
1435 }
1436 binaryOp(right,minus<double>());
1437 return (*this);
1438 }
1439
1440 Data&
1441 Data::operator-=(const boost::python::object& right)
1442 {
1443 Data tmp(right,getFunctionSpace(),false);
1444 binaryOp(tmp,minus<double>());
1445 return (*this);
1446 }
1447
1448 Data&
1449 Data::operator*=(const Data& right)
1450 {
1451 if (isProtected()) {
1452 throw DataException("Error - attempt to update protected Data object.");
1453 }
1454 binaryOp(right,multiplies<double>());
1455 return (*this);
1456 }
1457
1458 Data&
1459 Data::operator*=(const boost::python::object& right)
1460 {
1461 Data tmp(right,getFunctionSpace(),false);
1462 binaryOp(tmp,multiplies<double>());
1463 return (*this);
1464 }
1465
1466 Data&
1467 Data::operator/=(const Data& right)
1468 {
1469 if (isProtected()) {
1470 throw DataException("Error - attempt to update protected Data object.");
1471 }
1472 binaryOp(right,divides<double>());
1473 return (*this);
1474 }
1475
1476 Data&
1477 Data::operator/=(const boost::python::object& right)
1478 {
1479 Data tmp(right,getFunctionSpace(),false);
1480 binaryOp(tmp,divides<double>());
1481 return (*this);
1482 }
1483
1484 Data
1485 Data::rpowO(const boost::python::object& left) const
1486 {
1487 Data left_d(left,*this);
1488 return left_d.powD(*this);
1489 }
1490
1491 Data
1492 Data::powO(const boost::python::object& right) const
1493 {
1494 Data tmp(right,getFunctionSpace(),false);
1495 return powD(tmp);
1496 }
1497
1498 Data
1499 Data::powD(const Data& right) const
1500 {
1501 Data result;
1502 if (getDataPointRank()<right.getDataPointRank()) {
1503 result.copy(right);
1504 result.binaryOp(*this,escript::rpow);
1505 } else {
1506 result.copy(*this);
1507 result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1508 }
1509 return result;
1510 }
1511
1512
1513 //
1514 // NOTE: It is essential to specify the namespace this operator belongs to
1515 Data
1516 escript::operator+(const Data& left, const Data& right)
1517 {
1518 Data result;
1519 //
1520 // perform a deep copy
1521 if (left.getDataPointRank()<right.getDataPointRank()) {
1522 result.copy(right);
1523 result+=left;
1524 } else {
1525 result.copy(left);
1526 result+=right;
1527 }
1528 return result;
1529 }
1530
1531 //
1532 // NOTE: It is essential to specify the namespace this operator belongs to
1533 Data
1534 escript::operator-(const Data& left, const Data& right)
1535 {
1536 Data result;
1537 //
1538 // perform a deep copy
1539 if (left.getDataPointRank()<right.getDataPointRank()) {
1540 result=right.neg();
1541 result+=left;
1542 } else {
1543 result.copy(left);
1544 result-=right;
1545 }
1546 return result;
1547 }
1548
1549 //
1550 // NOTE: It is essential to specify the namespace this operator belongs to
1551 Data
1552 escript::operator*(const Data& left, const Data& right)
1553 {
1554 Data result;
1555 //
1556 // perform a deep copy
1557 if (left.getDataPointRank()<right.getDataPointRank()) {
1558 result.copy(right);
1559 result*=left;
1560 } else {
1561 result.copy(left);
1562 result*=right;
1563 }
1564 return result;
1565 }
1566
1567 //
1568 // NOTE: It is essential to specify the namespace this operator belongs to
1569 Data
1570 escript::operator/(const Data& left, const Data& right)
1571 {
1572 Data result;
1573 //
1574 // perform a deep copy
1575 if (left.getDataPointRank()<right.getDataPointRank()) {
1576 result=right.oneOver();
1577 result*=left;
1578 } else {
1579 result.copy(left);
1580 result/=right;
1581 }
1582 return result;
1583 }
1584
1585 //
1586 // NOTE: It is essential to specify the namespace this operator belongs to
1587 Data
1588 escript::operator+(const Data& left, const boost::python::object& right)
1589 {
1590 return left+Data(right,left.getFunctionSpace(),false);
1591 }
1592
1593 //
1594 // NOTE: It is essential to specify the namespace this operator belongs to
1595 Data
1596 escript::operator-(const Data& left, const boost::python::object& right)
1597 {
1598 return left-Data(right,left.getFunctionSpace(),false);
1599 }
1600
1601 //
1602 // NOTE: It is essential to specify the namespace this operator belongs to
1603 Data
1604 escript::operator*(const Data& left, const boost::python::object& right)
1605 {
1606 return left*Data(right,left.getFunctionSpace(),false);
1607 }
1608
1609 //
1610 // NOTE: It is essential to specify the namespace this operator belongs to
1611 Data
1612 escript::operator/(const Data& left, const boost::python::object& right)
1613 {
1614 return left/Data(right,left.getFunctionSpace(),false);
1615 }
1616
1617 //
1618 // NOTE: It is essential to specify the namespace this operator belongs to
1619 Data
1620 escript::operator+(const boost::python::object& left, const Data& right)
1621 {
1622 return Data(left,right.getFunctionSpace(),false)+right;
1623 }
1624
1625 //
1626 // NOTE: It is essential to specify the namespace this operator belongs to
1627 Data
1628 escript::operator-(const boost::python::object& left, const Data& right)
1629 {
1630 return Data(left,right.getFunctionSpace(),false)-right;
1631 }
1632
1633 //
1634 // NOTE: It is essential to specify the namespace this operator belongs to
1635 Data
1636 escript::operator*(const boost::python::object& left, const Data& right)
1637 {
1638 return Data(left,right.getFunctionSpace(),false)*right;
1639 }
1640
1641 //
1642 // NOTE: It is essential to specify the namespace this operator belongs to
1643 Data
1644 escript::operator/(const boost::python::object& left, const Data& right)
1645 {
1646 return Data(left,right.getFunctionSpace(),false)/right;
1647 }
1648
1649 //
1650 //bool escript::operator==(const Data& left, const Data& right)
1651 //{
1652 // /*
1653 // NB: this operator does very little at this point, and isn't to
1654 // be relied on. Requires further implementation.
1655 // */
1656 //
1657 // bool ret;
1658 //
1659 // if (left.isEmpty()) {
1660 // if(!right.isEmpty()) {
1661 // ret = false;
1662 // } else {
1663 // ret = true;
1664 // }
1665 // }
1666 //
1667 // if (left.isConstant()) {
1668 // if(!right.isConstant()) {
1669 // ret = false;
1670 // } else {
1671 // ret = true;
1672 // }
1673 // }
1674 //
1675 // if (left.isTagged()) {
1676 // if(!right.isTagged()) {
1677 // ret = false;
1678 // } else {
1679 // ret = true;
1680 // }
1681 // }
1682 //
1683 // if (left.isExpanded()) {
1684 // if(!right.isExpanded()) {
1685 // ret = false;
1686 // } else {
1687 // ret = true;
1688 // }
1689 // }
1690 //
1691 // return ret;
1692 //}
1693
1694 /* TODO */
1695 /* global reduction */
1696 Data
1697 Data::getItem(const boost::python::object& key) const
1698 {
1699 const DataArrayView& view=getPointDataView();
1700
1701 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1702
1703 if (slice_region.size()!=view.getRank()) {
1704 throw DataException("Error - slice size does not match Data rank.");
1705 }
1706
1707 return getSlice(slice_region);
1708 }
1709
1710 /* TODO */
1711 /* global reduction */
1712 Data
1713 Data::getSlice(const DataArrayView::RegionType& region) const
1714 {
1715 return Data(*this,region);
1716 }
1717
1718 /* TODO */
1719 /* global reduction */
1720 void
1721 Data::setItemO(const boost::python::object& key,
1722 const boost::python::object& value)
1723 {
1724 Data tempData(value,getFunctionSpace());
1725 setItemD(key,tempData);
1726 }
1727
1728 /* TODO */
1729 /* global reduction */
1730 void
1731 Data::setItemD(const boost::python::object& key,
1732 const Data& value)
1733 {
1734 const DataArrayView& view=getPointDataView();
1735
1736 DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1737 if (slice_region.size()!=view.getRank()) {
1738 throw DataException("Error - slice size does not match Data rank.");
1739 }
1740 if (getFunctionSpace()!=value.getFunctionSpace()) {
1741 setSlice(Data(value,getFunctionSpace()),slice_region);
1742 } else {
1743 setSlice(value,slice_region);
1744 }
1745 }
1746
1747 /* TODO */
1748 /* global reduction */
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
2573 int
2574 Data::get_MPISize() const
2575 {
2576 int error, size;
2577 #ifdef PASO_MPI
2578 error = MPI_Comm_size( get_MPIComm(), &size );
2579 #else
2580 size = 1;
2581 #endif
2582 return size;
2583 }
2584
2585 int
2586 Data::get_MPIRank() const
2587 {
2588 int error, rank;
2589 #ifdef PASO_MPI
2590 error = MPI_Comm_rank( get_MPIComm(), &rank );
2591 #else
2592 rank = 0;
2593 #endif
2594 return rank;
2595 }
2596
2597 MPI_Comm
2598 Data::get_MPIComm() const
2599 {
2600 #ifdef PASO_MPI
2601 return MPI_COMM_WORLD;
2602 #else
2603 return -1;
2604 #endif
2605 }
2606

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26