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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2716 - (hide annotations)
Tue Oct 13 04:58:08 2009 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 92507 byte(s)
Small mod - process C_Tensor??Operation on ExpandedData as a whole 
sample rather than a datapoint at a time.
This is just to bring Non-Lazy into line with Lazy.

1 jgs 94
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 jgs 121 /** \file Data.h */
16 jgs 94
17     #ifndef DATA_H
18     #define DATA_H
19 woo409 757 #include "system_dep.h"
20 jgs 94
21 jfenwick 2271 #include "DataTypes.h"
22 jgs 474 #include "DataAbstract.h"
23     #include "DataAlgorithm.h"
24     #include "FunctionSpace.h"
25     #include "BinaryOp.h"
26     #include "UnaryOp.h"
27     #include "DataException.h"
28 jgs 94
29 jfenwick 2271
30 jgs 94 extern "C" {
31 jgs 474 #include "DataC.h"
32 jfenwick 2271 //#include <omp.h>
33 jgs 94 }
34    
35 ksteube 1312 #include "esysmpi.h"
36 jgs 94 #include <string>
37     #include <algorithm>
38 jfenwick 1693 #include <sstream>
39 jgs 94
40     #include <boost/shared_ptr.hpp>
41     #include <boost/python/object.hpp>
42     #include <boost/python/tuple.hpp>
43    
44 jfenwick 2271 #include "BufferGroup.h"
45    
46 jgs 121 namespace escript {
47    
48     //
49 jgs 149 // Forward declaration for various implementations of Data.
50 jgs 121 class DataConstant;
51     class DataTagged;
52     class DataExpanded;
53 jfenwick 2271 class DataLazy;
54 jgs 121
55 jgs 94 /**
56     \brief
57 jfenwick 1802 Data represents a collection of datapoints.
58 jgs 94
59     Description:
60 jfenwick 1802 Internally, the datapoints are actually stored by a DataAbstract object.
61     The specific instance of DataAbstract used may vary over the lifetime
62     of the Data object.
63     Some methods on this class return references (eg getShape()).
64     These references should not be used after an operation which changes the underlying DataAbstract object.
65     Doing so will lead to invalid memory access.
66     This should not affect any methods exposed via boost::python.
67 jgs 94 */
68     class Data {
69    
70     public:
71    
72 jgs 110 // These typedefs allow function names to be cast to pointers
73     // to functions of the appropriate type when calling unaryOp etc.
74 jgs 94 typedef double (*UnaryDFunPtr)(double);
75     typedef double (*BinaryDFunPtr)(double,double);
76    
77 gross 854
78 jgs 94 /**
79 jgs 102 Constructors.
80     */
81    
82     /**
83 jgs 94 \brief
84     Default constructor.
85     Creates a DataEmpty object.
86     */
87 woo409 757 ESCRIPT_DLL_API
88 jgs 94 Data();
89    
90     /**
91     \brief
92     Copy constructor.
93     WARNING: Only performs a shallow copy.
94     */
95 woo409 757 ESCRIPT_DLL_API
96 jgs 94 Data(const Data& inData);
97    
98     /**
99     \brief
100     Constructor from another Data object. If "what" is different from the
101 jgs 102 function space of inData the inData are tried to be interpolated to what,
102 jgs 94 otherwise a shallow copy of inData is returned.
103     */
104 woo409 757 ESCRIPT_DLL_API
105 jgs 94 Data(const Data& inData,
106     const FunctionSpace& what);
107    
108     /**
109 jfenwick 1796 \brief Copy Data from an existing vector
110     */
111 jgs 94
112 woo409 757 ESCRIPT_DLL_API
113 jfenwick 1796 Data(const DataTypes::ValueType& value,
114     const DataTypes::ShapeType& shape,
115     const FunctionSpace& what=FunctionSpace(),
116     bool expanded=false);
117 jgs 94
118     /**
119     \brief
120 jfenwick 2459 Constructor which creates a Data with points having the specified shape.
121 jgs 94
122     \param value - Input - Single value applied to all Data.
123     \param dataPointShape - Input - The shape of each data point.
124     \param what - Input - A description of what this data represents.
125     \param expanded - Input - Flag, if true fill the entire container with
126     the given value. Otherwise a more efficient storage
127     mechanism will be used.
128     */
129 woo409 757 ESCRIPT_DLL_API
130 jgs 94 Data(double value,
131 jfenwick 1796 const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
132 jgs 94 const FunctionSpace& what=FunctionSpace(),
133     bool expanded=false);
134    
135     /**
136     \brief
137     Constructor which performs a deep copy of a region from another Data object.
138    
139     \param inData - Input - Input Data object.
140     \param region - Input - Region to copy.
141     */
142 woo409 757 ESCRIPT_DLL_API
143 jgs 94 Data(const Data& inData,
144 jfenwick 1796 const DataTypes::RegionType& region);
145 jgs 94
146     /**
147     \brief
148 jfenwick 2458 Constructor which copies data from any object that can be treated like a python array/sequence.
149 jgs 94
150     \param value - Input - Input data.
151     \param what - Input - A description of what this data represents.
152     \param expanded - Input - Flag, if true fill the entire container with
153     the value. Otherwise a more efficient storage
154     mechanism will be used.
155     */
156 woo409 757 ESCRIPT_DLL_API
157 jgs 94 Data(const boost::python::object& value,
158     const FunctionSpace& what=FunctionSpace(),
159     bool expanded=false);
160    
161     /**
162     \brief
163     Constructor which creates a DataConstant.
164 jfenwick 2458 Copies data from any object that can be treated like a python array/sequence.
165     All other parameters are copied from other.
166 jgs 94
167     \param value - Input - Input data.
168     \param other - Input - contains all other parameters.
169     */
170 woo409 757 ESCRIPT_DLL_API
171 jgs 94 Data(const boost::python::object& value,
172     const Data& other);
173    
174     /**
175     \brief
176     Constructor which creates a DataConstant of "shape" with constant value.
177     */
178 woo409 757 ESCRIPT_DLL_API
179 matt 1327 Data(double value,
180     const boost::python::tuple& shape=boost::python::make_tuple(),
181 jgs 94 const FunctionSpace& what=FunctionSpace(),
182     bool expanded=false);
183 jfenwick 1796
184    
185    
186 jgs 151 /**
187 jfenwick 1796 \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
188     Once you have passed the pointer, do not delete it.
189     */
190     ESCRIPT_DLL_API
191 jfenwick 2005 explicit Data(DataAbstract* underlyingdata);
192 jfenwick 1796
193 jfenwick 2005 /**
194     \brief Create a Data based on the supplied DataAbstract
195     */
196     ESCRIPT_DLL_API
197     explicit Data(DataAbstract_ptr underlyingdata);
198 jfenwick 1796
199     /**
200 jgs 151 \brief
201     Destructor
202     */
203 woo409 757 ESCRIPT_DLL_API
204 jgs 151 ~Data();
205 jgs 94
206     /**
207 jfenwick 1799 \brief Make this object a deep copy of "other".
208 jgs 94 */
209 woo409 757 ESCRIPT_DLL_API
210 jgs 94 void
211 jgs 102 copy(const Data& other);
212 jgs 94
213     /**
214 jfenwick 1799 \brief Return a pointer to a deep copy of this object.
215     */
216     ESCRIPT_DLL_API
217 jfenwick 2105 Data
218 jfenwick 1799 copySelf();
219    
220    
221 jfenwick 2005 /**
222     \brief produce a delayed evaluation version of this Data.
223     */
224     ESCRIPT_DLL_API
225     Data
226     delay();
227 jfenwick 1799
228 jfenwick 2005 /**
229     \brief convert the current data into lazy data.
230     */
231     ESCRIPT_DLL_API
232     void
233     delaySelf();
234 jfenwick 1799
235 jfenwick 2005
236 jfenwick 1799 /**
237 jgs 102 Member access methods.
238 jgs 94 */
239    
240     /**
241     \brief
242 matt 1327 switches on update protection
243 gross 783
244     */
245     ESCRIPT_DLL_API
246     void
247     setProtection();
248    
249     /**
250     \brief
251 jfenwick 2465 Returns true, if the data object is protected against update
252 gross 783
253     */
254     ESCRIPT_DLL_API
255     bool
256     isProtected() const;
257 gross 1034
258 jfenwick 2271
259 jfenwick 2455 /**
260     \brief
261 jfenwick 2465 Return the value of a data point as a python tuple.
262 jfenwick 2455 */
263 woo409 757 ESCRIPT_DLL_API
264 jfenwick 2271 const boost::python::object
265     getValueOfDataPointAsTuple(int dataPointNo);
266 jgs 117
267     /**
268     \brief
269 gross 1034 sets the values of a data-point from a python object on this process
270 jgs 121 */
271 woo409 757 ESCRIPT_DLL_API
272 gross 921 void
273 gross 1034 setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
274 jgs 121
275     /**
276     \brief
277 jfenwick 2458 sets the values of a data-point from a array-like object on this process
278 jgs 121 */
279 woo409 757 ESCRIPT_DLL_API
280 jgs 149 void
281 jfenwick 2271 setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
282 jgs 149
283     /**
284     \brief
285 gross 922 sets the values of a data-point on this process
286     */
287     ESCRIPT_DLL_API
288     void
289     setValueOfDataPoint(int dataPointNo, const double);
290    
291     /**
292 jfenwick 2455 \brief Return a data point across all processors as a python tuple.
293 gross 921 */
294     ESCRIPT_DLL_API
295 jfenwick 2271 const boost::python::object
296     getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
297    
298 gross 921 /**
299     \brief
300 jgs 149 Return the tag number associated with the given data-point.
301    
302     */
303 woo409 757 ESCRIPT_DLL_API
304 jgs 149 int
305     getTagNumber(int dpno);
306    
307     /**
308     \brief
309 jgs 94 Return the C wrapper for the Data object.
310     */
311 woo409 757 ESCRIPT_DLL_API
312 jgs 102 escriptDataC
313     getDataC();
314 jgs 94
315 jfenwick 1796
316    
317 jgs 94 /**
318     \brief
319     Return the C wrapper for the Data object - const version.
320     */
321 woo409 757 ESCRIPT_DLL_API
322 jgs 102 escriptDataC
323     getDataC() const;
324 jgs 94
325     /**
326 jfenwick 2271 \brief How much space is required to evaulate a sample of the Data.
327     */
328     ESCRIPT_DLL_API
329     size_t
330     getSampleBufferSize() const;
331    
332    
333    
334     /**
335 jgs 94 \brief
336 jfenwick 1693 Write the data as a string. For large amounts of data, a summary is printed.
337 jgs 94 */
338 woo409 757 ESCRIPT_DLL_API
339 jgs 102 std::string
340 jfenwick 1693 toString() const;
341 jgs 94
342     /**
343     \brief
344 jgs 102 Whatever the current Data type make this into a DataExpanded.
345 jgs 94 */
346 woo409 757 ESCRIPT_DLL_API
347 jgs 94 void
348     expand();
349    
350     /**
351     \brief
352 jgs 102 If possible convert this Data to DataTagged. This will only allow
353 jgs 94 Constant data to be converted to tagged. An attempt to convert
354     Expanded data to tagged will throw an exception.
355     */
356 woo409 757 ESCRIPT_DLL_API
357 jgs 94 void
358     tag();
359    
360     /**
361 jfenwick 2005 \brief If this data is lazy, then convert it to ready data.
362     What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
363     */
364     ESCRIPT_DLL_API
365     void
366     resolve();
367    
368    
369     /**
370 jfenwick 2271 \brief Ensures data is ready for write access.
371     This means that the data will be resolved if lazy and will be copied if shared with another Data object.
372     \warning This method should only be called in single threaded sections of code. (It modifies m_data).
373     Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
374     Doing so might introduce additional sharing.
375     */
376     ESCRIPT_DLL_API
377     void
378     requireWrite();
379    
380     /**
381 jgs 94 \brief
382     Return true if this Data is expanded.
383 jfenwick 2271 \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
384 jgs 94 */
385 woo409 757 ESCRIPT_DLL_API
386 jgs 94 bool
387     isExpanded() const;
388    
389     /**
390     \brief
391 jfenwick 2271 Return true if this Data is expanded or resolves to expanded.
392     That is, if it has a separate value for each datapoint in the sample.
393     */
394     ESCRIPT_DLL_API
395     bool
396     actsExpanded() const;
397    
398    
399     /**
400     \brief
401 jgs 94 Return true if this Data is tagged.
402     */
403 woo409 757 ESCRIPT_DLL_API
404 jgs 94 bool
405     isTagged() const;
406    
407     /**
408     \brief
409 jgs 102 Return true if this Data is constant.
410 jgs 94 */
411 woo409 757 ESCRIPT_DLL_API
412 jgs 94 bool
413 jgs 102 isConstant() const;
414 jgs 94
415     /**
416 jfenwick 2005 \brief Return true if this Data is lazy.
417     */
418     ESCRIPT_DLL_API
419     bool
420     isLazy() const;
421    
422     /**
423     \brief Return true if this data is ready.
424     */
425     ESCRIPT_DLL_API
426     bool
427     isReady() const;
428    
429     /**
430 jgs 94 \brief
431 jfenwick 1803 Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
432     contains datapoints.
433 jgs 94 */
434 woo409 757 ESCRIPT_DLL_API
435 jgs 94 bool
436 jgs 102 isEmpty() const;
437 jgs 94
438     /**
439     \brief
440     Return the function space.
441     */
442 woo409 757 ESCRIPT_DLL_API
443 jgs 94 inline
444     const FunctionSpace&
445     getFunctionSpace() const
446     {
447     return m_data->getFunctionSpace();
448     }
449    
450     /**
451     \brief
452     Return a copy of the function space.
453     */
454 woo409 757 ESCRIPT_DLL_API
455 jgs 94 const FunctionSpace
456     getCopyOfFunctionSpace() const;
457    
458     /**
459     \brief
460     Return the domain.
461     */
462 woo409 757 ESCRIPT_DLL_API
463 jgs 94 inline
464 jfenwick 1872 // const AbstractDomain&
465     const_Domain_ptr
466 jgs 94 getDomain() const
467     {
468     return getFunctionSpace().getDomain();
469     }
470    
471 jfenwick 1872
472 jgs 94 /**
473     \brief
474 jfenwick 1872 Return the domain.
475     TODO: For internal use only. This should be removed.
476     */
477     ESCRIPT_DLL_API
478     inline
479     // const AbstractDomain&
480     Domain_ptr
481     getDomainPython() const
482     {
483     return getFunctionSpace().getDomainPython();
484     }
485    
486     /**
487     \brief
488 jgs 94 Return a copy of the domain.
489     */
490 woo409 757 ESCRIPT_DLL_API
491 jgs 94 const AbstractDomain
492     getCopyOfDomain() const;
493    
494     /**
495     \brief
496     Return the rank of the point data.
497     */
498 woo409 757 ESCRIPT_DLL_API
499 jgs 94 inline
500 jfenwick 1952 unsigned int
501 jgs 94 getDataPointRank() const
502     {
503 jfenwick 1796 return m_data->getRank();
504 jgs 94 }
505    
506     /**
507     \brief
508 gross 921 Return the number of data points
509     */
510     ESCRIPT_DLL_API
511     inline
512     int
513     getNumDataPoints() const
514     {
515     return getNumSamples() * getNumDataPointsPerSample();
516     }
517     /**
518     \brief
519 jgs 102 Return the number of samples.
520 jgs 94 */
521 woo409 757 ESCRIPT_DLL_API
522 jgs 94 inline
523     int
524 jgs 102 getNumSamples() const
525 jgs 94 {
526 jgs 102 return m_data->getNumSamples();
527 jgs 94 }
528    
529     /**
530     \brief
531 jgs 102 Return the number of data points per sample.
532 jgs 94 */
533 woo409 757 ESCRIPT_DLL_API
534 jgs 102 inline
535 jgs 94 int
536 jgs 102 getNumDataPointsPerSample() const
537     {
538     return m_data->getNumDPPSample();
539     }
540 jfenwick 1796
541    
542 gross 950 /**
543 jfenwick 1796 \brief
544     Return the number of values in the shape for this object.
545     */
546     ESCRIPT_DLL_API
547     int
548     getNoValues() const
549     {
550     return m_data->getNoValues();
551     }
552    
553    
554     /**
555 gross 950 \brief
556 matt 1327 dumps the object into a netCDF file
557 gross 950 */
558     ESCRIPT_DLL_API
559     void
560 gross 1141 dump(const std::string fileName) const;
561 jfenwick 2005
562 jfenwick 2459 /**
563     \brief returns the values of the object as a list of tuples (one for each datapoint).
564 jfenwick 2271
565 jfenwick 2459 \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
566     If false, the result is a list of scalars [1, 2, ...]
567     */
568     ESCRIPT_DLL_API
569     const boost::python::object
570 jfenwick 2465 toListOfTuples(bool scalarastuple=true);
571 jfenwick 2271
572 jfenwick 2459
573 jfenwick 2271 /**
574     \brief
575     Return the sample data for the given sample no. This is not the
576     preferred interface but is provided for use by C code.
577     The bufferg parameter is only required for LazyData.
578     \param sampleNo - Input - the given sample no.
579 jfenwick 2519 \param bufferg - A buffer to compute (and store) sample data in will be selected from this group.
580 jfenwick 2271 \return pointer to the sample data.
581     */
582     ESCRIPT_DLL_API
583     inline
584     const DataAbstract::ValueType::value_type*
585     getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg=0);
586    
587    
588 jgs 94 /**
589     \brief
590     Return the sample data for the given sample no. This is not the
591     preferred interface but is provided for use by C code.
592     \param sampleNo - Input - the given sample no.
593 jfenwick 2271 \return pointer to the sample data.
594 jgs 94 */
595 woo409 757 ESCRIPT_DLL_API
596 jgs 102 inline
597 jgs 94 DataAbstract::ValueType::value_type*
598 jfenwick 2271 getSampleDataRW(DataAbstract::ValueType::size_type sampleNo);
599 jgs 94
600 jfenwick 2271
601 jgs 94 /**
602     \brief
603     Return the sample data for the given tag. If an attempt is made to
604     access data that isn't tagged an exception will be thrown.
605     \param tag - Input - the tag key.
606     */
607 woo409 757 ESCRIPT_DLL_API
608 jgs 102 inline
609 jgs 94 DataAbstract::ValueType::value_type*
610 jgs 102 getSampleDataByTag(int tag)
611     {
612     return m_data->getSampleDataByTag(tag);
613     }
614 jgs 94
615     /**
616     \brief
617 jfenwick 2459 Return a reference into the DataVector which points to the specified data point.
618 jgs 94 \param sampleNo - Input -
619     \param dataPointNo - Input -
620     */
621 woo409 757 ESCRIPT_DLL_API
622 jfenwick 1796 DataTypes::ValueType::const_reference
623 jfenwick 2271 getDataPointRO(int sampleNo, int dataPointNo);
624 jfenwick 1796
625 jfenwick 2459 /**
626     \brief
627     Return a reference into the DataVector which points to the specified data point.
628     \param sampleNo - Input -
629     \param dataPointNo - Input -
630     */
631 jfenwick 1796 ESCRIPT_DLL_API
632     DataTypes::ValueType::reference
633 jfenwick 2271 getDataPointRW(int sampleNo, int dataPointNo);
634 jfenwick 1796
635    
636    
637     /**
638     \brief
639     Return the offset for the given sample and point within the sample
640     */
641     ESCRIPT_DLL_API
642 jgs 94 inline
643 jfenwick 1796 DataTypes::ValueType::size_type
644     getDataOffset(int sampleNo,
645 jgs 94 int dataPointNo)
646     {
647 jfenwick 2005 return m_data->getPointOffset(sampleNo,dataPointNo);
648 jgs 94 }
649    
650     /**
651     \brief
652     Return a reference to the data point shape.
653     */
654 woo409 757 ESCRIPT_DLL_API
655 jfenwick 1796 inline
656     const DataTypes::ShapeType&
657     getDataPointShape() const
658     {
659     return m_data->getShape();
660     }
661 jgs 94
662     /**
663     \brief
664 jgs 102 Return the data point shape as a tuple of integers.
665 jgs 94 */
666 woo409 757 ESCRIPT_DLL_API
667 jgs 121 const boost::python::tuple
668 jgs 94 getShapeTuple() const;
669    
670     /**
671     \brief
672     Return the size of the data point. It is the product of the
673     data point shape dimensions.
674     */
675 woo409 757 ESCRIPT_DLL_API
676 jgs 94 int
677     getDataPointSize() const;
678    
679     /**
680     \brief
681 jgs 102 Return the number of doubles stored for this Data.
682 jgs 94 */
683 woo409 757 ESCRIPT_DLL_API
684 jfenwick 1796 DataTypes::ValueType::size_type
685 jgs 94 getLength() const;
686    
687 jfenwick 2646 /**
688     \brief Return true if this object contains no samples.
689     This is not the same as isEmpty()
690     */
691     ESCRIPT_DLL_API
692     bool
693     hasNoSamples() const
694     {
695     return getLength()==0;
696     }
697 gross 1044
698 jgs 94 /**
699     \brief
700 gross 1044 Assign the given value to the tag assocciated with name. Implicitly converts this
701     object to type DataTagged. Throws an exception if this object
702     cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
703 jfenwick 2519 \param name - Input - name of tag.
704 gross 1044 \param value - Input - Value to associate with given key.
705     */
706     ESCRIPT_DLL_API
707     void
708     setTaggedValueByName(std::string name,
709     const boost::python::object& value);
710    
711     /**
712     \brief
713 jgs 102 Assign the given value to the tag. Implicitly converts this
714 gross 1358 object to type DataTagged if it is constant.
715    
716 jgs 102 \param tagKey - Input - Integer key.
717     \param value - Input - Value to associate with given key.
718 jgs 544 ==>*
719 jgs 94 */
720 woo409 757 ESCRIPT_DLL_API
721 jgs 102 void
722     setTaggedValue(int tagKey,
723     const boost::python::object& value);
724    
725     /**
726     \brief
727     Assign the given value to the tag. Implicitly converts this
728 gross 1358 object to type DataTagged if it is constant.
729    
730 jgs 102 \param tagKey - Input - Integer key.
731 jfenwick 1796 \param pointshape - Input - The shape of the value parameter
732 jgs 102 \param value - Input - Value to associate with given key.
733 jfenwick 1796 \param dataOffset - Input - Offset of the begining of the point within the value parameter
734 jgs 102 */
735 woo409 757 ESCRIPT_DLL_API
736 jgs 102 void
737 jgs 121 setTaggedValueFromCPP(int tagKey,
738 jfenwick 1796 const DataTypes::ShapeType& pointshape,
739     const DataTypes::ValueType& value,
740     int dataOffset=0);
741 jgs 102
742 jfenwick 1796
743    
744 jgs 102 /**
745     \brief
746     Copy other Data object into this Data object where mask is positive.
747     */
748 woo409 757 ESCRIPT_DLL_API
749 jgs 102 void
750     copyWithMask(const Data& other,
751     const Data& mask);
752    
753     /**
754     Data object operation methods and operators.
755     */
756    
757     /**
758     \brief
759 gross 1118 set all values to zero
760     *
761     */
762     ESCRIPT_DLL_API
763     void
764     setToZero();
765    
766     /**
767     \brief
768 jgs 102 Interpolates this onto the given functionspace and returns
769     the result as a Data object.
770 jgs 123 *
771 jgs 102 */
772 woo409 757 ESCRIPT_DLL_API
773 jgs 94 Data
774     interpolate(const FunctionSpace& functionspace) const;
775 jfenwick 2628
776    
777     ESCRIPT_DLL_API
778     Data
779     interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
780 gross 2668 double undef, Data& B, double Bmin, double Bstep,bool check_boundaries);
781 jfenwick 2628
782 jfenwick 2646 ESCRIPT_DLL_API
783     Data
784     interpolateFromTable1D(const WrappedArray& table, double Amin, double Astep,
785 gross 2668 double undef,bool check_boundaries);
786 jfenwick 2628
787 jfenwick 2646
788    
789    
790 jfenwick 2628 ESCRIPT_DLL_API
791     Data
792 jfenwick 2646 interpolateFromTable2DP(boost::python::object table, double Amin, double Astep,
793 gross 2668 Data& B, double Bmin, double Bstep, double undef,bool check_boundaries);
794 jfenwick 2628
795 jfenwick 2646 ESCRIPT_DLL_API
796     Data
797     interpolateFromTable1DP(boost::python::object table, double Amin, double Astep,
798 gross 2668 double undef,bool check_boundaries);
799 jfenwick 2646
800 jgs 94 /**
801     \brief
802     Calculates the gradient of the data at the data points of functionspace.
803     If functionspace is not present the function space of Function(getDomain()) is used.
804 jgs 123 *
805 jgs 94 */
806 woo409 757 ESCRIPT_DLL_API
807 jgs 94 Data
808     gradOn(const FunctionSpace& functionspace) const;
809    
810 woo409 757 ESCRIPT_DLL_API
811 jgs 94 Data
812     grad() const;
813    
814     /**
815 jfenwick 2455 \brief
816     Calculate the integral over the function space domain as a python tuple.
817 jgs 94 */
818 woo409 757 ESCRIPT_DLL_API
819 jfenwick 2271 boost::python::object
820 jfenwick 2455 integrateToTuple_const() const;
821 jgs 94
822 jfenwick 2005
823 jfenwick 2455 /**
824     \brief
825     Calculate the integral over the function space domain as a python tuple.
826     */
827 jfenwick 2271 ESCRIPT_DLL_API
828     boost::python::object
829     integrateToTuple();
830    
831    
832    
833 jgs 94 /**
834     \brief
835 gross 854 Returns 1./ Data object
836     *
837     */
838     ESCRIPT_DLL_API
839     Data
840     oneOver() const;
841     /**
842     \brief
843 jgs 94 Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
844 jgs 123 *
845 jgs 94 */
846 woo409 757 ESCRIPT_DLL_API
847 jgs 94 Data
848 gross 698 wherePositive() const;
849 jgs 94
850     /**
851     \brief
852 jgs 102 Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.
853 jgs 123 *
854 jgs 94 */
855 woo409 757 ESCRIPT_DLL_API
856 jgs 94 Data
857 gross 698 whereNegative() const;
858 jgs 102
859     /**
860     \brief
861     Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.
862 jgs 123 *
863 jgs 102 */
864 woo409 757 ESCRIPT_DLL_API
865 jgs 102 Data
866 gross 698 whereNonNegative() const;
867 jgs 94
868     /**
869     \brief
870 jgs 102 Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.
871 jgs 123 *
872 jgs 94 */
873 woo409 757 ESCRIPT_DLL_API
874 jgs 94 Data
875 gross 698 whereNonPositive() const;
876 jgs 94
877     /**
878     \brief
879 jgs 102 Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.
880 jgs 123 *
881 jgs 94 */
882 woo409 757 ESCRIPT_DLL_API
883 jgs 94 Data
884 jgs 571 whereZero(double tol=0.0) const;
885 jgs 94
886     /**
887     \brief
888 jgs 102 Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.
889 jgs 123 *
890 jgs 94 */
891 woo409 757 ESCRIPT_DLL_API
892 jgs 94 Data
893 jgs 571 whereNonZero(double tol=0.0) const;
894 jgs 102
895     /**
896     \brief
897     Return the maximum absolute value of this Data object.
898 jfenwick 2005
899     The method is not const because lazy data needs to be expanded before Lsup can be computed.
900     The _const form can be used when the Data object is const, however this will only work for
901     Data which is not Lazy.
902    
903 jfenwick 1857 For Data which contain no samples (or tagged Data for which no tags in use have a value)
904     zero is returned.
905 jgs 94 */
906 woo409 757 ESCRIPT_DLL_API
907 jgs 102 double
908 jfenwick 2005 Lsup();
909 jgs 94
910 jfenwick 2005 ESCRIPT_DLL_API
911     double
912     Lsup_const() const;
913    
914    
915 jgs 94 /**
916     \brief
917 jgs 102 Return the maximum value of this Data object.
918 jfenwick 2005
919     The method is not const because lazy data needs to be expanded before sup can be computed.
920     The _const form can be used when the Data object is const, however this will only work for
921     Data which is not Lazy.
922    
923 jfenwick 1857 For Data which contain no samples (or tagged Data for which no tags in use have a value)
924     a large negative value is returned.
925 jgs 94 */
926 woo409 757 ESCRIPT_DLL_API
927 jgs 102 double
928 jfenwick 2005 sup();
929 jgs 94
930 jfenwick 2005 ESCRIPT_DLL_API
931     double
932     sup_const() const;
933    
934    
935 jgs 94 /**
936     \brief
937 jgs 102 Return the minimum value of this Data object.
938 jfenwick 2005
939     The method is not const because lazy data needs to be expanded before inf can be computed.
940     The _const form can be used when the Data object is const, however this will only work for
941     Data which is not Lazy.
942    
943 jfenwick 1857 For Data which contain no samples (or tagged Data for which no tags in use have a value)
944     a large positive value is returned.
945 jgs 94 */
946 woo409 757 ESCRIPT_DLL_API
947 jgs 94 double
948 jfenwick 2005 inf();
949 jgs 94
950 jfenwick 2005 ESCRIPT_DLL_API
951     double
952     inf_const() const;
953    
954    
955    
956 jgs 94 /**
957     \brief
958 jgs 102 Return the absolute value of each data point of this Data object.
959 jgs 123 *
960 jgs 94 */
961 woo409 757 ESCRIPT_DLL_API
962 jgs 102 Data
963     abs() const;
964 jgs 94
965     /**
966     \brief
967 jgs 102 Return the maximum value of each data point of this Data object.
968 jgs 123 *
969 jgs 94 */
970 woo409 757 ESCRIPT_DLL_API
971 jgs 102 Data
972     maxval() const;
973 jgs 94
974     /**
975     \brief
976 jgs 102 Return the minimum value of each data point of this Data object.
977 jgs 123 *
978 jgs 94 */
979 woo409 757 ESCRIPT_DLL_API
980 jgs 94 Data
981 jgs 102 minval() const;
982 jgs 94
983     /**
984     \brief
985 jgs 121 Return the (sample number, data-point number) of the data point with
986 jfenwick 2476 the minimum component value in this Data object.
987     \note If you are working in python, please consider using Locator
988     instead of manually manipulating process and point IDs.
989 jgs 121 */
990 woo409 757 ESCRIPT_DLL_API
991 jgs 121 const boost::python::tuple
992 gross 921 minGlobalDataPoint() const;
993 jgs 121
994 jfenwick 2476 /**
995     \brief
996     Return the (sample number, data-point number) of the data point with
997     the minimum component value in this Data object.
998     \note If you are working in python, please consider using Locator
999     instead of manually manipulating process and point IDs.
1000     */
1001 woo409 757 ESCRIPT_DLL_API
1002 jfenwick 2476 const boost::python::tuple
1003     maxGlobalDataPoint() const;
1004    
1005    
1006    
1007 jgs 121 /**
1008     \brief
1009 jgs 102 Return the sign of each data point of this Data object.
1010     -1 for negative values, zero for zero values, 1 for positive values.
1011 jgs 123 *
1012 jgs 94 */
1013 woo409 757 ESCRIPT_DLL_API
1014 jgs 102 Data
1015     sign() const;
1016 jgs 94
1017     /**
1018 jgs 123 \brief
1019 ksteube 775 Return the symmetric part of a matrix which is half the matrix plus its transpose.
1020     *
1021     */
1022     ESCRIPT_DLL_API
1023     Data
1024     symmetric() const;
1025    
1026     /**
1027     \brief
1028     Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
1029     *
1030     */
1031     ESCRIPT_DLL_API
1032     Data
1033     nonsymmetric() const;
1034    
1035     /**
1036     \brief
1037     Return the trace of a matrix
1038     *
1039     */
1040     ESCRIPT_DLL_API
1041     Data
1042 gross 800 trace(int axis_offset) const;
1043 ksteube 775
1044     /**
1045     \brief
1046     Transpose each data point of this Data object around the given axis.
1047     *
1048     */
1049     ESCRIPT_DLL_API
1050     Data
1051     transpose(int axis_offset) const;
1052    
1053     /**
1054     \brief
1055 gross 576 Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1056     Currently this function is restricted to rank 2, square shape, and dimension 3.
1057     *
1058     */
1059 woo409 757 ESCRIPT_DLL_API
1060 gross 576 Data
1061     eigenvalues() const;
1062    
1063     /**
1064     \brief
1065     Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1066 matt 1327 the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1067     tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1068 gross 576 first non-zero entry is positive.
1069     Currently this function is restricted to rank 2, square shape, and dimension 3
1070     *
1071     */
1072 woo409 757 ESCRIPT_DLL_API
1073 gross 576 const boost::python::tuple
1074     eigenvalues_and_eigenvectors(const double tol=1.e-12) const;
1075    
1076     /**
1077     \brief
1078 gross 804 swaps the components axis0 and axis1
1079 jgs 123 *
1080 jgs 102 */
1081 woo409 757 ESCRIPT_DLL_API
1082 jgs 102 Data
1083 gross 804 swapaxes(const int axis0, const int axis1) const;
1084 jgs 102
1085     /**
1086 jgs 123 \brief
1087 ksteube 876 Return the error function erf of each data point of this Data object.
1088     *
1089     */
1090     ESCRIPT_DLL_API
1091     Data
1092     erf() const;
1093    
1094     /**
1095     \brief
1096 jgs 123 Return the sin of each data point of this Data object.
1097     *
1098 jgs 102 */
1099 woo409 757 ESCRIPT_DLL_API
1100 jgs 102 Data
1101 jgs 123 sin() const;
1102    
1103     /**
1104     \brief
1105     Return the cos of each data point of this Data object.
1106     *
1107     */
1108 woo409 757 ESCRIPT_DLL_API
1109 jgs 123 Data
1110     cos() const;
1111    
1112     /**
1113     \brief
1114     Return the tan of each data point of this Data object.
1115     *
1116     */
1117 woo409 757 ESCRIPT_DLL_API
1118 jgs 123 Data
1119     tan() const;
1120    
1121     /**
1122     \brief
1123 jgs 150 Return the asin of each data point of this Data object.
1124     *
1125     */
1126 woo409 757 ESCRIPT_DLL_API
1127 jgs 150 Data
1128     asin() const;
1129    
1130     /**
1131     \brief
1132     Return the acos of each data point of this Data object.
1133     *
1134     */
1135 woo409 757 ESCRIPT_DLL_API
1136 jgs 150 Data
1137     acos() const;
1138    
1139     /**
1140     \brief
1141     Return the atan of each data point of this Data object.
1142     *
1143     */
1144 woo409 757 ESCRIPT_DLL_API
1145 jgs 150 Data
1146     atan() const;
1147    
1148     /**
1149     \brief
1150     Return the sinh of each data point of this Data object.
1151     *
1152     */
1153 woo409 757 ESCRIPT_DLL_API
1154 jgs 150 Data
1155     sinh() const;
1156    
1157     /**
1158     \brief
1159     Return the cosh of each data point of this Data object.
1160     *
1161     */
1162 woo409 757 ESCRIPT_DLL_API
1163 jgs 150 Data
1164     cosh() const;
1165    
1166     /**
1167     \brief
1168     Return the tanh of each data point of this Data object.
1169     *
1170     */
1171 woo409 757 ESCRIPT_DLL_API
1172 jgs 150 Data
1173     tanh() const;
1174    
1175     /**
1176     \brief
1177     Return the asinh of each data point of this Data object.
1178     *
1179     */
1180 woo409 757 ESCRIPT_DLL_API
1181 jgs 150 Data
1182     asinh() const;
1183    
1184     /**
1185     \brief
1186     Return the acosh of each data point of this Data object.
1187     *
1188     */
1189 woo409 757 ESCRIPT_DLL_API
1190 jgs 150 Data
1191     acosh() const;
1192    
1193     /**
1194     \brief
1195     Return the atanh of each data point of this Data object.
1196     *
1197     */
1198 woo409 757 ESCRIPT_DLL_API
1199 jgs 150 Data
1200     atanh() const;
1201    
1202     /**
1203     \brief
1204 jgs 123 Return the log to base 10 of each data point of this Data object.
1205     *
1206     */
1207 woo409 757 ESCRIPT_DLL_API
1208 jgs 123 Data
1209 gross 286 log10() const;
1210 jgs 123
1211     /**
1212     \brief
1213     Return the natural log of each data point of this Data object.
1214     *
1215     */
1216 woo409 757 ESCRIPT_DLL_API
1217 jgs 123 Data
1218 gross 286 log() const;
1219 jgs 123
1220     /**
1221     \brief
1222     Return the exponential function of each data point of this Data object.
1223     *
1224     */
1225 woo409 757 ESCRIPT_DLL_API
1226 jgs 123 Data
1227 jgs 102 exp() const;
1228    
1229     /**
1230 jgs 123 \brief
1231     Return the square root of each data point of this Data object.
1232     *
1233 jgs 102 */
1234 woo409 757 ESCRIPT_DLL_API
1235 jgs 102 Data
1236     sqrt() const;
1237    
1238     /**
1239 jgs 123 \brief
1240     Return the negation of each data point of this Data object.
1241     *
1242 jgs 121 */
1243 woo409 757 ESCRIPT_DLL_API
1244 jgs 121 Data
1245     neg() const;
1246    
1247     /**
1248 jgs 123 \brief
1249     Return the identity of each data point of this Data object.
1250     Simply returns this object unmodified.
1251     *
1252 jgs 121 */
1253 woo409 757 ESCRIPT_DLL_API
1254 jgs 121 Data
1255     pos() const;
1256    
1257     /**
1258 jgs 94 \brief
1259 jgs 102 Return the given power of each data point of this Data object.
1260 jgs 121
1261     \param right Input - the power to raise the object to.
1262 jgs 123 *
1263 jgs 102 */
1264 woo409 757 ESCRIPT_DLL_API
1265 jgs 102 Data
1266     powD(const Data& right) const;
1267    
1268 jgs 121 /**
1269 jgs 123 \brief
1270     Return the given power of each data point of this boost python object.
1271 matt 1327
1272 jgs 123 \param right Input - the power to raise the object to.
1273     *
1274 jgs 121 */
1275 woo409 757 ESCRIPT_DLL_API
1276 jgs 102 Data
1277     powO(const boost::python::object& right) const;
1278    
1279     /**
1280 jgs 123 \brief
1281 gross 699 Return the given power of each data point of this boost python object.
1282 matt 1327
1283 gross 699 \param left Input - the bases
1284     *
1285     */
1286    
1287 woo409 757 ESCRIPT_DLL_API
1288 gross 699 Data
1289     rpowO(const boost::python::object& left) const;
1290    
1291     /**
1292     \brief
1293 jgs 123 writes the object to a file in the DX file format
1294 jgs 104 */
1295 woo409 757 ESCRIPT_DLL_API
1296 jgs 104 void
1297     saveDX(std::string fileName) const;
1298    
1299     /**
1300 jgs 123 \brief
1301     writes the object to a file in the VTK file format
1302 jgs 110 */
1303 woo409 757 ESCRIPT_DLL_API
1304 jgs 110 void
1305     saveVTK(std::string fileName) const;
1306    
1307 jfenwick 2635
1308    
1309 jgs 110 /**
1310 jgs 102 \brief
1311     Overloaded operator +=
1312     \param right - Input - The right hand side.
1313 jgs 123 *
1314 jgs 102 */
1315 woo409 757 ESCRIPT_DLL_API
1316 jgs 102 Data& operator+=(const Data& right);
1317 woo409 757 ESCRIPT_DLL_API
1318 jgs 102 Data& operator+=(const boost::python::object& right);
1319    
1320 ksteube 1312 ESCRIPT_DLL_API
1321     Data& operator=(const Data& other);
1322    
1323 jgs 102 /**
1324     \brief
1325     Overloaded operator -=
1326     \param right - Input - The right hand side.
1327 jgs 123 *
1328 jgs 102 */
1329 woo409 757 ESCRIPT_DLL_API
1330 jgs 102 Data& operator-=(const Data& right);
1331 woo409 757 ESCRIPT_DLL_API
1332 jgs 102 Data& operator-=(const boost::python::object& right);
1333    
1334     /**
1335     \brief
1336     Overloaded operator *=
1337     \param right - Input - The right hand side.
1338 jgs 123 *
1339 jgs 102 */
1340 woo409 757 ESCRIPT_DLL_API
1341 jgs 102 Data& operator*=(const Data& right);
1342 woo409 757 ESCRIPT_DLL_API
1343 jgs 102 Data& operator*=(const boost::python::object& right);
1344    
1345     /**
1346     \brief
1347     Overloaded operator /=
1348     \param right - Input - The right hand side.
1349 jgs 123 *
1350 jgs 102 */
1351 woo409 757 ESCRIPT_DLL_API
1352 jgs 102 Data& operator/=(const Data& right);
1353 woo409 757 ESCRIPT_DLL_API
1354 jgs 102 Data& operator/=(const boost::python::object& right);
1355    
1356     /**
1357     \brief
1358 jgs 94 Returns true if this can be interpolated to functionspace.
1359     */
1360 woo409 757 ESCRIPT_DLL_API
1361 jgs 94 bool
1362     probeInterpolation(const FunctionSpace& functionspace) const;
1363    
1364     /**
1365 jgs 102 Data object slicing methods.
1366     */
1367    
1368     /**
1369 jgs 94 \brief
1370 jgs 102 Returns a slice from this Data object.
1371    
1372     /description
1373     Implements the [] get operator in python.
1374     Calls getSlice.
1375    
1376     \param key - Input - python slice tuple specifying
1377     slice to return.
1378 jgs 94 */
1379 woo409 757 ESCRIPT_DLL_API
1380 jgs 102 Data
1381     getItem(const boost::python::object& key) const;
1382    
1383     /**
1384     \brief
1385     Copies slice from value into this Data object.
1386    
1387     Implements the [] set operator in python.
1388     Calls setSlice.
1389    
1390     \param key - Input - python slice tuple specifying
1391     slice to copy from value.
1392     \param value - Input - Data object to copy from.
1393     */
1394 woo409 757 ESCRIPT_DLL_API
1395 jgs 94 void
1396 jgs 102 setItemD(const boost::python::object& key,
1397     const Data& value);
1398 jgs 94
1399 woo409 757 ESCRIPT_DLL_API
1400 jgs 102 void
1401     setItemO(const boost::python::object& key,
1402     const boost::python::object& value);
1403    
1404     // These following public methods should be treated as private.
1405    
1406 jgs 94 /**
1407     \brief
1408 jgs 102 Perform the given unary operation on every element of every data point in
1409     this Data object.
1410 jgs 94 */
1411 jgs 102 template <class UnaryFunction>
1412 woo409 757 ESCRIPT_DLL_API
1413 jgs 102 inline
1414 jgs 94 void
1415 matt 1334 unaryOp2(UnaryFunction operation);
1416 jgs 102
1417     /**
1418     \brief
1419     Return a Data object containing the specified slice of
1420     this Data object.
1421     \param region - Input - Region to copy.
1422 jgs 123 *
1423 jgs 94 */
1424 woo409 757 ESCRIPT_DLL_API
1425 jgs 102 Data
1426 jfenwick 1796 getSlice(const DataTypes::RegionType& region) const;
1427 jgs 94
1428 jgs 102 /**
1429     \brief
1430     Copy the specified slice from the given value into this
1431     Data object.
1432     \param value - Input - Data to copy from.
1433     \param region - Input - Region to copy.
1434 jgs 123 *
1435 jgs 102 */
1436 woo409 757 ESCRIPT_DLL_API
1437 jgs 102 void
1438     setSlice(const Data& value,
1439 jfenwick 1796 const DataTypes::RegionType& region);
1440 jgs 102
1441 jgs 119 /**
1442     \brief
1443 matt 1327 print the data values to stdout. Used for debugging
1444 bcumming 751 */
1445 bcumming 782 ESCRIPT_DLL_API
1446 matt 1327 void
1447 matt 1332 print(void);
1448 bcumming 751
1449 bcumming 782 /**
1450     \brief
1451     return the MPI rank number of the local data
1452 matt 1332 MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1453     is returned
1454 bcumming 782 */
1455     ESCRIPT_DLL_API
1456 matt 1332 int
1457     get_MPIRank(void) const;
1458 bcumming 782
1459     /**
1460     \brief
1461     return the MPI rank number of the local data
1462 matt 1332 MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1463     is returned
1464 bcumming 782 */
1465     ESCRIPT_DLL_API
1466 matt 1332 int
1467     get_MPISize(void) const;
1468 bcumming 782
1469     /**
1470     \brief
1471     return the MPI rank number of the local data
1472 matt 1332 MPI_COMM_WORLD is assumed and returned.
1473 bcumming 782 */
1474     ESCRIPT_DLL_API
1475 matt 1332 MPI_Comm
1476     get_MPIComm(void) const;
1477 ksteube 813
1478     /**
1479     \brief
1480     return the object produced by the factory, which is a DataConstant or DataExpanded
1481 jfenwick 1796 TODO Ownership of this object should be explained in doco.
1482 ksteube 813 */
1483     ESCRIPT_DLL_API
1484 matt 1332 DataAbstract*
1485     borrowData(void) const;
1486 ksteube 813
1487 jfenwick 2005 ESCRIPT_DLL_API
1488     DataAbstract_ptr
1489     borrowDataPtr(void) const;
1490 jfenwick 1796
1491 jfenwick 2005 ESCRIPT_DLL_API
1492     DataReady_ptr
1493     borrowReadyPtr(void) const;
1494    
1495    
1496    
1497 jfenwick 1796 /**
1498     \brief
1499     Return a pointer to the beginning of the datapoint at the specified offset.
1500     TODO Eventually these should be inlined.
1501     \param i - position(offset) in the underlying datastructure
1502     */
1503 jfenwick 2271
1504 jfenwick 1796 ESCRIPT_DLL_API
1505     DataTypes::ValueType::const_reference
1506 jfenwick 2271 getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1507 jfenwick 1796
1508    
1509     ESCRIPT_DLL_API
1510     DataTypes::ValueType::reference
1511 jfenwick 2271 getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1512 jfenwick 1796
1513 jfenwick 2271
1514    
1515     /**
1516     \brief Create a buffer for use by getSample
1517     Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1518     Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1519    
1520     In multi-threaded sections, this needs to be called on each thread.
1521    
1522     \return A BufferGroup* if Data is lazy, NULL otherwise.
1523     \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1524     */
1525     ESCRIPT_DLL_API
1526     BufferGroup*
1527     allocSampleBuffer() const;
1528    
1529     /**
1530     \brief Free a buffer allocated with allocSampleBuffer.
1531     \param buffer Input - pointer to the buffer to deallocate.
1532     */
1533     ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1534    
1535 jgs 102 protected:
1536    
1537 jgs 94 private:
1538    
1539 jfenwick 2005 double
1540     LsupWorker() const;
1541    
1542     double
1543     supWorker() const;
1544    
1545     double
1546     infWorker() const;
1547    
1548 jfenwick 2271 boost::python::object
1549 jfenwick 2005 integrateWorker() const;
1550    
1551 jfenwick 2476 void
1552     calc_minGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1553    
1554     void
1555     calc_maxGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1556    
1557    
1558 jgs 94 /**
1559     \brief
1560 jgs 102 Check *this and the right operand are compatible. Throws
1561     an exception if they aren't.
1562     \param right - Input - The right hand side.
1563     */
1564     inline
1565     void
1566     operandCheck(const Data& right) const
1567     {
1568     return m_data->operandCheck(*(right.m_data.get()));
1569     }
1570    
1571     /**
1572     \brief
1573     Perform the specified reduction algorithm on every element of every data point in
1574 jgs 113 this Data object according to the given function and return the single value result.
1575 jgs 102 */
1576 jgs 147 template <class BinaryFunction>
1577 jgs 102 inline
1578     double
1579 jgs 147 algorithm(BinaryFunction operation,
1580     double initial_value) const;
1581 jgs 102
1582 jgs 113 /**
1583     \brief
1584     Reduce each data-point in this Data object using the given operation. Return a Data
1585     object with the same number of data-points, but with each data-point containing only
1586     one value - the result of the reduction operation on the corresponding data-point in
1587     this Data object
1588     */
1589 jgs 147 template <class BinaryFunction>
1590 jgs 106 inline
1591     Data
1592 jgs 147 dp_algorithm(BinaryFunction operation,
1593     double initial_value) const;
1594 jgs 106
1595 jgs 102 /**
1596     \brief
1597     Perform the given binary operation on all of the data's elements.
1598     The underlying type of the right hand side (right) determines the final
1599     type of *this after the operation. For example if the right hand side
1600     is expanded *this will be expanded if necessary.
1601     RHS is a Data object.
1602     */
1603     template <class BinaryFunction>
1604     inline
1605     void
1606     binaryOp(const Data& right,
1607     BinaryFunction operation);
1608    
1609     /**
1610     \brief
1611     Convert the data type of the RHS to match this.
1612     \param right - Input - data type to match.
1613     */
1614     void
1615     typeMatchLeft(Data& right) const;
1616    
1617     /**
1618     \brief
1619     Convert the data type of this to match the RHS.
1620     \param right - Input - data type to match.
1621     */
1622     void
1623     typeMatchRight(const Data& right);
1624    
1625     /**
1626     \brief
1627 jgs 94 Construct a Data object of the appropriate type.
1628     */
1629 jfenwick 1796
1630 jgs 94 void
1631 jfenwick 1796 initialise(const DataTypes::ValueType& value,
1632     const DataTypes::ShapeType& shape,
1633 jgs 94 const FunctionSpace& what,
1634     bool expanded);
1635    
1636 jfenwick 1796 void
1637 jfenwick 2271 initialise(const WrappedArray& value,
1638 jfenwick 1796 const FunctionSpace& what,
1639     bool expanded);
1640    
1641 jgs 94 //
1642 gross 783 // flag to protect the data object against any update
1643     bool m_protected;
1644 jfenwick 2271 mutable bool m_shared;
1645     bool m_lazy;
1646 gross 783
1647     //
1648 jgs 102 // pointer to the actual data object
1649 jfenwick 1872 // boost::shared_ptr<DataAbstract> m_data;
1650     DataAbstract_ptr m_data;
1651 jgs 94
1652 jfenwick 2271 // If possible please use getReadyPtr instead.
1653     // But see warning below.
1654 jfenwick 2005 const DataReady*
1655     getReady() const;
1656    
1657     DataReady*
1658     getReady();
1659    
1660 jfenwick 2271
1661     // Be wary of using this for local operations since it (temporarily) increases reference count.
1662     // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1663     // getReady() instead
1664 jfenwick 2005 DataReady_ptr
1665     getReadyPtr();
1666    
1667     const_DataReady_ptr
1668     getReadyPtr() const;
1669    
1670    
1671 jfenwick 2271 /**
1672     \brief Update the Data's shared flag
1673     This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1674     For internal use only.
1675     */
1676     void updateShareStatus(bool nowshared) const
1677     {
1678     m_shared=nowshared; // m_shared is mutable
1679     }
1680    
1681     // In the isShared() method below:
1682     // A problem would occur if m_data (the address pointed to) were being modified
1683     // while the call m_data->is_shared is being executed.
1684     //
1685     // Q: So why do I think this code can be thread safe/correct?
1686     // A: We need to make some assumptions.
1687     // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1688     // 2. We assume that no constructions or assignments which will share previously unshared
1689     // will occur while this call is executing. This is consistent with the way Data:: and C are written.
1690     //
1691     // This means that the only transition we need to consider, is when a previously shared object is
1692     // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1693     // In those cases the m_shared flag changes to false after m_data has completed changing.
1694     // For any threads executing before the flag switches they will assume the object is still shared.
1695     bool isShared() const
1696     {
1697     return m_shared;
1698     /* if (m_shared) return true;
1699     if (m_data->isShared())
1700     {
1701     updateShareStatus(true);
1702     return true;
1703     }
1704     return false;*/
1705     }
1706    
1707     void forceResolve()
1708     {
1709     if (isLazy())
1710     {
1711     #ifdef _OPENMP
1712     if (omp_in_parallel())
1713     { // Yes this is throwing an exception out of an omp thread which is forbidden.
1714     throw DataException("Please do not call forceResolve() in a parallel region.");
1715     }
1716     #endif
1717     resolve();
1718     }
1719     }
1720    
1721     /**
1722     \brief if another object is sharing out member data make a copy to work with instead.
1723     This code should only be called from single threaded sections of code.
1724     */
1725     void exclusiveWrite()
1726     {
1727     #ifdef _OPENMP
1728     if (omp_in_parallel())
1729     {
1730     // *((int*)0)=17;
1731     throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1732     }
1733     #endif
1734     forceResolve();
1735     if (isShared())
1736     {
1737     DataAbstract* t=m_data->deepCopy();
1738     set_m_data(DataAbstract_ptr(t));
1739     }
1740     }
1741    
1742     /**
1743     \brief checks if caller can have exclusive write to the object
1744     */
1745     void checkExclusiveWrite()
1746     {
1747     if (isLazy() || isShared())
1748     {
1749     throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1750     }
1751     }
1752    
1753     /**
1754     \brief Modify the data abstract hosted by this Data object
1755     For internal use only.
1756     Passing a pointer to null is permitted (do this in the destructor)
1757     \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1758     */
1759     void set_m_data(DataAbstract_ptr p);
1760    
1761     friend class DataAbstract; // To allow calls to updateShareStatus
1762    
1763 jgs 94 };
1764    
1765 jfenwick 2005 } // end namespace escript
1766 jgs 94
1767 jfenwick 1796
1768 jfenwick 2005 // No, this is not supposed to be at the top of the file
1769     // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1770     // so that I can dynamic cast between them below.
1771     #include "DataReady.h"
1772 jfenwick 2271 #include "DataLazy.h"
1773 jfenwick 2005
1774     namespace escript
1775     {
1776    
1777     inline
1778     const DataReady*
1779     Data::getReady() const
1780     {
1781     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1782     EsysAssert((dr!=0), "Error - casting to DataReady.");
1783     return dr;
1784     }
1785    
1786     inline
1787     DataReady*
1788     Data::getReady()
1789     {
1790     DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1791     EsysAssert((dr!=0), "Error - casting to DataReady.");
1792     return dr;
1793     }
1794    
1795 jfenwick 2271 // Be wary of using this for local operations since it (temporarily) increases reference count.
1796     // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1797     // getReady() instead
1798 jfenwick 2005 inline
1799     DataReady_ptr
1800     Data::getReadyPtr()
1801     {
1802     DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1803     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1804     return dr;
1805     }
1806    
1807    
1808     inline
1809     const_DataReady_ptr
1810     Data::getReadyPtr() const
1811     {
1812     const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1813     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1814     return dr;
1815     }
1816    
1817     inline
1818     DataAbstract::ValueType::value_type*
1819 jfenwick 2271 Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1820 jfenwick 2005 {
1821     if (isLazy())
1822     {
1823 jfenwick 2271 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1824 jfenwick 2005 }
1825 jfenwick 2271 return getReady()->getSampleDataRW(sampleNo);
1826 jfenwick 2005 }
1827    
1828 jfenwick 2271 inline
1829     const DataAbstract::ValueType::value_type*
1830     Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1831     {
1832     DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1833     if (l!=0)
1834     {
1835     size_t offset=0;
1836     if (bufferg==NULL)
1837     {
1838     throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1839     }
1840     const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1841     return &((*res)[offset]);
1842     }
1843     return getReady()->getSampleDataRO(sampleNo);
1844     }
1845    
1846    
1847    
1848 jgs 102 /**
1849 ksteube 1748 Modify a filename for MPI parallel output to multiple files
1850     */
1851     char *Escript_MPI_appendRankToFileName(const char *, int, int);
1852    
1853     /**
1854 jgs 102 Binary Data object operators.
1855     */
1856 matt 1327 inline double rpow(double x,double y)
1857 gross 854 {
1858     return pow(y,x);
1859 gross 1028 }
1860 jgs 94
1861     /**
1862     \brief
1863     Operator+
1864     Takes two Data objects.
1865     */
1866 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
1867 jgs 94
1868     /**
1869     \brief
1870     Operator-
1871     Takes two Data objects.
1872     */
1873 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1874 jgs 94
1875     /**
1876     \brief
1877     Operator*
1878     Takes two Data objects.
1879     */
1880 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1881 jgs 94
1882     /**
1883     \brief
1884     Operator/
1885     Takes two Data objects.
1886     */
1887 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1888 jgs 94
1889     /**
1890     \brief
1891     Operator+
1892     Takes LHS Data object and RHS python::object.
1893     python::object must be convertable to Data type.
1894     */
1895 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1896 jgs 94
1897     /**
1898     \brief
1899     Operator-
1900     Takes LHS Data object and RHS python::object.
1901     python::object must be convertable to Data type.
1902     */
1903 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1904 jgs 94
1905     /**
1906     \brief
1907     Operator*
1908     Takes LHS Data object and RHS python::object.
1909     python::object must be convertable to Data type.
1910     */
1911 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1912 jgs 94
1913     /**
1914     \brief
1915     Operator/
1916     Takes LHS Data object and RHS python::object.
1917     python::object must be convertable to Data type.
1918     */
1919 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1920 jgs 94
1921     /**
1922     \brief
1923     Operator+
1924     Takes LHS python::object and RHS Data object.
1925     python::object must be convertable to Data type.
1926     */
1927 woo409 757 ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1928 jgs 94
1929     /**
1930     \brief
1931     Operator-
1932     Takes LHS python::object and RHS Data object.
1933     python::object must be convertable to Data type.
1934     */
1935 woo409 757 ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1936 jgs 94
1937     /**
1938     \brief
1939     Operator*
1940     Takes LHS python::object and RHS Data object.
1941     python::object must be convertable to Data type.
1942     */
1943 woo409 757 ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1944 jgs 94
1945     /**
1946     \brief
1947     Operator/
1948     Takes LHS python::object and RHS Data object.
1949     python::object must be convertable to Data type.
1950     */
1951 woo409 757 ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1952 jgs 94
1953 ksteube 1312
1954    
1955 jgs 94 /**
1956     \brief
1957     Output operator
1958     */
1959 woo409 757 ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1960 jgs 94
1961     /**
1962     \brief
1963 ksteube 813 Compute a tensor product of two Data objects
1964 jfenwick 2519 \param arg_0 - Input - Data object
1965     \param arg_1 - Input - Data object
1966 ksteube 813 \param axis_offset - Input - axis offset
1967     \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1968     */
1969     ESCRIPT_DLL_API
1970     Data
1971 jfenwick 2519 C_GeneralTensorProduct(Data& arg_0,
1972     Data& arg_1,
1973 ksteube 813 int axis_offset=0,
1974     int transpose=0);
1975    
1976 jgs 102 /**
1977     \brief
1978     Perform the given binary operation with this and right as operands.
1979     Right is a Data object.
1980     */
1981 jgs 94 template <class BinaryFunction>
1982     inline
1983     void
1984     Data::binaryOp(const Data& right,
1985     BinaryFunction operation)
1986     {
1987     //
1988     // if this has a rank of zero promote it to the rank of the RHS
1989 jfenwick 1796 if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1990 gross 854 throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1991 jgs 94 }
1992 jfenwick 2005
1993     if (isLazy() || right.isLazy())
1994     {
1995     throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1996     }
1997 jgs 94 //
1998     // initially make the temporary a shallow copy
1999     Data tempRight(right);
2000 jfenwick 2005
2001 jgs 94 if (getFunctionSpace()!=right.getFunctionSpace()) {
2002     if (right.probeInterpolation(getFunctionSpace())) {
2003     //
2004 matt 1327 // an interpolation is required so create a new Data
2005 jgs 94 tempRight=Data(right,this->getFunctionSpace());
2006     } else if (probeInterpolation(right.getFunctionSpace())) {
2007     //
2008     // interpolate onto the RHS function space
2009     Data tempLeft(*this,right.getFunctionSpace());
2010 jfenwick 2271 // m_data=tempLeft.m_data;
2011     set_m_data(tempLeft.m_data);
2012 jgs 94 }
2013     }
2014     operandCheck(tempRight);
2015     //
2016     // ensure this has the right type for the RHS
2017 jgs 102 typeMatchRight(tempRight);
2018 jgs 94 //
2019     // Need to cast to the concrete types so that the correct binaryOp
2020     // is called.
2021     if (isExpanded()) {
2022     //
2023     // Expanded data will be done in parallel, the right hand side can be
2024     // of any data type
2025     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2026     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2027 jfenwick 2005 escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2028 jgs 94 } else if (isTagged()) {
2029     //
2030     // Tagged data is operated on serially, the right hand side can be
2031     // either DataConstant or DataTagged
2032     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2033     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2034     if (right.isTagged()) {
2035     DataTagged* rightC=dynamic_cast<DataTagged*>(tempRight.m_data.get());
2036     EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
2037     escript::binaryOp(*leftC,*rightC,operation);
2038     } else {
2039     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2040     EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
2041     escript::binaryOp(*leftC,*rightC,operation);
2042     }
2043 jgs 102 } else if (isConstant()) {
2044 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2045     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2046     EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2047     escript::binaryOp(*leftC,*rightC,operation);
2048     }
2049     }
2050    
2051 jgs 102 /**
2052     \brief
2053     Perform the given Data object reduction algorithm on this and return the result.
2054     Given operation combines each element of each data point, thus argument
2055     object (*this) is a rank n Data object, and returned object is a scalar.
2056     Calls escript::algorithm.
2057     */
2058 jgs 147 template <class BinaryFunction>
2059 jgs 94 inline
2060     double
2061 jgs 147 Data::algorithm(BinaryFunction operation, double initial_value) const
2062 jgs 94 {
2063     if (isExpanded()) {
2064     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2065     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2066 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2067 jgs 102 } else if (isTagged()) {
2068 jgs 94 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2069     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2070 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2071 jgs 102 } else if (isConstant()) {
2072 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2073     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2074 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2075 jfenwick 1803 } else if (isEmpty()) {
2076     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2077 jfenwick 2005 } else if (isLazy()) {
2078     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2079     } else {
2080     throw DataException("Error - Data encapsulates an unknown type.");
2081 jgs 94 }
2082     }
2083    
2084 jgs 102 /**
2085     \brief
2086     Perform the given data point reduction algorithm on data and return the result.
2087     Given operation combines each element within each data point into a scalar,
2088 matt 1327 thus argument object is a rank n Data object, and returned object is a
2089 jgs 102 rank 0 Data object.
2090     Calls escript::dp_algorithm.
2091     */
2092 jgs 147 template <class BinaryFunction>
2093 jgs 94 inline
2094     Data
2095 jgs 147 Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2096 jgs 94 {
2097 jfenwick 1803 if (isEmpty()) {
2098     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2099     }
2100     else if (isExpanded()) {
2101 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2102 jgs 106 DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2103 jgs 102 DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2104     EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2105     EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2106 jgs 147 escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2107 jgs 559 return result;
2108 jfenwick 1803 }
2109     else if (isTagged()) {
2110 jgs 106 DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2111 jgs 102 EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2112 jfenwick 1796 DataTypes::ValueType defval(1);
2113     defval[0]=0;
2114     DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2115 jgs 147 escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2116 jfenwick 1796 return Data(resultT); // note: the Data object now owns the resultT pointer
2117 jfenwick 1803 }
2118     else if (isConstant()) {
2119 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2120 jgs 106 DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2121 jgs 102 DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2122     EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2123     EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2124 jgs 147 escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2125 jgs 559 return result;
2126 jfenwick 2005 } else if (isLazy()) {
2127     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2128     } else {
2129     throw DataException("Error - Data encapsulates an unknown type.");
2130 jgs 102 }
2131 jgs 94 }
2132    
2133 trankine 1430 /**
2134     \brief
2135     Compute a tensor operation with two Data objects
2136 jfenwick 2519 \param arg_0 - Input - Data object
2137     \param arg_1 - Input - Data object
2138 trankine 1430 \param operation - Input - Binary op functor
2139     */
2140 matt 1327 template <typename BinaryFunction>
2141 trankine 1430 inline
2142 matt 1327 Data
2143     C_TensorBinaryOperation(Data const &arg_0,
2144 matt 1332 Data const &arg_1,
2145     BinaryFunction operation)
2146 matt 1327 {
2147 jfenwick 1803 if (arg_0.isEmpty() || arg_1.isEmpty())
2148     {
2149     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2150     }
2151 jfenwick 2005 if (arg_0.isLazy() || arg_1.isLazy())
2152     {
2153     throw DataException("Error - Operations not permitted on lazy data.");
2154     }
2155 matt 1327 // Interpolate if necessary and find an appropriate function space
2156     Data arg_0_Z, arg_1_Z;
2157     if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2158     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2159     arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2160     arg_1_Z = Data(arg_1);
2161     }
2162     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2163     arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2164     arg_0_Z =Data(arg_0);
2165     }
2166     else {
2167     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2168     }
2169     } else {
2170     arg_0_Z = Data(arg_0);
2171     arg_1_Z = Data(arg_1);
2172     }
2173     // Get rank and shape of inputs
2174     int rank0 = arg_0_Z.getDataPointRank();
2175     int rank1 = arg_1_Z.getDataPointRank();
2176 jfenwick 1796 DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2177     DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2178 matt 1327 int size0 = arg_0_Z.getDataPointSize();
2179     int size1 = arg_1_Z.getDataPointSize();
2180     // Declare output Data object
2181     Data res;
2182    
2183     if (shape0 == shape1) {
2184     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2185 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2186 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2187     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2188     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2189 jfenwick 1796
2190 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2191     }
2192     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2193    
2194     // Prepare the DataConstant input
2195     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2196    
2197     // Borrow DataTagged input from Data object
2198     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2199    
2200     // Prepare a DataTagged output 2
2201 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2202 matt 1327 res.tag();
2203     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2204    
2205     // Prepare offset into DataConstant
2206     int offset_0 = tmp_0->getPointOffset(0,0);
2207 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2208 jfenwick 1796
2209 matt 1327 // Get the pointers to the actual data
2210 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2211     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2212 jfenwick 1796
2213 matt 1327 // Compute a result for the default
2214     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2215     // Compute a result for each tag
2216     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2217     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2218     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2219 jfenwick 1796 tmp_2->addTag(i->first);
2220 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2221     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2222 jfenwick 1796
2223 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2224 matt 1327 }
2225    
2226     }
2227     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2228     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2229     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2230     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2231     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2232    
2233     int sampleNo_1,dataPointNo_1;
2234     int numSamples_1 = arg_1_Z.getNumSamples();
2235     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2236     int offset_0 = tmp_0->getPointOffset(0,0);
2237 jfenwick 2271 res.requireWrite();
2238 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2239     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2240 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2241     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2242     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2243 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2244     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2245     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2246 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2247     }
2248 matt 1327 }
2249    
2250     }
2251     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2252     // Borrow DataTagged input from Data object
2253     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2254    
2255     // Prepare the DataConstant input
2256     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2257    
2258     // Prepare a DataTagged output 2
2259 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2260 matt 1327 res.tag();
2261     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2262    
2263     // Prepare offset into DataConstant
2264     int offset_1 = tmp_1->getPointOffset(0,0);
2265 jfenwick 2271
2266     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2267 matt 1327 // Get the pointers to the actual data
2268 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2269     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2270 matt 1327 // Compute a result for the default
2271     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2272     // Compute a result for each tag
2273     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2274     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2275     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2276 jfenwick 1796 tmp_2->addTag(i->first);
2277 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2278     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2279 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2280 matt 1327 }
2281    
2282     }
2283     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2284     // Borrow DataTagged input from Data object
2285     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2286    
2287     // Borrow DataTagged input from Data object
2288     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2289    
2290     // Prepare a DataTagged output 2
2291     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2292 matt 1332 res.tag(); // DataTagged output
2293 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2294    
2295     // Get the pointers to the actual data
2296 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2297     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2298     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2299 jfenwick 1796
2300 matt 1327 // Compute a result for the default
2301     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2302     // Merge the tags
2303     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2304     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2305     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2306     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2307 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2308 matt 1327 }
2309     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2310 jfenwick 1796 tmp_2->addTag(i->first);
2311 matt 1327 }
2312     // Compute a result for each tag
2313     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2314     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2315 jfenwick 1796
2316 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2317     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2318     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2319 jfenwick 1796
2320 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2321 matt 1327 }
2322    
2323     }
2324     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2325     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2326     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2327     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2328     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2329     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2330    
2331     int sampleNo_0,dataPointNo_0;
2332     int numSamples_0 = arg_0_Z.getNumSamples();
2333     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2334 jfenwick 2271 res.requireWrite();
2335 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2336     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2337 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2338 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2339 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2340     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2341     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2342 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2343     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2344 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2345     }
2346 matt 1327 }
2347    
2348     }
2349     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2350     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2351     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2352     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2353     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2354    
2355     int sampleNo_0,dataPointNo_0;
2356     int numSamples_0 = arg_0_Z.getNumSamples();
2357     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2358     int offset_1 = tmp_1->getPointOffset(0,0);
2359 jfenwick 2271 res.requireWrite();
2360 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2361     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2362 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2363     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2364     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2365 jfenwick 1796
2366 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2367     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2368     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2369 jfenwick 1796
2370    
2371 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2372     }
2373 matt 1327 }
2374    
2375     }
2376     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2377     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2378     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2379     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2380     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2381     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2382    
2383     int sampleNo_0,dataPointNo_0;
2384     int numSamples_0 = arg_0_Z.getNumSamples();
2385     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2386 jfenwick 2271 res.requireWrite();
2387 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2388     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2389 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2390 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2391 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2392     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2393     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2394 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2395     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2396 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2397     }
2398 matt 1327 }
2399    
2400     }
2401     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2402     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2403     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2404     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2405     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2406     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2407    
2408     int sampleNo_0,dataPointNo_0;
2409     int numSamples_0 = arg_0_Z.getNumSamples();
2410     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2411 jfenwick 2271 res.requireWrite();
2412 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2413     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2414 jfenwick 2716 dataPointNo_0=0;
2415     // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2416 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2417     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2418     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2419 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2420     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2421     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2422 jfenwick 2716 tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2423     // }
2424 matt 1327 }
2425    
2426     }
2427     else {
2428     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2429     }
2430    
2431     } else if (0 == rank0) {
2432     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2433 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataConstant output
2434 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2435     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2436     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2437 matt 1327 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2438     }
2439     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2440    
2441     // Prepare the DataConstant input
2442     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2443    
2444     // Borrow DataTagged input from Data object
2445     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2446    
2447     // Prepare a DataTagged output 2
2448 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataTagged output
2449 matt 1327 res.tag();
2450     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2451    
2452     // Prepare offset into DataConstant
2453     int offset_0 = tmp_0->getPointOffset(0,0);
2454 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2455 jfenwick 1796
2456 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2457     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2458    
2459 matt 1327 // Compute a result for the default
2460     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2461     // Compute a result for each tag
2462     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2463     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2464     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2465 jfenwick 1796 tmp_2->addTag(i->first);
2466 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2467     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2468 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2469 matt 1327 }
2470    
2471     }
2472     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2473    
2474     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2475     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2476     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2477     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2478    
2479     int sampleNo_1,dataPointNo_1;
2480     int numSamples_1 = arg_1_Z.getNumSamples();
2481     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2482     int offset_0 = tmp_0->getPointOffset(0,0);
2483 jfenwick 2271 res.requireWrite();
2484 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2485     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2486 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2487     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2488     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2489 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2490     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2491     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2492 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2493 matt 1327
2494 matt 1332 }
2495 matt 1327 }
2496    
2497     }
2498     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2499    
2500     // Borrow DataTagged input from Data object
2501     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2502    
2503     // Prepare the DataConstant input
2504     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2505    
2506     // Prepare a DataTagged output 2
2507 matt 1332 res = Data(0.0, shape1, arg_0_Z.getFunctionSpace()); // DataTagged output
2508 matt 1327 res.tag();
2509     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2510    
2511     // Prepare offset into DataConstant
2512     int offset_1 = tmp_1->getPointOffset(0,0);
2513 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2514 jfenwick 1796
2515     // Get the pointers to the actual data
2516 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2517     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2518 jfenwick 1796
2519    
2520 matt 1327 // Compute a result for the default
2521     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2522     // Compute a result for each tag
2523     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2524     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2525     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2526 jfenwick 1796 tmp_2->addTag(i->first);
2527 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2528     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2529 jfenwick 1796
2530 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2531 matt 1327 }
2532    
2533     }
2534     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2535    
2536     // Borrow DataTagged input from Data object
2537     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2538    
2539     // Borrow DataTagged input from Data object
2540     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2541    
2542     // Prepare a DataTagged output 2
2543     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2544 matt 1332 res.tag(); // DataTagged output
2545 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2546    
2547     // Get the pointers to the actual data
2548 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2549     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2550     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2551 jfenwick 1796
2552 matt 1327 // Compute a result for the default
2553     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2554     // Merge the tags
2555     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2556     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2557     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2558     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2559 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2560 matt 1327 }
2561     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2562 jfenwick 1796 tmp_2->addTag(i->first);
2563 matt 1327 }
2564     // Compute a result for each tag
2565     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2566     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2567 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2568     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2569     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2570 jfenwick 1796
2571 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2572 matt 1327 }
2573    
2574     }
2575     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2576    
2577     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2578     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2579     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2580     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2581     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2582    
2583     int sampleNo_0,dataPointNo_0;
2584     int numSamples_0 = arg_0_Z.getNumSamples();
2585     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2586 jfenwick 2271 res.requireWrite();
2587 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2588     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2589 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2590 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2591 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2592     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2593     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2594 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2595     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2596 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2597     }
2598 matt 1327 }
2599    
2600     }
2601     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2602     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2603     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2604     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2605     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2606    
2607     int sampleNo_0,dataPointNo_0;
2608     int numSamples_0 = arg_0_Z.getNumSamples();
2609     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2610     int offset_1 = tmp_1->getPointOffset(0,0);
2611 jfenwick 2271 res.requireWrite();
2612 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2613     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2614 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2615     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2616     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2617 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2618     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2619     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2620 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2621     }
2622 matt 1327 }
2623    
2624    
2625     }
2626     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2627    
2628     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2629     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2630     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2631     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2632     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2633    
2634     int sampleNo_0,dataPointNo_0;
2635     int numSamples_0 = arg_0_Z.getNumSamples();
2636     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2637 jfenwick 2271 res.requireWrite();
2638 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2639     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2640 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2641 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2642 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2643     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2644     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2645 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2646     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2647 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2648     }
2649 matt 1327 }
2650    
2651     }
2652     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2653    
2654     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2655     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2656     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2657     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2658     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2659    
2660     int sampleNo_0,dataPointNo_0;
2661     int numSamples_0 = arg_0_Z.getNumSamples();
2662     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2663 jfenwick 2271 res.requireWrite();
2664 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2665     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2666 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2667     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2668     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2669     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2670 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2671     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2672     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2673 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2674     }
2675 matt 1327 }
2676    
2677     }
2678     else {
2679     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2680     }
2681    
2682     } else if (0 == rank1) {
2683     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2684 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2685 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2686     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2687     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2688 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2689     }
2690     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2691    
2692     // Prepare the DataConstant input
2693     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2694    
2695     // Borrow DataTagged input from Data object
2696     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2697    
2698     // Prepare a DataTagged output 2
2699 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2700 matt 1327 res.tag();
2701     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2702    
2703     // Prepare offset into DataConstant
2704     int offset_0 = tmp_0->getPointOffset(0,0);
2705 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2706    
2707 jfenwick 1796 //Get the pointers to the actual data
2708 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2709     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2710 jfenwick 1796
2711 matt 1327 // Compute a result for the default
2712     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2713     // Compute a result for each tag
2714     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2715     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2716     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2717 jfenwick 1796 tmp_2->addTag(i->first);
2718 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2719     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2720 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2721 matt 1327 }
2722     }
2723     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2724    
2725     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2726     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2727     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2728     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2729    
2730     int sampleNo_1,dataPointNo_1;
2731     int numSamples_1 = arg_1_Z.getNumSamples();
2732     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2733     int offset_0 = tmp_0->getPointOffset(0,0);
2734 jfenwick 2271 res.requireWrite();
2735 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2736     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2737 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2738     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2739     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2740 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2741     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2742     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2743 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2744     }
2745 matt 1327 }
2746    
2747     }
2748     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2749    
2750     // Borrow DataTagged input from Data object
2751     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2752    
2753     // Prepare the DataConstant input
2754     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2755    
2756     // Prepare a DataTagged output 2
2757 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2758 matt 1327 res.tag();
2759     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2760    
2761     // Prepare offset into DataConstant
2762     int offset_1 = tmp_1->getPointOffset(0,0);
2763 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2764 matt 1327 // Get the pointers to the actual data
2765 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2766     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2767 matt 1327 // Compute a result for the default
2768     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2769     // Compute a result for each tag
2770     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2771     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2772     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2773 jfenwick 1796 tmp_2->addTag(i->first);
2774 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2775     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2776 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2777 matt 1327 }
2778    
2779     }
2780     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2781    
2782     // Borrow DataTagged input from Data object
2783     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2784    
2785     // Borrow DataTagged input from Data object
2786     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2787    
2788     // Prepare a DataTagged output 2
2789     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2790 matt 1332 res.tag(); // DataTagged output
2791 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2792    
2793     // Get the pointers to the actual data
2794 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2795     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2796     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2797 jfenwick 1796
2798 matt 1327 // Compute a result for the default
2799     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2800     // Merge the tags
2801     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2802     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2803     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2804     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2805 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2806 matt 1327 }
2807     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2808 jfenwick 1796 tmp_2->addTag(i->first);
2809 matt 1327 }
2810     // Compute a result for each tag
2811     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2812     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2813 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2814     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2815     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2816 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2817 matt 1327 }
2818    
2819     }
2820     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2821    
2822     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2823     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2824     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2825     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2826     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2827    
2828     int sampleNo_0,dataPointNo_0;
2829     int numSamples_0 = arg_0_Z.getNumSamples();
2830     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2831 jfenwick 2271 res.requireWrite();
2832 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2833     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2834 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2835 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2836 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2837     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2838     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2839 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2840     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2841 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2842     }
2843 matt 1327 }
2844    
2845     }
2846     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2847     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2848     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2849     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2850     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2851    
2852     int sampleNo_0,dataPointNo_0;
2853     int numSamples_0 = arg_0_Z.getNumSamples();
2854     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2855     int offset_1 = tmp_1->getPointOffset(0,0);
2856 jfenwick 2271 res.requireWrite();
2857 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2858     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2859 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2860     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2861     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2862 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2863     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2864     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2865 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2866     }
2867 matt 1327 }
2868    
2869    
2870     }
2871     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2872    
2873     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2874     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2875     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2876     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2877     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2878    
2879     int sampleNo_0,dataPointNo_0;
2880     int numSamples_0 = arg_0_Z.getNumSamples();
2881     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2882 jfenwick 2271 res.requireWrite();
2883 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2884     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2885 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2886 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2887 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2888     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2889     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2890 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2891     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2892 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2893     }
2894 matt 1327 }
2895    
2896     }
2897     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2898    
2899     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2900     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2901     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2902     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2903     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2904    
2905     int sampleNo_0,dataPointNo_0;
2906     int numSamples_0 = arg_0_Z.getNumSamples();
2907     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2908 jfenwick 2271 res.requireWrite();
2909 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2910     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2911 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2912     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2913     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2914     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2915 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2916     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2917     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2918 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2919     }
2920 matt 1327 }
2921    
2922     }
2923     else {
2924     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2925     }
2926    
2927     } else {
2928     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2929     }
2930    
2931     return res;
2932 jgs 94 }
2933 matt 1327
2934 matt 1334 template <typename UnaryFunction>
2935     Data
2936     C_TensorUnaryOperation(Data const &arg_0,
2937     UnaryFunction operation)
2938     {
2939 jfenwick 1803 if (arg_0.isEmpty()) // do this before we attempt to interpolate
2940     {
2941     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2942     }
2943 jfenwick 2005 if (arg_0.isLazy())
2944     {
2945     throw DataException("Error - Operations not permitted on lazy data.");
2946     }
2947 matt 1334 // Interpolate if necessary and find an appropriate function space
2948     Data arg_0_Z = Data(arg_0);
2949    
2950     // Get rank and shape of inputs
2951 jfenwick 1796 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2952 matt 1334 int size0 = arg_0_Z.getDataPointSize();
2953    
2954     // Declare output Data object
2955     Data res;
2956    
2957     if (arg_0_Z.isConstant()) {
2958     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataConstant output
2959 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2960     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2961 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2962     }
2963     else if (arg_0_Z.isTagged()) {
2964    
2965     // Borrow DataTagged input from Data object
2966     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2967    
2968     // Prepare a DataTagged output 2
2969     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2970     res.tag();
2971     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2972    
2973     // Get the pointers to the actual data
2974 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2975     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2976 matt 1334 // Compute a result for the default
2977     tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2978     // Compute a result for each tag
2979     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2980     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2981     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2982 jfenwick 1796 tmp_2->addTag(i->first);
2983 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2984     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2985 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2986     }
2987    
2988     }
2989     else if (arg_0_Z.isExpanded()) {
2990    
2991     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2992     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2993     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2994    
2995     int sampleNo_0,dataPointNo_0;
2996     int numSamples_0 = arg_0_Z.getNumSamples();
2997     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2998     #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2999     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
3000 jfenwick 2716 dataPointNo_0=0;
3001     // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3002 matt 1334 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3003     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3004 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3005     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3006 jfenwick 2716 tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3007     // }
3008 matt 1334 }
3009     }
3010     else {
3011     throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3012     }
3013    
3014     return res;
3015 matt 1327 }
3016 matt 1334
3017     }
3018 jgs 94 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26