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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2271 - (hide annotations)
Mon Feb 16 05:08:29 2009 UTC (11 years ago) by jfenwick
File MIME type: text/plain
File size: 93041 byte(s)
Merging version 2269 to trunk

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26