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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2246 - (show annotations)
Thu Feb 5 06:04:55 2009 UTC (11 years, 2 months ago) by jfenwick
File size: 85997 byte(s)
Made some changes to the way copyWithMask works.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26