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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2635 - (show annotations)
Thu Aug 27 04:54:41 2009 UTC (10 years, 7 months ago) by jfenwick
File size: 89284 byte(s)
A bunch of changes related to saveDataCSV.
[Not completed or unit tested yet]

Added saveDataCSV to util.py
AbstractDomain (and MeshAdapter) have a commonFunctionSpace method to 
take a group of FunctionSpaces and return something they can all be interpolated to.

Added pointToStream() in DataTypes to help print points.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26