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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2282 - (hide annotations)
Thu Feb 19 06:30:19 2009 UTC (11 years, 1 month ago) by jfenwick
File size: 81407 byte(s)
setTaggedValueByName now throws if the tag name does not exist.
The documentation says this is what happens - now the code matches the doco.
This resolves issue 238
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 ksteube 1312 const
1019 jfenwick 2271 boost::python::numeric :: array
1020 ksteube 1312 Data:: getValueOfDataPoint(int dataPointNo)
1021 jgs 121 {
1022 jfenwick 2271 return boost::python::numeric::array(getValueOfDataPointAsTuple(dataPointNo));
1023     }
1024 jfenwick 2005
1025 jfenwick 2271 const boost::python::object
1026     Data::getValueOfDataPointAsTuple(int dataPointNo)
1027     {
1028     forceResolve();
1029     if (getNumDataPointsPerSample()>0) {
1030 gross 921 int sampleNo = dataPointNo/getNumDataPointsPerSample();
1031     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1032     //
1033     // Check a valid sample number has been supplied
1034 trankine 924 if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1035 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid sampleNo.");
1036 gross 921 }
1037 ksteube 1312
1038 gross 921 //
1039     // Check a valid data point number has been supplied
1040 trankine 924 if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1041 jfenwick 2271 throw DataException("Error - Data::getValueOfDataPointAsTuple: invalid dataPointNoInSample.");
1042 gross 921 }
1043     // TODO: global error handling
1044 jfenwick 1796 DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1045 jfenwick 2271 return pointToTuple(getDataPointShape(),&(getDataAtOffsetRO(offset)));
1046     }
1047     else
1048     {
1049     // The old method would return an empty numArray of the given shape
1050     // I'm going to throw an exception because if we have zero points per sample we have problems
1051     throw DataException("Error - need at least 1 datapoint per sample.");
1052     }
1053 ksteube 1312
1054 jgs 117 }
1055 jfenwick 1796
1056 jfenwick 2271
1057 gross 921 void
1058 ksteube 1312 Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
1059 jgs 121 {
1060 gross 1034 // this will throw if the value cannot be represented
1061 jfenwick 2271 setValueOfDataPointToArray(dataPointNo,py_object);
1062 gross 1034 }
1063    
1064     void
1065 jfenwick 2271 Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::object& obj)
1066 gross 1034 {
1067 gross 921 if (isProtected()) {
1068     throw DataException("Error - attempt to update protected Data object.");
1069     }
1070 jfenwick 2271 forceResolve();
1071    
1072     WrappedArray w(obj);
1073 gross 921 //
1074     // check rank
1075 jfenwick 2271 if (static_cast<unsigned int>(w.getRank())<getDataPointRank())
1076 gross 921 throw DataException("Rank of numarray does not match Data object rank");
1077 bcumming 790
1078 jgs 121 //
1079 gross 921 // check shape of num_array
1080 phornby 1953 for (unsigned int i=0; i<getDataPointRank(); i++) {
1081 jfenwick 2271 if (w.getShape()[i]!=getDataPointShape()[i])
1082 gross 921 throw DataException("Shape of numarray does not match Data object rank");
1083 jgs 121 }
1084     //
1085 gross 921 // make sure data is expanded:
1086 gross 1859 //
1087 gross 921 if (!isExpanded()) {
1088     expand();
1089 jgs 121 }
1090 gross 921 if (getNumDataPointsPerSample()>0) {
1091     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1092     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1093 jfenwick 2271 m_data->copyToDataPoint(sampleNo, dataPointNoInSample,w);
1094 gross 921 } else {
1095 jfenwick 2271 m_data->copyToDataPoint(-1, 0,w);
1096 gross 921 }
1097     }
1098 jgs 121
1099 gross 922 void
1100     Data::setValueOfDataPoint(int dataPointNo, const double value)
1101     {
1102     if (isProtected()) {
1103     throw DataException("Error - attempt to update protected Data object.");
1104     }
1105     //
1106     // make sure data is expanded:
1107 jfenwick 2271 forceResolve();
1108 gross 922 if (!isExpanded()) {
1109     expand();
1110     }
1111     if (getNumDataPointsPerSample()>0) {
1112     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1113     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1114     m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
1115     } else {
1116     m_data->copyToDataPoint(-1, 0,value);
1117     }
1118     }
1119    
1120 ksteube 1312 const
1121 gross 921 boost::python::numeric::array
1122 ksteube 1312 Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
1123 gross 921 {
1124 jfenwick 2271 return boost::python::numeric::array(getValueOfGlobalDataPointAsTuple(procNo, dataPointNo));
1125     }
1126 jgs 121
1127 jfenwick 2271 const
1128     boost::python::object
1129     Data::getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo)
1130     {
1131     // This could be lazier than it is now
1132     forceResolve();
1133 jgs 121
1134 jfenwick 2271 // copy datapoint into a buffer
1135     // broadcast buffer to all nodes
1136     // convert buffer to tuple
1137     // return tuple
1138 jgs 121
1139 jfenwick 2271 const DataTypes::ShapeType& dataPointShape = getDataPointShape();
1140     size_t length=DataTypes::noValues(dataPointShape);
1141 jgs 121
1142 gross 921 // added for the MPI communication
1143     double *tmpData = new double[length];
1144 bcumming 790
1145 jfenwick 2271 // updated for the MPI case
1146     if( get_MPIRank()==procNo ){
1147     if (getNumDataPointsPerSample()>0) {
1148     int sampleNo = dataPointNo/getNumDataPointsPerSample();
1149     int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
1150     //
1151     // Check a valid sample number has been supplied
1152     if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
1153     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
1154 bcumming 790 }
1155    
1156 jfenwick 2271 //
1157     // Check a valid data point number has been supplied
1158     if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
1159     throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
1160 bcumming 790 }
1161 jfenwick 2271 // TODO: global error handling
1162     DataTypes::ValueType::size_type offset=getDataOffset(sampleNo, dataPointNoInSample);
1163 bcumming 790
1164 jfenwick 2271 memcpy(tmpData,&(getDataAtOffsetRO(offset)),length*sizeof(double));
1165     }
1166     }
1167     #ifdef PASO_MPI
1168     // broadcast the data to all other processes
1169     MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
1170     #endif
1171    
1172     boost::python::tuple t=pointToTuple(dataPointShape,tmpData);
1173     delete [] tmpData;
1174 jgs 121 //
1175     // return the loaded array
1176 jfenwick 2271 return t;
1177    
1178 jgs 121 }
1179    
1180 gross 921
1181 jfenwick 2271 boost::python::object
1182     Data::integrateToTuple_const() const
1183 jfenwick 2005 {
1184     if (isLazy())
1185     {
1186     throw DataException("Error - cannot integrate for constant lazy data.");
1187     }
1188     return integrateWorker();
1189     }
1190 gross 921
1191 jfenwick 2271 boost::python::object
1192     Data::integrateToTuple()
1193 jgs 94 {
1194 jfenwick 2005 if (isLazy())
1195     {
1196     expand();
1197     }
1198     return integrateWorker();
1199 jfenwick 2271
1200 jfenwick 2005 }
1201    
1202    
1203 jfenwick 2271 boost::python::object
1204     Data::integrate_const() const
1205     {
1206     if (isLazy())
1207     {
1208     throw DataException("Error - cannot integrate for constant lazy data.");
1209     }
1210     return boost::python::numeric::array(integrateWorker());
1211     }
1212 jfenwick 2005
1213 jfenwick 2271 boost::python::object
1214     Data::integrate()
1215     {
1216     if (isLazy())
1217     {
1218     expand();
1219     }
1220     return boost::python::numeric::array(integrateWorker());
1221     }
1222    
1223    
1224    
1225     boost::python::object
1226 jfenwick 2005 Data::integrateWorker() const
1227     {
1228 jfenwick 1796 DataTypes::ShapeType shape = getDataPointShape();
1229 ksteube 1312 int dataPointSize = getDataPointSize();
1230 jgs 94
1231     //
1232     // calculate the integral values
1233 ksteube 1312 vector<double> integrals(dataPointSize);
1234     vector<double> integrals_local(dataPointSize);
1235 jfenwick 2089 const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1236     if (dom==0)
1237     {
1238     throw DataException("Can not integrate over non-continuous domains.");
1239     }
1240 ksteube 1312 #ifdef PASO_MPI
1241 jfenwick 2089 dom->setToIntegrals(integrals_local,*this);
1242 ksteube 1312 // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1243     double *tmp = new double[dataPointSize];
1244     double *tmp_local = new double[dataPointSize];
1245     for (int i=0; i<dataPointSize; i++) { tmp_local[i] = integrals_local[i]; }
1246     MPI_Allreduce( &tmp_local[0], &tmp[0], dataPointSize, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
1247     for (int i=0; i<dataPointSize; i++) { integrals[i] = tmp[i]; }
1248 jfenwick 2271 tuple result=pointToTuple(shape,tmp);
1249 ksteube 1312 delete[] tmp;
1250     delete[] tmp_local;
1251     #else
1252 jfenwick 2089 dom->setToIntegrals(integrals,*this);
1253 jfenwick 2271 /* double *tmp = new double[dataPointSize];
1254     for (int i=0; i<dataPointSize; i++) { tmp[i]=integrals[i]; }*/
1255     tuple result=pointToTuple(shape,integrals);
1256     // delete tmp;
1257 ksteube 1312 #endif
1258 jgs 94
1259    
1260 jfenwick 2271 return result;
1261 jgs 94 }
1262    
1263     Data
1264     Data::sin() const
1265     {
1266 jfenwick 2146 MAKELAZYOP(SIN)
1267 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sin);
1268 jgs 94 }
1269    
1270     Data
1271     Data::cos() const
1272     {
1273 jfenwick 2146 MAKELAZYOP(COS)
1274 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cos);
1275 jgs 94 }
1276    
1277     Data
1278     Data::tan() const
1279     {
1280 jfenwick 2146 MAKELAZYOP(TAN)
1281 matt 1349 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tan);
1282 jgs 94 }
1283    
1284     Data
1285 jgs 150 Data::asin() const
1286     {
1287 jfenwick 2146 MAKELAZYOP(ASIN)
1288 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::asin);
1289 jgs 150 }
1290    
1291     Data
1292     Data::acos() const
1293     {
1294 jfenwick 2146 MAKELAZYOP(ACOS)
1295 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::acos);
1296 jgs 150 }
1297    
1298 phornby 1026
1299 jgs 150 Data
1300     Data::atan() const
1301     {
1302 jfenwick 2146 MAKELAZYOP(ATAN)
1303 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::atan);
1304 jgs 150 }
1305    
1306     Data
1307     Data::sinh() const
1308     {
1309 jfenwick 2146 MAKELAZYOP(SINH)
1310 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sinh);
1311 jgs 150 }
1312    
1313     Data
1314     Data::cosh() const
1315     {
1316 jfenwick 2146 MAKELAZYOP(COSH)
1317 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::cosh);
1318 jgs 150 }
1319    
1320     Data
1321     Data::tanh() const
1322     {
1323 jfenwick 2146 MAKELAZYOP(TANH)
1324 jfenwick 2005 return C_TensorUnaryOperation<double (*)(double)>(*this, ::tanh);
1325 jgs 150 }
1326    
1327 phornby 1026
1328 jgs 150 Data
1329 phornby 1026 Data::erf() const
1330     {
1331 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1332 gross 1028 throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1333     #else
1334 jfenwick 2146 MAKELAZYOP(ERF)
1335 matt 1334 return C_TensorUnaryOperation(*this, ::erf);
1336 phornby 1026 #endif
1337     }
1338    
1339     Data
1340 jgs 150 Data::asinh() const
1341     {
1342 jfenwick 2146 MAKELAZYOP(ASINH)
1343 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1344 matt 1334 return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1345 phornby 1032 #else
1346 matt 1334 return C_TensorUnaryOperation(*this, ::asinh);
1347 phornby 1032 #endif
1348 jgs 150 }
1349    
1350     Data
1351     Data::acosh() const
1352     {
1353 jfenwick 2146 MAKELAZYOP(ACOSH)
1354 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1355 matt 1334 return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1356 phornby 1032 #else
1357 matt 1334 return C_TensorUnaryOperation(*this, ::acosh);
1358 phornby 1032 #endif
1359 jgs 150 }
1360    
1361     Data
1362     Data::atanh() const
1363     {
1364 jfenwick 2146 MAKELAZYOP(ATANH)
1365 phornby 2048 #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1366 matt 1334 return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1367 phornby 1032 #else
1368 matt 1334 return C_TensorUnaryOperation(*this, ::atanh);
1369 phornby 1032 #endif
1370 jgs 150 }
1371    
1372     Data
1373 gross 286 Data::log10() const
1374 jfenwick 2146 {
1375     MAKELAZYOP(LOG10)
1376 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log10);
1377 jgs 94 }
1378    
1379     Data
1380 gross 286 Data::log() const
1381 jgs 94 {
1382 jfenwick 2146 MAKELAZYOP(LOG)
1383 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::log);
1384 jgs 94 }
1385    
1386 jgs 106 Data
1387     Data::sign() const
1388 jgs 94 {
1389 jfenwick 2146 MAKELAZYOP(SIGN)
1390 matt 1334 return C_TensorUnaryOperation(*this, escript::fsign);
1391 jgs 94 }
1392    
1393 jgs 106 Data
1394     Data::abs() const
1395 jgs 94 {
1396 jfenwick 2146 MAKELAZYOP(ABS)
1397 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::fabs);
1398 jgs 94 }
1399    
1400 jgs 106 Data
1401     Data::neg() const
1402 jgs 94 {
1403 jfenwick 2146 MAKELAZYOP(NEG)
1404 matt 1334 return C_TensorUnaryOperation(*this, negate<double>());
1405 jgs 94 }
1406    
1407 jgs 102 Data
1408 jgs 106 Data::pos() const
1409 jgs 94 {
1410 jfenwick 2005 // not doing lazy check here is deliberate.
1411     // since a deep copy of lazy data should be cheap, I'll just let it happen now
1412 jgs 148 Data result;
1413     // perform a deep copy
1414     result.copy(*this);
1415     return result;
1416 jgs 102 }
1417    
1418     Data
1419 jgs 106 Data::exp() const
1420 jfenwick 2146 {
1421     MAKELAZYOP(EXP)
1422 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::exp);
1423 jgs 102 }
1424    
1425     Data
1426 jgs 106 Data::sqrt() const
1427 jgs 102 {
1428 jfenwick 2146 MAKELAZYOP(SQRT)
1429 matt 1350 return C_TensorUnaryOperation<double (*)(double)>(*this, ::sqrt);
1430 jgs 102 }
1431    
1432 jgs 106 double
1433 jfenwick 2005 Data::Lsup_const() const
1434 jgs 102 {
1435 jfenwick 2005 if (isLazy())
1436     {
1437     throw DataException("Error - cannot compute Lsup for constant lazy data.");
1438     }
1439     return LsupWorker();
1440     }
1441    
1442     double
1443     Data::Lsup()
1444     {
1445     if (isLazy())
1446     {
1447 jfenwick 2085 resolve();
1448 jfenwick 2005 }
1449     return LsupWorker();
1450     }
1451    
1452     double
1453     Data::sup_const() const
1454     {
1455     if (isLazy())
1456     {
1457     throw DataException("Error - cannot compute sup for constant lazy data.");
1458     }
1459     return supWorker();
1460     }
1461    
1462     double
1463     Data::sup()
1464     {
1465     if (isLazy())
1466     {
1467 jfenwick 2085 resolve();
1468 jfenwick 2005 }
1469     return supWorker();
1470     }
1471    
1472     double
1473     Data::inf_const() const
1474     {
1475     if (isLazy())
1476     {
1477     throw DataException("Error - cannot compute inf for constant lazy data.");
1478     }
1479     return infWorker();
1480     }
1481    
1482     double
1483     Data::inf()
1484     {
1485     if (isLazy())
1486     {
1487 jfenwick 2085 resolve();
1488 jfenwick 2005 }
1489     return infWorker();
1490     }
1491    
1492     double
1493     Data::LsupWorker() const
1494     {
1495 phornby 1455 double localValue;
1496 jgs 106 //
1497     // set the initial absolute maximum value to zero
1498 bcumming 751
1499 jgs 147 AbsMax abs_max_func;
1500 bcumming 751 localValue = algorithm(abs_max_func,0);
1501     #ifdef PASO_MPI
1502 phornby 1455 double globalValue;
1503 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1504     return globalValue;
1505     #else
1506     return localValue;
1507     #endif
1508 jgs 117 }
1509    
1510     double
1511 jfenwick 2005 Data::supWorker() const
1512 jgs 102 {
1513 phornby 1455 double localValue;
1514 jgs 106 //
1515     // set the initial maximum value to min possible double
1516 jgs 147 FMax fmax_func;
1517 bcumming 751 localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1518     #ifdef PASO_MPI
1519 phornby 1455 double globalValue;
1520 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1521     return globalValue;
1522     #else
1523     return localValue;
1524     #endif
1525 jgs 102 }
1526    
1527 jgs 106 double
1528 jfenwick 2005 Data::infWorker() const
1529 jgs 102 {
1530 phornby 1455 double localValue;
1531 jgs 106 //
1532     // set the initial minimum value to max possible double
1533 jgs 147 FMin fmin_func;
1534 bcumming 751 localValue = algorithm(fmin_func,numeric_limits<double>::max());
1535     #ifdef PASO_MPI
1536 phornby 1455 double globalValue;
1537 bcumming 751 MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1538     return globalValue;
1539     #else
1540     return localValue;
1541     #endif
1542 jgs 102 }
1543    
1544 bcumming 751 /* TODO */
1545     /* global reduction */
1546 jgs 102 Data
1547 jgs 106 Data::maxval() const
1548 jgs 102 {
1549 jfenwick 2037 if (isLazy())
1550     {
1551     Data temp(*this); // to get around the fact that you can't resolve a const Data
1552     temp.resolve();
1553     return temp.maxval();
1554     }
1555 jgs 113 //
1556     // set the initial maximum value to min possible double
1557 jgs 147 FMax fmax_func;
1558     return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1559 jgs 102 }
1560    
1561     Data
1562 jgs 106 Data::minval() const
1563 jgs 102 {
1564 jfenwick 2037 if (isLazy())
1565     {
1566     Data temp(*this); // to get around the fact that you can't resolve a const Data
1567     temp.resolve();
1568     return temp.minval();
1569     }
1570 jgs 113 //
1571     // set the initial minimum value to max possible double
1572 jgs 147 FMin fmin_func;
1573     return dp_algorithm(fmin_func,numeric_limits<double>::max());
1574 jgs 102 }
1575    
1576 jgs 123 Data
1577 gross 804 Data::swapaxes(const int axis0, const int axis1) const
1578 jgs 123 {
1579 jfenwick 2085 if (isLazy())
1580     {
1581     Data temp(*this);
1582     temp.resolve();
1583     return temp.swapaxes(axis0,axis1);
1584     }
1585 gross 804 int axis0_tmp,axis1_tmp;
1586 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1587     DataTypes::ShapeType ev_shape;
1588 gross 800 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1589     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1590     int rank=getDataPointRank();
1591 gross 804 if (rank<2) {
1592     throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1593 gross 800 }
1594 gross 804 if (axis0<0 || axis0>rank-1) {
1595     throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1596     }
1597     if (axis1<0 || axis1>rank-1) {
1598     throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1599     }
1600     if (axis0 == axis1) {
1601     throw DataException("Error - Data::swapaxes: axis indices must be different.");
1602     }
1603     if (axis0 > axis1) {
1604     axis0_tmp=axis1;
1605     axis1_tmp=axis0;
1606     } else {
1607     axis0_tmp=axis0;
1608     axis1_tmp=axis1;
1609     }
1610 gross 800 for (int i=0; i<rank; i++) {
1611 gross 804 if (i == axis0_tmp) {
1612 ksteube 1312 ev_shape.push_back(s[axis1_tmp]);
1613 gross 804 } else if (i == axis1_tmp) {
1614 ksteube 1312 ev_shape.push_back(s[axis0_tmp]);
1615 gross 800 } else {
1616 ksteube 1312 ev_shape.push_back(s[i]);
1617 gross 800 }
1618     }
1619     Data ev(0.,ev_shape,getFunctionSpace());
1620     ev.typeMatchRight(*this);
1621 gross 804 m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1622 gross 800 return ev;
1623    
1624 jgs 123 }
1625    
1626     Data
1627 ksteube 775 Data::symmetric() const
1628 jgs 123 {
1629 ksteube 775 // check input
1630 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1631 ksteube 775 if (getDataPointRank()==2) {
1632 ksteube 1312 if(s[0] != s[1])
1633 ksteube 775 throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1634     }
1635     else if (getDataPointRank()==4) {
1636     if(!(s[0] == s[2] && s[1] == s[3]))
1637     throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1638     }
1639     else {
1640     throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1641     }
1642 jfenwick 2146 MAKELAZYOP(SYM)
1643 ksteube 775 Data ev(0.,getDataPointShape(),getFunctionSpace());
1644     ev.typeMatchRight(*this);
1645     m_data->symmetric(ev.m_data.get());
1646     return ev;
1647     }
1648    
1649     Data
1650     Data::nonsymmetric() const
1651     {
1652 jfenwick 2146 MAKELAZYOP(NSYM)
1653 ksteube 775 // check input
1654 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1655 ksteube 775 if (getDataPointRank()==2) {
1656 ksteube 1312 if(s[0] != s[1])
1657 ksteube 775 throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1658 jfenwick 1796 DataTypes::ShapeType ev_shape;
1659 ksteube 775 ev_shape.push_back(s[0]);
1660     ev_shape.push_back(s[1]);
1661     Data ev(0.,ev_shape,getFunctionSpace());
1662     ev.typeMatchRight(*this);
1663     m_data->nonsymmetric(ev.m_data.get());
1664     return ev;
1665     }
1666     else if (getDataPointRank()==4) {
1667     if(!(s[0] == s[2] && s[1] == s[3]))
1668     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1669 jfenwick 1796 DataTypes::ShapeType ev_shape;
1670 ksteube 775 ev_shape.push_back(s[0]);
1671     ev_shape.push_back(s[1]);
1672     ev_shape.push_back(s[2]);
1673     ev_shape.push_back(s[3]);
1674     Data ev(0.,ev_shape,getFunctionSpace());
1675     ev.typeMatchRight(*this);
1676     m_data->nonsymmetric(ev.m_data.get());
1677     return ev;
1678     }
1679     else {
1680     throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1681     }
1682     }
1683    
1684     Data
1685 gross 800 Data::trace(int axis_offset) const
1686 jfenwick 2084 {
1687 jfenwick 2146 MAKELAZYOPOFF(TRACE,axis_offset)
1688 jfenwick 2200 if ((axis_offset<0) || (axis_offset>getDataPointRank()))
1689     {
1690     throw DataException("Error - Data::trace, axis_offset must be between 0 and rank-2 inclusive.");
1691     }
1692 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1693 ksteube 775 if (getDataPointRank()==2) {
1694 jfenwick 1796 DataTypes::ShapeType ev_shape;
1695 ksteube 775 Data ev(0.,ev_shape,getFunctionSpace());
1696     ev.typeMatchRight(*this);
1697 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1698 ksteube 775 return ev;
1699     }
1700     if (getDataPointRank()==3) {
1701 jfenwick 1796 DataTypes::ShapeType ev_shape;
1702 ksteube 775 if (axis_offset==0) {
1703     int s2=s[2];
1704     ev_shape.push_back(s2);
1705     }
1706     else if (axis_offset==1) {
1707     int s0=s[0];
1708     ev_shape.push_back(s0);
1709     }
1710     Data ev(0.,ev_shape,getFunctionSpace());
1711     ev.typeMatchRight(*this);
1712 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1713 ksteube 775 return ev;
1714     }
1715     if (getDataPointRank()==4) {
1716 jfenwick 1796 DataTypes::ShapeType ev_shape;
1717 ksteube 775 if (axis_offset==0) {
1718     ev_shape.push_back(s[2]);
1719     ev_shape.push_back(s[3]);
1720     }
1721     else if (axis_offset==1) {
1722     ev_shape.push_back(s[0]);
1723     ev_shape.push_back(s[3]);
1724     }
1725     else if (axis_offset==2) {
1726     ev_shape.push_back(s[0]);
1727     ev_shape.push_back(s[1]);
1728     }
1729     Data ev(0.,ev_shape,getFunctionSpace());
1730     ev.typeMatchRight(*this);
1731 gross 800 m_data->trace(ev.m_data.get(), axis_offset);
1732 ksteube 775 return ev;
1733     }
1734     else {
1735 gross 800 throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1736 ksteube 775 }
1737     }
1738    
1739     Data
1740     Data::transpose(int axis_offset) const
1741 jfenwick 2005 {
1742 jfenwick 2146 MAKELAZYOPOFF(TRANS,axis_offset)
1743 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1744     DataTypes::ShapeType ev_shape;
1745 ksteube 775 // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1746     // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1747     int rank=getDataPointRank();
1748     if (axis_offset<0 || axis_offset>rank) {
1749     throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1750     }
1751     for (int i=0; i<rank; i++) {
1752     int index = (axis_offset+i)%rank;
1753     ev_shape.push_back(s[index]); // Append to new shape
1754     }
1755     Data ev(0.,ev_shape,getFunctionSpace());
1756     ev.typeMatchRight(*this);
1757     m_data->transpose(ev.m_data.get(), axis_offset);
1758     return ev;
1759 jgs 123 }
1760    
1761 gross 576 Data
1762     Data::eigenvalues() const
1763     {
1764 jfenwick 2005 if (isLazy())
1765     {
1766     Data temp(*this); // to get around the fact that you can't resolve a const Data
1767     temp.resolve();
1768     return temp.eigenvalues();
1769     }
1770 gross 576 // check input
1771 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1772 ksteube 1312 if (getDataPointRank()!=2)
1773 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1774 ksteube 1312 if(s[0] != s[1])
1775 gross 576 throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1776     // create return
1777 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1778 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1779     ev.typeMatchRight(*this);
1780     m_data->eigenvalues(ev.m_data.get());
1781     return ev;
1782     }
1783    
1784 jgs 121 const boost::python::tuple
1785 gross 576 Data::eigenvalues_and_eigenvectors(const double tol) const
1786     {
1787 jfenwick 2005 if (isLazy())
1788     {
1789     Data temp(*this); // to get around the fact that you can't resolve a const Data
1790     temp.resolve();
1791     return temp.eigenvalues_and_eigenvectors(tol);
1792     }
1793 jfenwick 1796 DataTypes::ShapeType s=getDataPointShape();
1794 ksteube 1312 if (getDataPointRank()!=2)
1795 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1796 ksteube 1312 if(s[0] != s[1])
1797 gross 576 throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1798     // create return
1799 jfenwick 1796 DataTypes::ShapeType ev_shape(1,s[0]);
1800 gross 576 Data ev(0.,ev_shape,getFunctionSpace());
1801     ev.typeMatchRight(*this);
1802 jfenwick 1796 DataTypes::ShapeType V_shape(2,s[0]);
1803 gross 576 Data V(0.,V_shape,getFunctionSpace());
1804     V.typeMatchRight(*this);
1805     m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1806     return make_tuple(boost::python::object(ev),boost::python::object(V));
1807     }
1808    
1809     const boost::python::tuple
1810 gross 921 Data::minGlobalDataPoint() const
1811 jgs 121 {
1812 gross 921 // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1813 jgs 148 // abort (for unknown reasons) if there are openmp directives with it in the
1814     // surrounding function
1815    
1816     int DataPointNo;
1817 gross 921 int ProcNo;
1818     calc_minGlobalDataPoint(ProcNo,DataPointNo);
1819     return make_tuple(ProcNo,DataPointNo);
1820 jgs 148 }
1821    
1822     void
1823 gross 921 Data::calc_minGlobalDataPoint(int& ProcNo,
1824     int& DataPointNo) const
1825 jgs 148 {
1826 jfenwick 2005 if (isLazy())
1827     {
1828     Data temp(*this); // to get around the fact that you can't resolve a const Data
1829     temp.resolve();
1830     return temp.calc_minGlobalDataPoint(ProcNo,DataPointNo);
1831     }
1832 jgs 148 int i,j;
1833     int lowi=0,lowj=0;
1834     double min=numeric_limits<double>::max();
1835    
1836 jgs 121 Data temp=minval();
1837    
1838     int numSamples=temp.getNumSamples();
1839     int numDPPSample=temp.getNumDataPointsPerSample();
1840    
1841 jgs 148 double next,local_min;
1842 jfenwick 1946 int local_lowi=0,local_lowj=0;
1843 jgs 121
1844 caltinay 2081 #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1845 jgs 148 {
1846     local_min=min;
1847     #pragma omp for private(i,j) schedule(static)
1848     for (i=0; i<numSamples; i++) {
1849     for (j=0; j<numDPPSample; j++) {
1850 jfenwick 2271 next=temp.getDataAtOffsetRO(temp.getDataOffset(i,j));
1851 jgs 148 if (next<local_min) {
1852     local_min=next;
1853     local_lowi=i;
1854     local_lowj=j;
1855     }
1856 jgs 121 }
1857     }
1858 jgs 148 #pragma omp critical
1859     if (local_min<min) {
1860     min=local_min;
1861     lowi=local_lowi;
1862     lowj=local_lowj;
1863     }
1864 jgs 121 }
1865    
1866 bcumming 782 #ifdef PASO_MPI
1867     // determine the processor on which the minimum occurs
1868 jfenwick 2271 next = temp.getDataPointRO(lowi,lowj);
1869 bcumming 782 int lowProc = 0;
1870     double *globalMins = new double[get_MPISize()+1];
1871 caltinay 2081 int error;
1872     error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1873 ksteube 1312
1874 bcumming 782 if( get_MPIRank()==0 ){
1875     next = globalMins[lowProc];
1876     for( i=1; i<get_MPISize(); i++ )
1877     if( next>globalMins[i] ){
1878     lowProc = i;
1879     next = globalMins[i];
1880     }
1881     }
1882     MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1883    
1884     delete [] globalMins;
1885     ProcNo = lowProc;
1886 bcumming 790 #else
1887     ProcNo = 0;
1888 bcumming 782 #endif
1889 gross 921 DataPointNo = lowj + lowi * numDPPSample;
1890 jgs 121 }
1891    
1892 jgs 104 void
1893     Data::saveDX(std::string fileName) const
1894     {
1895 jfenwick 1803 if (isEmpty())
1896     {
1897     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1898     }
1899 jfenwick 2005 if (isLazy())
1900     {
1901     Data temp(*this); // to get around the fact that you can't resolve a const Data
1902     temp.resolve();
1903     temp.saveDX(fileName);
1904     return;
1905     }
1906 jgs 153 boost::python::dict args;
1907     args["data"]=boost::python::object(this);
1908 jfenwick 1872 getDomain()->saveDX(fileName,args);
1909 jgs 104 return;
1910     }
1911    
1912 jgs 110 void
1913     Data::saveVTK(std::string fileName) const
1914     {
1915 jfenwick 1803 if (isEmpty())
1916     {
1917     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
1918     }
1919 jfenwick 2005 if (isLazy())
1920     {
1921     Data temp(*this); // to get around the fact that you can't resolve a const Data
1922     temp.resolve();
1923     temp.saveVTK(fileName);
1924     return;
1925     }
1926 jgs 153 boost::python::dict args;
1927     args["data"]=boost::python::object(this);
1928 jfenwick 1872 getDomain()->saveVTK(fileName,args);
1929 jgs 110 return;
1930     }
1931    
1932 jgs 102 Data&
1933     Data::operator+=(const Data& right)
1934     {
1935 gross 783 if (isProtected()) {
1936     throw DataException("Error - attempt to update protected Data object.");
1937     }
1938 jfenwick 2146 MAKELAZYBINSELF(right,ADD) // for lazy + is equivalent to +=
1939 jfenwick 2271 exclusiveWrite(); // Since Lazy data does not modify its leaves we only need to worry here
1940 jfenwick 2146 binaryOp(right,plus<double>());
1941     return (*this);
1942 jgs 94 }
1943    
1944 jgs 102 Data&
1945     Data::operator+=(const boost::python::object& right)
1946 jgs 94 {
1947 jfenwick 2005 if (isProtected()) {
1948     throw DataException("Error - attempt to update protected Data object.");
1949     }
1950 gross 854 Data tmp(right,getFunctionSpace(),false);
1951 jfenwick 2271 (*this)+=tmp;
1952     return *this;
1953 jgs 94 }
1954 jfenwick 2005
1955     // Hmmm, operator= makes a deep copy but the copy constructor does not?
1956 ksteube 1312 Data&
1957     Data::operator=(const Data& other)
1958     {
1959 jfenwick 2271 #if defined ASSIGNMENT_MEANS_DEEPCOPY
1960     // This should not be used.
1961     // Just leaving this here until I have completed transition
1962 ksteube 1312 copy(other);
1963 jfenwick 2271 #else
1964     m_protected=false; // since any changes should be caught by exclusiveWrite();
1965     // m_data=other.m_data;
1966     set_m_data(other.m_data);
1967     #endif
1968 ksteube 1312 return (*this);
1969     }
1970 jgs 94
1971 jgs 102 Data&
1972     Data::operator-=(const Data& right)
1973 jgs 94 {
1974 gross 783 if (isProtected()) {
1975     throw DataException("Error - attempt to update protected Data object.");
1976     }
1977 jfenwick 2146 MAKELAZYBINSELF(right,SUB)
1978 jfenwick 2271 exclusiveWrite();
1979 jfenwick 2146 binaryOp(right,minus<double>());
1980     return (*this);
1981 jgs 94 }
1982    
1983 jgs 102 Data&
1984     Data::operator-=(const boost::python::object& right)
1985 jgs 94 {
1986 jfenwick 2005 if (isProtected()) {
1987     throw DataException("Error - attempt to update protected Data object.");
1988     }
1989 gross 854 Data tmp(right,getFunctionSpace(),false);
1990 jfenwick 2271 (*this)-=tmp;
1991 jfenwick 2146 return (*this);
1992 jgs 94 }
1993    
1994 jgs 102 Data&
1995     Data::operator*=(const Data& right)
1996 jgs 94 {
1997 gross 783 if (isProtected()) {
1998     throw DataException("Error - attempt to update protected Data object.");
1999     }
2000 jfenwick 2146 MAKELAZYBINSELF(right,MUL)
2001 jfenwick 2271 exclusiveWrite();
2002 jfenwick 2146 binaryOp(right,multiplies<double>());
2003     return (*this);
2004 jgs 94 }
2005    
2006 jgs 102 Data&
2007     Data::operator*=(const boost::python::object& right)
2008 jfenwick 2005 {
2009     if (isProtected()) {
2010     throw DataException("Error - attempt to update protected Data object.");
2011     }
2012 gross 854 Data tmp(right,getFunctionSpace(),false);
2013 jfenwick 2271 (*this)*=tmp;
2014 jfenwick 2146 return (*this);
2015 jgs 94 }
2016    
2017 jgs 102 Data&
2018     Data::operator/=(const Data& right)
2019 jgs 94 {
2020 gross 783 if (isProtected()) {
2021     throw DataException("Error - attempt to update protected Data object.");
2022     }
2023 jfenwick 2146 MAKELAZYBINSELF(right,DIV)
2024 jfenwick 2271 exclusiveWrite();
2025 jfenwick 2146 binaryOp(right,divides<double>());
2026     return (*this);
2027 jgs 94 }
2028    
2029 jgs 102 Data&
2030     Data::operator/=(const boost::python::object& right)
2031 jgs 94 {
2032 jfenwick 2005 if (isProtected()) {
2033     throw DataException("Error - attempt to update protected Data object.");
2034     }
2035 gross 854 Data tmp(right,getFunctionSpace(),false);
2036 jfenwick 2271 (*this)/=tmp;
2037 jfenwick 2146 return (*this);
2038 jgs 94 }
2039    
2040 jgs 102 Data
2041 gross 699 Data::rpowO(const boost::python::object& left) const
2042     {
2043     Data left_d(left,*this);
2044     return left_d.powD(*this);
2045     }
2046    
2047     Data
2048 jgs 102 Data::powO(const boost::python::object& right) const
2049 jgs 94 {
2050 gross 854 Data tmp(right,getFunctionSpace(),false);
2051     return powD(tmp);
2052 jgs 94 }
2053    
2054 jgs 102 Data
2055     Data::powD(const Data& right) const
2056 jgs 94 {
2057 jfenwick 2146 MAKELAZYBIN(right,POW)
2058 matt 1350 return C_TensorBinaryOperation<double (*)(double, double)>(*this, right, ::pow);
2059 jgs 94 }
2060    
2061     //
2062 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2063 jgs 102 Data
2064     escript::operator+(const Data& left, const Data& right)
2065 jgs 94 {
2066 jfenwick 2146 MAKELAZYBIN2(left,right,ADD)
2067 matt 1327 return C_TensorBinaryOperation(left, right, plus<double>());
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 Data& right)
2074 jgs 94 {
2075 jfenwick 2146 MAKELAZYBIN2(left,right,SUB)
2076 matt 1327 return C_TensorBinaryOperation(left, right, minus<double>());
2077 jgs 94 }
2078    
2079     //
2080 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2081 jgs 102 Data
2082     escript::operator*(const Data& left, const Data& right)
2083 jgs 94 {
2084 jfenwick 2146 MAKELAZYBIN2(left,right,MUL)
2085 matt 1327 return C_TensorBinaryOperation(left, right, multiplies<double>());
2086 jgs 94 }
2087    
2088     //
2089 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2090 jgs 102 Data
2091     escript::operator/(const Data& left, const Data& right)
2092 jgs 94 {
2093 jfenwick 2146 MAKELAZYBIN2(left,right,DIV)
2094 matt 1327 return C_TensorBinaryOperation(left, right, divides<double>());
2095 jgs 94 }
2096    
2097     //
2098 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2099 jgs 102 Data
2100     escript::operator+(const Data& left, const boost::python::object& right)
2101 jgs 94 {
2102 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2103     MAKELAZYBIN2(left,tmp,ADD)
2104     return left+tmp;
2105 jgs 94 }
2106    
2107     //
2108 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2109 jgs 102 Data
2110     escript::operator-(const Data& left, const boost::python::object& right)
2111 jgs 94 {
2112 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2113     MAKELAZYBIN2(left,tmp,SUB)
2114     return left-tmp;
2115 jgs 94 }
2116    
2117     //
2118 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2119 jgs 102 Data
2120     escript::operator*(const Data& left, const boost::python::object& right)
2121 jgs 94 {
2122 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2123     MAKELAZYBIN2(left,tmp,MUL)
2124     return left*tmp;
2125 jgs 94 }
2126    
2127     //
2128 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2129 jgs 102 Data
2130     escript::operator/(const Data& left, const boost::python::object& right)
2131 jgs 94 {
2132 jfenwick 2146 Data tmp(right,left.getFunctionSpace(),false);
2133     MAKELAZYBIN2(left,tmp,DIV)
2134     return left/tmp;
2135 jgs 94 }
2136    
2137     //
2138 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2139 jgs 102 Data
2140     escript::operator+(const boost::python::object& left, const Data& right)
2141 jgs 94 {
2142 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2143     MAKELAZYBIN2(tmp,right,ADD)
2144     return tmp+right;
2145 jgs 94 }
2146    
2147     //
2148 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2149 jgs 102 Data
2150     escript::operator-(const boost::python::object& left, const Data& right)
2151 jgs 94 {
2152 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2153     MAKELAZYBIN2(tmp,right,SUB)
2154     return tmp-right;
2155 jgs 94 }
2156    
2157     //
2158 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2159 jgs 102 Data
2160     escript::operator*(const boost::python::object& left, const Data& right)
2161 jgs 94 {
2162 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2163     MAKELAZYBIN2(tmp,right,MUL)
2164     return tmp*right;
2165 jgs 94 }
2166    
2167     //
2168 jgs 123 // NOTE: It is essential to specify the namespace this operator belongs to
2169 jgs 102 Data
2170     escript::operator/(const boost::python::object& left, const Data& right)
2171 jgs 94 {
2172 jfenwick 2146 Data tmp(left,right.getFunctionSpace(),false);
2173     MAKELAZYBIN2(tmp,right,DIV)
2174     return tmp/right;
2175 jgs 94 }
2176    
2177 jgs 102
2178 bcumming 751 /* TODO */
2179     /* global reduction */
2180 jgs 102 Data
2181 ksteube 1312 Data::getItem(const boost::python::object& key) const
2182 jgs 94 {
2183    
2184 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2185 jgs 94
2186 jfenwick 1796 if (slice_region.size()!=getDataPointRank()) {
2187 jgs 102 throw DataException("Error - slice size does not match Data rank.");
2188 jgs 94 }
2189    
2190 jgs 102 return getSlice(slice_region);
2191 jgs 94 }
2192    
2193 bcumming 751 /* TODO */
2194     /* global reduction */
2195 jgs 94 Data
2196 jfenwick 1796 Data::getSlice(const DataTypes::RegionType& region) const
2197 jgs 94 {
2198 jgs 102 return Data(*this,region);
2199 jgs 94 }
2200    
2201 bcumming 751 /* TODO */
2202     /* global reduction */
2203 jgs 94 void
2204 jgs 102 Data::setItemO(const boost::python::object& key,
2205     const boost::python::object& value)
2206 jgs 94 {
2207 jgs 102 Data tempData(value,getFunctionSpace());
2208     setItemD(key,tempData);
2209     }
2210    
2211     void
2212     Data::setItemD(const boost::python::object& key,
2213     const Data& value)
2214     {
2215 jfenwick 1796 DataTypes::RegionType slice_region=DataTypes::getSliceRegion(getDataPointShape(),key);
2216     if (slice_region.size()!=getDataPointRank()) {
2217 jgs 94 throw DataException("Error - slice size does not match Data rank.");
2218     }
2219 jfenwick 2271 exclusiveWrite();
2220 jgs 108 if (getFunctionSpace()!=value.getFunctionSpace()) {
2221     setSlice(Data(value,getFunctionSpace()),slice_region);
2222     } else {
2223     setSlice(value,slice_region);
2224     }
2225 jgs 94 }
2226    
2227     void
2228 jgs 102 Data::setSlice(const Data& value,
2229 jfenwick 1796 const DataTypes::RegionType& region)
2230 jgs 94 {
2231 gross 783 if (isProtected()) {
2232     throw DataException("Error - attempt to update protected Data object.");
2233     }
2234 jfenwick 2271 forceResolve();
2235     exclusiveWrite(); // In case someone finds a way to call this without going through setItemD
2236 jgs 102 Data tempValue(value);
2237     typeMatchLeft(tempValue);
2238     typeMatchRight(tempValue);
2239 jfenwick 2005 getReady()->setSlice(tempValue.m_data.get(),region);
2240 jgs 102 }
2241    
2242     void
2243     Data::typeMatchLeft(Data& right) const
2244     {
2245 jfenwick 2005 if (right.isLazy() && !isLazy())
2246     {
2247     right.resolve();
2248     }
2249 jgs 102 if (isExpanded()){
2250     right.expand();
2251     } else if (isTagged()) {
2252     if (right.isConstant()) {
2253     right.tag();
2254     }
2255     }
2256     }
2257    
2258     void
2259     Data::typeMatchRight(const Data& right)
2260     {
2261 jfenwick 2005 if (isLazy() && !right.isLazy())
2262     {
2263     resolve();
2264     }
2265 jgs 94 if (isTagged()) {
2266     if (right.isExpanded()) {
2267     expand();
2268     }
2269     } else if (isConstant()) {
2270     if (right.isExpanded()) {
2271     expand();
2272     } else if (right.isTagged()) {
2273     tag();
2274     }
2275     }
2276     }
2277    
2278 jfenwick 2282 // The normal TaggedValue adds the tag if it is not already present
2279     // This form does not. It throws instead.
2280     // This is because the names are maintained by the domain and cannot be added
2281     // without knowing the tag number to map it to.
2282 jgs 94 void
2283 gross 1044 Data::setTaggedValueByName(std::string name,
2284 ksteube 1312 const boost::python::object& value)
2285 gross 1044 {
2286 jfenwick 1872 if (getFunctionSpace().getDomain()->isValidTagName(name)) {
2287 jfenwick 2271 forceResolve();
2288     exclusiveWrite();
2289 jfenwick 1872 int tagKey=getFunctionSpace().getDomain()->getTag(name);
2290 gross 1044 setTaggedValue(tagKey,value);
2291     }
2292 jfenwick 2282 else
2293     { // The
2294     throw DataException("Error - unknown tag in setTaggedValueByName.");
2295     }
2296 gross 1044 }
2297 jfenwick 2271
2298 gross 1044 void
2299 jgs 94 Data::setTaggedValue(int tagKey,
2300     const boost::python::object& value)
2301     {
2302 gross 783 if (isProtected()) {
2303     throw DataException("Error - attempt to update protected Data object.");
2304     }
2305 jgs 94 //
2306     // Ensure underlying data object is of type DataTagged
2307 jfenwick 2271 forceResolve();
2308     exclusiveWrite();
2309 gross 1358 if (isConstant()) tag();
2310 jfenwick 2271 WrappedArray w(value);
2311 jgs 94
2312 jfenwick 1796 DataVector temp_data2;
2313 jfenwick 2271 temp_data2.copyFromArray(w,1);
2314 jfenwick 1796
2315 jfenwick 2271 m_data->setTaggedValue(tagKey,w.getShape(), temp_data2);
2316 jgs 94 }
2317    
2318 jfenwick 1796
2319 jgs 110 void
2320 jgs 121 Data::setTaggedValueFromCPP(int tagKey,
2321 jfenwick 1796 const DataTypes::ShapeType& pointshape,
2322     const DataTypes::ValueType& value,
2323     int dataOffset)
2324 jgs 121 {
2325 gross 783 if (isProtected()) {
2326     throw DataException("Error - attempt to update protected Data object.");
2327     }
2328 jgs 121 //
2329     // Ensure underlying data object is of type DataTagged
2330 jfenwick 2271 forceResolve();
2331 gross 1358 if (isConstant()) tag();
2332 jfenwick 2271 exclusiveWrite();
2333 jgs 121 //
2334     // Call DataAbstract::setTaggedValue
2335 jfenwick 1796 m_data->setTaggedValue(tagKey,pointshape, value, dataOffset);
2336 jgs 121 }
2337    
2338 jgs 149 int
2339     Data::getTagNumber(int dpno)
2340     {
2341 jfenwick 1803 if (isEmpty())
2342     {
2343     throw DataException("Error - operation not permitted on instances of DataEmpty.");
2344     }
2345 gross 1409 return getFunctionSpace().getTagFromDataPointNo(dpno);
2346 jgs 149 }
2347    
2348 jgs 119
2349 jgs 94 ostream& escript::operator<<(ostream& o, const Data& data)
2350     {
2351     o << data.toString();
2352     return o;
2353     }
2354 bcumming 782
2355 ksteube 813 Data
2356     escript::C_GeneralTensorProduct(Data& arg_0,
2357     Data& arg_1,
2358     int axis_offset,
2359     int transpose)
2360     {
2361 gross 826 // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2362 ksteube 813 // SM is the product of the last axis_offset entries in arg_0.getShape().
2363    
2364 jfenwick 2005 // deal with any lazy data
2365 jfenwick 2066 // if (arg_0.isLazy()) {arg_0.resolve();}
2366     // if (arg_1.isLazy()) {arg_1.resolve();}
2367 jfenwick 2199 if (arg_0.isLazy() || arg_1.isLazy() || (AUTOLAZYON && (arg_0.isExpanded() || arg_1.isExpanded())))
2368 jfenwick 2066 {
2369     DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2370     return Data(c);
2371     }
2372 jfenwick 2005
2373 ksteube 813 // Interpolate if necessary and find an appropriate function space
2374 gross 826 Data arg_0_Z, arg_1_Z;
2375 ksteube 813 if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2376     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2377 gross 826 arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2378     arg_1_Z = Data(arg_1);
2379 ksteube 813 }
2380     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2381 gross 826 arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2382     arg_0_Z =Data(arg_0);
2383 ksteube 813 }
2384     else {
2385     throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2386     }
2387 gross 826 } else {
2388     arg_0_Z = Data(arg_0);
2389     arg_1_Z = Data(arg_1);
2390 ksteube 813 }
2391     // Get rank and shape of inputs
2392 gross 826 int rank0 = arg_0_Z.getDataPointRank();
2393     int rank1 = arg_1_Z.getDataPointRank();
2394 jfenwick 1796 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2395     const DataTypes::ShapeType& shape1 = arg_1_Z.getDataPointShape();
2396 ksteube 813
2397     // Prepare for the loops of the product and verify compatibility of shapes
2398     int start0=0, start1=0;
2399     if (transpose == 0) {}
2400     else if (transpose == 1) { start0 = axis_offset; }
2401     else if (transpose == 2) { start1 = rank1-axis_offset; }
2402     else { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2403    
2404 jfenwick 1796
2405 ksteube 813 // Adjust the shapes for transpose
2406 jfenwick 1796 DataTypes::ShapeType tmpShape0(rank0); // pre-sizing the vectors rather
2407     DataTypes::ShapeType tmpShape1(rank1); // than using push_back
2408     for (int i=0; i<rank0; i++) { tmpShape0[i]=shape0[(i+start0)%rank0]; }
2409     for (int i=0; i<rank1; i++) { tmpShape1[i]=shape1[(i+start1)%rank1]; }
2410 ksteube 813
2411     #if 0
2412     // For debugging: show shape after transpose
2413     char tmp[100];
2414     std::string shapeStr;
2415     shapeStr = "(";
2416     for (int i=0; i<rank0; i++) { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2417     shapeStr += ")";
2418     cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2419     shapeStr = "(";
2420     for (int i=0; i<rank1; i++) { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2421     shapeStr += ")";
2422     cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2423     #endif
2424    
2425     // Prepare for the loops of the product
2426     int SL=1, SM=1, SR=1;
2427     for (int i=0; i<rank0-axis_offset; i++) {
2428     SL *= tmpShape0[i];
2429     }
2430     for (int i=rank0-axis_offset; i<rank0; i++) {
2431     if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2432     throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2433     }
2434     SM *= tmpShape0[i];
2435     }
2436     for (int i=axis_offset; i<rank1; i++) {
2437     SR *= tmpShape1[i];
2438     }
2439    
2440 jfenwick 1796 // Define the shape of the output (rank of shape is the sum of the loop ranges below)
2441     DataTypes::ShapeType shape2(rank0+rank1-2*axis_offset);
2442     { // block to limit the scope of out_index
2443     int out_index=0;
2444     for (int i=0; i<rank0-axis_offset; i++, ++out_index) { shape2[out_index]=tmpShape0[i]; } // First part of arg_0_Z
2445     for (int i=axis_offset; i<rank1; i++, ++out_index) { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2446     }
2447 ksteube 813
2448 jfenwick 2086 if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2449     {
2450     ostringstream os;
2451     os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2452     throw DataException(os.str());
2453     }
2454    
2455 ksteube 813 // Declare output Data object
2456 gross 826 Data res;
2457 ksteube 813
2458 gross 826 if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2459     res = Data(0.0, shape2, arg_1_Z.getFunctionSpace()); // DataConstant output
2460 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2461     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2462     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2463 ksteube 813 matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2464     }