/[escript]/branches/numpy/escript/src/Data.cpp
ViewVC logotype

Annotation of /branches/numpy/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2453 - (hide annotations)
Tue Jun 2 22:54:20 2009 UTC (9 years, 9 months ago) by jfenwick
File size: 80688 byte(s)
Removing NONUMARRAY defines

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