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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26