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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2653 - (show annotations)
Tue Sep 8 04:26:30 2009 UTC (10 years, 4 months ago) by jfenwick
File size: 90989 byte(s)
Fix bug in maxGlobalDataPoint and minGlobalDataPoint.
They now give the correct answers and the datapoint ids returned are globally
correct.

Removed some #defines from before COW
Removed hasNoSamples() - I don't trust myself to use that properly let alone anybody else.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26