/[escript]/branches/schroedinger/escript/src/Data.h
ViewVC logotype

Annotation of /branches/schroedinger/escript/src/Data.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1943 - (hide annotations)
Wed Oct 29 04:05:14 2008 UTC (10 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 90197 byte(s)
Branch commit
Added interpolation.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26