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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2783 - (hide annotations)
Thu Nov 26 05:07:33 2009 UTC (10 years, 3 months ago) by lgao
File MIME type: text/plain
File size: 92080 byte(s)
process C_TensorBinaryOperation at the level of a whole sample rather than 
each datapoint when an Expanded data is operating with a Constant data. This
could improve the efficiency of the non-lazy version escipt. 


1 jgs 94
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 ksteube 1312
14 ksteube 1811
15 jgs 121 /** \file Data.h */
16 jgs 94
17     #ifndef DATA_H
18     #define DATA_H
19 woo409 757 #include "system_dep.h"
20 jgs 94
21 jfenwick 2271 #include "DataTypes.h"
22 jgs 474 #include "DataAbstract.h"
23     #include "DataAlgorithm.h"
24     #include "FunctionSpace.h"
25     #include "BinaryOp.h"
26     #include "UnaryOp.h"
27     #include "DataException.h"
28 jgs 94
29 jfenwick 2271
30 jgs 94 extern "C" {
31 jgs 474 #include "DataC.h"
32 jfenwick 2271 //#include <omp.h>
33 jgs 94 }
34    
35 jfenwick 2771 #ifdef _OPENMP
36     #include <omp.h>
37     #endif
38    
39 ksteube 1312 #include "esysmpi.h"
40 jgs 94 #include <string>
41     #include <algorithm>
42 jfenwick 1693 #include <sstream>
43 jgs 94
44     #include <boost/shared_ptr.hpp>
45     #include <boost/python/object.hpp>
46     #include <boost/python/tuple.hpp>
47    
48 jgs 121 namespace escript {
49    
50     //
51 jgs 149 // Forward declaration for various implementations of Data.
52 jgs 121 class DataConstant;
53     class DataTagged;
54     class DataExpanded;
55 jfenwick 2271 class DataLazy;
56 jgs 121
57 jgs 94 /**
58     \brief
59 jfenwick 1802 Data represents a collection of datapoints.
60 jgs 94
61     Description:
62 jfenwick 1802 Internally, the datapoints are actually stored by a DataAbstract object.
63     The specific instance of DataAbstract used may vary over the lifetime
64     of the Data object.
65     Some methods on this class return references (eg getShape()).
66     These references should not be used after an operation which changes the underlying DataAbstract object.
67     Doing so will lead to invalid memory access.
68     This should not affect any methods exposed via boost::python.
69 jgs 94 */
70     class Data {
71    
72     public:
73    
74 jgs 110 // These typedefs allow function names to be cast to pointers
75     // to functions of the appropriate type when calling unaryOp etc.
76 jgs 94 typedef double (*UnaryDFunPtr)(double);
77     typedef double (*BinaryDFunPtr)(double,double);
78    
79 gross 854
80 jgs 94 /**
81 jgs 102 Constructors.
82     */
83    
84     /**
85 jgs 94 \brief
86     Default constructor.
87     Creates a DataEmpty object.
88     */
89 woo409 757 ESCRIPT_DLL_API
90 jgs 94 Data();
91    
92     /**
93     \brief
94     Copy constructor.
95     WARNING: Only performs a shallow copy.
96     */
97 woo409 757 ESCRIPT_DLL_API
98 jgs 94 Data(const Data& inData);
99    
100     /**
101     \brief
102     Constructor from another Data object. If "what" is different from the
103 jgs 102 function space of inData the inData are tried to be interpolated to what,
104 jgs 94 otherwise a shallow copy of inData is returned.
105     */
106 woo409 757 ESCRIPT_DLL_API
107 jgs 94 Data(const Data& inData,
108     const FunctionSpace& what);
109    
110     /**
111 jfenwick 1796 \brief Copy Data from an existing vector
112     */
113 jgs 94
114 woo409 757 ESCRIPT_DLL_API
115 jfenwick 1796 Data(const DataTypes::ValueType& value,
116     const DataTypes::ShapeType& shape,
117     const FunctionSpace& what=FunctionSpace(),
118     bool expanded=false);
119 jgs 94
120     /**
121     \brief
122 jfenwick 2459 Constructor which creates a Data with points having the specified shape.
123 jgs 94
124     \param value - Input - Single value applied to all Data.
125     \param dataPointShape - Input - The shape of each data point.
126     \param what - Input - A description of what this data represents.
127     \param expanded - Input - Flag, if true fill the entire container with
128     the given value. Otherwise a more efficient storage
129     mechanism will be used.
130     */
131 woo409 757 ESCRIPT_DLL_API
132 jgs 94 Data(double value,
133 jfenwick 1796 const DataTypes::ShapeType& dataPointShape=DataTypes::ShapeType(),
134 jgs 94 const FunctionSpace& what=FunctionSpace(),
135     bool expanded=false);
136    
137     /**
138     \brief
139     Constructor which performs a deep copy of a region from another Data object.
140    
141     \param inData - Input - Input Data object.
142     \param region - Input - Region to copy.
143     */
144 woo409 757 ESCRIPT_DLL_API
145 jgs 94 Data(const Data& inData,
146 jfenwick 1796 const DataTypes::RegionType& region);
147 jgs 94
148     /**
149     \brief
150 jfenwick 2458 Constructor which copies data from any object that can be treated like a python array/sequence.
151 jgs 94
152     \param value - Input - Input data.
153     \param what - Input - A description of what this data represents.
154     \param expanded - Input - Flag, if true fill the entire container with
155     the value. Otherwise a more efficient storage
156     mechanism will be used.
157     */
158 woo409 757 ESCRIPT_DLL_API
159 jgs 94 Data(const boost::python::object& value,
160     const FunctionSpace& what=FunctionSpace(),
161     bool expanded=false);
162    
163     /**
164     \brief
165     Constructor which creates a DataConstant.
166 jfenwick 2458 Copies data from any object that can be treated like a python array/sequence.
167     All other parameters are copied from other.
168 jgs 94
169     \param value - Input - Input data.
170     \param other - Input - contains all other parameters.
171     */
172 woo409 757 ESCRIPT_DLL_API
173 jgs 94 Data(const boost::python::object& value,
174     const Data& other);
175    
176     /**
177     \brief
178     Constructor which creates a DataConstant of "shape" with constant value.
179     */
180 woo409 757 ESCRIPT_DLL_API
181 matt 1327 Data(double value,
182     const boost::python::tuple& shape=boost::python::make_tuple(),
183 jgs 94 const FunctionSpace& what=FunctionSpace(),
184     bool expanded=false);
185 jfenwick 1796
186    
187    
188 jgs 151 /**
189 jfenwick 1796 \brief Create a Data using an existing DataAbstract. Warning: The new object assumes ownership of the pointer!
190     Once you have passed the pointer, do not delete it.
191     */
192     ESCRIPT_DLL_API
193 jfenwick 2005 explicit Data(DataAbstract* underlyingdata);
194 jfenwick 1796
195 jfenwick 2005 /**
196     \brief Create a Data based on the supplied DataAbstract
197     */
198     ESCRIPT_DLL_API
199     explicit Data(DataAbstract_ptr underlyingdata);
200 jfenwick 1796
201     /**
202 jgs 151 \brief
203     Destructor
204     */
205 woo409 757 ESCRIPT_DLL_API
206 jgs 151 ~Data();
207 jgs 94
208     /**
209 jfenwick 1799 \brief Make this object a deep copy of "other".
210 jgs 94 */
211 woo409 757 ESCRIPT_DLL_API
212 jgs 94 void
213 jgs 102 copy(const Data& other);
214 jgs 94
215     /**
216 jfenwick 1799 \brief Return a pointer to a deep copy of this object.
217     */
218     ESCRIPT_DLL_API
219 jfenwick 2105 Data
220 jfenwick 1799 copySelf();
221    
222    
223 jfenwick 2005 /**
224     \brief produce a delayed evaluation version of this Data.
225     */
226     ESCRIPT_DLL_API
227     Data
228     delay();
229 jfenwick 1799
230 jfenwick 2005 /**
231     \brief convert the current data into lazy data.
232     */
233     ESCRIPT_DLL_API
234     void
235     delaySelf();
236 jfenwick 1799
237 jfenwick 2005
238 jfenwick 1799 /**
239 jgs 102 Member access methods.
240 jgs 94 */
241    
242     /**
243     \brief
244 matt 1327 switches on update protection
245 gross 783
246     */
247     ESCRIPT_DLL_API
248     void
249     setProtection();
250    
251     /**
252     \brief
253 jfenwick 2465 Returns true, if the data object is protected against update
254 gross 783
255     */
256     ESCRIPT_DLL_API
257     bool
258     isProtected() const;
259 gross 1034
260 jfenwick 2271
261 jfenwick 2455 /**
262     \brief
263 jfenwick 2465 Return the value of a data point as a python tuple.
264 jfenwick 2455 */
265 woo409 757 ESCRIPT_DLL_API
266 jfenwick 2271 const boost::python::object
267     getValueOfDataPointAsTuple(int dataPointNo);
268 jgs 117
269     /**
270     \brief
271 gross 1034 sets the values of a data-point from a python object on this process
272 jgs 121 */
273 woo409 757 ESCRIPT_DLL_API
274 gross 921 void
275 gross 1034 setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object);
276 jgs 121
277     /**
278     \brief
279 jfenwick 2458 sets the values of a data-point from a array-like object on this process
280 jgs 121 */
281 woo409 757 ESCRIPT_DLL_API
282 jgs 149 void
283 jfenwick 2271 setValueOfDataPointToArray(int dataPointNo, const boost::python::object&);
284 jgs 149
285     /**
286     \brief
287 gross 922 sets the values of a data-point on this process
288     */
289     ESCRIPT_DLL_API
290     void
291     setValueOfDataPoint(int dataPointNo, const double);
292    
293     /**
294 jfenwick 2455 \brief Return a data point across all processors as a python tuple.
295 gross 921 */
296     ESCRIPT_DLL_API
297 jfenwick 2271 const boost::python::object
298     getValueOfGlobalDataPointAsTuple(int procNo, int dataPointNo);
299    
300 gross 921 /**
301     \brief
302 jgs 149 Return the tag number associated with the given data-point.
303    
304     */
305 woo409 757 ESCRIPT_DLL_API
306 jgs 149 int
307     getTagNumber(int dpno);
308    
309     /**
310     \brief
311 jgs 94 Return the C wrapper for the Data object.
312     */
313 woo409 757 ESCRIPT_DLL_API
314 jgs 102 escriptDataC
315     getDataC();
316 jgs 94
317 jfenwick 1796
318    
319 jgs 94 /**
320     \brief
321     Return the C wrapper for the Data object - const version.
322     */
323 woo409 757 ESCRIPT_DLL_API
324 jgs 102 escriptDataC
325     getDataC() const;
326 jgs 94
327 jfenwick 2271
328     /**
329 jgs 94 \brief
330 jfenwick 1693 Write the data as a string. For large amounts of data, a summary is printed.
331 jgs 94 */
332 woo409 757 ESCRIPT_DLL_API
333 jgs 102 std::string
334 jfenwick 1693 toString() const;
335 jgs 94
336     /**
337     \brief
338 jgs 102 Whatever the current Data type make this into a DataExpanded.
339 jgs 94 */
340 woo409 757 ESCRIPT_DLL_API
341 jgs 94 void
342     expand();
343    
344     /**
345     \brief
346 jgs 102 If possible convert this Data to DataTagged. This will only allow
347 jgs 94 Constant data to be converted to tagged. An attempt to convert
348     Expanded data to tagged will throw an exception.
349     */
350 woo409 757 ESCRIPT_DLL_API
351 jgs 94 void
352     tag();
353    
354     /**
355 jfenwick 2005 \brief If this data is lazy, then convert it to ready data.
356     What type of ready data depends on the expression. For example, Constant+Tagged==Tagged.
357     */
358     ESCRIPT_DLL_API
359     void
360     resolve();
361    
362    
363     /**
364 jfenwick 2271 \brief Ensures data is ready for write access.
365     This means that the data will be resolved if lazy and will be copied if shared with another Data object.
366     \warning This method should only be called in single threaded sections of code. (It modifies m_data).
367     Do not create any Data objects from this one between calling requireWrite and getSampleDataRW.
368     Doing so might introduce additional sharing.
369     */
370     ESCRIPT_DLL_API
371     void
372     requireWrite();
373    
374     /**
375 jgs 94 \brief
376     Return true if this Data is expanded.
377 jfenwick 2271 \note To determine if a sample will contain separate values for each datapoint. Use actsExpanded instead.
378 jgs 94 */
379 woo409 757 ESCRIPT_DLL_API
380 jgs 94 bool
381     isExpanded() const;
382    
383     /**
384     \brief
385 jfenwick 2271 Return true if this Data is expanded or resolves to expanded.
386     That is, if it has a separate value for each datapoint in the sample.
387     */
388     ESCRIPT_DLL_API
389     bool
390     actsExpanded() const;
391    
392    
393     /**
394     \brief
395 jgs 94 Return true if this Data is tagged.
396     */
397 woo409 757 ESCRIPT_DLL_API
398 jgs 94 bool
399     isTagged() const;
400    
401     /**
402     \brief
403 jgs 102 Return true if this Data is constant.
404 jgs 94 */
405 woo409 757 ESCRIPT_DLL_API
406 jgs 94 bool
407 jgs 102 isConstant() const;
408 jgs 94
409     /**
410 jfenwick 2005 \brief Return true if this Data is lazy.
411     */
412     ESCRIPT_DLL_API
413     bool
414     isLazy() const;
415    
416     /**
417     \brief Return true if this data is ready.
418     */
419     ESCRIPT_DLL_API
420     bool
421     isReady() const;
422    
423     /**
424 jgs 94 \brief
425 jfenwick 1803 Return true if this Data holds an instance of DataEmpty. This is _not_ the same as asking if the object
426     contains datapoints.
427 jgs 94 */
428 woo409 757 ESCRIPT_DLL_API
429 jgs 94 bool
430 jgs 102 isEmpty() const;
431 jgs 94
432     /**
433     \brief
434     Return the function space.
435     */
436 woo409 757 ESCRIPT_DLL_API
437 jgs 94 inline
438     const FunctionSpace&
439     getFunctionSpace() const
440     {
441     return m_data->getFunctionSpace();
442     }
443    
444     /**
445     \brief
446     Return a copy of the function space.
447     */
448 woo409 757 ESCRIPT_DLL_API
449 jgs 94 const FunctionSpace
450     getCopyOfFunctionSpace() const;
451    
452     /**
453     \brief
454     Return the domain.
455     */
456 woo409 757 ESCRIPT_DLL_API
457 jgs 94 inline
458 jfenwick 1872 // const AbstractDomain&
459     const_Domain_ptr
460 jgs 94 getDomain() const
461     {
462     return getFunctionSpace().getDomain();
463     }
464    
465 jfenwick 1872
466 jgs 94 /**
467     \brief
468 jfenwick 1872 Return the domain.
469     TODO: For internal use only. This should be removed.
470     */
471     ESCRIPT_DLL_API
472     inline
473     // const AbstractDomain&
474     Domain_ptr
475     getDomainPython() const
476     {
477     return getFunctionSpace().getDomainPython();
478     }
479    
480     /**
481     \brief
482 jgs 94 Return a copy of the domain.
483     */
484 woo409 757 ESCRIPT_DLL_API
485 jgs 94 const AbstractDomain
486     getCopyOfDomain() const;
487    
488     /**
489     \brief
490     Return the rank of the point data.
491     */
492 woo409 757 ESCRIPT_DLL_API
493 jgs 94 inline
494 jfenwick 1952 unsigned int
495 jgs 94 getDataPointRank() const
496     {
497 jfenwick 1796 return m_data->getRank();
498 jgs 94 }
499    
500     /**
501     \brief
502 gross 921 Return the number of data points
503     */
504     ESCRIPT_DLL_API
505     inline
506     int
507     getNumDataPoints() const
508     {
509     return getNumSamples() * getNumDataPointsPerSample();
510     }
511     /**
512     \brief
513 jgs 102 Return the number of samples.
514 jgs 94 */
515 woo409 757 ESCRIPT_DLL_API
516 jgs 94 inline
517     int
518 jgs 102 getNumSamples() const
519 jgs 94 {
520 jgs 102 return m_data->getNumSamples();
521 jgs 94 }
522    
523     /**
524     \brief
525 jgs 102 Return the number of data points per sample.
526 jgs 94 */
527 woo409 757 ESCRIPT_DLL_API
528 jgs 102 inline
529 jgs 94 int
530 jgs 102 getNumDataPointsPerSample() const
531     {
532     return m_data->getNumDPPSample();
533     }
534 jfenwick 1796
535    
536 gross 950 /**
537 jfenwick 1796 \brief
538     Return the number of values in the shape for this object.
539     */
540     ESCRIPT_DLL_API
541     int
542     getNoValues() const
543     {
544     return m_data->getNoValues();
545     }
546    
547    
548     /**
549 gross 950 \brief
550 matt 1327 dumps the object into a netCDF file
551 gross 950 */
552     ESCRIPT_DLL_API
553     void
554 gross 1141 dump(const std::string fileName) const;
555 jfenwick 2005
556 jfenwick 2459 /**
557     \brief returns the values of the object as a list of tuples (one for each datapoint).
558 jfenwick 2271
559 jfenwick 2459 \param scalarastuple If true, scalar data will produce single valued tuples [(1,) (2,) ...]
560     If false, the result is a list of scalars [1, 2, ...]
561     */
562     ESCRIPT_DLL_API
563     const boost::python::object
564 jfenwick 2465 toListOfTuples(bool scalarastuple=true);
565 jfenwick 2271
566 jfenwick 2459
567 jfenwick 2271 /**
568     \brief
569     Return the sample data for the given sample no. This is not the
570     preferred interface but is provided for use by C code.
571     The bufferg parameter is only required for LazyData.
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    
1760 jgs 94 };
1761    
1762 jfenwick 2005 } // end namespace escript
1763 jgs 94
1764 jfenwick 1796
1765 jfenwick 2005 // No, this is not supposed to be at the top of the file
1766     // DataAbstact needs to be declared first, then DataReady needs to be fully declared
1767     // so that I can dynamic cast between them below.
1768     #include "DataReady.h"
1769 jfenwick 2271 #include "DataLazy.h"
1770 jfenwick 2005
1771     namespace escript
1772     {
1773    
1774     inline
1775     const DataReady*
1776     Data::getReady() const
1777     {
1778     const DataReady* dr=dynamic_cast<const DataReady*>(m_data.get());
1779     EsysAssert((dr!=0), "Error - casting to DataReady.");
1780     return dr;
1781     }
1782    
1783     inline
1784     DataReady*
1785     Data::getReady()
1786     {
1787     DataReady* dr=dynamic_cast<DataReady*>(m_data.get());
1788     EsysAssert((dr!=0), "Error - casting to DataReady.");
1789     return dr;
1790     }
1791    
1792 jfenwick 2271 // Be wary of using this for local operations since it (temporarily) increases reference count.
1793     // If you are just using this to call a method on DataReady instead of DataAbstract consider using
1794     // getReady() instead
1795 jfenwick 2005 inline
1796     DataReady_ptr
1797     Data::getReadyPtr()
1798     {
1799     DataReady_ptr dr=boost::dynamic_pointer_cast<DataReady>(m_data);
1800     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1801     return dr;
1802     }
1803    
1804    
1805     inline
1806     const_DataReady_ptr
1807     Data::getReadyPtr() const
1808     {
1809     const_DataReady_ptr dr=boost::dynamic_pointer_cast<const DataReady>(m_data);
1810     EsysAssert((dr.get()!=0), "Error - casting to DataReady.");
1811     return dr;
1812     }
1813    
1814     inline
1815     DataAbstract::ValueType::value_type*
1816 jfenwick 2271 Data::getSampleDataRW(DataAbstract::ValueType::size_type sampleNo)
1817 jfenwick 2005 {
1818     if (isLazy())
1819     {
1820 jfenwick 2271 throw DataException("Error, attempt to acquire RW access to lazy data. Please call requireWrite() first.");
1821 jfenwick 2005 }
1822 jfenwick 2271 return getReady()->getSampleDataRW(sampleNo);
1823 jfenwick 2005 }
1824    
1825 jfenwick 2271 inline
1826     const DataAbstract::ValueType::value_type*
1827 jfenwick 2770 Data::getSampleDataRO(DataAbstract::ValueType::size_type sampleNo)
1828 jfenwick 2271 {
1829     DataLazy* l=dynamic_cast<DataLazy*>(m_data.get());
1830     if (l!=0)
1831     {
1832     size_t offset=0;
1833 jfenwick 2770 const DataTypes::ValueType* res=l->resolveSample(sampleNo,offset);
1834 jfenwick 2271 return &((*res)[offset]);
1835     }
1836     return getReady()->getSampleDataRO(sampleNo);
1837     }
1838    
1839    
1840    
1841 jgs 102 /**
1842 ksteube 1748 Modify a filename for MPI parallel output to multiple files
1843     */
1844     char *Escript_MPI_appendRankToFileName(const char *, int, int);
1845    
1846     /**
1847 jgs 102 Binary Data object operators.
1848     */
1849 matt 1327 inline double rpow(double x,double y)
1850 gross 854 {
1851     return pow(y,x);
1852 gross 1028 }
1853 jgs 94
1854     /**
1855     \brief
1856     Operator+
1857     Takes two Data objects.
1858     */
1859 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const Data& right);
1860 jgs 94
1861     /**
1862     \brief
1863     Operator-
1864     Takes two Data objects.
1865     */
1866 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const Data& right);
1867 jgs 94
1868     /**
1869     \brief
1870     Operator*
1871     Takes two Data objects.
1872     */
1873 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const Data& right);
1874 jgs 94
1875     /**
1876     \brief
1877     Operator/
1878     Takes two Data objects.
1879     */
1880 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const Data& right);
1881 jgs 94
1882     /**
1883     \brief
1884     Operator+
1885     Takes LHS Data object and RHS python::object.
1886     python::object must be convertable to Data type.
1887     */
1888 woo409 757 ESCRIPT_DLL_API Data operator+(const Data& left, const boost::python::object& right);
1889 jgs 94
1890     /**
1891     \brief
1892     Operator-
1893     Takes LHS Data object and RHS python::object.
1894     python::object must be convertable to Data type.
1895     */
1896 woo409 757 ESCRIPT_DLL_API Data operator-(const Data& left, const boost::python::object& right);
1897 jgs 94
1898     /**
1899     \brief
1900     Operator*
1901     Takes LHS Data object and RHS python::object.
1902     python::object must be convertable to Data type.
1903     */
1904 woo409 757 ESCRIPT_DLL_API Data operator*(const Data& left, const boost::python::object& right);
1905 jgs 94
1906     /**
1907     \brief
1908     Operator/
1909     Takes LHS Data object and RHS python::object.
1910     python::object must be convertable to Data type.
1911     */
1912 woo409 757 ESCRIPT_DLL_API Data operator/(const Data& left, const boost::python::object& right);
1913 jgs 94
1914     /**
1915     \brief
1916     Operator+
1917     Takes LHS python::object and RHS Data object.
1918     python::object must be convertable to Data type.
1919     */
1920 woo409 757 ESCRIPT_DLL_API Data operator+(const boost::python::object& left, const Data& right);
1921 jgs 94
1922     /**
1923     \brief
1924     Operator-
1925     Takes LHS python::object and RHS Data object.
1926     python::object must be convertable to Data type.
1927     */
1928 woo409 757 ESCRIPT_DLL_API Data operator-(const boost::python::object& left, const Data& right);
1929 jgs 94
1930     /**
1931     \brief
1932     Operator*
1933     Takes LHS python::object and RHS Data object.
1934     python::object must be convertable to Data type.
1935     */
1936 woo409 757 ESCRIPT_DLL_API Data operator*(const boost::python::object& left, const Data& right);
1937 jgs 94
1938     /**
1939     \brief
1940     Operator/
1941     Takes LHS python::object and RHS Data object.
1942     python::object must be convertable to Data type.
1943     */
1944 woo409 757 ESCRIPT_DLL_API Data operator/(const boost::python::object& left, const Data& right);
1945 jgs 94
1946 ksteube 1312
1947    
1948 jgs 94 /**
1949     \brief
1950     Output operator
1951     */
1952 woo409 757 ESCRIPT_DLL_API std::ostream& operator<<(std::ostream& o, const Data& data);
1953 jgs 94
1954     /**
1955     \brief
1956 ksteube 813 Compute a tensor product of two Data objects
1957 jfenwick 2519 \param arg_0 - Input - Data object
1958     \param arg_1 - Input - Data object
1959 ksteube 813 \param axis_offset - Input - axis offset
1960     \param transpose - Input - 0: transpose neither, 1: transpose arg0, 2: transpose arg1
1961     */
1962     ESCRIPT_DLL_API
1963     Data
1964 jfenwick 2519 C_GeneralTensorProduct(Data& arg_0,
1965     Data& arg_1,
1966 ksteube 813 int axis_offset=0,
1967     int transpose=0);
1968    
1969 jgs 102 /**
1970     \brief
1971     Perform the given binary operation with this and right as operands.
1972     Right is a Data object.
1973     */
1974 jgs 94 template <class BinaryFunction>
1975     inline
1976     void
1977     Data::binaryOp(const Data& right,
1978     BinaryFunction operation)
1979     {
1980     //
1981     // if this has a rank of zero promote it to the rank of the RHS
1982 jfenwick 1796 if (getDataPointRank()==0 && right.getDataPointRank()!=0) {
1983 gross 854 throw DataException("Error - attempt to update rank zero object with object with rank bigger than zero.");
1984 jgs 94 }
1985 jfenwick 2005
1986     if (isLazy() || right.isLazy())
1987     {
1988     throw DataException("Programmer error - attempt to call binaryOp with Lazy Data.");
1989     }
1990 jgs 94 //
1991     // initially make the temporary a shallow copy
1992     Data tempRight(right);
1993 jfenwick 2005
1994 jgs 94 if (getFunctionSpace()!=right.getFunctionSpace()) {
1995     if (right.probeInterpolation(getFunctionSpace())) {
1996     //
1997 matt 1327 // an interpolation is required so create a new Data
1998 jgs 94 tempRight=Data(right,this->getFunctionSpace());
1999     } else if (probeInterpolation(right.getFunctionSpace())) {
2000     //
2001     // interpolate onto the RHS function space
2002     Data tempLeft(*this,right.getFunctionSpace());
2003 jfenwick 2271 // m_data=tempLeft.m_data;
2004     set_m_data(tempLeft.m_data);
2005 jgs 94 }
2006     }
2007     operandCheck(tempRight);
2008     //
2009     // ensure this has the right type for the RHS
2010 jgs 102 typeMatchRight(tempRight);
2011 jgs 94 //
2012     // Need to cast to the concrete types so that the correct binaryOp
2013     // is called.
2014     if (isExpanded()) {
2015     //
2016     // Expanded data will be done in parallel, the right hand side can be
2017     // of any data type
2018     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2019     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2020 jfenwick 2005 escript::binaryOp(*leftC,*(tempRight.getReady()),operation);
2021 jgs 94 } else if (isTagged()) {
2022     //
2023     // Tagged data is operated on serially, the right hand side can be
2024     // either DataConstant or DataTagged
2025     DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2026     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2027     if (right.isTagged()) {
2028     DataTagged* rightC=dynamic_cast<DataTagged*>(tempRight.m_data.get());
2029     EsysAssert((rightC!=0), "Programming error - casting to DataTagged.");
2030     escript::binaryOp(*leftC,*rightC,operation);
2031     } else {
2032     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2033     EsysAssert((rightC!=0), "Programming error - casting to DataConstant.");
2034     escript::binaryOp(*leftC,*rightC,operation);
2035     }
2036 jgs 102 } else if (isConstant()) {
2037 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2038     DataConstant* rightC=dynamic_cast<DataConstant*>(tempRight.m_data.get());
2039     EsysAssert((leftC!=0 && rightC!=0), "Programming error - casting to DataConstant.");
2040     escript::binaryOp(*leftC,*rightC,operation);
2041     }
2042     }
2043    
2044 jgs 102 /**
2045     \brief
2046     Perform the given Data object reduction algorithm on this and return the result.
2047     Given operation combines each element of each data point, thus argument
2048     object (*this) is a rank n Data object, and returned object is a scalar.
2049     Calls escript::algorithm.
2050     */
2051 jgs 147 template <class BinaryFunction>
2052 jgs 94 inline
2053     double
2054 jgs 147 Data::algorithm(BinaryFunction operation, double initial_value) const
2055 jgs 94 {
2056     if (isExpanded()) {
2057     DataExpanded* leftC=dynamic_cast<DataExpanded*>(m_data.get());
2058     EsysAssert((leftC!=0), "Programming error - casting to DataExpanded.");
2059 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2060 jgs 102 } else if (isTagged()) {
2061 jgs 94 DataTagged* leftC=dynamic_cast<DataTagged*>(m_data.get());
2062     EsysAssert((leftC!=0), "Programming error - casting to DataTagged.");
2063 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2064 jgs 102 } else if (isConstant()) {
2065 jgs 94 DataConstant* leftC=dynamic_cast<DataConstant*>(m_data.get());
2066     EsysAssert((leftC!=0), "Programming error - casting to DataConstant.");
2067 jgs 147 return escript::algorithm(*leftC,operation,initial_value);
2068 jfenwick 1803 } else if (isEmpty()) {
2069     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2070 jfenwick 2005 } else if (isLazy()) {
2071     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2072     } else {
2073     throw DataException("Error - Data encapsulates an unknown type.");
2074 jgs 94 }
2075     }
2076    
2077 jgs 102 /**
2078     \brief
2079     Perform the given data point reduction algorithm on data and return the result.
2080     Given operation combines each element within each data point into a scalar,
2081 matt 1327 thus argument object is a rank n Data object, and returned object is a
2082 jgs 102 rank 0 Data object.
2083     Calls escript::dp_algorithm.
2084     */
2085 jgs 147 template <class BinaryFunction>
2086 jgs 94 inline
2087     Data
2088 jgs 147 Data::dp_algorithm(BinaryFunction operation, double initial_value) const
2089 jgs 94 {
2090 jfenwick 1803 if (isEmpty()) {
2091     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2092     }
2093     else if (isExpanded()) {
2094 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2095 jgs 106 DataExpanded* dataE=dynamic_cast<DataExpanded*>(m_data.get());
2096 jgs 102 DataExpanded* resultE=dynamic_cast<DataExpanded*>(result.m_data.get());
2097     EsysAssert((dataE!=0), "Programming error - casting data to DataExpanded.");
2098     EsysAssert((resultE!=0), "Programming error - casting result to DataExpanded.");
2099 jgs 147 escript::dp_algorithm(*dataE,*resultE,operation,initial_value);
2100 jgs 559 return result;
2101 jfenwick 1803 }
2102     else if (isTagged()) {
2103 jgs 106 DataTagged* dataT=dynamic_cast<DataTagged*>(m_data.get());
2104 jgs 102 EsysAssert((dataT!=0), "Programming error - casting data to DataTagged.");
2105 jfenwick 1796 DataTypes::ValueType defval(1);
2106     defval[0]=0;
2107     DataTagged* resultT=new DataTagged(getFunctionSpace(), DataTypes::scalarShape, defval, dataT);
2108 jgs 147 escript::dp_algorithm(*dataT,*resultT,operation,initial_value);
2109 jfenwick 1796 return Data(resultT); // note: the Data object now owns the resultT pointer
2110 jfenwick 1803 }
2111     else if (isConstant()) {
2112 jfenwick 1796 Data result(0,DataTypes::ShapeType(),getFunctionSpace(),isExpanded());
2113 jgs 106 DataConstant* dataC=dynamic_cast<DataConstant*>(m_data.get());
2114 jgs 102 DataConstant* resultC=dynamic_cast<DataConstant*>(result.m_data.get());
2115     EsysAssert((dataC!=0), "Programming error - casting data to DataConstant.");
2116     EsysAssert((resultC!=0), "Programming error - casting result to DataConstant.");
2117 jgs 147 escript::dp_algorithm(*dataC,*resultC,operation,initial_value);
2118 jgs 559 return result;
2119 jfenwick 2005 } else if (isLazy()) {
2120     throw DataException("Error - Operations not permitted on instances of DataLazy.");
2121     } else {
2122     throw DataException("Error - Data encapsulates an unknown type.");
2123 jgs 102 }
2124 jgs 94 }
2125    
2126 trankine 1430 /**
2127     \brief
2128     Compute a tensor operation with two Data objects
2129 jfenwick 2519 \param arg_0 - Input - Data object
2130     \param arg_1 - Input - Data object
2131 trankine 1430 \param operation - Input - Binary op functor
2132     */
2133 matt 1327 template <typename BinaryFunction>
2134 trankine 1430 inline
2135 matt 1327 Data
2136     C_TensorBinaryOperation(Data const &arg_0,
2137 matt 1332 Data const &arg_1,
2138     BinaryFunction operation)
2139 matt 1327 {
2140 jfenwick 1803 if (arg_0.isEmpty() || arg_1.isEmpty())
2141     {
2142     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2143     }
2144 jfenwick 2005 if (arg_0.isLazy() || arg_1.isLazy())
2145     {
2146     throw DataException("Error - Operations not permitted on lazy data.");
2147     }
2148 matt 1327 // Interpolate if necessary and find an appropriate function space
2149     Data arg_0_Z, arg_1_Z;
2150     if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2151     if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2152     arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2153     arg_1_Z = Data(arg_1);
2154     }
2155     else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2156     arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2157     arg_0_Z =Data(arg_0);
2158     }
2159     else {
2160     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible function spaces.");
2161     }
2162     } else {
2163     arg_0_Z = Data(arg_0);
2164     arg_1_Z = Data(arg_1);
2165     }
2166     // Get rank and shape of inputs
2167     int rank0 = arg_0_Z.getDataPointRank();
2168     int rank1 = arg_1_Z.getDataPointRank();
2169 jfenwick 1796 DataTypes::ShapeType shape0 = arg_0_Z.getDataPointShape();
2170     DataTypes::ShapeType shape1 = arg_1_Z.getDataPointShape();
2171 matt 1327 int size0 = arg_0_Z.getDataPointSize();
2172     int size1 = arg_1_Z.getDataPointSize();
2173     // Declare output Data object
2174     Data res;
2175    
2176     if (shape0 == shape1) {
2177     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2178 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2179 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2180     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2181     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2182 jfenwick 1796
2183 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2184     }
2185     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2186    
2187     // Prepare the DataConstant input
2188     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2189    
2190     // Borrow DataTagged input from Data object
2191     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2192    
2193     // Prepare a DataTagged output 2
2194 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2195 matt 1327 res.tag();
2196     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2197    
2198     // Prepare offset into DataConstant
2199     int offset_0 = tmp_0->getPointOffset(0,0);
2200 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2201 jfenwick 1796
2202 matt 1327 // Get the pointers to the actual data
2203 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2204     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2205 jfenwick 1796
2206 matt 1327 // Compute a result for the default
2207     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2208     // Compute a result for each tag
2209     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2210     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2211     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2212 jfenwick 1796 tmp_2->addTag(i->first);
2213 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2214     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2215 jfenwick 1796
2216 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2217 matt 1327 }
2218    
2219     }
2220     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2221     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2222     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2223     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2224     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2225    
2226     int sampleNo_1,dataPointNo_1;
2227     int numSamples_1 = arg_1_Z.getNumSamples();
2228     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2229     int offset_0 = tmp_0->getPointOffset(0,0);
2230 jfenwick 2271 res.requireWrite();
2231 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2232     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2233 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2234     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2235     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2236 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2237     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2238     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2239 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2240     }
2241 matt 1327 }
2242    
2243     }
2244     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2245     // Borrow DataTagged input from Data object
2246     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2247    
2248     // Prepare the DataConstant input
2249     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2250    
2251     // Prepare a DataTagged output 2
2252 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2253 matt 1327 res.tag();
2254     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2255    
2256     // Prepare offset into DataConstant
2257     int offset_1 = tmp_1->getPointOffset(0,0);
2258 jfenwick 2271
2259     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2260 matt 1327 // Get the pointers to the actual data
2261 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2262     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2263 matt 1327 // Compute a result for the default
2264     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2265     // Compute a result for each tag
2266     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2267     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2268     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2269 jfenwick 1796 tmp_2->addTag(i->first);
2270 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2271     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2272 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2273 matt 1327 }
2274    
2275     }
2276     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2277     // Borrow DataTagged input from Data object
2278     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2279    
2280     // Borrow DataTagged input from Data object
2281     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2282    
2283     // Prepare a DataTagged output 2
2284     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2285 matt 1332 res.tag(); // DataTagged output
2286 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2287    
2288     // Get the pointers to the actual data
2289 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2290     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2291     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2292 jfenwick 1796
2293 matt 1327 // Compute a result for the default
2294     tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2295     // Merge the tags
2296     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2297     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2298     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2299     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2300 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2301 matt 1327 }
2302     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2303 jfenwick 1796 tmp_2->addTag(i->first);
2304 matt 1327 }
2305     // Compute a result for each tag
2306     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2307     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2308 jfenwick 1796
2309 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2310     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2311     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2312 jfenwick 1796
2313 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2314 matt 1327 }
2315    
2316     }
2317     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2318     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2319     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2320     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2321     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2322     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2323    
2324     int sampleNo_0,dataPointNo_0;
2325     int numSamples_0 = arg_0_Z.getNumSamples();
2326     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2327 jfenwick 2271 res.requireWrite();
2328 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2329     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2330 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2331 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2332 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2333     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2334     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2335 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2336     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2337 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2338     }
2339 matt 1327 }
2340    
2341     }
2342     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2343     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2344     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2345     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2346     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2347    
2348     int sampleNo_0,dataPointNo_0;
2349     int numSamples_0 = arg_0_Z.getNumSamples();
2350     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2351     int offset_1 = tmp_1->getPointOffset(0,0);
2352 jfenwick 2271 res.requireWrite();
2353 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2354     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2355 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2356     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2357     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2358 jfenwick 1796
2359 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2360     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2361     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2362 jfenwick 1796
2363    
2364 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2365     }
2366 matt 1327 }
2367    
2368     }
2369     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2370     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2371     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2372     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2373     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2374     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2375    
2376     int sampleNo_0,dataPointNo_0;
2377     int numSamples_0 = arg_0_Z.getNumSamples();
2378     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2379 jfenwick 2271 res.requireWrite();
2380 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2381     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2382 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2383 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2384 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2385     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2386     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2387 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2388     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2389 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1, ptr_2, operation);
2390     }
2391 matt 1327 }
2392    
2393     }
2394     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2395     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2396     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2397     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2398     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2399     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2400    
2401     int sampleNo_0,dataPointNo_0;
2402     int numSamples_0 = arg_0_Z.getNumSamples();
2403     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2404 jfenwick 2271 res.requireWrite();
2405 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2406     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2407 jfenwick 2716 dataPointNo_0=0;
2408     // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2409 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2410     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2411     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2412 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2413     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2414     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2415 jfenwick 2716 tensor_binary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_1, ptr_2, operation);
2416     // }
2417 matt 1327 }
2418    
2419     }
2420     else {
2421     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2422     }
2423    
2424     } else if (0 == rank0) {
2425     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2426 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataConstant output
2427 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2428     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2429     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2430 matt 1327 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2431     }
2432     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2433    
2434     // Prepare the DataConstant input
2435     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2436    
2437     // Borrow DataTagged input from Data object
2438     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2439    
2440     // Prepare a DataTagged output 2
2441 matt 1332 res = Data(0.0, shape1, arg_1_Z.getFunctionSpace()); // DataTagged output
2442 matt 1327 res.tag();
2443     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2444    
2445     // Prepare offset into DataConstant
2446     int offset_0 = tmp_0->getPointOffset(0,0);
2447 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2448 jfenwick 1796
2449 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2450     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2451    
2452 matt 1327 // Compute a result for the default
2453     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2454     // Compute a result for each tag
2455     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2456     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2457     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2458 jfenwick 1796 tmp_2->addTag(i->first);
2459 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2460     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2461 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2462 matt 1327 }
2463    
2464     }
2465     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2466    
2467     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2468     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2469     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2470     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2471    
2472     int sampleNo_1,dataPointNo_1;
2473     int numSamples_1 = arg_1_Z.getNumSamples();
2474     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2475     int offset_0 = tmp_0->getPointOffset(0,0);
2476 lgao 2783 const double *ptr_src = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2477     double ptr_0 = ptr_src[0];
2478     int size = size1*numDataPointsPerSample_1;
2479 jfenwick 2271 res.requireWrite();
2480 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2481     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2482 lgao 2783 // for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2483     int offset_1 = tmp_1->getPointOffset(sampleNo_1,0);
2484     int offset_2 = tmp_2->getPointOffset(sampleNo_1,0);
2485     // const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2486 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2487     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2488 lgao 2783 tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
2489 matt 1327
2490 lgao 2783 // }
2491 matt 1327 }
2492    
2493     }
2494     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2495    
2496     // Borrow DataTagged input from Data object
2497     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2498    
2499     // Prepare the DataConstant input
2500     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2501    
2502     // Prepare a DataTagged output 2
2503 matt 1332 res = Data(0.0, shape1, arg_0_Z.getFunctionSpace()); // DataTagged output
2504 matt 1327 res.tag();
2505     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2506    
2507     // Prepare offset into DataConstant
2508     int offset_1 = tmp_1->getPointOffset(0,0);
2509 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2510 jfenwick 1796
2511     // Get the pointers to the actual data
2512 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2513     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2514 jfenwick 1796
2515    
2516 matt 1327 // Compute a result for the default
2517     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2518     // Compute a result for each tag
2519     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2520     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2521     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2522 jfenwick 1796 tmp_2->addTag(i->first);
2523 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2524     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2525 jfenwick 1796
2526 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2527 matt 1327 }
2528    
2529     }
2530     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2531    
2532     // Borrow DataTagged input from Data object
2533     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2534    
2535     // Borrow DataTagged input from Data object
2536     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2537    
2538     // Prepare a DataTagged output 2
2539     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace());
2540 matt 1332 res.tag(); // DataTagged output
2541 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2542    
2543     // Get the pointers to the actual data
2544 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2545     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2546     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2547 jfenwick 1796
2548 matt 1327 // Compute a result for the default
2549     tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2550     // Merge the tags
2551     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2552     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2553     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2554     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2555 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2556 matt 1327 }
2557     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2558 jfenwick 1796 tmp_2->addTag(i->first);
2559 matt 1327 }
2560     // Compute a result for each tag
2561     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2562     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2563 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2564     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2565     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2566 jfenwick 1796
2567 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2568 matt 1327 }
2569    
2570     }
2571     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2572    
2573     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2574     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2575     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2576     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2577     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2578    
2579     int sampleNo_0,dataPointNo_0;
2580     int numSamples_0 = arg_0_Z.getNumSamples();
2581     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2582 jfenwick 2271 res.requireWrite();
2583 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2584     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2585 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2586 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2587 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2588     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2589     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2590 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2591     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2592 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2593     }
2594 matt 1327 }
2595    
2596     }
2597     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2598     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2599     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2600     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2601     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2602    
2603     int sampleNo_0,dataPointNo_0;
2604     int numSamples_0 = arg_0_Z.getNumSamples();
2605     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2606     int offset_1 = tmp_1->getPointOffset(0,0);
2607 jfenwick 2271 res.requireWrite();
2608 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2609     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2610 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2611     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2612     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2613 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2614     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2615     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2616 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2617     }
2618 matt 1327 }
2619    
2620    
2621     }
2622     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2623    
2624     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2625     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2626     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2627     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2628     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2629    
2630     int sampleNo_0,dataPointNo_0;
2631     int numSamples_0 = arg_0_Z.getNumSamples();
2632     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2633 jfenwick 2271 res.requireWrite();
2634 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2635     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2636 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2637 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2638 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2639     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2640     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2641 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2642     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2643 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2644     }
2645 matt 1327 }
2646    
2647     }
2648     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2649    
2650     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2651     res = Data(0.0, shape1, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2652     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2653     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2654     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2655    
2656     int sampleNo_0,dataPointNo_0;
2657     int numSamples_0 = arg_0_Z.getNumSamples();
2658     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2659 jfenwick 2271 res.requireWrite();
2660 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2661     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2662 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2663     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2664     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2665     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2666 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2667     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2668     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2669 matt 1332 tensor_binary_operation(size1, ptr_0[0], ptr_1, ptr_2, operation);
2670     }
2671 matt 1327 }
2672    
2673     }
2674     else {
2675     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2676     }
2677    
2678     } else if (0 == rank1) {
2679     if (arg_0_Z.isConstant() && arg_1_Z.isConstant()) {
2680 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataConstant output
2681 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2682     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(0));
2683     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2684 matt 1327 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2685     }
2686     else if (arg_0_Z.isConstant() && arg_1_Z.isTagged()) {
2687    
2688     // Prepare the DataConstant input
2689     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2690    
2691     // Borrow DataTagged input from Data object
2692     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2693    
2694     // Prepare a DataTagged output 2
2695 matt 1332 res = Data(0.0, shape0, arg_1_Z.getFunctionSpace()); // DataTagged output
2696 matt 1327 res.tag();
2697     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2698    
2699     // Prepare offset into DataConstant
2700     int offset_0 = tmp_0->getPointOffset(0,0);
2701 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2702    
2703 jfenwick 1796 //Get the pointers to the actual data
2704 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2705     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2706 jfenwick 1796
2707 matt 1327 // Compute a result for the default
2708     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2709     // Compute a result for each tag
2710     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2711     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2712     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2713 jfenwick 1796 tmp_2->addTag(i->first);
2714 jfenwick 2271 const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2715     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2716 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2717 matt 1327 }
2718     }
2719     else if (arg_0_Z.isConstant() && arg_1_Z.isExpanded()) {
2720    
2721     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2722     DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2723     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2724     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2725    
2726     int sampleNo_1,dataPointNo_1;
2727     int numSamples_1 = arg_1_Z.getNumSamples();
2728     int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2729     int offset_0 = tmp_0->getPointOffset(0,0);
2730 jfenwick 2271 res.requireWrite();
2731 matt 1327 #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2732     for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2733 matt 1332 for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2734     int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2735     int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2736 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2737     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2738     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2739 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2740     }
2741 matt 1327 }
2742    
2743     }
2744     else if (arg_0_Z.isTagged() && arg_1_Z.isConstant()) {
2745    
2746     // Borrow DataTagged input from Data object
2747     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2748    
2749     // Prepare the DataConstant input
2750     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2751    
2752     // Prepare a DataTagged output 2
2753 matt 1332 res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2754 matt 1327 res.tag();
2755     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2756    
2757     // Prepare offset into DataConstant
2758     int offset_1 = tmp_1->getPointOffset(0,0);
2759 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2760 matt 1327 // Get the pointers to the actual data
2761 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2762     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2763 matt 1327 // Compute a result for the default
2764     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2765     // Compute a result for each tag
2766     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2767     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2768     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2769 jfenwick 1796 tmp_2->addTag(i->first);
2770 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2771     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2772 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2773 matt 1327 }
2774    
2775     }
2776     else if (arg_0_Z.isTagged() && arg_1_Z.isTagged()) {
2777    
2778     // Borrow DataTagged input from Data object
2779     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2780    
2781     // Borrow DataTagged input from Data object
2782     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2783    
2784     // Prepare a DataTagged output 2
2785     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace());
2786 matt 1332 res.tag(); // DataTagged output
2787 matt 1327 DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2788    
2789     // Get the pointers to the actual data
2790 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2791     const double *ptr_1 = &(tmp_1->getDefaultValueRO(0));
2792     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2793 jfenwick 1796
2794 matt 1327 // Compute a result for the default
2795     tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2796     // Merge the tags
2797     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2798     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2799     const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2800     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2801 jfenwick 1796 tmp_2->addTag(i->first); // use tmp_2 to get correct shape
2802 matt 1327 }
2803     for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2804 jfenwick 1796 tmp_2->addTag(i->first);
2805 matt 1327 }
2806     // Compute a result for each tag
2807     const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2808     for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2809 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2810     const double *ptr_1 = &(tmp_1->getDataByTagRO(i->first,0));
2811     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2812 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2813 matt 1327 }
2814    
2815     }
2816     else if (arg_0_Z.isTagged() && arg_1_Z.isExpanded()) {
2817    
2818     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2819     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2820     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2821     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2822     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2823    
2824     int sampleNo_0,dataPointNo_0;
2825     int numSamples_0 = arg_0_Z.getNumSamples();
2826     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2827 jfenwick 2271 res.requireWrite();
2828 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2829     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2830 matt 1332 int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2831 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2832 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2833     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2834     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2835 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2836     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2837 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2838     }
2839 matt 1327 }
2840    
2841     }
2842     else if (arg_0_Z.isExpanded() && arg_1_Z.isConstant()) {
2843     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2844     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2845     DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2846     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2847    
2848     int sampleNo_0,dataPointNo_0;
2849     int numSamples_0 = arg_0_Z.getNumSamples();
2850     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2851     int offset_1 = tmp_1->getPointOffset(0,0);
2852 lgao 2783 const double *ptr_src = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2853     double ptr_1 = ptr_src[0];
2854     int size = size0 * numDataPointsPerSample_0;
2855 jfenwick 2271 res.requireWrite();
2856 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2857     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2858 lgao 2783 // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2859     int offset_0 = tmp_0->getPointOffset(sampleNo_0,0);
2860     int offset_2 = tmp_2->getPointOffset(sampleNo_0,0);
2861 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2862 lgao 2783 // const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2863 jfenwick 2271 double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2864 lgao 2783 tensor_binary_operation(size, ptr_0, ptr_1, ptr_2, operation);
2865     // }
2866 matt 1327 }
2867    
2868    
2869     }
2870     else if (arg_0_Z.isExpanded() && arg_1_Z.isTagged()) {
2871    
2872     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2873     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2874     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2875     DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2876     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2877    
2878     int sampleNo_0,dataPointNo_0;
2879     int numSamples_0 = arg_0_Z.getNumSamples();
2880     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2881 jfenwick 2271 res.requireWrite();
2882 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2883     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2884 matt 1332 int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2885 jfenwick 2271 const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2886 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2887     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2888     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2889 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2890     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2891 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2892     }
2893 matt 1327 }
2894    
2895     }
2896     else if (arg_0_Z.isExpanded() && arg_1_Z.isExpanded()) {
2897    
2898     // After finding a common function space above the two inputs have the same numSamples and num DPPS
2899     res = Data(0.0, shape0, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2900     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2901     DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2902     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2903    
2904     int sampleNo_0,dataPointNo_0;
2905     int numSamples_0 = arg_0_Z.getNumSamples();
2906     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2907 jfenwick 2271 res.requireWrite();
2908 matt 1327 #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2909     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2910 matt 1332 for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2911     int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2912     int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2913     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2914 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
2915     const double *ptr_1 = &(arg_1_Z.getDataAtOffsetRO(offset_1));
2916     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
2917 matt 1332 tensor_binary_operation(size0, ptr_0, ptr_1[0], ptr_2, operation);
2918     }
2919 matt 1327 }
2920    
2921     }
2922     else {
2923     throw DataException("Error - C_TensorBinaryOperation: unknown combination of inputs");
2924     }
2925    
2926     } else {
2927     throw DataException("Error - C_TensorBinaryOperation: arguments have incompatible shapes");
2928     }
2929    
2930     return res;
2931 jgs 94 }
2932 matt 1327
2933 matt 1334 template <typename UnaryFunction>
2934     Data
2935     C_TensorUnaryOperation(Data const &arg_0,
2936     UnaryFunction operation)
2937     {
2938 jfenwick 1803 if (arg_0.isEmpty()) // do this before we attempt to interpolate
2939     {
2940     throw DataException("Error - Operations not permitted on instances of DataEmpty.");
2941     }
2942 jfenwick 2005 if (arg_0.isLazy())
2943     {
2944     throw DataException("Error - Operations not permitted on lazy data.");
2945     }
2946 matt 1334 // Interpolate if necessary and find an appropriate function space
2947     Data arg_0_Z = Data(arg_0);
2948    
2949     // Get rank and shape of inputs
2950 jfenwick 1796 const DataTypes::ShapeType& shape0 = arg_0_Z.getDataPointShape();
2951 matt 1334 int size0 = arg_0_Z.getDataPointSize();
2952    
2953     // Declare output Data object
2954     Data res;
2955    
2956     if (arg_0_Z.isConstant()) {
2957     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataConstant output
2958 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(0));
2959     double *ptr_2 = &(res.getDataAtOffsetRW(0));
2960 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2961     }
2962     else if (arg_0_Z.isTagged()) {
2963    
2964     // Borrow DataTagged input from Data object
2965     DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2966    
2967     // Prepare a DataTagged output 2
2968     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace()); // DataTagged output
2969     res.tag();
2970     DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2971    
2972     // Get the pointers to the actual data
2973 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDefaultValueRO(0));
2974     double *ptr_2 = &(tmp_2->getDefaultValueRW(0));
2975 matt 1334 // Compute a result for the default
2976     tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2977     // Compute a result for each tag
2978     const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2979     DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2980     for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2981 jfenwick 1796 tmp_2->addTag(i->first);
2982 jfenwick 2271 const double *ptr_0 = &(tmp_0->getDataByTagRO(i->first,0));
2983     double *ptr_2 = &(tmp_2->getDataByTagRW(i->first,0));
2984 matt 1334 tensor_unary_operation(size0, ptr_0, ptr_2, operation);
2985     }
2986    
2987     }
2988     else if (arg_0_Z.isExpanded()) {
2989    
2990     res = Data(0.0, shape0, arg_0_Z.getFunctionSpace(),true); // DataExpanded output
2991     DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2992     DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2993    
2994     int sampleNo_0,dataPointNo_0;
2995     int numSamples_0 = arg_0_Z.getNumSamples();
2996     int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2997     #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2998     for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2999 jfenwick 2716 dataPointNo_0=0;
3000     // for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
3001 matt 1334 int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
3002     int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
3003 jfenwick 2271 const double *ptr_0 = &(arg_0_Z.getDataAtOffsetRO(offset_0));
3004     double *ptr_2 = &(res.getDataAtOffsetRW(offset_2));
3005 jfenwick 2716 tensor_unary_operation(size0*numDataPointsPerSample_0, ptr_0, ptr_2, operation);
3006     // }
3007 matt 1334 }
3008     }
3009     else {
3010     throw DataException("Error - C_TensorUnaryOperation: unknown combination of inputs");
3011     }
3012    
3013     return res;
3014 matt 1327 }
3015 matt 1334
3016     }
3017 jgs 94 #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26