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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (10 years ago) by jfenwick
File MIME type: text/plain
File size: 92296 byte(s)
Don't panic.
Updating copyright stamps

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26