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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26