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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26