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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26