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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2644 - (hide annotations)
Wed Sep 2 04:14:03 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 91882 byte(s)
Add unit tests for saveDataCSV which should be ready for use now.
Keyword args are now output in sorted order.

Removed some unneeded functions.
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 gross 1044
688    
689 jgs 94 /**
690     \brief
691 gross 1044 Assign the given value to the tag assocciated with name. Implicitly converts this
692     object to type DataTagged. Throws an exception if this object
693     cannot be converted to a DataTagged object or name cannot be mapped onto a tag key.
694 jfenwick 2519 \param name - Input - name of tag.
695 gross 1044 \param value - Input - Value to associate with given key.
696     */
697     ESCRIPT_DLL_API
698     void
699     setTaggedValueByName(std::string name,
700     const boost::python::object& value);
701    
702     /**
703     \brief
704 jgs 102 Assign the given value to the tag. Implicitly converts this
705 gross 1358 object to type DataTagged if it is constant.
706    
707 jgs 102 \param tagKey - Input - Integer key.
708     \param value - Input - Value to associate with given key.
709 jgs 544 ==>*
710 jgs 94 */
711 woo409 757 ESCRIPT_DLL_API
712 jgs 102 void
713     setTaggedValue(int tagKey,
714     const boost::python::object& value);
715    
716     /**
717     \brief
718     Assign the given value to the tag. Implicitly converts this
719 gross 1358 object to type DataTagged if it is constant.
720    
721 jgs 102 \param tagKey - Input - Integer key.
722 jfenwick 1796 \param pointshape - Input - The shape of the value parameter
723 jgs 102 \param value - Input - Value to associate with given key.
724 jfenwick 1796 \param dataOffset - Input - Offset of the begining of the point within the value parameter
725 jgs 102 */
726 woo409 757 ESCRIPT_DLL_API
727 jgs 102 void
728 jgs 121 setTaggedValueFromCPP(int tagKey,
729 jfenwick 1796 const DataTypes::ShapeType& pointshape,
730     const DataTypes::ValueType& value,
731     int dataOffset=0);
732 jgs 102
733 jfenwick 1796
734    
735 jgs 102 /**
736     \brief
737     Copy other Data object into this Data object where mask is positive.
738     */
739 woo409 757 ESCRIPT_DLL_API
740 jgs 102 void
741     copyWithMask(const Data& other,
742     const Data& mask);
743    
744     /**
745     Data object operation methods and operators.
746     */
747    
748     /**
749     \brief
750 gross 1118 set all values to zero
751     *
752     */
753     ESCRIPT_DLL_API
754     void
755     setToZero();
756    
757     /**
758     \brief
759 jgs 102 Interpolates this onto the given functionspace and returns
760     the result as a Data object.
761 jgs 123 *
762 jgs 102 */
763 woo409 757 ESCRIPT_DLL_API
764 jgs 94 Data
765     interpolate(const FunctionSpace& functionspace) const;
766 jfenwick 2628
767    
768     ESCRIPT_DLL_API
769     Data
770     interpolateFromTable2D(const WrappedArray& table, double Amin, double Astep,
771     double undef, Data& B, double Bmin, double Bstep);
772    
773    
774     // This signature needs to be tuned
775     ESCRIPT_DLL_API
776     Data
777     interpolateFromTable(boost::python::object table, double Amin, double Astep,
778     double undef, Data& B, double Bmin, double Bstep);
779    
780 jgs 94 /**
781     \brief
782     Calculates the gradient of the data at the data points of functionspace.
783     If functionspace is not present the function space of Function(getDomain()) is used.
784 jgs 123 *
785 jgs 94 */
786 woo409 757 ESCRIPT_DLL_API
787 jgs 94 Data
788     gradOn(const FunctionSpace& functionspace) const;
789    
790 woo409 757 ESCRIPT_DLL_API
791 jgs 94 Data
792     grad() const;
793    
794     /**
795 jfenwick 2455 \brief
796     Calculate the integral over the function space domain as a python tuple.
797 jgs 94 */
798 woo409 757 ESCRIPT_DLL_API
799 jfenwick 2271 boost::python::object
800 jfenwick 2455 integrateToTuple_const() const;
801 jgs 94
802 jfenwick 2005
803 jfenwick 2455 /**
804     \brief
805     Calculate the integral over the function space domain as a python tuple.
806     */
807 jfenwick 2271 ESCRIPT_DLL_API
808     boost::python::object
809     integrateToTuple();
810    
811    
812    
813 jgs 94 /**
814     \brief
815 gross 854 Returns 1./ Data object
816     *
817     */
818     ESCRIPT_DLL_API
819     Data
820     oneOver() const;
821     /**
822     \brief
823 jgs 94 Return a Data with a 1 for +ive values and a 0 for 0 or -ive values.
824 jgs 123 *
825 jgs 94 */
826 woo409 757 ESCRIPT_DLL_API
827 jgs 94 Data
828 gross 698 wherePositive() const;
829 jgs 94
830     /**
831     \brief
832 jgs 102 Return a Data with a 1 for -ive values and a 0 for +ive or 0 values.
833 jgs 123 *
834 jgs 94 */
835 woo409 757 ESCRIPT_DLL_API
836 jgs 94 Data
837 gross 698 whereNegative() const;
838 jgs 102
839     /**
840     \brief
841     Return a Data with a 1 for +ive or 0 values and a 0 for -ive values.
842 jgs 123 *
843 jgs 102 */
844 woo409 757 ESCRIPT_DLL_API
845 jgs 102 Data
846 gross 698 whereNonNegative() const;
847 jgs 94
848     /**
849     \brief
850 jgs 102 Return a Data with a 1 for -ive or 0 values and a 0 for +ive values.
851 jgs 123 *
852 jgs 94 */
853 woo409 757 ESCRIPT_DLL_API
854 jgs 94 Data
855 gross 698 whereNonPositive() const;
856 jgs 94
857     /**
858     \brief
859 jgs 102 Return a Data with a 1 for 0 values and a 0 for +ive or -ive values.
860 jgs 123 *
861 jgs 94 */
862 woo409 757 ESCRIPT_DLL_API
863 jgs 94 Data
864 jgs 571 whereZero(double tol=0.0) const;
865 jgs 94
866     /**
867     \brief
868 jgs 102 Return a Data with a 0 for 0 values and a 1 for +ive or -ive values.
869 jgs 123 *
870 jgs 94 */
871 woo409 757 ESCRIPT_DLL_API
872 jgs 94 Data
873 jgs 571 whereNonZero(double tol=0.0) const;
874 jgs 102
875     /**
876     \brief
877     Return the maximum absolute value of this Data object.
878 jfenwick 2005
879     The method is not const because lazy data needs to be expanded before Lsup can be computed.
880     The _const form can be used when the Data object is const, however this will only work for
881     Data which is not Lazy.
882    
883 jfenwick 1857 For Data which contain no samples (or tagged Data for which no tags in use have a value)
884     zero is returned.
885 jgs 94 */
886 woo409 757 ESCRIPT_DLL_API
887 jgs 102 double
888 jfenwick 2005 Lsup();
889 jgs 94
890 jfenwick 2005 ESCRIPT_DLL_API
891     double
892     Lsup_const() const;
893    
894    
895 jgs 94 /**
896     \brief
897 jgs 102 Return the maximum value of this Data object.
898 jfenwick 2005
899     The method is not const because lazy data needs to be expanded before sup 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     a large negative value is returned.
905 jgs 94 */
906 woo409 757 ESCRIPT_DLL_API
907 jgs 102 double
908 jfenwick 2005 sup();
909 jgs 94
910 jfenwick 2005 ESCRIPT_DLL_API
911     double
912     sup_const() const;
913    
914    
915 jgs 94 /**
916     \brief
917 jgs 102 Return the minimum value of this Data object.
918 jfenwick 2005
919     The method is not const because lazy data needs to be expanded before inf 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 positive value is returned.
925 jgs 94 */
926 woo409 757 ESCRIPT_DLL_API
927 jgs 94 double
928 jfenwick 2005 inf();
929 jgs 94
930 jfenwick 2005 ESCRIPT_DLL_API
931     double
932     inf_const() const;
933    
934    
935    
936 jgs 94 /**
937     \brief
938 jgs 102 Return the absolute value of each data point of this Data object.
939 jgs 123 *
940 jgs 94 */
941 woo409 757 ESCRIPT_DLL_API
942 jgs 102 Data
943     abs() const;
944 jgs 94
945     /**
946     \brief
947 jgs 102 Return the maximum value of each data point of this Data object.
948 jgs 123 *
949 jgs 94 */
950 woo409 757 ESCRIPT_DLL_API
951 jgs 102 Data
952     maxval() const;
953 jgs 94
954     /**
955     \brief
956 jgs 102 Return the minimum value of each data point of this Data object.
957 jgs 123 *
958 jgs 94 */
959 woo409 757 ESCRIPT_DLL_API
960 jgs 94 Data
961 jgs 102 minval() const;
962 jgs 94
963     /**
964     \brief
965 jgs 121 Return the (sample number, data-point number) of the data point with
966 jfenwick 2476 the minimum component value in this Data object.
967     \note If you are working in python, please consider using Locator
968     instead of manually manipulating process and point IDs.
969 jgs 121 */
970 woo409 757 ESCRIPT_DLL_API
971 jgs 121 const boost::python::tuple
972 gross 921 minGlobalDataPoint() const;
973 jgs 121
974 jfenwick 2476 /**
975     \brief
976     Return the (sample number, data-point number) of the data point with
977     the minimum component value in this Data object.
978     \note If you are working in python, please consider using Locator
979     instead of manually manipulating process and point IDs.
980     */
981 woo409 757 ESCRIPT_DLL_API
982 jfenwick 2476 const boost::python::tuple
983     maxGlobalDataPoint() const;
984    
985    
986    
987 jgs 121 /**
988     \brief
989 jgs 102 Return the sign of each data point of this Data object.
990     -1 for negative values, zero for zero values, 1 for positive values.
991 jgs 123 *
992 jgs 94 */
993 woo409 757 ESCRIPT_DLL_API
994 jgs 102 Data
995     sign() const;
996 jgs 94
997     /**
998 jgs 123 \brief
999 ksteube 775 Return the symmetric part of a matrix which is half the matrix plus its transpose.
1000     *
1001     */
1002     ESCRIPT_DLL_API
1003     Data
1004     symmetric() const;
1005    
1006     /**
1007     \brief
1008     Return the nonsymmetric part of a matrix which is half the matrix minus its transpose.
1009     *
1010     */
1011     ESCRIPT_DLL_API
1012     Data
1013     nonsymmetric() const;
1014    
1015     /**
1016     \brief
1017     Return the trace of a matrix
1018     *
1019     */
1020     ESCRIPT_DLL_API
1021     Data
1022 gross 800 trace(int axis_offset) const;
1023 ksteube 775
1024     /**
1025     \brief
1026     Transpose each data point of this Data object around the given axis.
1027     *
1028     */
1029     ESCRIPT_DLL_API
1030     Data
1031     transpose(int axis_offset) const;
1032    
1033     /**
1034     \brief
1035 gross 576 Return the eigenvalues of the symmetric part at each data point of this Data object in increasing values.
1036     Currently this function is restricted to rank 2, square shape, and dimension 3.
1037     *
1038     */
1039 woo409 757 ESCRIPT_DLL_API
1040 gross 576 Data
1041     eigenvalues() const;
1042    
1043     /**
1044     \brief
1045     Return the eigenvalues and corresponding eigenvcetors of the symmetric part at each data point of this Data object.
1046 matt 1327 the eigenvalues are ordered in increasing size where eigenvalues with relative difference less than
1047     tol are treated as equal. The eigenvectors are orthogonal, normalized and the sclaed such that the
1048 gross 576 first non-zero entry is positive.
1049     Currently this function is restricted to rank 2, square shape, and dimension 3
1050     *
1051     */
1052 woo409 757 ESCRIPT_DLL_API
1053 gross 576 const boost::python::tuple
1054     eigenvalues_and_eigenvectors(const double tol=1.e-12) const;
1055    
1056     /**
1057     \brief
1058 gross 804 swaps the components axis0 and axis1
1059 jgs 123 *
1060 jgs 102 */
1061 woo409 757 ESCRIPT_DLL_API
1062 jgs 102 Data
1063 gross 804 swapaxes(const int axis0, const int axis1) const;
1064 jgs 102
1065     /**
1066 jgs 123 \brief
1067 ksteube 876 Return the error function erf of each data point of this Data object.
1068     *
1069     */
1070     ESCRIPT_DLL_API
1071     Data
1072     erf() const;
1073    
1074     /**
1075     \brief
1076 jgs 123 Return the sin of each data point of this Data object.
1077     *
1078 jgs 102 */
1079 woo409 757 ESCRIPT_DLL_API
1080 jgs 102 Data
1081 jgs 123 sin() const;
1082    
1083     /**
1084     \brief
1085     Return the cos of each data point of this Data object.
1086     *
1087     */
1088 woo409 757 ESCRIPT_DLL_API
1089 jgs 123 Data
1090     cos() const;
1091    
1092     /**
1093     \brief
1094     Return the tan of each data point of this Data object.
1095     *
1096     */
1097 woo409 757 ESCRIPT_DLL_API
1098 jgs 123 Data
1099     tan() const;
1100    
1101     /**
1102     \brief
1103 jgs 150 Return the asin of each data point of this Data object.
1104     *
1105     */
1106 woo409 757 ESCRIPT_DLL_API
1107 jgs 150 Data
1108     asin() const;
1109    
1110     /**
1111     \brief
1112     Return the acos of each data point of this Data object.
1113     *
1114     */
1115 woo409 757 ESCRIPT_DLL_API
1116 jgs 150 Data
1117     acos() const;
1118    
1119     /**
1120     \brief
1121     Return the atan of each data point of this Data object.
1122     *
1123     */
1124 woo409 757 ESCRIPT_DLL_API
1125 jgs 150 Data
1126     atan() const;
1127    
1128     /**
1129     \brief
1130     Return the sinh of each data point of this Data object.
1131     *
1132     */
1133 woo409 757 ESCRIPT_DLL_API
1134 jgs 150 Data
1135     sinh() const;
1136    
1137     /**
1138     \brief
1139     Return the cosh of each data point of this Data object.
1140     *
1141     */
1142 woo409 757 ESCRIPT_DLL_API
1143 jgs 150 Data
1144     cosh() const;
1145    
1146     /**
1147     \brief
1148     Return the tanh of each data point of this Data object.
1149     *
1150     */
1151 woo409 757 ESCRIPT_DLL_API
1152 jgs 150 Data
1153     tanh() const;
1154    
1155     /**
1156     \brief
1157     Return the asinh of each data point of this Data object.
1158     *
1159     */
1160 woo409 757 ESCRIPT_DLL_API
1161 jgs 150 Data
1162     asinh() const;
1163    
1164     /**
1165     \brief
1166     Return the acosh of each data point of this Data object.
1167     *
1168     */
1169 woo409 757 ESCRIPT_DLL_API
1170 jgs 150 Data
1171     acosh() const;
1172    
1173     /**
1174     \brief
1175     Return the atanh of each data point of this Data object.
1176     *
1177     */
1178 woo409 757 ESCRIPT_DLL_API
1179 jgs 150 Data
1180     atanh() const;
1181    
1182     /**
1183     \brief
1184 jgs 123 Return the log to base 10 of each data point of this Data object.
1185     *
1186     */
1187 woo409 757 ESCRIPT_DLL_API
1188 jgs 123 Data
1189 gross 286 log10() const;
1190 jgs 123
1191     /**
1192     \brief
1193     Return the natural log of each data point of this Data object.
1194     *
1195     */
1196 woo409 757 ESCRIPT_DLL_API
1197 jgs 123 Data
1198 gross 286 log() const;
1199 jgs 123
1200     /**
1201     \brief
1202     Return the exponential function of each data point of this Data object.
1203     *
1204     */
1205 woo409 757 ESCRIPT_DLL_API
1206 jgs 123 Data
1207 jgs 102 exp() const;
1208    
1209     /**
1210 jgs 123 \brief
1211     Return the square root of each data point of this Data object.
1212     *
1213 jgs 102 */
1214 woo409 757 ESCRIPT_DLL_API
1215 jgs 102 Data
1216     sqrt() const;
1217    
1218     /**
1219 jgs 123 \brief
1220     Return the negation of each data point of this Data object.
1221     *
1222 jgs 121 */
1223 woo409 757 ESCRIPT_DLL_API
1224 jgs 121 Data
1225     neg() const;
1226    
1227     /**
1228 jgs 123 \brief
1229     Return the identity of each data point of this Data object.
1230     Simply returns this object unmodified.
1231     *
1232 jgs 121 */
1233 woo409 757 ESCRIPT_DLL_API
1234 jgs 121 Data
1235     pos() const;
1236    
1237     /**
1238 jgs 94 \brief
1239 jgs 102 Return the given power of each data point of this Data object.
1240 jgs 121
1241     \param right Input - the power to raise the object to.
1242 jgs 123 *
1243 jgs 102 */
1244 woo409 757 ESCRIPT_DLL_API
1245 jgs 102 Data
1246     powD(const Data& right) const;
1247    
1248 jgs 121 /**
1249 jgs 123 \brief
1250     Return the given power of each data point of this boost python object.
1251 matt 1327
1252 jgs 123 \param right Input - the power to raise the object to.
1253     *
1254 jgs 121 */
1255 woo409 757 ESCRIPT_DLL_API
1256 jgs 102 Data
1257     powO(const boost::python::object& right) const;
1258    
1259     /**
1260 jgs 123 \brief
1261 gross 699 Return the given power of each data point of this boost python object.
1262 matt 1327
1263 gross 699 \param left Input - the bases
1264     *
1265     */
1266    
1267 woo409 757 ESCRIPT_DLL_API
1268 gross 699 Data
1269     rpowO(const boost::python::object& left) const;
1270    
1271     /**
1272     \brief
1273 jgs 123 writes the object to a file in the DX file format
1274 jgs 104 */
1275 woo409 757 ESCRIPT_DLL_API
1276 jgs 104 void
1277     saveDX(std::string fileName) const;
1278    
1279     /**
1280 jgs 123 \brief
1281     writes the object to a file in the VTK file format
1282 jgs 110 */
1283 woo409 757 ESCRIPT_DLL_API
1284 jgs 110 void
1285     saveVTK(std::string fileName) const;
1286    
1287 jfenwick 2635
1288    
1289 jgs 110 /**
1290 jgs 102 \brief
1291     Overloaded operator +=
1292     \param right - Input - The right hand side.
1293 jgs 123 *
1294 jgs 102 */
1295 woo409 757 ESCRIPT_DLL_API
1296 jgs 102 Data& operator+=(const Data& right);
1297 woo409 757 ESCRIPT_DLL_API
1298 jgs 102 Data& operator+=(const boost::python::object& right);
1299    
1300 ksteube 1312 ESCRIPT_DLL_API
1301     Data& operator=(const Data& other);
1302    
1303 jgs 102 /**
1304     \brief
1305     Overloaded operator -=
1306     \param right - Input - The right hand side.
1307 jgs 123 *
1308 jgs 102 */
1309 woo409 757 ESCRIPT_DLL_API
1310 jgs 102 Data& operator-=(const Data& right);
1311 woo409 757 ESCRIPT_DLL_API
1312 jgs 102 Data& operator-=(const boost::python::object& right);
1313    
1314     /**
1315     \brief
1316     Overloaded operator *=
1317     \param right - Input - The right hand side.
1318 jgs 123 *
1319 jgs 102 */
1320 woo409 757 ESCRIPT_DLL_API
1321 jgs 102 Data& operator*=(const Data& right);
1322 woo409 757 ESCRIPT_DLL_API
1323 jgs 102 Data& operator*=(const boost::python::object& right);
1324    
1325     /**
1326     \brief
1327     Overloaded operator /=
1328     \param right - Input - The right hand side.
1329 jgs 123 *
1330 jgs 102 */
1331 woo409 757 ESCRIPT_DLL_API
1332 jgs 102 Data& operator/=(const Data& right);
1333 woo409 757 ESCRIPT_DLL_API
1334 jgs 102 Data& operator/=(const boost::python::object& right);
1335    
1336     /**
1337     \brief
1338 jgs 94 Returns true if this can be interpolated to functionspace.
1339     */
1340 woo409 757 ESCRIPT_DLL_API
1341 jgs 94 bool
1342     probeInterpolation(const FunctionSpace& functionspace) const;
1343    
1344     /**
1345 jgs 102 Data object slicing methods.
1346     */
1347    
1348     /**
1349 jgs 94 \brief
1350 jgs 102 Returns a slice from this Data object.
1351    
1352     /description
1353     Implements the [] get operator in python.
1354     Calls getSlice.
1355    
1356     \param key - Input - python slice tuple specifying
1357     slice to return.
1358 jgs 94 */
1359 woo409 757 ESCRIPT_DLL_API
1360 jgs 102 Data
1361     getItem(const boost::python::object& key) const;
1362    
1363     /**
1364     \brief
1365     Copies slice from value into this Data object.
1366    
1367     Implements the [] set operator in python.
1368     Calls setSlice.
1369    
1370     \param key - Input - python slice tuple specifying
1371     slice to copy from value.
1372     \param value - Input - Data object to copy from.
1373     */
1374 woo409 757 ESCRIPT_DLL_API
1375 jgs 94 void
1376 jgs 102 setItemD(const boost::python::object& key,
1377     const Data& value);
1378 jgs 94
1379 woo409 757 ESCRIPT_DLL_API
1380 jgs 102 void
1381     setItemO(const boost::python::object& key,
1382     const boost::python::object& value);
1383    
1384     // These following public methods should be treated as private.
1385    
1386 jgs 94 /**
1387     \brief
1388 jgs 102 Perform the given unary operation on every element of every data point in
1389     this Data object.
1390 jgs 94 */
1391 jgs 102 template <class UnaryFunction>
1392 woo409 757 ESCRIPT_DLL_API
1393 jgs 102 inline
1394 jgs 94 void
1395 matt 1334 unaryOp2(UnaryFunction operation);
1396 jgs 102
1397     /**
1398     \brief
1399     Return a Data object containing the specified slice of
1400     this Data object.
1401     \param region - Input - Region to copy.
1402 jgs 123 *
1403 jgs 94 */
1404 woo409 757 ESCRIPT_DLL_API
1405 jgs 102 Data
1406 jfenwick 1796 getSlice(const DataTypes::RegionType& region) const;
1407 jgs 94
1408 jgs 102 /**
1409     \brief
1410     Copy the specified slice from the given value into this
1411     Data object.
1412     \param value - Input - Data to copy from.
1413     \param region - Input - Region to copy.
1414 jgs 123 *
1415 jgs 102 */
1416 woo409 757 ESCRIPT_DLL_API
1417 jgs 102 void
1418     setSlice(const Data& value,
1419 jfenwick 1796 const DataTypes::RegionType& region);
1420 jgs 102
1421 jgs 119 /**
1422     \brief
1423 matt 1327 print the data values to stdout. Used for debugging
1424 bcumming 751 */
1425 bcumming 782 ESCRIPT_DLL_API
1426 matt 1327 void
1427 matt 1332 print(void);
1428 bcumming 751
1429 bcumming 782 /**
1430     \brief
1431     return the MPI rank number of the local data
1432 matt 1332 MPI_COMM_WORLD is assumed and the result of MPI_Comm_size()
1433     is returned
1434 bcumming 782 */
1435     ESCRIPT_DLL_API
1436 matt 1332 int
1437     get_MPIRank(void) const;
1438 bcumming 782
1439     /**
1440     \brief
1441     return the MPI rank number of the local data
1442 matt 1332 MPI_COMM_WORLD is assumed and the result of MPI_Comm_rank()
1443     is returned
1444 bcumming 782 */
1445     ESCRIPT_DLL_API
1446 matt 1332 int
1447     get_MPISize(void) const;
1448 bcumming 782
1449     /**
1450     \brief
1451     return the MPI rank number of the local data
1452 matt 1332 MPI_COMM_WORLD is assumed and returned.
1453 bcumming 782 */
1454     ESCRIPT_DLL_API
1455 matt 1332 MPI_Comm
1456     get_MPIComm(void) const;
1457 ksteube 813
1458     /**
1459     \brief
1460     return the object produced by the factory, which is a DataConstant or DataExpanded
1461 jfenwick 1796 TODO Ownership of this object should be explained in doco.
1462 ksteube 813 */
1463     ESCRIPT_DLL_API
1464 matt 1332 DataAbstract*
1465     borrowData(void) const;
1466 ksteube 813
1467 jfenwick 2005 ESCRIPT_DLL_API
1468     DataAbstract_ptr
1469     borrowDataPtr(void) const;
1470 jfenwick 1796
1471 jfenwick 2005 ESCRIPT_DLL_API
1472     DataReady_ptr
1473     borrowReadyPtr(void) const;
1474    
1475    
1476    
1477 jfenwick 1796 /**
1478     \brief
1479     Return a pointer to the beginning of the datapoint at the specified offset.
1480     TODO Eventually these should be inlined.
1481     \param i - position(offset) in the underlying datastructure
1482     */
1483 jfenwick 2271
1484 jfenwick 1796 ESCRIPT_DLL_API
1485     DataTypes::ValueType::const_reference
1486 jfenwick 2271 getDataAtOffsetRO(DataTypes::ValueType::size_type i);
1487 jfenwick 1796
1488    
1489     ESCRIPT_DLL_API
1490     DataTypes::ValueType::reference
1491 jfenwick 2271 getDataAtOffsetRW(DataTypes::ValueType::size_type i);
1492 jfenwick 1796
1493 jfenwick 2271
1494    
1495     /**
1496     \brief Create a buffer for use by getSample
1497     Allocates a DataVector large enough for DataLazy::resolveSample to operate on for the current Data.
1498     Do not use this buffer for other Data instances (unless you are sure they will be the same size).
1499    
1500     In multi-threaded sections, this needs to be called on each thread.
1501    
1502     \return A BufferGroup* if Data is lazy, NULL otherwise.
1503     \warning This pointer must be deallocated using freeSampleBuffer to avoid cross library memory issues.
1504     */
1505     ESCRIPT_DLL_API
1506     BufferGroup*
1507     allocSampleBuffer() const;
1508    
1509     /**
1510     \brief Free a buffer allocated with allocSampleBuffer.
1511     \param buffer Input - pointer to the buffer to deallocate.
1512     */
1513     ESCRIPT_DLL_API void freeSampleBuffer(BufferGroup* buffer);
1514    
1515 jgs 102 protected:
1516    
1517 jgs 94 private:
1518    
1519 jfenwick 2005 double
1520     LsupWorker() const;
1521    
1522     double
1523     supWorker() const;
1524    
1525     double
1526     infWorker() const;
1527    
1528 jfenwick 2271 boost::python::object
1529 jfenwick 2005 integrateWorker() const;
1530    
1531 jfenwick 2476 void
1532     calc_minGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1533    
1534     void
1535     calc_maxGlobalDataPoint(int& ProcNo, int& DataPointNo) const;
1536    
1537    
1538 jgs 94 /**
1539     \brief
1540 jgs 102 Check *this and the right operand are compatible. Throws
1541     an exception if they aren't.
1542     \param right - Input - The right hand side.
1543     */
1544     inline
1545     void
1546     operandCheck(const Data& right) const
1547     {
1548     return m_data->operandCheck(*(right.m_data.get()));
1549     }
1550    
1551     /**
1552     \brief
1553     Perform the specified reduction algorithm on every element of every data point in
1554 jgs 113 this Data object according to the given function and return the single value result.
1555 jgs 102 */
1556 jgs 147 template <class BinaryFunction>
1557 jgs 102 inline
1558     double
1559 jgs 147 algorithm(BinaryFunction operation,
1560     double initial_value) const;
1561 jgs 102
1562 jgs 113 /**
1563     \brief
1564     Reduce each data-point in this Data object using the given operation. Return a Data
1565     object with the same number of data-points, but with each data-point containing only
1566     one value - the result of the reduction operation on the corresponding data-point in
1567     this Data object
1568     */
1569 jgs 147 template <class BinaryFunction>
1570 jgs 106 inline
1571     Data
1572 jgs 147 dp_algorithm(BinaryFunction operation,
1573     double initial_value) const;
1574 jgs 106
1575 jgs 102 /**
1576     \brief
1577     Perform the given binary operation on all of the data's elements.
1578     The underlying type of the right hand side (right) determines the final
1579     type of *this after the operation. For example if the right hand side
1580     is expanded *this will be expanded if necessary.
1581     RHS is a Data object.
1582     */
1583     template <class BinaryFunction>
1584     inline
1585     void
1586     binaryOp(const Data& right,
1587     BinaryFunction operation);
1588    
1589     /**
1590     \brief
1591     Convert the data type of the RHS to match this.
1592     \param right - Input - data type to match.
1593     */
1594     void
1595     typeMatchLeft(Data& right) const;
1596    
1597     /**
1598     \brief
1599     Convert the data type of this to match the RHS.
1600     \param right - Input - data type to match.
1601     */
1602     void
1603     typeMatchRight(const Data& right);
1604    
1605     /**
1606     \brief
1607 jgs 94 Construct a Data object of the appropriate type.
1608     */
1609 jfenwick 1796
1610 jgs 94 void
1611 jfenwick 1796 initialise(const DataTypes::ValueType& value,
1612     const DataTypes::ShapeType& shape,
1613 jgs 94 const FunctionSpace& what,
1614     bool expanded);
1615    
1616 jfenwick 1796 void
1617 jfenwick 2271 initialise(const WrappedArray& value,
1618 jfenwick 1796 const FunctionSpace& what,
1619     bool expanded);
1620    
1621 jgs 94 //
1622 gross 783 // flag to protect the data object against any update
1623     bool m_protected;
1624 jfenwick 2271 mutable bool m_shared;
1625     bool m_lazy;
1626 gross 783
1627     //
1628 jgs 102 // pointer to the actual data object
1629 jfenwick 1872 // boost::shared_ptr<DataAbstract> m_data;
1630     DataAbstract_ptr m_data;
1631 jgs 94
1632 jfenwick 2271 // If possible please use getReadyPtr instead.
1633     // But see warning below.
1634 jfenwick 2005 const DataReady*
1635     getReady() const;
1636    
1637     DataReady*
1638     getReady();
1639    
1640 jfenwick 2271
1641     // Be wary of using this for local operations since it (temporarily) increases reference count.
1642     // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1643     // getReady() instead
1644 jfenwick 2005 DataReady_ptr
1645     getReadyPtr();
1646    
1647     const_DataReady_ptr
1648     getReadyPtr() const;
1649    
1650    
1651 jfenwick 2271 /**
1652     \brief Update the Data's shared flag
1653     This indicates that the DataAbstract used by this object is now shared (or no longer shared).
1654     For internal use only.
1655     */
1656     void updateShareStatus(bool nowshared) const
1657     {
1658     m_shared=nowshared; // m_shared is mutable
1659     }
1660    
1661     // In the isShared() method below:
1662     // A problem would occur if m_data (the address pointed to) were being modified
1663     // while the call m_data->is_shared is being executed.
1664     //
1665     // Q: So why do I think this code can be thread safe/correct?
1666     // A: We need to make some assumptions.
1667     // 1. We assume it is acceptable to return true under some conditions when we aren't shared.
1668     // 2. We assume that no constructions or assignments which will share previously unshared
1669     // will occur while this call is executing. This is consistent with the way Data:: and C are written.
1670     //
1671     // This means that the only transition we need to consider, is when a previously shared object is
1672     // not shared anymore. ie. the other objects have been destroyed or a deep copy has been made.
1673     // In those cases the m_shared flag changes to false after m_data has completed changing.
1674     // For any threads executing before the flag switches they will assume the object is still shared.
1675     bool isShared() const
1676     {
1677     return m_shared;
1678     /* if (m_shared) return true;
1679     if (m_data->isShared())
1680     {
1681     updateShareStatus(true);
1682     return true;
1683     }
1684     return false;*/
1685     }
1686    
1687     void forceResolve()
1688     {
1689     if (isLazy())
1690     {
1691     #ifdef _OPENMP
1692     if (omp_in_parallel())
1693     { // Yes this is throwing an exception out of an omp thread which is forbidden.
1694     throw DataException("Please do not call forceResolve() in a parallel region.");
1695     }
1696     #endif
1697     resolve();
1698     }
1699     }
1700    
1701     /**
1702     \brief if another object is sharing out member data make a copy to work with instead.
1703     This code should only be called from single threaded sections of code.
1704     */
1705     void exclusiveWrite()
1706     {
1707     #ifdef _OPENMP
1708     if (omp_in_parallel())
1709     {
1710     // *((int*)0)=17;
1711     throw DataException("Programming error. Please do not run exclusiveWrite() in multi-threaded sections.");
1712     }
1713     #endif
1714     forceResolve();
1715     if (isShared())
1716     {
1717     DataAbstract* t=m_data->deepCopy();
1718     set_m_data(DataAbstract_ptr(t));
1719     }
1720     }
1721    
1722     /**
1723     \brief checks if caller can have exclusive write to the object
1724     */
1725     void checkExclusiveWrite()
1726     {
1727     if (isLazy() || isShared())
1728     {
1729     throw DataException("Programming error. ExclusiveWrite required - please call requireWrite()");
1730     }
1731     }
1732    
1733     /**
1734     \brief Modify the data abstract hosted by this Data object
1735     For internal use only.
1736     Passing a pointer to null is permitted (do this in the destructor)
1737     \warning Only to be called in single threaded code or inside a single/critical section. This method needs to be atomic.
1738     */
1739     void set_m_data(DataAbstract_ptr p);
1740    
1741     friend class DataAbstract; // To allow calls to updateShareStatus
1742    
1743 jgs 94 };
1744    
1745 jfenwick 2005 } // end namespace escript
1746 jgs 94
1747 jfenwick 1796
1748 jfenwick 2005 // No, this is not supposed to be at the top of the file
1749     // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1750     // so that I can dynamic cast between them below.
1751     #include "DataReady.h"
1752 jfenwick 2271 #include "DataLazy.h"
1753 jfenwick 2005
1754     namespace escript
1755     {
1756    
1757     inline
1758     const DataReady*
1759     Data::getReady() const
1760     {
1761     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1762     EsysAssert((dr!=0), "Error - casting to DataReady.");
1763     return dr;
1764     }
1765    
1766     inline
1767     DataReady*
1768     Data::getReady()
1769     {
1770     DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1771     EsysAssert((dr!=0), "Error - casting to DataReady.");
1772     return dr;
1773     }
1774    
1775 jfenwick 2271 // Be wary of using this for local operations since it (temporarily) increases reference count.
1776     // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1777     // getReady() instead
1778 jfenwick 2005 inline
1779     DataReady_ptr
1780     Data::getReadyPtr()
1781     {
1782     DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1783     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1784     return dr;
1785     }
1786    
1787    
1788     inline
1789     const_DataReady_ptr
1790     Data::getReadyPtr() const
1791     {
1792     const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1793     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1794     return dr;
1795     }
1796    
1797     inline
1798     DataAbstract::ValueType::value_type*
1799 jfenwick 2271 Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1800 jfenwick 2005 {
1801     if (isLazy())
1802     {
1803 jfenwick 2271 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1804 jfenwick 2005 }
1805 jfenwick 2271 return getReady()->getSampleDataRW(sampleNo);
1806 jfenwick 2005 }
1807    
1808 jfenwick 2271 inline
1809     const DataAbstract::ValueType::value_type*
1810     Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo, BufferGroup* bufferg)
1811     {
1812     DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1813     if (l!=0)
1814     {
1815     size_t offset=0;
1816     if (bufferg==NULL)
1817     {
1818     throw DataException("Error, attempt to getSampleDataRO for lazy Data with buffer==NULL");
1819     }
1820     const DataTypes::ValueType* res=l->resolveSample(*bufferg,sampleNo,offset);
1821     return &((*res)[offset]);
1822     }
1823     return getReady()->getSampleDataRO(sampleNo);
1824     }
1825    
1826    
1827    
1828 jgs 102 /**
1829 ksteube 1748 Modify a filename for MPI parallel output to multiple files
1830     */
1831     char *Escript_MPI_appendRankToFileName(const char *, int, int);
1832    
1833     /**
1834 jgs 102 Binary Data object operators.
1835     */
1836 matt 1327 inline double rpow(double x,double y)
1837 gross 854 {
1838     return pow(y,x);
1839 gross 1028 }
1840 jgs 94
1841     /**
1842     \brief
1843     Operator+
1844     Takes two Data objects.
1845     */
1846 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
1847 jgs 94
1848     /**
1849     \brief
1850     Operator-
1851     Takes two Data objects.
1852     */
1853 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1854 jgs 94
1855     /**
1856     \brief
1857     Operator*
1858     Takes two Data objects.
1859     */
1860 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1861 jgs 94
1862     /**
1863     \brief
1864     Operator/
1865     Takes two Data objects.
1866     */
1867 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1868 jgs 94
1869     /**
1870     \brief
1871     Operator+
1872     Takes LHS Data object and RHS python::object.
1873     python::object must be convertable to Data type.
1874     */
1875 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1876 jgs 94
1877     /**
1878     \brief
1879     Operator-
1880     Takes LHS Data object and RHS python::object.
1881     python::object must be convertable to Data type.
1882     */
1883 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1884 jgs 94
1885     /**
1886     \brief
1887     Operator*
1888     Takes LHS Data object and RHS python::object.
1889     python::object must be convertable to Data type.
1890     */
1891 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1892 jgs 94
1893     /**
1894     \brief
1895     Operator/
1896     Takes LHS Data object and RHS python::object.
1897     python::object must be convertable to Data type.
1898     */
1899 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1900 jgs 94
1901     /**
1902     \brief
1903     Operator+
1904     Takes LHS python::object and RHS Data object.
1905     python::object must be convertable to Data type.
1906     */
1907 woo409 757 ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1908 jgs 94
1909     /**
1910     \brief
1911     Operator-
1912     Takes LHS python::object and RHS Data object.
1913     python::object must be convertable to Data type.
1914     */
1915 woo409 757 ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1916 jgs 94
1917     /**
1918     \brief
1919     Operator*
1920     Takes LHS python::object and RHS Data object.
1921     python::object must be convertable to Data type.
1922     */
1923 woo409 757 ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1924 jgs 94
1925     /**
1926     \brief
1927     Operator/
1928     Takes LHS python::object and RHS Data object.
1929     python::object must be convertable to Data type.
1930     */
1931 woo409 757 ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1932 jgs 94
1933 ksteube 1312
1934    
1935 jgs 94 /**
1936     \brief
1937     Output operator
1938     */
1939 woo409 757 ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1940 jgs 94
1941     /**
1942     \brief
1943 ksteube 813 Compute a tensor product of two Data objects
1944 jfenwick 2519 \param arg_0 - Input - Data object
1945     \param arg_1 - Input - Data object
1946 ksteube 813 \param axis_offset - Input - axis offset
1947     \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1948     */
1949     ESCRIPT_DLL_API
1950     Data
1951 jfenwick 2519 C_GeneralTensorProduct(Data& arg_0,
1952     Data& arg_1,
1953 ksteube 813 int axis_offset=0,
1954     int transpose=0);
1955    
1956 jgs 102 /**
1957     \brief
1958     Perform the given binary operation with this and right as operands.
1959     Right is a Data object.
1960     */
1961 jgs 94 template <class BinaryFunction>
1962     inline
1963     void
1964     Data::binaryOp(const Data& right,
1965     BinaryFunction operation)
1966     {
1967     //
1968     // if this has a rank of zero promote it to the rank of the RHS
1969 jfenwick 1796 if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1970 gross 854 throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1971 jgs 94 }
1972 jfenwick 2005
1973     if (isLazy() || right.isLazy())
1974     {
1975     throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1976     }
1977 jgs 94 //
1978     // initially make the temporary a shallow copy
1979     Data tempRight(right);
1980 jfenwick 2005
1981 jgs 94 if (getFunctionSpace()!=right.getFunctionSpace()) {
1982     if (right.probeInterpolation(getFunctionSpace())) {
1983     //
1984 matt 1327 // an interpolation is required so create a new Data
1985 jgs 94 tempRight=Data(right,this->getFunctionSpace());
1986     } else if (probeInterpolation(right.getFunctionSpace())) {
1987     //
1988     // interpolate onto the RHS function space
1989     Data tempLeft(*this,right.getFunctionSpace());
1990 jfenwick 2271 // m_data=tempLeft.m_data;
1991     set_m_data(tempLeft.m_data);
1992 jgs 94 }
1993     }
1994     operandCheck(tempRight);
1995     //
1996     // ensure this has the right type for the RHS
1997 jgs 102 typeMatchRight(tempRight);
1998 jgs 94 //
1999     // Need to cast to the concrete types so that the correct binaryOp
2000     // is called.
2001     if (isExpanded()) {
2002     //
2003     // Expanded data will be done in parallel, the right hand side can be
2004     // of any data type
2005     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2006     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2007 jfenwick 2005 escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2008 jgs 94 } else if (isTagged()) {
2009     //
2010     // Tagged data is operated on serially, the right hand side can be
2011     // either DataConstant or DataTagged
2012     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2013     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2014     if (right.isTagged()) {
2015     DataTagged* rightC=dynamic_cast<DataTagged*>(tempRight.m_data.get());
2016     EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
2017     escript::binaryOp(*leftC,*rightC,operation);
2018     } else {
2019     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2020     EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
2021     escript::binaryOp(*leftC,*rightC,operation);
2022     }
2023 jgs 102 } else if (isConstant()) {
2024 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2025     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2026     EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2027     escript::binaryOp(*leftC,*rightC,operation);
2028     }
2029     }
2030    
2031 jgs 102 /**
2032     \brief
2033     Perform the given Data object reduction algorithm on this and return the result.
2034     Given operation combines each element of each data point, thus argument
2035     object (*this) is a rank n Data object, and returned object is a scalar.
2036     Calls escript::algorithm.
2037     */
2038 jgs 147 template <class BinaryFunction>
2039 jgs 94 inline
2040     double
2041 jgs 147 Data::algorithm(BinaryFunction operation, double initial_value) const
2042 jgs 94 {
2043     if (isExpanded()) {
2044     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2045     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2046 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2047 jgs 102 } else if (isTagged()) {
2048 jgs 94 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2049     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2050 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2051 jgs 102 } else if (isConstant()) {
2052 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2053     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2054 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2055 jfenwick 1803 } else if (isEmpty()) {
2056     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2057 jfenwick 2005 } else if (isLazy()) {
2058     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2059     } else {
2060     throw DataException("Error - Data encapsulates an unknown type.");
2061 jgs 94 }
2062     }
2063    
2064 jgs 102 /**
2065     \brief
2066     Perform the given data point reduction algorithm on data and return the result.
2067     Given operation combines each element within each data point into a scalar,
2068 matt 1327 thus argument object is a rank n Data object, and returned object is a
2069 jgs 102 rank 0 Data object.
2070     Calls escript::dp_algorithm.
2071     */
2072 jgs 147 template <class BinaryFunction>
2073 jgs 94 inline
2074     Data
2075 jgs 147 Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2076 jgs 94 {
2077 jfenwick 1803 if (isEmpty()) {
2078     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2079     }
2080     else if (isExpanded()) {
2081 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2082 jgs 106 DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2083 jgs 102 DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2084     EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2085     EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2086 jgs 147 escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2087 jgs 559 return result;
2088 jfenwick 1803 }
2089     else if (isTagged()) {
2090 jgs 106 DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2091 jgs 102 EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2092 jfenwick 1796 DataTypes::ValueType defval(1);
2093     defval[0]=0;
2094     DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2095 jgs 147 escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2096 jfenwick 1796 return Data(resultT); // note: the Data object now owns the resultT pointer
2097 jfenwick 1803 }
2098     else if (isConstant()) {
2099 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2100 jgs 106 DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2101 jgs 102 DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2102     EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2103     EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2104 jgs 147 escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2105 jgs 559 return result;
2106 jfenwick 2005 } else if (isLazy()) {
2107     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2108     } else {
2109     throw DataException("Error - Data encapsulates an unknown type.");
2110 jgs 102 }
2111 jgs 94 }
2112    
2113 trankine 1430 /**
2114     \brief
2115     Compute a tensor operation with two Data objects
2116 jfenwick 2519 \param arg_0 - Input - Data object
2117     \param arg_1 - Input - Data object
2118 trankine 1430 \param operation - Input - Binary op functor
2119     */
2120 matt 1327 template <typename BinaryFunction>
2121 trankine 1430 inline
2122 matt 1327 Data
2123     C_TensorBinaryOperation(Data const &arg_0,
2124 matt 1332 Data const &arg_1,
2125     BinaryFunction operation)
2126 matt 1327 {
2127 jfenwick 1803 if (arg_0.isEmpty() || arg_1.isEmpty())
2128     {
2129     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2130     }
2131 jfenwick 2005 if (arg_0.isLazy() || arg_1.isLazy())
2132     {
2133     throw DataException("Error - Operations not permitted on lazy data.");
2134     }
2135 matt 1327 // Interpolate if necessary and find an appropriate function space
2136     Data arg_0_Z, arg_1_Z;
2137     if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2138     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2139     arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2140     arg_1_Z = Data(arg_1);
2141     }
2142     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2143     arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2144     arg_0_Z =Data(arg_0);
2145     }
2146     else {
2147     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2148     }
2149     } else {
2150     arg_0_Z = Data(arg_0);
2151     arg_1_Z = Data(arg_1);
2152     }
2153     // Get rank and shape of inputs
2154     int rank0 = arg_0_Z.getDataPointRank();
2155     int rank1 = arg_1_Z.getDataPointRank();
2156 jfenwick 1796 DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2157     DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2158 matt 1327 int size0 = arg_0_Z.getDataPointSize();
2159     int size1 = arg_1_Z.getDataPointSize();
2160     // Declare output Data object
2161     Data res;
2162    
2163     if (shape0 == shape1) {
2164     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2165 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2166 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2167     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2168     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2169 jfenwick 1796
2170 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2171     }
2172     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2173    
2174     // Prepare the DataConstant input
2175     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2176    
2177     // Borrow DataTagged input from Data object
2178     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2179    
2180     // Prepare a DataTagged output 2
2181 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2182 matt 1327 res.tag();
2183     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2184    
2185     // Prepare offset into DataConstant
2186     int offset_0 = tmp_0->getPointOffset(0,0);
2187 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2188 jfenwick 1796
2189 matt 1327 // Get the pointers to the actual data
2190 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2191     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2192 jfenwick 1796
2193 matt 1327 // Compute a result for the default
2194     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2195     // Compute a result for each tag
2196     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2197     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2198     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2199 jfenwick 1796 tmp_2->addTag(i->first);
2200 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2201     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2202 jfenwick 1796
2203 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2204 matt 1327 }
2205    
2206     }
2207     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2208     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2209     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2210     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2211     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2212    
2213     int sampleNo_1,dataPointNo_1;
2214     int numSamples_1 = arg_1_Z.getNumSamples();
2215     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2216     int offset_0 = tmp_0->getPointOffset(0,0);
2217 jfenwick 2271 res.requireWrite();
2218 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2219     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2220 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2221     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2222     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2223 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2224     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2225     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2226 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2227     }
2228 matt 1327 }
2229    
2230     }
2231     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2232     // Borrow DataTagged input from Data object
2233     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2234    
2235     // Prepare the DataConstant input
2236     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2237    
2238     // Prepare a DataTagged output 2
2239 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2240 matt 1327 res.tag();
2241     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2242    
2243     // Prepare offset into DataConstant
2244     int offset_1 = tmp_1->getPointOffset(0,0);
2245 jfenwick 2271
2246     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2247 matt 1327 // Get the pointers to the actual data
2248 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2249     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2250 matt 1327 // Compute a result for the default
2251     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2252     // Compute a result for each tag
2253     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2254     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2255     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2256 jfenwick 1796 tmp_2->addTag(i->first);
2257 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2258     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2259 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2260 matt 1327 }
2261    
2262     }
2263     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2264     // Borrow DataTagged input from Data object
2265     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2266    
2267     // Borrow DataTagged input from Data object
2268     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2269    
2270     // Prepare a DataTagged output 2
2271     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2272 matt 1332 res.tag(); // DataTagged output
2273 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2274    
2275     // Get the pointers to the actual data
2276 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2277     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2278     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2279 jfenwick 1796
2280 matt 1327 // Compute a result for the default
2281     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2282     // Merge the tags
2283     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2284     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2285     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2286     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2287 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2288 matt 1327 }
2289     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2290 jfenwick 1796 tmp_2->addTag(i->first);
2291 matt 1327 }
2292     // Compute a result for each tag
2293     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2294     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2295 jfenwick 1796
2296 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2297     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2298     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2299 jfenwick 1796
2300 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2301 matt 1327 }
2302    
2303     }
2304     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2305     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2306     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2307     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2308     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2309     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2310    
2311     int sampleNo_0,dataPointNo_0;
2312     int numSamples_0 = arg_0_Z.getNumSamples();
2313     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2314 jfenwick 2271 res.requireWrite();
2315 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2316     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2317 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2318 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2319 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2320     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2321     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2322 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2323     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2324 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2325     }
2326 matt 1327 }
2327    
2328     }
2329     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2330     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2331     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2332     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2333     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2334    
2335     int sampleNo_0,dataPointNo_0;
2336     int numSamples_0 = arg_0_Z.getNumSamples();
2337     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2338     int offset_1 = tmp_1->getPointOffset(0,0);
2339 jfenwick 2271 res.requireWrite();
2340 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2341     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2342 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2343     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2344     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2345 jfenwick 1796
2346 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2347     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2348     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2349 jfenwick 1796
2350    
2351 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2352     }
2353 matt 1327 }
2354    
2355     }
2356     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2357     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2358     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2359     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2360     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2361     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2362    
2363     int sampleNo_0,dataPointNo_0;
2364     int numSamples_0 = arg_0_Z.getNumSamples();
2365     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2366 jfenwick 2271 res.requireWrite();
2367 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2368     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2369 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2370 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2371 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2372     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2373     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2374 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2375     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2376 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2377     }
2378 matt 1327 }
2379    
2380     }
2381     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2382     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2383     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2384     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2385     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2386     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2387    
2388     int sampleNo_0,dataPointNo_0;
2389     int numSamples_0 = arg_0_Z.getNumSamples();
2390     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2391 jfenwick 2271 res.requireWrite();
2392 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2393     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2394 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2395     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2396     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2397     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2398 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2399     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2400     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2401 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2402     }
2403 matt 1327 }
2404    
2405     }
2406     else {
2407     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2408     }
2409    
2410     } else if (0 == rank0) {
2411     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2412 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataConstant output
2413 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2414     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2415     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2416 matt 1327 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2417     }
2418     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2419    
2420     // Prepare the DataConstant input
2421     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2422    
2423     // Borrow DataTagged input from Data object
2424     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2425    
2426     // Prepare a DataTagged output 2
2427 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataTagged output
2428 matt 1327 res.tag();
2429     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2430    
2431     // Prepare offset into DataConstant
2432     int offset_0 = tmp_0->getPointOffset(0,0);
2433 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2434 jfenwick 1796
2435 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2436     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2437    
2438 matt 1327 // Compute a result for the default
2439     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2440     // Compute a result for each tag
2441     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2442     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2443     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2444 jfenwick 1796 tmp_2->addTag(i->first);
2445 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2446     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2447 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2448 matt 1327 }
2449    
2450     }
2451     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2452    
2453     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2454     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2455     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2456     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2457    
2458     int sampleNo_1,dataPointNo_1;
2459     int numSamples_1 = arg_1_Z.getNumSamples();
2460     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2461     int offset_0 = tmp_0->getPointOffset(0,0);
2462 jfenwick 2271 res.requireWrite();
2463 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2464     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2465 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2466     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2467     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2468 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2469     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2470     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2471 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2472 matt 1327
2473 matt 1332 }
2474 matt 1327 }
2475    
2476     }
2477     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2478    
2479     // Borrow DataTagged input from Data object
2480     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2481    
2482     // Prepare the DataConstant input
2483     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2484    
2485     // Prepare a DataTagged output 2
2486 matt 1332 res = Data(0.0, shape1, arg_0_Z.getFunctionSpace()); // DataTagged output
2487 matt 1327 res.tag();
2488     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2489    
2490     // Prepare offset into DataConstant
2491     int offset_1 = tmp_1->getPointOffset(0,0);
2492 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2493 jfenwick 1796
2494     // Get the pointers to the actual data
2495 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2496     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2497 jfenwick 1796
2498    
2499 matt 1327 // Compute a result for the default
2500     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2501     // Compute a result for each tag
2502     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2503     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2504     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2505 jfenwick 1796 tmp_2->addTag(i->first);
2506 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2507     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2508 jfenwick 1796
2509 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2510 matt 1327 }
2511    
2512     }
2513     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2514    
2515     // Borrow DataTagged input from Data object
2516     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2517    
2518     // Borrow DataTagged input from Data object
2519     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2520    
2521     // Prepare a DataTagged output 2
2522     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2523 matt 1332 res.tag(); // DataTagged output
2524 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2525    
2526     // Get the pointers to the actual data
2527 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2528     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2529     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2530 jfenwick 1796
2531 matt 1327 // Compute a result for the default
2532     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2533     // Merge the tags
2534     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2535     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2536     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2537     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2538 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2539 matt 1327 }
2540     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2541 jfenwick 1796 tmp_2->addTag(i->first);
2542 matt 1327 }
2543     // Compute a result for each tag
2544     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2545     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2546 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2547     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2548     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2549 jfenwick 1796
2550 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2551 matt 1327 }
2552    
2553     }
2554     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2555    
2556     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2557     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2558     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2559     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2560     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2561    
2562     int sampleNo_0,dataPointNo_0;
2563     int numSamples_0 = arg_0_Z.getNumSamples();
2564     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2565 jfenwick 2271 res.requireWrite();
2566 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2567     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2568 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2569 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2570 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2571     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2572     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2573 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2574     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2575 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2576     }
2577 matt 1327 }
2578    
2579     }
2580     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2581     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2582     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2583     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2584     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2585    
2586     int sampleNo_0,dataPointNo_0;
2587     int numSamples_0 = arg_0_Z.getNumSamples();
2588     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2589     int offset_1 = tmp_1->getPointOffset(0,0);
2590 jfenwick 2271 res.requireWrite();
2591 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2592     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2593 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2594     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2595     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2596 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2597     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2598     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2599 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2600     }
2601 matt 1327 }
2602    
2603    
2604     }
2605     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2606    
2607     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2608     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2609     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2610     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2611     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2612    
2613     int sampleNo_0,dataPointNo_0;
2614     int numSamples_0 = arg_0_Z.getNumSamples();
2615     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2616 jfenwick 2271 res.requireWrite();
2617 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2618     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2619 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2620 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2621 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2622     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2623     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2624 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2625     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2626 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2627     }
2628 matt 1327 }
2629    
2630     }
2631     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2632    
2633     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2634     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2635     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2636     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2637     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2638    
2639     int sampleNo_0,dataPointNo_0;
2640     int numSamples_0 = arg_0_Z.getNumSamples();
2641     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2642 jfenwick 2271 res.requireWrite();
2643 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2644     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2645 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2646     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2647     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2648     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2649 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2650     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2651     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2652 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2653     }
2654 matt 1327 }
2655    
2656     }
2657     else {
2658     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2659     }
2660    
2661     } else if (0 == rank1) {
2662     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2663 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2664 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2665     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2666     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2667 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2668     }
2669     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2670    
2671     // Prepare the DataConstant input
2672     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2673    
2674     // Borrow DataTagged input from Data object
2675     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2676    
2677     // Prepare a DataTagged output 2
2678 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2679 matt 1327 res.tag();
2680     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2681    
2682     // Prepare offset into DataConstant
2683     int offset_0 = tmp_0->getPointOffset(0,0);
2684 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2685    
2686 jfenwick 1796 //Get the pointers to the actual data
2687 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2688     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2689 jfenwick 1796
2690 matt 1327 // Compute a result for the default
2691     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2692     // Compute a result for each tag
2693     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2694     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2695     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2696 jfenwick 1796 tmp_2->addTag(i->first);
2697 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2698     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2699 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2700 matt 1327 }
2701     }
2702     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2703    
2704     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2705     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2706     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2707     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2708    
2709     int sampleNo_1,dataPointNo_1;
2710     int numSamples_1 = arg_1_Z.getNumSamples();
2711     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2712     int offset_0 = tmp_0->getPointOffset(0,0);
2713 jfenwick 2271 res.requireWrite();
2714 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2715     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2716 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2717     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2718     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2719 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2720     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2721     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2722 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2723     }
2724 matt 1327 }
2725    
2726     }
2727     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2728    
2729     // Borrow DataTagged input from Data object
2730     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2731    
2732     // Prepare the DataConstant input
2733     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2734    
2735     // Prepare a DataTagged output 2
2736 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2737 matt 1327 res.tag();
2738     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2739    
2740     // Prepare offset into DataConstant
2741     int offset_1 = tmp_1->getPointOffset(0,0);
2742 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2743 matt 1327 // Get the pointers to the actual data
2744 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2745     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2746 matt 1327 // Compute a result for the default
2747     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2748     // Compute a result for each tag
2749     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2750     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2751     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2752 jfenwick 1796 tmp_2->addTag(i->first);
2753 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2754     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2755 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2756 matt 1327 }
2757    
2758     }
2759     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2760    
2761     // Borrow DataTagged input from Data object
2762     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2763    
2764     // Borrow DataTagged input from Data object
2765     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2766    
2767     // Prepare a DataTagged output 2
2768     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2769 matt 1332 res.tag(); // DataTagged output
2770 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2771    
2772     // Get the pointers to the actual data
2773 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2774     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2775     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2776 jfenwick 1796
2777 matt 1327 // Compute a result for the default
2778     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2779     // Merge the tags
2780     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2781     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2782     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2783     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2784 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2785 matt 1327 }
2786     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2787 jfenwick 1796 tmp_2->addTag(i->first);
2788 matt 1327 }
2789     // Compute a result for each tag
2790     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2791     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2792 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2793     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2794     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2795 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2796 matt 1327 }
2797    
2798     }
2799     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2800    
2801     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2802     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2803     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2804     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2805     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2806    
2807     int sampleNo_0,dataPointNo_0;
2808     int numSamples_0 = arg_0_Z.getNumSamples();
2809     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2810 jfenwick 2271 res.requireWrite();
2811 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2812     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2813 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2814 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2815 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2816     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2817     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2818 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2819     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2820 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2821     }
2822 matt 1327 }
2823    
2824     }
2825     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2826     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2827     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2828     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2829     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2830    
2831     int sampleNo_0,dataPointNo_0;
2832     int numSamples_0 = arg_0_Z.getNumSamples();
2833     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2834     int offset_1 = tmp_1->getPointOffset(0,0);
2835 jfenwick 2271 res.requireWrite();
2836 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2839     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2840     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2841 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2842     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2843     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2844 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2845     }
2846 matt 1327 }
2847    
2848    
2849     }
2850     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2851    
2852     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2853     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2856     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857    
2858     int sampleNo_0,dataPointNo_0;
2859     int numSamples_0 = arg_0_Z.getNumSamples();
2860     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2861 jfenwick 2271 res.requireWrite();
2862 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2863     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2864 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2865 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2866 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2867     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2868     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2869 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2870     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2871 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2872     }
2873 matt 1327 }
2874    
2875     }
2876     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2877    
2878     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2879     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2880     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2881     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2882     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2883    
2884     int sampleNo_0,dataPointNo_0;
2885     int numSamples_0 = arg_0_Z.getNumSamples();
2886     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2887 jfenwick 2271 res.requireWrite();
2888 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2889     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2890 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2891     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2892     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2893     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2894 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2895     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2896     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2897 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2898     }
2899 matt 1327 }
2900    
2901     }
2902     else {
2903     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2904     }
2905    
2906     } else {
2907     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2908     }
2909    
2910     return res;
2911 jgs 94 }
2912 matt 1327
2913 matt 1334 template <typename UnaryFunction>
2914     Data
2915     C_TensorUnaryOperation(Data const &arg_0,
2916     UnaryFunction operation)
2917     {
2918 jfenwick 1803 if (arg_0.isEmpty()) // do this before we attempt to interpolate
2919     {
2920     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2921     }
2922 jfenwick 2005 if (arg_0.isLazy())
2923     {
2924     throw DataException("Error - Operations not permitted on lazy data.");
2925     }
2926 matt 1334 // Interpolate if necessary and find an appropriate function space
2927     Data arg_0_Z = Data(arg_0);
2928    
2929     // Get rank and shape of inputs
2930 jfenwick 1796 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2931 matt 1334 int size0 = arg_0_Z.getDataPointSize();
2932    
2933     // Declare output Data object
2934     Data res;
2935    
2936     if (arg_0_Z.isConstant()) {
2937     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataConstant output
2938 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2939     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2940 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2941     }
2942     else if (arg_0_Z.isTagged()) {
2943    
2944     // Borrow DataTagged input from Data object
2945     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2946    
2947     // Prepare a DataTagged output 2
2948     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2949     res.tag();
2950     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2951    
2952     // Get the pointers to the actual data
2953 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2954     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2955 matt 1334 // Compute a result for the default
2956     tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2957     // Compute a result for each tag
2958     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2959     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2960     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2961 jfenwick 1796 tmp_2->addTag(i->first);
2962 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2963     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2964 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2965     }
2966    
2967     }
2968     else if (arg_0_Z.isExpanded()) {
2969    
2970     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2971     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2972     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2973    
2974     int sampleNo_0,dataPointNo_0;
2975     int numSamples_0 = arg_0_Z.getNumSamples();
2976     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2977     #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2978     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2979     for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2980     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2981     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2982 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2983     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2984 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2985     }
2986     }
2987     }
2988     else {
2989     throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
2990     }
2991    
2992     return res;
2993 matt 1327 }
2994 matt 1334
2995     }
2996 jgs 94 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26