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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2771 - (hide annotations)
Wed Nov 25 01:38:49 2009 UTC (10 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 91818 byte(s)
Fix

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26