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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3390 - (hide annotations)
Thu Dec 2 00:34:37 2010 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 93400 byte(s)
RandomData added

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26