/[escript]/branches/arrayview_from_1695_trunk/escript/src/Data.cpp
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1752 - (show annotations)
Fri Sep 5 06:26:29 2008 UTC (10 years, 8 months ago) by jfenwick
File size: 83743 byte(s)
Branch commit

Now passes all c++ unit tests but fails some python ones.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26