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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26