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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1693 - (hide annotations)
Fri Aug 8 04:22:58 2008 UTC (11 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 83936 byte(s)
Modified toString() on Data objects.
str(x) or print x will now print a summary of the data rather than all 
the points if the output would take more than 80 lines.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26