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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1803 - (show annotations)
Wed Sep 24 06:20:29 2008 UTC (10 years, 10 months ago) by jfenwick
File size: 73760 byte(s)
All about making DataEmpty instances throw.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Exposed getDim from AbstractDomain to python to fix bug.

Added isEmpty member to DataAbstract to allow it to throw is queries are 
made about a DataEmpty instance.


Added exceptions to DataAbstract, DataEmpty and Data to prevent calls 
being made against DataEmpty objects.
The following still work as expected on DataEmpty instances

copy, getDomain, getFunctionSpace, isEmpty, isExpanded, isProtected, 
isTagged, setprotection.

You can also call interpolate, however it should throw if you try to 
change FunctionSpaces.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26