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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2005 - (show annotations)
Mon Nov 10 01:21:39 2008 UTC (11 years, 4 months ago) by jfenwick
File size: 85415 byte(s)
Bringing all changes across from schroedinger.
(Note this does not mean development is done, just that it will happen
on the trunk for now).
If anyone notices any problems please contact me.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26