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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2735 - (hide annotations)
Mon Nov 2 02:03:24 2009 UTC (10 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 92869 byte(s)
Fixed bug where calling minval() reintroduced laziness after it had been removed.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26