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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2721 - (hide annotations)
Fri Oct 16 05:40:12 2009 UTC (10 years, 5 months ago) by jfenwick
File size: 96572 byte(s)
minval and maxval are now lazy operations (they weren't before).
Whether or not Lsup, sup and inf resolve their arguments before computing answers is controlled by the escriptParam 'RESOLVE_COLLECTIVE'.
Note: integrate() still forces a resolve.

Added some unit tests for operations which weren't tested before.
Added deepcopy implementations for lazy operations which got missed somehow.
1 jgs 480
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 474 #include "Data.h"
16 jgs 94
17 jgs 480 #include "DataExpanded.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20     #include "DataEmpty.h"
21 jfenwick 2005 #include "DataLazy.h"
22 jgs 480 #include "FunctionSpaceFactory.h"
23     #include "AbstractContinuousDomain.h"
24     #include "UnaryFuncs.h"
25 jfenwick 1796 #include "FunctionSpaceException.h"
26 jfenwick 1897 #include "EscriptParams.h"
27 jfenwick 1796
28 ksteube 1312 extern "C" {
29 phornby 2078 #include "esysUtils/blocktimer.h"
30 ksteube 1312 }
31 jgs 480
32 jgs 119 #include <fstream>
33 jgs 94 #include <algorithm>
34     #include <vector>
35     #include <functional>
36 jfenwick 2086 #include <sstream> // so we can throw messages about ranks
37 jgs 94
38 jgs 153 #include <boost/python/dict.hpp>
39 jgs 94 #include <boost/python/extract.hpp>
40     #include <boost/python/long.hpp>
41 jfenwick 2271 #include "WrappedArray.h"
42 jgs 94
43     using namespace std;
44     using namespace boost::python;
45     using namespace boost;
46     using namespace escript;
47    
48 jfenwick 2005 // 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 jfenwick 2271 // #define forceResolve() if (isLazy()) {#resolve();}
51    
52 jfenwick 2146 #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 jfenwick 2005
64 jfenwick 2496 #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 jfenwick 2146 #define MAKELAZYBINSELF(R,X) if (isLazy() || R.isLazy() || (AUTOLAZYON && (isExpanded() || R.isExpanded()))) \
71     {\
72     DataLazy* c=new DataLazy(m_data,R.borrowDataPtr(),X);\
73 jfenwick 2271 /* m_data=c->getPtr();*/ set_m_data(c->getPtr());\
74 jfenwick 2146 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 jfenwick 2721 #define CHECK_DO_CRES escriptParams.getRESOLVE_COLLECTIVE()
91    
92 jfenwick 2271 namespace
93     {
94    
95 jfenwick 2459 template <class ARR>
96     inline
97     boost::python::tuple
98     pointToTuple1(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
99     {
100     using namespace boost::python;
101     using boost::python::tuple;
102     using boost::python::list;
103    
104     list l;
105     unsigned int dim=shape[0];
106     for (size_t i=0;i<dim;++i)
107     {
108     l.append(v[i+offset]);
109     }
110     return tuple(l);
111     }
112    
113     template <class ARR>
114     inline
115     boost::python::tuple
116     pointToTuple2(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
117     {
118     using namespace boost::python;
119     using boost::python::tuple;
120     using boost::python::list;
121    
122     unsigned int shape0=shape[0];
123     unsigned int shape1=shape[1];
124     list lj;
125     for (size_t j=0;j<shape0;++j)
126     {
127     list li;
128     for (size_t i=0;i<shape1;++i)
129     {
130     li.append(v[offset+DataTypes::getRelIndex(shape,j,i)]);
131     }
132     lj.append(tuple(li));
133     }
134     return tuple(lj);
135     }
136    
137     template <class ARR>
138     inline
139     boost::python::tuple
140     pointToTuple3(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
141     {
142     using namespace boost::python;
143     using boost::python::tuple;
144     using boost::python::list;
145    
146     unsigned int shape0=shape[0];
147     unsigned int shape1=shape[1];
148     unsigned int shape2=shape[2];
149    
150     list lk;
151     for (size_t k=0;k<shape0;++k)
152     {
153     list lj;
154     for (size_t j=0;j<shape1;++j)
155     {
156     list li;
157     for (size_t i=0;i<shape2;++i)
158     {
159 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,k,j,i)]);
160 jfenwick 2459 }
161     lj.append(tuple(li));
162     }
163     lk.append(tuple(lj));
164     }
165     return tuple(lk);
166     }
167    
168     template <class ARR>
169     inline
170     boost::python::tuple
171     pointToTuple4(const DataTypes::ShapeType& shape, ARR v, unsigned long offset)
172     {
173     using namespace boost::python;
174     using boost::python::tuple;
175     using boost::python::list;
176    
177     unsigned int shape0=shape[0];
178     unsigned int shape1=shape[1];
179     unsigned int shape2=shape[2];
180     unsigned int shape3=shape[3];
181    
182     list ll;
183     for (size_t l=0;l<shape0;++l)
184     {
185     list lk;
186     for (size_t k=0;k<shape1;++k)
187     {
188     list lj;
189     for (size_t j=0;j<shape2;++j)
190     {
191     list li;
192     for (size_t i=0;i<shape3;++i)
193     {
194 jfenwick 2482 li.append(v[offset+DataTypes::getRelIndex(shape,l,k,j,i)]);
195 jfenwick 2459 }
196     lj.append(tuple(li));
197     }
198     lk.append(tuple(lj));
199     }
200     ll.append(tuple(lk));
201     }
202     return tuple(ll);
203     }
204    
205    
206 jfenwick 2271 // This should be safer once the DataC RO changes have been brought in
207     template <class ARR>
208     boost::python::tuple
209     pointToTuple( const DataTypes::ShapeType& shape,ARR v)
210     {
211     using namespace boost::python;
212     using boost::python::list;
213     int rank=shape.size();
214     if (rank==0)
215     {
216     return make_tuple(v[0]);
217     }
218     else if (rank==1)
219     {
220 jfenwick 2459 return pointToTuple1(shape,v,0);
221 jfenwick 2271 }
222     else if (rank==2)
223     {
224 jfenwick 2459 return pointToTuple2(shape,v,0);
225 jfenwick 2271 }
226     else if (rank==3)
227     {
228 jfenwick 2459 return pointToTuple3(shape,v,0);
229 jfenwick 2271 }
230     else if (rank==4)
231     {
232 jfenwick 2459 return pointToTuple4(shape,v,0);
233 jfenwick 2271 }
234     else
235     throw DataException("Unknown rank in pointToTuple.");
236     }
237    
238     } // anonymous namespace
239    
240 jgs 94 Data::Data()
241 jfenwick 2271 : m_shared(false), m_lazy(false)
242 jgs 94 {
243     //
244     // Default data is type DataEmpty
245     DataAbstract* temp=new DataEmpty();
246 jfenwick 2271 // m_data=temp->getPtr();
247     set_m_data(temp->getPtr());
248 gross 783 m_protected=false;
249 jgs 94 }
250    
251     Data::Data(double value,
252     const tuple& shape,
253     const FunctionSpace& what,
254     bool expanded)
255 jfenwick 2271 : m_shared(false), m_lazy(false)
256 jgs 94 {
257 jfenwick 1796 DataTypes::ShapeType dataPointShape;
258 jgs 94 for (int i = 0; i < shape.attr("__len__")(); ++i) {
259     dataPointShape.push_back(extract<const int>(shape[i]));
260     }
261 matt 1319
262 jfenwick 1796 int len = DataTypes::noValues(dataPointShape);
263 matt 1319 DataVector temp_data(len,value,len);
264 jfenwick 1796 initialise(temp_data, dataPointShape, what, expanded);
265 gross 783 m_protected=false;
266 jgs 94 }
267    
268     Data::Data(double value,
269 jfenwick 1796 const DataTypes::ShapeType& dataPointShape,
270 jgs 94 const FunctionSpace& what,
271     bool expanded)
272 jfenwick 2271 : m_shared(false), m_lazy(false)
273 jgs 94 {
274 jfenwick 1796 int len = DataTypes::noValues(dataPointShape);
275 matt 1319 DataVector temp_data(len,value,len);
276 jfenwick 1796 initialise(temp_data, dataPointShape, what, expanded);
277 gross 783 m_protected=false;
278 jgs 94 }
279    
280 jgs 102 Data::Data(const Data& inData)
281 jfenwick 2271 : m_shared(false), m_lazy(false)
282 jgs 94 {
283 jfenwick 2271 set_m_data(inData.m_data);
284 gross 783 m_protected=inData.isProtected();
285 jgs 94 }
286    
287 jfenwick 1796
288 jgs 94 Data::Data(const Data& inData,
289 jfenwick 1796 const DataTypes::RegionType& region)
290 jfenwick 2271 : m_shared(false), m_lazy(false)
291 jgs 94 {
292 jfenwick 2005 DataAbstract_ptr dat=inData.m_data;
293     if (inData.isLazy())
294     {
295     dat=inData.m_data->resolve();
296     }
297     else
298     {
299     dat=inData.m_data;
300     }
301 jgs 94 //
302 jgs 102 // Create Data which is a slice of another Data
303 jfenwick 2005 DataAbstract* tmp = dat->getSlice(region);
304 jfenwick 2271 set_m_data(DataAbstract_ptr(tmp));
305 gross 783 m_protected=false;
306 jfenwick 2271
307 jgs 94 }
308    
309     Data::Data(const Data& inData,
310     const FunctionSpace& functionspace)
311 jfenwick 2271 : m_shared(false), m_lazy(false)
312 jgs 94 {
313 jfenwick 1803 if (inData.isEmpty())
314     {
315     throw DataException("Error - will not interpolate for instances of DataEmpty.");
316     }
317 jgs 94 if (inData.getFunctionSpace()==functionspace) {
318 jfenwick 2271 set_m_data(inData.m_data);
319 jfenwick 2005 }
320     else
321     {
322    
323     if (inData.isConstant()) { // for a constant function, we just need to use the new function space
324     if (!inData.probeInterpolation(functionspace))
325     { // Even though this is constant, we still need to check whether interpolation is allowed
326 jfenwick 2087 throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
327 jfenwick 2005 }
328     // if the data is not lazy, this will just be a cast to DataReady
329     DataReady_ptr dr=inData.m_data->resolve();
330 jfenwick 2271 DataConstant* dc=new DataConstant(functionspace,inData.m_data->getShape(),dr->getVectorRO());
331     // m_data=DataAbstract_ptr(dc);
332     set_m_data(DataAbstract_ptr(dc));
333 jgs 94 } else {
334 jfenwick 2005 Data tmp(0,inData.getDataPointShape(),functionspace,true);
335     // Note: Must use a reference or pointer to a derived object
336     // in order to get polymorphic behaviour. Shouldn't really
337     // be able to create an instance of AbstractDomain but that was done
338     // as a boost:python work around which may no longer be required.
339     /*const AbstractDomain& inDataDomain=inData.getDomain();*/
340     const_Domain_ptr inDataDomain=inData.getDomain();
341     if (inDataDomain==functionspace.getDomain()) {
342     inDataDomain->interpolateOnDomain(tmp,inData);
343     } else {
344     inDataDomain->interpolateACross(tmp,inData);
345     }
346 jfenwick 2271 // m_data=tmp.m_data;
347     set_m_data(tmp.m_data);
348 jgs 94 }
349     }
350 gross 783 m_protected=false;
351 jgs 94 }
352    
353 jfenwick 1796 Data::Data(DataAbstract* underlyingdata)
354 jfenwick 2271 : m_shared(false), m_lazy(false)
355 jgs 94 {
356 jfenwick 2271 set_m_data(underlyingdata->getPtr());
357 jfenwick 1796 m_protected=false;
358 jgs 94 }
359    
360 jfenwick 2005 Data::Data(DataAbstract_ptr underlyingdata)
361 jfenwick 2271 : m_shared(false), m_lazy(false)
362 jfenwick 2005 {
363 jfenwick 2271 set_m_data(underlyingdata);
364 jfenwick 2005 m_protected=false;
365     }
366    
367 jfenwick 1796 Data::Data(const DataTypes::ValueType& value,
368     const DataTypes::ShapeType& shape,
369     const FunctionSpace& what,
370     bool expanded)
371 jfenwick 2271 : m_shared(false), m_lazy(false)
372 jfenwick 1796 {
373     initialise(value,shape,what,expanded);
374     m_protected=false;
375 jgs 94 }
376    
377 jfenwick 1796
378 jgs 94 Data::Data(const object& value,
379     const FunctionSpace& what,
380     bool expanded)
381 jfenwick 2271 : m_shared(false), m_lazy(false)
382 jgs 94 {
383 jfenwick 2271 WrappedArray w(value);
384     initialise(w,what,expanded);
385 gross 783 m_protected=false;
386 jgs 94 }
387    
388 matt 1319
389 jgs 94 Data::Data(const object& value,
390     const Data& other)
391 jfenwick 2271 : m_shared(false), m_lazy(false)
392 jgs 94 {
393 jfenwick 2271 WrappedArray w(value);
394 matt 1319
395 jfenwick 2458 // extract the shape of the array
396 jfenwick 2271 const DataTypes::ShapeType& tempShape=w.getShape();
397     if (w.getRank()==0) {
398 matt 1319
399    
400 jfenwick 1796 // get the space for the data vector
401     int len1 = DataTypes::noValues(tempShape);
402     DataVector temp_data(len1, 0.0, len1);
403 jfenwick 2271 temp_data.copyFromArray(w,1);
404 jfenwick 1796
405     int len = DataTypes::noValues(other.getDataPointShape());
406    
407 jfenwick 2271 DataVector temp2_data(len, temp_data[0], len);
408 jfenwick 1796 DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
409 jfenwick 2271 // m_data=DataAbstract_ptr(t);
410     set_m_data(DataAbstract_ptr(t));
411 jfenwick 1796
412 jgs 94 } else {
413     //
414     // Create a DataConstant with the same sample shape as other
415 jfenwick 2271 DataConstant* t=new DataConstant(w,other.getFunctionSpace());
416     // m_data=DataAbstract_ptr(t);
417     set_m_data(DataAbstract_ptr(t));
418 jgs 94 }
419 gross 783 m_protected=false;
420 jgs 94 }
421    
422 jgs 151 Data::~Data()
423     {
424 jfenwick 2271 set_m_data(DataAbstract_ptr());
425 jgs 151 }
426    
427 jfenwick 1796
428 jfenwick 2271 // only call in thread safe contexts.
429     // This method should be atomic
430     void Data::set_m_data(DataAbstract_ptr p)
431     {
432     if (m_data.get()!=0) // release old ownership
433     {
434     m_data->removeOwner(this);
435     }
436     if (p.get()!=0)
437     {
438     m_data=p;
439     m_data->addOwner(this);
440     m_shared=m_data->isShared();
441     m_lazy=m_data->isLazy();
442     }
443     }
444 jfenwick 1796
445 jfenwick 2271 void Data::initialise(const WrappedArray& value,
446 jfenwick 1796 const FunctionSpace& what,
447     bool expanded)
448     {
449     //
450     // Construct a Data object of the appropriate type.
451     // Construct the object first as there seems to be a bug which causes
452     // undefined behaviour if an exception is thrown during construction
453     // within the shared_ptr constructor.
454     if (expanded) {
455     DataAbstract* temp=new DataExpanded(value, what);
456 jfenwick 2271 // m_data=temp->getPtr();
457     set_m_data(temp->getPtr());
458 jfenwick 1796 } else {
459     DataAbstract* temp=new DataConstant(value, what);
460 jfenwick 2271 // m_data=temp->getPtr();
461     set_m_data(temp->getPtr());
462 jfenwick 1796 }
463     }
464    
465    
466     void
467     Data::initialise(const DataTypes::ValueType& value,
468     const DataTypes::ShapeType& shape,
469     const FunctionSpace& what,
470     bool expanded)
471     {
472     //
473     // Construct a Data object of the appropriate type.
474     // Construct the object first as there seems to be a bug which causes
475     // undefined behaviour if an exception is thrown during construction
476     // within the shared_ptr constructor.
477     if (expanded) {
478     DataAbstract* temp=new DataExpanded(what, shape, value);
479 jfenwick 2271 // m_data=temp->getPtr();
480     set_m_data(temp->getPtr());
481 jfenwick 1796 } else {
482     DataAbstract* temp=new DataConstant(what, shape, value);
483 jfenwick 2271 // m_data=temp->getPtr();
484     set_m_data(temp->getPtr());
485 jfenwick 1796 }
486     }
487    
488    
489 jgs 94 escriptDataC
490     Data::getDataC()
491     {
492     escriptDataC temp;
493     temp.m_dataPtr=(void*)this;
494     return temp;
495     }
496    
497     escriptDataC
498     Data::getDataC() const
499     {
500     escriptDataC temp;
501     temp.m_dataPtr=(void*)this;
502     return temp;
503     }
504    
505 jfenwick 2271 size_t
506     Data::getSampleBufferSize() const
507     {
508     return m_data->getSampleBufferSize();
509     }
510    
511    
512 jgs 121 const boost::python::tuple
513 jgs 94 Data::getShapeTuple() const
514     {
515 jfenwick 1796 const DataTypes::ShapeType& shape=getDataPointShape();
516 jgs 94 switch(getDataPointRank()) {
517     case 0:
518     return make_tuple();
519     case 1:
520     return make_tuple(long_(shape[0]));
521     case 2:
522     return make_tuple(long_(shape[0]),long_(shape[1]));
523     case 3:
524     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]));
525     case 4:
526     return make_tuple(long_(shape[0]),long_(shape[1]),long_(shape[2]),long_(shape[3]));
527     default:
528     throw DataException("Error - illegal Data rank.");
529     }
530     }
531 jfenwick 1799
532    
533     // The different name is needed because boost has trouble with overloaded functions.
534     // It can't work out what type the function is based soley on its name.
535     // There are ways to fix this involving creating function pointer variables for each form
536     // but there doesn't seem to be a need given that the methods have the same name from the python point of view
537 jfenwick 2105 Data
538 jfenwick 1799 Data::copySelf()
539     {
540     DataAbstract* temp=m_data->deepCopy();
541 jfenwick 2105 return Data(temp);
542 jfenwick 1799 }
543    
544 jgs 94 void
545     Data::copy(const Data& other)
546     {
547 jfenwick 1799 DataAbstract* temp=other.m_data->deepCopy();
548 jfenwick 1872 DataAbstract_ptr p=temp->getPtr();
549 jfenwick 2271 // m_data=p;
550     set_m_data(p);
551 jgs 94 }
552    
553 gross 1118
554 jfenwick 2005 Data
555     Data::delay()
556     {
557 jfenwick 2721 if (!isLazy())
558     {
559     DataLazy* dl=new DataLazy(m_data);
560     return Data(dl);
561     }
562     return *this;
563 jfenwick 2005 }
564    
565 jgs 94 void
566 jfenwick 2005 Data::delaySelf()
567     {
568     if (!isLazy())
569     {
570 jfenwick 2271 // m_data=(new DataLazy(m_data))->getPtr();
571     set_m_data((new DataLazy(m_data))->getPtr());
572 jfenwick 2005 }
573     }
574    
575 jfenwick 2271
576     // For lazy data, it would seem that DataTagged will need to be treated differently since even after setting all tags
577     // to zero, all the tags from all the DataTags would be in the result.
578     // However since they all have the same value (0) whether they are there or not should not matter.
579     // So I have decided that for all types this method will create a constant 0.
580     // It can be promoted up as required.
581     // A possible efficiency concern might be expanded->constant->expanded which has an extra memory management
582     // but we can deal with that if it arises.
583     //
584 jfenwick 2005 void
585 gross 1118 Data::setToZero()
586     {
587 jfenwick 1803 if (isEmpty())
588 gross 1118 {
589 jfenwick 1803 throw DataException("Error - Operations not permitted on instances of DataEmpty.");
590     }
591 jfenwick 2271 if (isLazy())
592     {
593     DataTypes::ValueType v(getNoValues(),0);
594     DataConstant* dc=new DataConstant(getFunctionSpace(),getDataPointShape(),v);
595     DataLazy* dl=new DataLazy(dc->getPtr());
596     set_m_data(dl->getPtr());
597     }
598     else
599     {
600     exclusiveWrite();
601     m_data->setToZero();
602     }
603 gross 1118 }
604    
605 jfenwick 2271
606 gross 1118 void
607 jgs 94 Data::copyWithMask(const Data& other,
608     const Data& mask)
609     {
610 jfenwick 1855 // 1. Interpolate if required so all Datas use the same FS as this
611     // 2. Tag or Expand so that all Data's are the same type
612     // 3. Iterate over the data vectors copying values where mask is >0
613 jfenwick 1803 if (other.isEmpty() || mask.isEmpty())
614     {
615     throw DataException("Error - copyWithMask not permitted using instances of DataEmpty.");
616     }
617 jfenwick 1855 Data other2(other);
618     Data mask2(mask);
619 jfenwick 2005 other2.resolve();
620     mask2.resolve();
621     this->resolve();
622 jfenwick 1855 FunctionSpace myFS=getFunctionSpace();
623     FunctionSpace oFS=other2.getFunctionSpace();
624     FunctionSpace mFS=mask2.getFunctionSpace();
625     if (oFS!=myFS)
626     {
627     if (other2.probeInterpolation(myFS))
628     {
629     other2=other2.interpolate(myFS);
630     }
631     else
632     {
633     throw DataException("Error - copyWithMask: other FunctionSpace is not compatible with this one.");
634     }
635     }
636     if (mFS!=myFS)
637     {
638     if (mask2.probeInterpolation(myFS))
639     {
640     mask2=mask2.interpolate(myFS);
641     }
642     else
643     {
644     throw DataException("Error - copyWithMask: mask FunctionSpace is not compatible with this one.");
645     }
646     }
647     // Ensure that all args have the same type
648     if (this->isExpanded() || mask2.isExpanded() || other2.isExpanded())
649     {
650     this->expand();
651     other2.expand();
652     mask2.expand();
653     }
654     else if (this->isTagged() || mask2.isTagged() || other2.isTagged())
655     {
656     this->tag();
657     other2.tag();
658     mask2.tag();
659     }
660     else if (this->isConstant() && mask2.isConstant() && other2.isConstant())
661     {
662     }
663     else
664     {
665     throw DataException("Error - Unknown DataAbstract passed to copyWithMask.");
666     }
667 jfenwick 2246 unsigned int selfrank=getDataPointRank();
668     unsigned int otherrank=other2.getDataPointRank();
669     unsigned int maskrank=mask2.getDataPointRank();
670     if ((selfrank==0) && (otherrank>0 || maskrank>0))
671     {
672     // to get here we must be copying from a large object into a scalar
673     // I am not allowing this.
674     // If you are calling copyWithMask then you are considering keeping some existing values
675     // and so I'm going to assume that you don't want your data objects getting a new shape.
676     throw DataException("Attempt to copyWithMask into a scalar from an object or mask with rank>0.");
677     }
678 jfenwick 2271 exclusiveWrite();
679 jfenwick 1855 // Now we iterate over the elements
680 jfenwick 2271 DataVector& self=getReady()->getVectorRW();;
681     const DataVector& ovec=other2.getReadyPtr()->getVectorRO();
682     const DataVector& mvec=mask2.getReadyPtr()->getVectorRO();
683 jfenwick 2246
684     if ((selfrank>0) && (otherrank==0) &&(maskrank==0))
685 jfenwick 1855 {
686 jfenwick 2246 // Not allowing this combination.
687     // it is not clear what the rank of the target should be.
688     // Should it be filled with the scalar (rank stays the same);
689     // or should the target object be reshaped to be a scalar as well.
690     throw DataException("Attempt to copyWithMask from scalar mask and data into non-scalar target.");
691 jfenwick 1855 }
692 jfenwick 2246 if ((selfrank>0) && (otherrank>0) &&(maskrank==0))
693     {
694     if (mvec[0]>0) // copy whole object if scalar is >0
695     {
696     copy(other);
697     }
698     return;
699     }
700     if (isTagged()) // so all objects involved will also be tagged
701     {
702     // note the !
703     if (!((getDataPointShape()==mask2.getDataPointShape()) &&
704     ((other2.getDataPointShape()==mask2.getDataPointShape()) || (otherrank==0))))
705     {
706     throw DataException("copyWithMask, shape mismatch.");
707     }
708    
709     // We need to consider the possibility that tags are missing or in the wrong order
710     // My guiding assumption here is: All tagged Datas are assumed to have the default value for
711     // all tags which are not explicitly defined
712    
713     const DataTagged* mptr=dynamic_cast<const DataTagged*>(mask2.m_data.get());
714     const DataTagged* optr=dynamic_cast<const DataTagged*>(other2.m_data.get());
715     DataTagged* tptr=dynamic_cast<DataTagged*>(m_data.get());
716    
717     // first, add any tags missing from other or mask
718     const DataTagged::DataMapType& olookup=optr->getTagLookup();
719     const DataTagged::DataMapType& mlookup=mptr->getTagLookup();
720 jfenwick 2260 const DataTagged::DataMapType& tlookup=tptr->getTagLookup();
721 jfenwick 2246 DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
722     for (i=olookup.begin();i!=olookup.end();i++)
723     {
724     tptr->addTag(i->first);
725     }
726     for (i=mlookup.begin();i!=mlookup.end();i++) {
727     tptr->addTag(i->first);
728     }
729     // now we know that *this has all the required tags but they aren't guaranteed to be in
730     // the same order
731    
732     // There are two possibilities: 1. all objects have the same rank. 2. other is a scalar
733     if ((selfrank==otherrank) && (otherrank==maskrank))
734     {
735 jfenwick 2260 for (i=tlookup.begin();i!=tlookup.end();i++)
736 jfenwick 2246 {
737     // get the target offset
738     DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
739     DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
740     DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
741 jfenwick 2260 for (int j=0;j<getDataPointSize();++j)
742 jfenwick 2246 {
743 jfenwick 2260 if (mvec[j+moff]>0)
744 jfenwick 2246 {
745 jfenwick 2260 self[j+toff]=ovec[j+ooff];
746 jfenwick 2246 }
747     }
748     }
749 jfenwick 2260 // now for the default value
750     for (int j=0;j<getDataPointSize();++j)
751     {
752     if (mvec[j+mptr->getDefaultOffset()]>0)
753     {
754     self[j+tptr->getDefaultOffset()]=ovec[j+optr->getDefaultOffset()];
755     }
756     }
757 jfenwick 2246 }
758     else // other is a scalar
759     {
760 jfenwick 2260 for (i=tlookup.begin();i!=tlookup.end();i++)
761 jfenwick 2246 {
762     // get the target offset
763     DataTypes::ValueType::size_type toff=tptr->getOffsetForTag(i->first);
764     DataTypes::ValueType::size_type moff=mptr->getOffsetForTag(i->first);
765     DataTypes::ValueType::size_type ooff=optr->getOffsetForTag(i->first);
766 jfenwick 2260 for (int j=0;j<getDataPointSize();++j)
767 jfenwick 2246 {
768 jfenwick 2260 if (mvec[j+moff]>0)
769 jfenwick 2246 {
770 jfenwick 2260 self[j+toff]=ovec[ooff];
771 jfenwick 2246 }
772     }
773     }
774 jfenwick 2260 // now for the default value
775     for (int j=0;j<getDataPointSize();++j)
776     {
777     if (mvec[j+mptr->getDefaultOffset()]>0)
778     {
779     self[j+tptr->getDefaultOffset()]=ovec[0];
780     }
781     }
782 jfenwick 2246 }
783    
784     return; // ugly
785     }
786     // mixed scalar and non-scalar operation
787     if ((selfrank>0) && (otherrank==0) && (mask2.getDataPointShape()==getDataPointShape()))
788     {
789     size_t num_points=self.size();
790     // OPENMP 3.0 allows unsigned loop vars.
791     #if defined(_OPENMP) && (_OPENMP < 200805)
792     long i;
793     #else
794     size_t i;
795 jfenwick 2260 #endif
796     size_t psize=getDataPointSize();
797 jfenwick 2246 #pragma omp parallel for private(i) schedule(static)
798     for (i=0;i<num_points;++i)
799     {
800     if (mvec[i]>0)
801     {
802 jfenwick 2260 self[i]=ovec[i/psize]; // since this is expanded there is one scalar
803     } // dest point
804 jfenwick 2246 }
805     return; // ugly!
806     }
807     // tagged data is already taken care of so we only need to worry about shapes
808     // special cases with scalars are already dealt with so all we need to worry about is shape
809     if ((getDataPointShape()!=other2.getDataPointShape()) || getDataPointShape()!=mask2.getDataPointShape())
810     {
811     ostringstream oss;
812     oss <<"Error - size mismatch in arguments to copyWithMask.";
813     oss << "\nself_shape=" << DataTypes::shapeToString(getDataPointShape());
814     oss << " other2_shape=" << DataTypes::shapeToString(other2.getDataPointShape());
815     oss << " mask2_shape=" << DataTypes::shapeToString(mask2.getDataPointShape());
816     throw DataException(oss.str());
817     }
818 jfenwick 1855 size_t num_points=self.size();
819 phornby 1921
820     // OPENMP 3.0 allows unsigned loop vars.
821     #if defined(_OPENMP) && (_OPENMP < 200805)
822     long i;
823     #else
824 phornby 1918 size_t i;
825 phornby 1921 #endif
826 jfenwick 1855 #pragma omp parallel for private(i) schedule(static)
827     for (i=0;i<num_points;++i)
828     {
829     if (mvec[i]>0)
830     {
831     self[i]=ovec[i];
832     }
833     }
834     }
835 jgs 94
836     bool
837     Data::isExpanded() const
838     {
839     DataExpanded* temp=dynamic_cast<DataExpanded*>(m_data.get());
840     return (temp!=0);
841     }
842    
843     bool
844 jfenwick 2271 Data::actsExpanded() const
845     {
846     return m_data->actsExpanded();
847    
848     }
849    
850    
851     bool
852 jgs 94 Data::isTagged() const
853     {
854     DataTagged* temp=dynamic_cast<DataTagged*>(m_data.get());
855     return (temp!=0);
856     }
857    
858     bool
859     Data::isEmpty() const
860     {
861     DataEmpty* temp=dynamic_cast<DataEmpty*>(m_data.get());
862     return (temp!=0);
863     }
864    
865     bool
866     Data::isConstant() const
867     {
868     DataConstant* temp=dynamic_cast<DataConstant*>(m_data.get());
869     return (temp!=0);
870     }
871    
872 jfenwick 2005 bool
873     Data::isLazy() const
874     {
875 jfenwick 2271 return m_lazy; // not asking m_data because we need to be able to ask this while m_data is changing
876 jfenwick 2005 }
877    
878     // at the moment this is synonymous with !isLazy() but that could change
879     bool
880     Data::isReady() const
881     {
882     return (dynamic_cast<DataReady*>(m_data.get())!=0);
883     }
884    
885    
886 jgs 94 void
887 ksteube 1312 Data::setProtection()
888     {
889 gross 783 m_protected=true;
890     }
891    
892     bool
893 ksteube 1312 Data::isProtected() const
894     {
895 gross 783 return m_protected;
896     }
897    
898    
899    
900     void
901 jgs 94 Data::expand()
902     {
903     if (isConstant()) {
904     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
905     DataAbstract* temp=new DataExpanded(*tempDataConst);
906 jfenwick 2271 // m_data=temp->getPtr();
907     set_m_data(temp->getPtr());
908 jgs 94 } else if (isTagged()) {
909     DataTagged* tempDataTag=dynamic_cast<DataTagged*>(m_data.get());
910     DataAbstract* temp=new DataExpanded(*tempDataTag);
911 jfenwick 2271 // m_data=temp->getPtr();
912     set_m_data(temp->getPtr());
913 jgs 94 } else if (isExpanded()) {
914     //
915     // do nothing
916     } else if (isEmpty()) {
917     throw DataException("Error - Expansion of DataEmpty not possible.");
918 jfenwick 2005 } else if (isLazy()) {
919     resolve();
920     expand(); // resolve might not give us expanded data
921 jgs 94 } else {
922     throw DataException("Error - Expansion not implemented for this Data type.");
923     }
924     }
925    
926     void
927     Data::tag()
928     {
929     if (isConstant()) {
930     DataConstant* tempDataConst=dynamic_cast<DataConstant*>(m_data.get());
931     DataAbstract* temp=new DataTagged(*tempDataConst);
932 jfenwick 2271 // m_data=temp->getPtr();
933     set_m_data(temp->getPtr());
934 jgs 94 } else if (isTagged()) {
935     // do nothing
936     } else if (isExpanded()) {
937     throw DataException("Error - Creating tag data from DataExpanded not possible.");
938     } else if (isEmpty()) {
939     throw DataException("Error - Creating tag data from DataEmpty not possible.");
940 jfenwick 2005 } else if (isLazy()) {
941     DataAbstract_ptr res=m_data->resolve();
942     if (m_data->isExpanded())
943     {
944     throw DataException("Error - data would resolve to DataExpanded, tagging is not possible.");
945     }
946 jfenwick 2271 // m_data=res;
947     set_m_data(res);
948 jfenwick 2005 tag();
949 jgs 94 } else {
950     throw DataException("Error - Tagging not implemented for this Data type.");
951     }
952     }
953    
954 jfenwick 2005 void
955     Data::resolve()
956     {
957     if (isLazy())
958     {
959 jfenwick 2271 // m_data=m_data->resolve();
960     set_m_data(m_data->resolve());
961 jfenwick 2005 }
962     }
963    
964 jfenwick 2271 void
965     Data::requireWrite()
966     {
967     resolve();
968     exclusiveWrite();
969     }
970 jfenwick 2005
971 gross 854 Data
972     Data::oneOver() const
973 jgs 102 {
974 jfenwick 2146 MAKELAZYOP(RECIP)
975 matt 1334 return C_TensorUnaryOperation(*this, bind1st(divides<double>(),1.));
976 jgs 102 }
977    
978 jgs 94 Data
979 gross 698 Data::wherePositive() const
980 jgs 94 {
981 jfenwick 2146 MAKELAZYOP(GZ)
982 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater<double>(),0.0));
983 jgs 94 }
984    
985     Data
986 gross 698 Data::whereNegative() const
987 jgs 102 {
988 jfenwick 2146 MAKELAZYOP(LZ)
989 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less<double>(),0.0));
990 jgs 102 }
991    
992     Data
993 gross 698 Data::whereNonNegative() const
994 jgs 94 {
995 jfenwick 2146 MAKELAZYOP(GEZ)
996 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(greater_equal<double>(),0.0));
997 jgs 94 }
998    
999     Data
1000 gross 698 Data::whereNonPositive() const
1001 jgs 94 {
1002 jfenwick 2146 MAKELAZYOP(LEZ)
1003 matt 1334 return C_TensorUnaryOperation(*this, bind2nd(less_equal<double>(),0.0));
1004 jgs 94 }
1005    
1006     Data
1007 jgs 571 Data::whereZero(double tol) const
1008 jgs 94 {
1009 jfenwick 2147 // Data dataAbs=abs();
1010     // return C_TensorUnaryOperation(dataAbs, bind2nd(less_equal<double>(),tol));
1011     MAKELAZYOPOFF(EZ,tol)
1012     return C_TensorUnaryOperation(*this, bind2nd(AbsLTE(),tol));
1013    
1014 jgs 94 }
1015    
1016     Data
1017 jgs 571 Data::whereNonZero(double tol) const
1018 jgs 102 {
1019 jfenwick 2147 // Data dataAbs=abs();
1020     // return C_TensorUnaryOperation(dataAbs, bind2nd(greater<double>(),tol));
1021     MAKELAZYOPOFF(NEZ,tol)
1022     return C_TensorUnaryOperation(*this, bind2nd(AbsGT(),tol));
1023    
1024 jgs 102 }
1025    
1026     Data
1027 jgs 94 Data::interpolate(const FunctionSpace& functionspace) const
1028     {
1029     return Data(*this,functionspace);
1030     }
1031    
1032     bool
1033     Data::probeInterpolation(const FunctionSpace& functionspace) const
1034     {
1035 jfenwick 2005 return getFunctionSpace().probeInterpolation(functionspace);
1036 jgs 94 }
1037    
1038     Data
1039     Data::gradOn(const FunctionSpace& functionspace) const
1040     {
1041 jfenwick 1803 if (isEmpty())
1042     {
1043     throw DataException("Error - operation not permitted on instances of DataEmpty.");
1044     }
1045 matt 1353 double blocktimer_start = blocktimer_time();
1046 jgs 94 if (functionspace.getDomain()!=getDomain())
1047     throw DataException("Error - gradient cannot be calculated on different domains.");
1048 jfenwick 1796 DataTypes::ShapeType grad_shape=getDataPointShape();
1049 jgs 94 grad_shape.push_back(functionspace.getDim());
1050     Data out(0.0,grad_shape,functionspace,true);
1051 jfenwick 1872 getDomain()->setToGradient(out,*this);
1052 matt 1353 blocktimer_increment("grad()", blocktimer_start);
1053 jgs 94 return out;
1054     }
1055    
1056     Data
1057     Data::grad() const
1058     {
1059 jfenwick 1803 if (isEmpty())
1060     {
1061     throw DataException("Error - operation not permitted on instances of DataEmpty.");
1062     }
1063 jfenwick 1872 return gradOn(escript::function(*getDomain()));
1064 jgs 94 }
1065    
1066     int
1067     Data::getDataPointSize() const
1068     {
1069 jfenwick 1796 return m_data->getNoValues();
1070 jgs 94 }
1071    
1072 jfenwick 2459
1073 jfenwick 1796 DataTypes::ValueType::size_type
1074 jgs 94 Data::getLength() const
1075     {
1076     return m_data->getLength();
1077     }
1078    
1079 jfenwick 2481
1080 jfenwick 2459 // There is no parallelism here ... elements need to be added in the correct order.
1081     // If we could presize the list and then fill in the elements it might work
1082     // This would need setting elements to be threadsafe.
1083     // Having mulitple C threads calling into one interpreter is aparently a no-no.
1084 jfenwick 2271 const boost::python::object
1085 jfenwick 2459 Data::toListOfTuples(bool scalarastuple)
1086     {
1087     using namespace boost::python;
1088 jfenwick 2461 using boost::python::list;
1089 jfenwick 2459 if (get_MPISize()>1)
1090     {
1091     throw DataException("::toListOfTuples is not available for MPI with more than one process.");
1092     }
1093     unsigned int rank=getDataPointRank();
1094     unsigned int size=getDataPointSize();
1095 jfenwick 2461
1096 jfenwick 2459 int npoints=getNumDataPoints();
1097     expand(); // This will also resolve if required
1098     const DataTypes::ValueType& vec=getReady()->getVectorRO();
1099 jfenwick 2461 boost::python::list temp;
1100     temp.append(object());
1101     boost::python::list res(temp*npoints);// presize the list by the "[None] * npoints" trick
1102 jfenwick 2459 if (rank==0)
1103     {
1104     long count;
1105     if (scalarastuple)
1106     {
1107     for (count=0;count<npoints;++count)
1108     {
1109 jfenwick 2461 res[count]=make_tuple(vec[count]);
1110 jfenwick 2459 }
1111     }
1112     else
1113     {
1114     for (count=0;count<npoints;++count)
1115     {
1116 jfenwick 2461 res[count]=vec[count];
1117 jfenwick 2459 }
1118     }
1119     }
1120     else if (rank==1)
1121     {
1122     size_t count;
1123     size_t offset=0;
1124     for (count=0;count<npoints;++count,offset+=size)
1125     {
1126 jfenwick 2461 res[count]=pointToTuple1(getDataPointShape(), vec, offset);
1127 jfenwick 2459 }
1128     }
1129     else if (rank==2)
1130     {
1131     size_t count;
1132     size_t offset=0;
1133     for (count=0;count<npoints;++count,offset+=size)
1134     {
1135 jfenwick 2461 res[count]=pointToTuple2(getDataPointShape(), vec, offset);
1136 jfenwick 2459 }
1137     }
1138     else if (rank==3)
1139     {
1140     size_t count;
1141     size_t offset=0;
1142     for (count=0;count<npoints;++count,offset+=size)
1143     {
1144 jfenwick 2461 res[count]=pointToTuple3(getDataPointShape(), vec, offset);
1145 jfenwick 2459 }
1146     }
1147     else if (rank==4)
1148     {
1149     size_t count;
1150     size_t offset=0;
1151     for (count=0;count<npoints;++count,offset+=size)
1152     {
1153 jfenwick 2461 res[count]=pointToTuple4(getDataPointShape(), vec, offset);
1154 jfenwick 2459 }
1155     }
1156     else
1157     {
1158     throw DataException("Unknown rank in ::toListOfTuples()");
1159     }
1160     return res;
1161     }
1162    
1163     const boost::python::object
1164 jfenwick 2271 Data::getValueOfDataPointAsTuple(int dataPointNo)
1165     {
1166     forceResolve();
1167     if (getNumDataPointsPerSample()>0) {
1168 gross 921 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1169     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1170     //
1171     // Check a valid sample number has been supplied
1172 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1173 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1174 gross 921 }
1175 ksteube 1312
1176 gross 921 //
1177     // Check a valid data point number has been supplied
1178 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1179 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1180 gross 921 }
1181     // TODO: global error handling
1182 jfenwick 1796 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1183 jfenwick 2271 return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1184     }
1185     else
1186     {
1187 jfenwick 2471 // The pre-numpy method would return an empty array of the given shape
1188 jfenwick 2271 // I'm going to throw an exception because if we have zero points per sample we have problems
1189     throw DataException("Error - need at least 1 datapoint per sample.");
1190     }
1191 ksteube 1312
1192 jgs 117 }
1193 jfenwick 1796
1194 jfenwick 2271
1195 gross 921 void
1196 ksteube 1312 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1197 jgs 121 {
1198 gross 1034 // this will throw if the value cannot be represented
1199 jfenwick 2271 setValueOfDataPointToArray(dataPointNo,py_object);
1200 gross 1034 }
1201    
1202     void
1203 jfenwick 2271 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1204 gross 1034 {
1205 gross 921 if (isProtected()) {
1206     throw DataException("Error - attempt to update protected Data object.");
1207     }
1208 jfenwick 2271
1209     WrappedArray w(obj);
1210 gross 921 //
1211     // check rank
1212 jfenwick 2271 if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1213 jfenwick 2458 throw DataException("Rank of array does not match Data object rank");
1214 bcumming 790
1215 jgs 121 //
1216 jfenwick 2458 // check shape of array
1217 phornby 1953 for (unsigned int i=0; i<getDataPointRank(); i++) {
1218 jfenwick 2271 if (w.getShape()[i]!=getDataPointShape()[i])
1219 jfenwick 2458 throw DataException("Shape of array does not match Data object rank");
1220 jgs 121 }
1221 jfenwick 2721
1222     exclusiveWrite();
1223    
1224 jgs 121 //
1225 gross 921 // make sure data is expanded:
1226 gross 1859 //
1227 gross 921 if (!isExpanded()) {
1228     expand();
1229 jgs 121 }
1230 gross 921 if (getNumDataPointsPerSample()>0) {
1231     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1232     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1233 jfenwick 2271 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1234 gross 921 } else {
1235 jfenwick 2271 m_data->copyToDataPoint(-1, 0,w);
1236 gross 921 }
1237     }
1238 jgs 121
1239 gross 922 void
1240     Data::setValueOfDataPoint(int dataPointNo, const double value)
1241     {
1242     if (isProtected()) {
1243     throw DataException("Error - attempt to update protected Data object.");
1244     }
1245     //
1246     // make sure data is expanded:
1247 jfenwick 2721 exclusiveWrite();
1248 gross 922 if (!isExpanded()) {
1249     expand();
1250     }
1251     if (getNumDataPointsPerSample()>0) {
1252     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1253     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1254     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1255     } else {
1256     m_data->copyToDataPoint(-1, 0,value);
1257     }
1258     }
1259    
1260 ksteube 1312 const
1261 jfenwick 2271 boost::python::object
1262     Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1263     {
1264     // This could be lazier than it is now
1265     forceResolve();
1266 jgs 121
1267 jfenwick 2271 // copy datapoint into a buffer
1268     // broadcast buffer to all nodes
1269     // convert buffer to tuple
1270     // return tuple
1271 jgs 121
1272 jfenwick 2271 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1273     size_t length=DataTypes::noValues(dataPointShape);
1274 jgs 121
1275 gross 921 // added for the MPI communication
1276     double *tmpData = new double[length];
1277 bcumming 790
1278 jfenwick 2271 // updated for the MPI case
1279     if( get_MPIRank()==procNo ){
1280     if (getNumDataPointsPerSample()>0) {
1281     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1282     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1283     //
1284     // Check a valid sample number has been supplied
1285     if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1286 jfenwick 2471 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid sampleNo.");
1287 bcumming 790 }
1288    
1289 jfenwick 2271 //
1290     // Check a valid data point number has been supplied
1291     if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1292 jfenwick 2471 throw DataException("Error - Data::getValueOfGlobalDataPointAsTuple: invalid dataPointNoInSample.");
1293 bcumming 790 }
1294 jfenwick 2271 // TODO: global error handling
1295     DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1296 bcumming 790
1297 jfenwick 2271 memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1298     }
1299     }
1300     #ifdef PASO_MPI
1301     // broadcast the data to all other processes
1302     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1303     #endif
1304    
1305     boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1306     delete [] tmpData;
1307 jgs 121 //
1308     // return the loaded array
1309 jfenwick 2271 return t;
1310    
1311 jgs 121 }
1312    
1313 gross 921
1314 jfenwick 2271 boost::python::object
1315     Data::integrateToTuple_const() const
1316 jfenwick 2005 {
1317     if (isLazy())
1318     {
1319     throw DataException("Error - cannot integrate for constant lazy data.");
1320     }
1321     return integrateWorker();
1322     }
1323 gross 921
1324 jfenwick 2271 boost::python::object
1325     Data::integrateToTuple()
1326 jgs 94 {
1327 jfenwick 2005 if (isLazy())
1328     {
1329 jfenwick 2721 expand(); // Can't do a non-resolving version of this without changing the domain object
1330     } // see the dom->setToIntegrals call. Not saying it can't be done, just not doing it yet.
1331 jfenwick 2005 return integrateWorker();
1332 jfenwick 2271
1333 jfenwick 2005 }
1334    
1335 jfenwick 2271 boost::python::object
1336 jfenwick 2005 Data::integrateWorker() const
1337     {
1338 jfenwick 1796 DataTypes::ShapeType shape = getDataPointShape();
1339 ksteube 1312 int dataPointSize = getDataPointSize();
1340 jgs 94
1341     //
1342     // calculate the integral values
1343 ksteube 1312 vector<double> integrals(dataPointSize);
1344     vector<double> integrals_local(dataPointSize);
1345 jfenwick 2089 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1346     if (dom==0)
1347     {
1348     throw DataException("Can not integrate over non-continuous domains.");
1349     }
1350 ksteube 1312 #ifdef PASO_MPI
1351 jfenwick 2089 dom->setToIntegrals(integrals_local,*this);
1352 ksteube 1312 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1353     double *tmp = new double[dataPointSize];
1354     double *tmp_local = new double[dataPointSize];
1355     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1356     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1357     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1358 jfenwick 2271 tuple result=pointToTuple(shape,tmp);
1359 ksteube 1312 delete[] tmp;
1360     delete[] tmp_local;
1361     #else
1362 jfenwick 2089 dom->setToIntegrals(integrals,*this);
1363 jfenwick 2271 /* double *tmp = new double[dataPointSize];
1364     for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1365     tuple result=pointToTuple(shape,integrals);
1366     // delete tmp;
1367 ksteube 1312 #endif
1368 jgs 94
1369    
1370 jfenwick 2271 return result;
1371 jgs 94 }
1372    
1373     Data
1374     Data::sin() const
1375     {
1376 jfenwick 2146 MAKELAZYOP(SIN)
1377 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1378 jgs 94 }
1379    
1380     Data
1381     Data::cos() const
1382     {
1383 jfenwick 2146 MAKELAZYOP(COS)
1384 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1385 jgs 94 }
1386    
1387     Data
1388     Data::tan() const
1389     {
1390 jfenwick 2146 MAKELAZYOP(TAN)
1391 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1392 jgs 94 }
1393    
1394     Data
1395 jgs 150 Data::asin() const
1396     {
1397 jfenwick 2146 MAKELAZYOP(ASIN)
1398 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1399 jgs 150 }
1400    
1401     Data
1402     Data::acos() const
1403     {
1404 jfenwick 2146 MAKELAZYOP(ACOS)
1405 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1406 jgs 150 }
1407    
1408 phornby 1026
1409 jgs 150 Data
1410     Data::atan() const
1411     {
1412 jfenwick 2146 MAKELAZYOP(ATAN)
1413 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1414 jgs 150 }
1415    
1416     Data
1417     Data::sinh() const
1418     {
1419 jfenwick 2146 MAKELAZYOP(SINH)
1420 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1421 jgs 150 }
1422    
1423     Data
1424     Data::cosh() const
1425     {
1426 jfenwick 2146 MAKELAZYOP(COSH)
1427 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1428 jgs 150 }
1429    
1430     Data
1431     Data::tanh() const
1432     {
1433 jfenwick 2146 MAKELAZYOP(TANH)
1434 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1435 jgs 150 }
1436    
1437 phornby 1026
1438 jgs 150 Data
1439 phornby 1026 Data::erf() const
1440     {
1441 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1442 gross 1028 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1443     #else
1444 jfenwick 2146 MAKELAZYOP(ERF)
1445 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
1446 phornby 1026 #endif
1447     }
1448    
1449     Data
1450 jgs 150 Data::asinh() const
1451     {
1452 jfenwick 2146 MAKELAZYOP(ASINH)
1453 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1454 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1455 phornby 1032 #else
1456 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
1457 phornby 1032 #endif
1458 jgs 150 }
1459    
1460     Data
1461     Data::acosh() const
1462     {
1463 jfenwick 2146 MAKELAZYOP(ACOSH)
1464 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1465 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1466 phornby 1032 #else
1467 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
1468 phornby 1032 #endif
1469 jgs 150 }
1470    
1471     Data
1472     Data::atanh() const
1473     {
1474 jfenwick 2146 MAKELAZYOP(ATANH)
1475 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1476 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1477 phornby 1032 #else
1478 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
1479 phornby 1032 #endif
1480 jgs 150 }
1481    
1482     Data
1483 gross 286 Data::log10() const
1484 jfenwick 2146 {
1485     MAKELAZYOP(LOG10)
1486 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1487 jgs 94 }
1488    
1489     Data
1490 gross 286 Data::log() const
1491 jgs 94 {
1492 jfenwick 2146 MAKELAZYOP(LOG)
1493 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1494 jgs 94 }
1495    
1496 jgs 106 Data
1497     Data::sign() const
1498 jgs 94 {
1499 jfenwick 2146 MAKELAZYOP(SIGN)
1500 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1501 jgs 94 }
1502    
1503 jgs 106 Data
1504     Data::abs() const
1505 jgs 94 {
1506 jfenwick 2146 MAKELAZYOP(ABS)
1507 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1508 jgs 94 }
1509    
1510 jgs 106 Data
1511     Data::neg() const
1512 jgs 94 {
1513 jfenwick 2146 MAKELAZYOP(NEG)
1514 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1515 jgs 94 }
1516    
1517 jgs 102 Data
1518 jgs 106 Data::pos() const
1519 jgs 94 {
1520 jfenwick 2005 // not doing lazy check here is deliberate.
1521     // since a deep copy of lazy data should be cheap, I'll just let it happen now
1522 jgs 148 Data result;
1523     // perform a deep copy
1524     result.copy(*this);
1525     return result;
1526 jgs 102 }
1527    
1528     Data
1529 jgs 106 Data::exp() const
1530 jfenwick 2146 {
1531     MAKELAZYOP(EXP)
1532 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1533 jgs 102 }
1534    
1535     Data
1536 jgs 106 Data::sqrt() const
1537 jgs 102 {
1538 jfenwick 2146 MAKELAZYOP(SQRT)
1539 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1540 jgs 102 }
1541    
1542 jgs 106 double
1543 jfenwick 2005 Data::Lsup_const() const
1544 jgs 102 {
1545 jfenwick 2005 if (isLazy())
1546     {
1547     throw DataException("Error - cannot compute Lsup for constant lazy data.");
1548     }
1549     return LsupWorker();
1550     }
1551    
1552     double
1553     Data::Lsup()
1554     {
1555     if (isLazy())
1556     {
1557 jfenwick 2721 if (!actsExpanded() || CHECK_DO_CRES)
1558     {
1559     resolve();
1560     }
1561     else
1562     {
1563     #ifdef PASO_MPI
1564     return lazyAlgWorker<AbsMax>(0,MPI_MAX);
1565     #else
1566     return lazyAlgWorker<AbsMax>(0,0);
1567     #endif
1568     }
1569 jfenwick 2005 }
1570     return LsupWorker();
1571     }
1572    
1573     double
1574     Data::sup_const() const
1575     {
1576     if (isLazy())
1577     {
1578     throw DataException("Error - cannot compute sup for constant lazy data.");
1579     }
1580     return supWorker();
1581     }
1582    
1583     double
1584     Data::sup()
1585     {
1586     if (isLazy())
1587     {
1588 jfenwick 2721 if (!actsExpanded() || CHECK_DO_CRES)
1589     {
1590     resolve();
1591     }
1592     else
1593     {
1594     #ifdef PASO_MPI
1595     return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, MPI_MAX);
1596     #else
1597     return lazyAlgWorker<FMax>(numeric_limits<double>::max()*-1, 0);
1598     #endif
1599     }
1600 jfenwick 2005 }
1601     return supWorker();
1602     }
1603    
1604     double
1605     Data::inf_const() const
1606     {
1607     if (isLazy())
1608     {
1609     throw DataException("Error - cannot compute inf for constant lazy data.");
1610     }
1611     return infWorker();
1612     }
1613    
1614     double
1615     Data::inf()
1616     {
1617     if (isLazy())
1618     {
1619 jfenwick 2721 if (!actsExpanded() || CHECK_DO_CRES)
1620     {
1621     resolve();
1622     }
1623     else
1624     {
1625     #ifdef PASO_MPI
1626     return lazyAlgWorker<FMin>(numeric_limits<double>::max(), MPI_MIN);
1627     #else
1628     return lazyAlgWorker<FMin>(numeric_limits<double>::max(), 0);
1629     #endif
1630     }
1631 jfenwick 2005 }
1632     return infWorker();
1633     }
1634    
1635 jfenwick 2721 template <class BinaryOp>
1636 jfenwick 2005 double
1637 jfenwick 2721 Data::lazyAlgWorker(double init, int mpiop_type)
1638     {
1639     if (!isLazy() || !m_data->actsExpanded())
1640     {
1641     throw DataException("Error - lazyAlgWorker can only be called on lazy(expanded) data.");
1642     }
1643     DataLazy* dl=dynamic_cast<DataLazy*>(m_data.get());
1644     EsysAssert((dl!=0), "Programming error - casting to DataLazy.");
1645     BufferGroup* bg=allocSampleBuffer();
1646     double val=init;
1647     int i=0;
1648     const size_t numsamples=getNumSamples();
1649     const size_t samplesize=getNoValues()*getNumDataPointsPerSample();
1650     BinaryOp operation;
1651     #pragma omp parallel private(i)
1652     {
1653     double localtot=init;
1654     #pragma omp for schedule(static) private(i)
1655     for (i=0;i<numsamples;++i)
1656     {
1657     size_t roffset=0;
1658     const DataTypes::ValueType* v=dl->resolveSample(*bg, i, roffset);
1659     // Now we have the sample, run operation on all points
1660     for (size_t j=0;j<samplesize;++j)
1661     {
1662     localtot=operation(localtot,(*v)[j+roffset]);
1663     }
1664     }
1665     #pragma omp critical
1666     val=operation(val,localtot);
1667     }
1668     freeSampleBuffer(bg);
1669     #ifdef PASO_MPI
1670     double globalValue;
1671     MPI_Allreduce( &val, &globalValue, 1, MPI_DOUBLE, mpiop_type, MPI_COMM_WORLD );
1672     return globalValue;
1673     #else
1674     return val;
1675     #endif
1676     }
1677    
1678     double
1679 jfenwick 2005 Data::LsupWorker() const
1680     {
1681 phornby 1455 double localValue;
1682 jgs 106 //
1683     // set the initial absolute maximum value to zero
1684 bcumming 751
1685 jgs 147 AbsMax abs_max_func;
1686 bcumming 751 localValue = algorithm(abs_max_func,0);
1687     #ifdef PASO_MPI
1688 phornby 1455 double globalValue;
1689 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1690     return globalValue;
1691     #else
1692     return localValue;
1693     #endif
1694 jgs 117 }
1695    
1696     double
1697 jfenwick 2005 Data::supWorker() const
1698 jgs 102 {
1699 phornby 1455 double localValue;
1700 jgs 106 //
1701     // set the initial maximum value to min possible double
1702 jgs 147 FMax fmax_func;
1703 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1704     #ifdef PASO_MPI
1705 phornby 1455 double globalValue;
1706 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1707     return globalValue;
1708     #else
1709     return localValue;
1710     #endif
1711 jgs 102 }
1712    
1713 jgs 106 double
1714 jfenwick 2005 Data::infWorker() const
1715 jgs 102 {
1716 phornby 1455 double localValue;
1717 jgs 106 //
1718     // set the initial minimum value to max possible double
1719 jgs 147 FMin fmin_func;
1720 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1721     #ifdef PASO_MPI
1722 phornby 1455 double globalValue;
1723 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1724     return globalValue;
1725     #else
1726     return localValue;
1727     #endif
1728 jgs 102 }
1729    
1730 bcumming 751 /* TODO */
1731     /* global reduction */
1732 jgs 102 Data
1733 jgs 106 Data::maxval() const
1734 jgs 102 {
1735 jfenwick 2721 MAKELAZYOP(MAXVAL)
1736 jgs 113 //
1737     // set the initial maximum value to min possible double
1738 jgs 147 FMax fmax_func;
1739     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1740 jgs 102 }
1741    
1742     Data
1743 jgs 106 Data::minval() const
1744 jgs 102 {
1745 jfenwick 2721 MAKELAZYOP(MINVAL)
1746 jgs 113 //
1747     // set the initial minimum value to max possible double
1748 jgs 147 FMin fmin_func;
1749     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1750 jgs 102 }
1751    
1752 jgs 123 Data
1753 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1754 jgs 123 {
1755 gross 804 int axis0_tmp,axis1_tmp;
1756 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1757     DataTypes::ShapeType ev_shape;
1758 gross 800 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1759     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1760     int rank=getDataPointRank();
1761 gross 804 if (rank<2) {
1762     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1763 gross 800 }
1764 gross 804 if (axis0<0 || axis0>rank-1) {
1765     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1766     }
1767     if (axis1<0 || axis1>rank-1) {
1768     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1769     }
1770     if (axis0 == axis1) {
1771     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1772     }
1773 jfenwick 2496 MAKELAZYOP2(SWAP,axis0,axis1)
1774     if (axis0 > axis1)
1775     {
1776     axis0_tmp=axis1;
1777     axis1_tmp=axis0;
1778 gross 804 }
1779 jfenwick 2496 else
1780     {
1781     axis0_tmp=axis0;
1782     axis1_tmp=axis1;
1783 gross 800 }
1784 jfenwick 2496 for (int i=0; i<rank; i++)
1785     {
1786     if (i == axis0_tmp)
1787     {
1788     ev_shape.push_back(s[axis1_tmp]);
1789     }
1790     else if (i == axis1_tmp)
1791     {
1792     ev_shape.push_back(s[axis0_tmp]);
1793     }
1794     else
1795     {
1796     ev_shape.push_back(s[i]);
1797     }
1798     }
1799 gross 800 Data ev(0.,ev_shape,getFunctionSpace());
1800     ev.typeMatchRight(*this);
1801 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1802 gross 800 return ev;
1803 jgs 123 }
1804    
1805     Data
1806 ksteube 775 Data::symmetric() const
1807 jgs 123 {
1808 ksteube 775 // check input
1809 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1810 ksteube 775 if (getDataPointRank()==2) {
1811 ksteube 1312 if(s[0] != s[1])
1812 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1813     }
1814     else if (getDataPointRank()==4) {
1815     if(!(s[0] == s[2] && s[1] == s[3]))
1816     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1817     }
1818     else {
1819     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1820     }
1821 jfenwick 2146 MAKELAZYOP(SYM)
1822 ksteube 775 Data ev(0.,getDataPointShape(),getFunctionSpace());
1823     ev.typeMatchRight(*this);
1824     m_data->symmetric(ev.m_data.get());
1825     return ev;
1826     }
1827    
1828     Data
1829     Data::nonsymmetric() const
1830     {
1831 jfenwick 2146 MAKELAZYOP(NSYM)
1832 ksteube 775 // check input
1833 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1834 ksteube 775 if (getDataPointRank()==2) {
1835 ksteube 1312 if(s[0] != s[1])
1836 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1837 jfenwick 1796 DataTypes::ShapeType ev_shape;
1838 ksteube 775 ev_shape.push_back(s[0]);
1839     ev_shape.push_back(s[1]);
1840     Data ev(0.,ev_shape,getFunctionSpace());
1841     ev.typeMatchRight(*this);
1842     m_data->nonsymmetric(ev.m_data.get());
1843     return ev;
1844     }
1845     else if (getDataPointRank()==4) {
1846     if(!(s[0] == s[2] && s[1] == s[3]))
1847     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1848 jfenwick 1796 DataTypes::ShapeType ev_shape;
1849 ksteube 775 ev_shape.push_back(s[0]);
1850     ev_shape.push_back(s[1]);
1851     ev_shape.push_back(s[2]);
1852     ev_shape.push_back(s[3]);
1853     Data ev(0.,ev_shape,getFunctionSpace());
1854     ev.typeMatchRight(*this);
1855     m_data->nonsymmetric(ev.m_data.get());
1856     return ev;
1857     }
1858     else {
1859     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1860     }
1861     }
1862    
1863     Data
1864 gross 800 Data::trace(int axis_offset) const
1865 jfenwick 2084 {
1866 jfenwick 2146 MAKELAZYOPOFF(TRACE,axis_offset)
1867 jfenwick 2200 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1868     {
1869     throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
1870     }
1871 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1872 ksteube 775 if (getDataPointRank()==2) {
1873 jfenwick 1796 DataTypes::ShapeType ev_shape;
1874 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1875     ev.typeMatchRight(*this);
1876 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1877 ksteube 775 return ev;
1878     }
1879     if (getDataPointRank()==3) {
1880 jfenwick 1796 DataTypes::ShapeType ev_shape;
1881 ksteube 775 if (axis_offset==0) {
1882     int s2=s[2];
1883     ev_shape.push_back(s2);
1884     }
1885     else if (axis_offset==1) {
1886     int s0=s[0];
1887     ev_shape.push_back(s0);
1888     }
1889     Data ev(0.,ev_shape,getFunctionSpace());
1890     ev.typeMatchRight(*this);
1891 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1892 ksteube 775 return ev;
1893     }
1894     if (getDataPointRank()==4) {
1895 jfenwick 1796 DataTypes::ShapeType ev_shape;
1896 ksteube 775 if (axis_offset==0) {
1897     ev_shape.push_back(s[2]);
1898     ev_shape.push_back(s[3]);
1899     }
1900     else if (axis_offset==1) {
1901     ev_shape.push_back(s[0]);
1902     ev_shape.push_back(s[3]);
1903     }
1904     else if (axis_offset==2) {
1905     ev_shape.push_back(s[0]);
1906     ev_shape.push_back(s[1]);
1907     }
1908     Data ev(0.,ev_shape,getFunctionSpace());
1909     ev.typeMatchRight(*this);
1910 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1911 ksteube 775 return ev;
1912     }
1913     else {
1914 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1915 ksteube 775 }
1916     }
1917    
1918     Data
1919     Data::transpose(int axis_offset) const
1920 jfenwick 2005 {
1921 jfenwick 2146 MAKELAZYOPOFF(TRANS,axis_offset)
1922 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1923     DataTypes::ShapeType ev_shape;
1924 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1925     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1926     int rank=getDataPointRank();
1927     if (axis_offset<0 || axis_offset>rank) {
1928     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1929     }
1930     for (int i=0; i<rank; i++) {
1931 jfenwick 2721
1932 ksteube 775 int index = (axis_offset+i)%rank;
1933     ev_shape.push_back(s[index]); // Append to new shape
1934     }
1935     Data ev(0.,ev_shape,getFunctionSpace());
1936     ev.typeMatchRight(*this);
1937     m_data->transpose(ev.m_data.get(), axis_offset);
1938     return ev;
1939 jgs 123 }
1940    
1941 gross 576 Data
1942     Data::eigenvalues() const
1943     {
1944 jfenwick 2005 if (isLazy())
1945     {
1946     Data temp(*this); // to get around the fact that you can't resolve a const Data
1947     temp.resolve();
1948     return temp.eigenvalues();
1949     }
1950 gross 576 // check input
1951 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1952 ksteube 1312 if (getDataPointRank()!=2)
1953 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1954 ksteube 1312 if(s[0] != s[1])
1955 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1956     // create return
1957 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1958 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1959     ev.typeMatchRight(*this);
1960     m_data->eigenvalues(ev.m_data.get());
1961     return ev;
1962     }
1963    
1964 jgs 121 const boost::python::tuple
1965 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1966     {
1967 jfenwick 2005 if (isLazy())
1968     {
1969     Data temp(*this); // to get around the fact that you can't resolve a const Data
1970     temp.resolve();
1971     return temp.eigenvalues_and_eigenvectors(tol);
1972     }
1973 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1974 ksteube 1312 if (getDataPointRank()!=2)
1975 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1976 ksteube 1312 if(s[0] != s[1])
1977 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1978     // create return
1979 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1980 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1981     ev.typeMatchRight(*this);
1982 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1983 gross 576 Data V(0.,V_shape,getFunctionSpace());
1984     V.typeMatchRight(*this);
1985     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1986     return make_tuple(boost::python::object(ev),boost::python::object(V));
1987     }
1988    
1989     const boost::python::tuple
1990 gross 921 Data::minGlobalDataPoint() const
1991 jgs 121 {
1992 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1993 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1994     // surrounding function
1995    
1996     int DataPointNo;
1997 gross 921 int ProcNo;
1998     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1999     return make_tuple(ProcNo,DataPointNo);
2000 jgs 148 }
2001    
2002     void
2003 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
2004     int& DataPointNo) const
2005 jgs 148 {
2006 jfenwick 2005 if (isLazy())
2007     {
2008     Data temp(*this); // to get around the fact that you can't resolve a const Data
2009     temp.resolve();
2010     return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
2011     }
2012 jgs 148 int i,j;
2013     int lowi=0,lowj=0;
2014     double min=numeric_limits<double>::max();
2015    
2016 jgs 121 Data temp=minval();
2017    
2018     int numSamples=temp.getNumSamples();
2019     int numDPPSample=temp.getNumDataPointsPerSample();
2020    
2021 jgs 148 double next,local_min;
2022 jfenwick 1946 int local_lowi=0,local_lowj=0;
2023 jgs 121
2024 caltinay 2081 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
2025 jgs 148 {
2026     local_min=min;
2027     #pragma omp for private(i,j) schedule(static)
2028     for (i=0; i<numSamples; i++) {
2029     for (j=0; j<numDPPSample; j++) {
2030 jfenwick 2271 next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2031 jgs 148 if (next<local_min) {
2032     local_min=next;
2033     local_lowi=i;
2034     local_lowj=j;
2035     }
2036 jgs 121 }
2037     }
2038 jgs 148 #pragma omp critical
2039 jfenwick 2476 if (local_min<min) { // If we found a smaller value than our sentinel
2040 jgs 148 min=local_min;
2041     lowi=local_lowi;
2042     lowj=local_lowj;
2043     }
2044 jgs 121 }
2045    
2046 bcumming 782 #ifdef PASO_MPI
2047 jfenwick 2476 // determine the processor on which the minimum occurs
2048     next = temp.getDataPointRO(lowi,lowj);
2049     int lowProc = 0;
2050     double *globalMins = new double[get_MPISize()+1];
2051     int error;
2052     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
2053 ksteube 1312
2054 jfenwick 2476 if( get_MPIRank()==0 ){
2055     next = globalMins[lowProc];
2056     for( i=1; i<get_MPISize(); i++ )
2057     if( next>globalMins[i] ){
2058     lowProc = i;
2059     next = globalMins[i];
2060     }
2061     }
2062 gross 2492 MPI_Bcast( &lowProc, 1, MPI_INT, 0, get_MPIComm() );
2063 jfenwick 2653 DataPointNo = lowj + lowi * numDPPSample;
2064     MPI_Bcast(&DataPointNo, 1, MPI_INT, lowProc, get_MPIComm() );
2065 jfenwick 2476 delete [] globalMins;
2066     ProcNo = lowProc;
2067 bcumming 790 #else
2068 jfenwick 2476 ProcNo = 0;
2069 jfenwick 2653 DataPointNo = lowj + lowi * numDPPSample;
2070 bcumming 782 #endif
2071 jgs 121 }
2072    
2073 jfenwick 2476
2074     const boost::python::tuple
2075     Data::maxGlobalDataPoint() const
2076     {
2077     int DataPointNo;
2078     int ProcNo;
2079     calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2080     return make_tuple(ProcNo,DataPointNo);
2081     }
2082    
2083 jgs 104 void
2084 jfenwick 2476 Data::calc_maxGlobalDataPoint(int& ProcNo,
2085     int& DataPointNo) const
2086     {
2087     if (isLazy())
2088     {
2089     Data temp(*this); // to get around the fact that you can't resolve a const Data
2090     temp.resolve();
2091     return temp.calc_maxGlobalDataPoint(ProcNo,DataPointNo);
2092     }
2093     int i,j;
2094     int highi=0,highj=0;
2095     //-------------
2096 jfenwick 2653 double max= -numeric_limits<double>::max();
2097 jfenwick 2476
2098     Data temp=maxval();
2099    
2100     int numSamples=temp.getNumSamples();
2101     int numDPPSample=temp.getNumDataPointsPerSample();
2102    
2103     double next,local_max;
2104     int local_highi=0,local_highj=0;
2105    
2106     #pragma omp parallel firstprivate(local_highi,local_highj) private(next,local_max)
2107     {
2108     local_max=max;
2109     #pragma omp for private(i,j) schedule(static)
2110     for (i=0; i<numSamples; i++) {
2111     for (j=0; j<numDPPSample; j++) {
2112     next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
2113     if (next>local_max) {
2114     local_max=next;
2115     local_highi=i;
2116     local_highj=j;
2117     }
2118     }
2119     }
2120     #pragma omp critical
2121     if (local_max>max) { // If we found a larger value than our sentinel
2122     max=local_max;
2123     highi=local_highi;
2124     highj=local_highj;
2125     }
2126     }
2127     #ifdef PASO_MPI
2128     // determine the processor on which the maximum occurs
2129     next = temp.getDataPointRO(highi,highj);
2130     int highProc = 0;
2131     double *globalMaxs = new double[get_MPISize()+1];
2132     int error;
2133     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMaxs, 1, MPI_DOUBLE, 0, get_MPIComm() );
2134     if( get_MPIRank()==0 ){
2135 jfenwick 2653 next = globalMaxs[highProc];
2136     for( i=1; i<get_MPISize(); i++ )
2137     {
2138     if( next<globalMaxs[i] )
2139     {
2140 jfenwick 2476 highProc = i;
2141     next = globalMaxs[i];
2142     }
2143 jfenwick 2653 }
2144 jfenwick 2476 }
2145 gross 2492 MPI_Bcast( &highProc, 1, MPI_INT, 0, get_MPIComm() );
2146 jfenwick 2653 DataPointNo = highj + highi * numDPPSample;
2147     MPI_Bcast(&DataPointNo, 1, MPI_INT, highProc, get_MPIComm() );
2148    
2149 jfenwick 2476 delete [] globalMaxs;
2150     ProcNo = highProc;
2151     #else
2152     ProcNo = 0;
2153 jfenwick 2653 DataPointNo = highj + highi * numDPPSample;
2154 jfenwick 2476 #endif
2155     }
2156    
2157     void
2158 jgs 104 Data::saveDX(std::string fileName) const
2159     {
2160 jfenwick 1803 if (isEmpty())
2161     {
2162     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2163     }
2164 jfenwick 2005 if (isLazy())
2165     {
2166     Data temp(*this); // to get around the fact that you can't resolve a const Data
2167     temp.resolve();
2168     temp.saveDX(fileName);
2169     return;
2170     }
2171 jgs 153 boost::python::dict args;
2172     args["data"]=boost::python::object(this);
2173 jfenwick 1872 getDomain()->saveDX(fileName,args);
2174 jgs 104 return;
2175     }
2176    
2177 jgs 110 void
2178     Data::saveVTK(std::string fileName) const
2179     {
2180 jfenwick 1803 if (isEmpty())
2181     {
2182     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2183     }
2184 jfenwick 2005 if (isLazy())
2185     {
2186     Data temp(*this); // to get around the fact that you can't resolve a const Data
2187     temp.resolve();
2188     temp.saveVTK(fileName);
2189     return;
2190     }
2191 jgs 153 boost::python::dict args;
2192     args["data"]=boost::python::object(this);
2193 gross 2421 getDomain()->saveVTK(fileName,args,"","");
2194 jgs 110 return;
2195     }
2196    
2197 jfenwick 2635
2198    
2199 jgs 102 Data&
2200     Data::operator+=(const Data& right)
2201     {
2202 gross 783 if (isProtected()) {
2203     throw DataException("Error - attempt to update protected Data object.");
2204     }
2205 jfenwick 2146 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
2206 jfenwick 2271 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
2207 jfenwick 2146 binaryOp(right,plus<double>());
2208     return (*this);
2209 jgs 94 }
2210    
2211 jgs 102 Data&
2212     Data::operator+=(const boost::python::object& right)
2213 jgs 94 {
2214 jfenwick 2005 if (isProtected()) {
2215     throw DataException("Error - attempt to update protected Data object.");
2216     }
2217 gross 854 Data tmp(right,getFunctionSpace(),false);
2218 jfenwick 2271 (*this)+=tmp;
2219     return *this;
2220 jgs 94 }
2221 jfenwick 2005
2222     // Hmmm, operator= makes a deep copy but the copy constructor does not?
2223 ksteube 1312 Data&
2224     Data::operator=(const Data& other)
2225     {
2226 jfenwick 2271 m_protected=false; // since any changes should be caught by exclusiveWrite();
2227     // m_data=other.m_data;
2228     set_m_data(other.m_data);
2229 ksteube 1312 return (*this);
2230     }
2231 jgs 94
2232 jgs 102 Data&
2233     Data::operator-=(const Data& right)
2234 jgs 94 {
2235 gross 783 if (isProtected()) {
2236     throw DataException("Error - attempt to update protected Data object.");
2237     }
2238 jfenwick 2146 MAKELAZYBINSELF(right,SUB)
2239 jfenwick 2271 exclusiveWrite();
2240 jfenwick 2146 binaryOp(right,minus<double>());
2241     return (*this);
2242 jgs 94 }
2243    
2244 jgs 102 Data&
2245     Data::operator-=(const boost::python::object& right)
2246 jgs 94 {
2247 jfenwick 2005 if (isProtected()) {
2248     throw DataException("Error - attempt to update protected Data object.");
2249     }
2250 gross 854 Data tmp(right,getFunctionSpace(),false);
2251 jfenwick 2271 (*this)-=tmp;
2252 jfenwick 2146 return (*this);
2253 jgs 94 }
2254    
2255 jgs 102 Data&
2256     Data::operator*=(const Data& right)
2257 jgs 94 {
2258 gross 783 if (isProtected()) {
2259     throw DataException("Error - attempt to update protected Data object.");
2260     }
2261 jfenwick 2146 MAKELAZYBINSELF(right,MUL)
2262 jfenwick 2271 exclusiveWrite();
2263 jfenwick 2146 binaryOp(right,multiplies<double>());
2264     return (*this);
2265 jgs 94 }
2266    
2267 jgs 102 Data&
2268     Data::operator*=(const boost::python::object& right)
2269 jfenwick 2005 {
2270     if (isProtected()) {
2271     throw DataException("Error - attempt to update protected Data object.");
2272     }
2273 gross 854 Data tmp(right,getFunctionSpace(),false);
2274 jfenwick 2271 (*this)*=tmp;
2275 jfenwick 2146 return (*this);
2276 jgs 94 }
2277    
2278 jgs 102 Data&
2279     Data::operator/=(const Data& right)
2280 jgs 94 {
2281 gross 783 if (isProtected()) {
2282     throw DataException("Error - attempt to update protected Data object.");
2283     }
2284 jfenwick 2146 MAKELAZYBINSELF(right,DIV)
2285 jfenwick 2271 exclusiveWrite();
2286 jfenwick 2146 binaryOp(right,divides<double>());
2287     return (*this);
2288 jgs 94 }
2289    
2290 jgs 102 Data&
2291     Data::operator/=(const boost::python::object& right)
2292 jgs 94 {
2293 jfenwick 2005 if (isProtected()) {
2294     throw DataException("Error - attempt to update protected Data object.");
2295     }
2296 gross 854 Data tmp(right,getFunctionSpace(),false);
2297 jfenwick 2271 (*this)/=tmp;
2298 jfenwick 2146 return (*this);
2299 jgs 94 }
2300    
2301 jgs 102 Data
2302 gross 699 Data::rpowO(const boost::python::object& left) const
2303     {
2304     Data left_d(left,*this);
2305     return left_d.powD(*this);
2306     }
2307    
2308     Data
2309 jgs 102 Data::powO(const boost::python::object& right) const
2310 jgs 94 {
2311 gross 854 Data tmp(right,getFunctionSpace(),false);
2312     return powD(tmp);
2313 jgs 94 }
2314    
2315 jgs 102 Data
2316     Data::powD(const Data& right) const
2317 jgs 94 {
2318 jfenwick 2146 MAKELAZYBIN(right,POW)
2319 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2320 jgs 94 }
2321    
2322     //
2323 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2324 jgs 102 Data
2325     escript::operator+(const Data& left, const Data& right)
2326 jgs 94 {
2327 jfenwick 2146 MAKELAZYBIN2(left,right,ADD)
2328 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
2329 jgs 94 }
2330    
2331     //
2332 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2333 jgs 102 Data
2334     escript::operator-(const Data& left, const Data& right)
2335 jgs 94 {
2336 jfenwick 2146 MAKELAZYBIN2(left,right,SUB)
2337 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2338 jgs 94 }
2339    
2340     //
2341 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2342 jgs 102 Data
2343     escript::operator*(const Data& left, const Data& right)
2344 jgs 94 {
2345 jfenwick 2146 MAKELAZYBIN2(left,right,MUL)
2346 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2347 jgs 94 }
2348    
2349     //
2350 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2351 jgs 102 Data
2352     escript::operator/(const Data& left, const Data& right)
2353 jgs 94 {
2354 jfenwick 2146 MAKELAZYBIN2(left,right,DIV)
2355 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2356 jgs 94 }
2357    
2358     //
2359 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2360 jgs 102 Data
2361     escript::operator+(const Data& left, const boost::python::object& right)
2362 jgs 94 {
2363 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2364     MAKELAZYBIN2(left,tmp,ADD)
2365     return left+tmp;
2366 jgs 94 }
2367    
2368     //
2369 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2370 jgs 102 Data
2371     escript::operator-(const Data& left, const boost::python::object& right)
2372 jgs 94 {
2373 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2374     MAKELAZYBIN2(left,tmp,SUB)
2375     return left-tmp;
2376 jgs 94 }
2377    
2378     //
2379 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2380 jgs 102 Data
2381     escript::operator*(const Data& left, const boost::python::object& right)
2382 jgs 94 {
2383 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2384     MAKELAZYBIN2(left,tmp,MUL)
2385     return left*tmp;
2386 jgs 94 }
2387    
2388     //
2389 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2390 jgs 102 Data
2391     escript::operator/(const Data& left, const boost::python::object& right)
2392 jgs 94 {
2393 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2394     MAKELAZYBIN2(left,tmp,DIV)
2395     return left/tmp;
2396 jgs 94 }
2397    
2398     //
2399 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2400 jgs 102 Data
2401     escript::operator+(const boost::python::object& left, const Data& right)
2402 jgs 94 {
2403 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2404     MAKELAZYBIN2(tmp,right,ADD)
2405     return tmp+right;
2406 jgs 94 }
2407    
2408     //
2409 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2410 jgs 102 Data
2411     escript::operator-(const boost::python::object& left, const Data& right)
2412 jgs 94 {
2413 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2414     MAKELAZYBIN2(tmp,right,SUB)
2415     return tmp-right;
2416 jgs 94 }
2417    
2418     //
2419 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2420 jgs 102 Data
2421     escript::operator*(const boost::python::object& left, const Data& right)
2422 jgs 94 {
2423 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2424     MAKELAZYBIN2(tmp,right,MUL)
2425     return tmp*right;
2426 jgs 94 }
2427    
2428     //
2429 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2430 jgs 102 Data
2431     escript::operator/(const boost::python::object& left, const Data& right)
2432 jgs 94 {
2433 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2434     MAKELAZYBIN2(tmp,right,DIV)
2435     return tmp/right;
2436 jgs 94 }
2437    
2438 jgs 102
2439 bcumming 751 /* TODO */
2440     /* global reduction */
2441 jgs 102 Data
2442 ksteube 1312 Data::getItem(const boost::python::object& key) const
2443 jgs 94 {
2444    
2445 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2446 jgs 94
2447 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2448 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2449 jgs 94 }
2450    
2451 jgs 102 return getSlice(slice_region);
2452 jgs 94 }
2453    
2454 bcumming 751 /* TODO */
2455     /* global reduction */
2456 jgs 94 Data
2457 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2458 jgs 94 {
2459 jgs 102 return Data(*this,region);
2460 jgs 94 }
2461    
2462 bcumming 751 /* TODO */
2463     /* global reduction */
2464 jgs 94 void
2465 jgs 102 Data::setItemO(const boost::python::object& key,
2466     const boost::python::object& value)
2467 jgs 94 {
2468 jgs 102 Data tempData(value,getFunctionSpace());
2469     setItemD(key,tempData);
2470