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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1748 - (hide annotations)
Wed Sep 3 06:10:39 2008 UTC (11 years, 1 month ago) by ksteube
Original Path: trunk/escript/src/Data.h
File MIME type: text/plain
File size: 83950 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26