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

Diff of /trunk/escript/src/Data.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/src/Data/Data.cpp revision 108 by jgs, Thu Jan 27 06:21:59 2005 UTC trunk/escript/src/Data.cpp revision 922 by gross, Fri Jan 5 04:23:05 2007 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
 /*=============================================================================  
2    
3   ******************************************************************************  /*
4   *                                                                            *   ************************************************************
5   *       COPYRIGHT ACcESS 2004 -  All Rights Reserved                         *   *          Copyright 2006 by ACcESS MNRF                   *
6   *                                                                            *   *                                                          *
7   * This software is the property of ACcESS.  No part of this code             *   *              http://www.access.edu.au                    *
8   * may be copied in any form or by any means without the expressed written    *   *       Primary Business: Queensland, Australia            *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *  Licensed under the Open Software License version 3.0    *
10   * by any unauthorised person is illegal unless that                          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
11   * person has a software license agreement with ACcESS.                       *   *                                                          *
12   *                                                                            *   ************************************************************
13   ******************************************************************************  */
14    #include "Data.h"
 ******************************************************************************/  
15    
16  #include "escript/Data/Data.h"  #include "DataExpanded.h"
17    #include "DataConstant.h"
18    #include "DataTagged.h"
19    #include "DataEmpty.h"
20    #include "DataArray.h"
21    #include "DataArrayView.h"
22    #include "DataProf.h"
23    #include "FunctionSpaceFactory.h"
24    #include "AbstractContinuousDomain.h"
25    #include "UnaryFuncs.h"
26    
27  #include <iostream>  #include <fstream>
28  #include <algorithm>  #include <algorithm>
29  #include <vector>  #include <vector>
 #include <exception>  
30  #include <functional>  #include <functional>
 #include <math.h>  
31    
32  #include <boost/python/str.hpp>  #include <boost/python/dict.hpp>
33  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
34  #include <boost/python/long.hpp>  #include <boost/python/long.hpp>
35    
 #include "escript/Data/DataException.h"  
   
 #include "escript/Data/DataExpanded.h"  
 #include "escript/Data/DataConstant.h"  
 #include "escript/Data/DataTagged.h"  
 #include "escript/Data/DataEmpty.h"  
 #include "escript/Data/DataArray.h"  
 #include "escript/Data/DataAlgorithm.h"  
 #include "escript/Data/FunctionSpaceFactory.h"  
 #include "escript/Data/AbstractContinuousDomain.h"  
 #include "escript/Data/UnaryFuncs.h"  
   
36  using namespace std;  using namespace std;
37  using namespace boost::python;  using namespace boost::python;
38  using namespace boost;  using namespace boost;
39  using namespace escript;  using namespace escript;
40    
41    #if defined DOPROF
42    //
43    // global table of profiling data for all Data objects
44    DataProf dataProfTable;
45    #endif
46    
47  Data::Data()  Data::Data()
48  {  {
49    //    //
# Line 52  Data::Data() Line 51  Data::Data()
51    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
52    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
53    m_data=temp_data;    m_data=temp_data;
54      m_protected=false;
55    #if defined DOPROF
56      // create entry in global profiling table for this object
57      profData = dataProfTable.newData();
58    #endif
59  }  }
60    
61  Data::Data(double value,  Data::Data(double value,
# Line 65  Data::Data(double value, Line 69  Data::Data(double value,
69    }    }
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
72      m_protected=false;
73    #if defined DOPROF
74      // create entry in global profiling table for this object
75      profData = dataProfTable.newData();
76    #endif
77  }  }
78    
79  Data::Data(double value,  Data::Data(double value,
# Line 75  Data::Data(double value, Line 84  Data::Data(double value,
84    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
85    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
86    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
87      m_protected=false;
88    #if defined DOPROF
89      // create entry in global profiling table for this object
90      profData = dataProfTable.newData();
91    #endif
92  }  }
93    
94  Data::Data(const Data& inData)  Data::Data(const Data& inData)
95  {  {
96    m_data=inData.m_data;    m_data=inData.m_data;
97      m_protected=inData.isProtected();
98    #if defined DOPROF
99      // create entry in global profiling table for this object
100      profData = dataProfTable.newData();
101    #endif
102  }  }
103    
104  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 90  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
110    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
111    m_data=temp_data;    m_data=temp_data;
112      m_protected=false;
113    #if defined DOPROF
114      // create entry in global profiling table for this object
115      profData = dataProfTable.newData();
116    #endif
117  }  }
118    
119  Data::Data(const Data& inData,  Data::Data(const Data& inData,
120             const FunctionSpace& functionspace)             const FunctionSpace& functionspace)
121  {  {
122    #if defined DOPROF
123      // create entry in global profiling table for this object
124      profData = dataProfTable.newData();
125    #endif
126    if (inData.getFunctionSpace()==functionspace) {    if (inData.getFunctionSpace()==functionspace) {
127      m_data=inData.m_data;      m_data=inData.m_data;
128    } else {    } else {
129        #if defined DOPROF
130        profData->interpolate++;
131        #endif
132      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
133      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
134      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
135      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
136      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
137      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
138      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
139        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 111  Data::Data(const Data& inData, Line 142  Data::Data(const Data& inData,
142      }      }
143      m_data=tmp.m_data;      m_data=tmp.m_data;
144    }    }
145      m_protected=false;
146  }  }
147    
148  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 122  Data::Data(const DataTagged::TagListType Line 154  Data::Data(const DataTagged::TagListType
154    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
155    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
156    m_data=temp_data;    m_data=temp_data;
157      m_protected=false;
158    if (expanded) {    if (expanded) {
159      expand();      expand();
160    }    }
161    #if defined DOPROF
162      // create entry in global profiling table for this object
163      profData = dataProfTable.newData();
164    #endif
165  }  }
166    
167  Data::Data(const numeric::array& value,  Data::Data(const numeric::array& value,
# Line 132  Data::Data(const numeric::array& value, Line 169  Data::Data(const numeric::array& value,
169             bool expanded)             bool expanded)
170  {  {
171    initialise(value,what,expanded);    initialise(value,what,expanded);
172      m_protected=false;
173    #if defined DOPROF
174      // create entry in global profiling table for this object
175      profData = dataProfTable.newData();
176    #endif
177  }  }
178    
179  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 139  Data::Data(const DataArrayView& value, Line 181  Data::Data(const DataArrayView& value,
181             bool expanded)             bool expanded)
182  {  {
183    initialise(value,what,expanded);    initialise(value,what,expanded);
184      m_protected=false;
185    #if defined DOPROF
186      // create entry in global profiling table for this object
187      profData = dataProfTable.newData();
188    #endif
189  }  }
190    
191  Data::Data(const object& value,  Data::Data(const object& value,
# Line 147  Data::Data(const object& value, Line 194  Data::Data(const object& value,
194  {  {
195    numeric::array asNumArray(value);    numeric::array asNumArray(value);
196    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
197      m_protected=false;
198    #if defined DOPROF
199      // create entry in global profiling table for this object
200      profData = dataProfTable.newData();
201    #endif
202  }  }
203    
204  Data::Data(const object& value,  Data::Data(const object& value,
# Line 167  Data::Data(const object& value, Line 219  Data::Data(const object& value,
219      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
220      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
221    }    }
222      m_protected=false;
223    #if defined DOPROF
224      // create entry in global profiling table for this object
225      profData = dataProfTable.newData();
226    #endif
227    }
228    
229    Data::~Data()
230    {
231    
232  }  }
233    
234  escriptDataC  escriptDataC
# Line 185  Data::getDataC() const Line 247  Data::getDataC() const
247    return temp;    return temp;
248  }  }
249    
250  tuple  const boost::python::tuple
251  Data::getShapeTuple() const  Data::getShapeTuple() const
252  {  {
253    const DataArrayView::ShapeType& shape=getDataPointShape();    const DataArrayView::ShapeType& shape=getDataPointShape();
# Line 303  Data::isConstant() const Line 365  Data::isConstant() const
365  }  }
366    
367  void  void
368    Data::setProtection()
369    {
370       m_protected=true;
371    }
372    
373    bool
374    Data::isProtected() const
375    {
376       return m_protected;
377    }
378    
379    
380    
381    void
382  Data::expand()  Data::expand()
383  {  {
384    if (isConstant()) {    if (isConstant()) {
# Line 344  Data::tag() Line 420  Data::tag()
420    }    }
421  }  }
422    
423  void  Data
424  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
425  {  {
426    m_data->reshapeDataPoint(shape);  #if defined DOPROF
427      profData->where++;
428    #endif
429      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
430  }  }
431    
432  Data  Data
433  Data::wherePositive() const  Data::wherePositive() const
434  {  {
435    #if defined DOPROF
436      profData->where++;
437    #endif
438    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater<double>(),0.0));
439  }  }
440    
441  Data  Data
442  Data::whereNegative() const  Data::whereNegative() const
443  {  {
444    #if defined DOPROF
445      profData->where++;
446    #endif
447    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less<double>(),0.0));
448  }  }
449    
450  Data  Data
451  Data::whereNonNegative() const  Data::whereNonNegative() const
452  {  {
453    #if defined DOPROF
454      profData->where++;
455    #endif
456    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(greater_equal<double>(),0.0));
457  }  }
458    
459  Data  Data
460  Data::whereNonPositive() const  Data::whereNonPositive() const
461  {  {
462    #if defined DOPROF
463      profData->where++;
464    #endif
465    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));    return escript::unaryOp(*this,bind2nd(less_equal<double>(),0.0));
466  }  }
467    
468  Data  Data
469  Data::whereZero() const  Data::whereZero(double tol) const
470  {  {
471    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));  #if defined DOPROF
472      profData->where++;
473    #endif
474      Data dataAbs=abs();
475      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
476  }  }
477    
478  Data  Data
479  Data::whereNonZero() const  Data::whereNonZero(double tol) const
480  {  {
481    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));  #if defined DOPROF
482      profData->where++;
483    #endif
484      Data dataAbs=abs();
485      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
486  }  }
487    
488  Data  Data
489  Data::interpolate(const FunctionSpace& functionspace) const  Data::interpolate(const FunctionSpace& functionspace) const
490  {  {
491    #if defined DOPROF
492      profData->interpolate++;
493    #endif
494    return Data(*this,functionspace);    return Data(*this,functionspace);
495  }  }
496    
# Line 410  Data::probeInterpolation(const FunctionS Line 512  Data::probeInterpolation(const FunctionS
512  Data  Data
513  Data::gradOn(const FunctionSpace& functionspace) const  Data::gradOn(const FunctionSpace& functionspace) const
514  {  {
515    #if defined DOPROF
516      profData->grad++;
517    #endif
518    if (functionspace.getDomain()!=getDomain())    if (functionspace.getDomain()!=getDomain())
519      throw DataException("Error - gradient cannot be calculated on different domains.");      throw DataException("Error - gradient cannot be calculated on different domains.");
520    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();    DataArrayView::ShapeType grad_shape=getPointDataView().getShape();
# Line 443  Data::getDataPointShape() const Line 548  Data::getDataPointShape() const
548    return getPointDataView().getShape();    return getPointDataView().getShape();
549  }  }
550    
551    
552    
553    void
554    Data::fillFromNumArray(const boost::python::numeric::array num_array)
555    {
556      if (isProtected()) {
557            throw DataException("Error - attempt to update protected Data object.");
558      }
559      //
560      // check rank
561      if (num_array.getrank()<getDataPointRank())
562          throw DataException("Rank of numarray does not match Data object rank");
563    
564      //
565      // check shape of num_array
566      for (int i=0; i<getDataPointRank(); i++) {
567        if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])
568           throw DataException("Shape of numarray does not match Data object rank");
569      }
570    
571      //
572      // make sure data is expanded:
573      if (!isExpanded()) {
574        expand();
575      }
576    
577      //
578      // and copy over
579      m_data->copyAll(num_array);
580    }
581    
582    const
583    boost::python::numeric::array
584    Data::convertToNumArray()
585    {
586      //
587      // determine the total number of data points
588      int numSamples = getNumSamples();
589      int numDataPointsPerSample = getNumDataPointsPerSample();
590      int numDataPoints = numSamples * numDataPointsPerSample;
591    
592      //
593      // determine the rank and shape of each data point
594      int dataPointRank = getDataPointRank();
595      DataArrayView::ShapeType dataPointShape = getDataPointShape();
596    
597      //
598      // create the numeric array to be returned
599      boost::python::numeric::array numArray(0.0);
600    
601      //
602      // the rank of the returned numeric array will be the rank of
603      // the data points, plus one. Where the rank of the array is n,
604      // the last n-1 dimensions will be equal to the shape of the
605      // data points, whilst the first dimension will be equal to the
606      // total number of data points. Thus the array will consist of
607      // a serial vector of the data points.
608      int arrayRank = dataPointRank + 1;
609      DataArrayView::ShapeType arrayShape;
610      arrayShape.push_back(numDataPoints);
611      for (int d=0; d<dataPointRank; d++) {
612         arrayShape.push_back(dataPointShape[d]);
613      }
614    
615      //
616      // resize the numeric array to the shape just calculated
617      if (arrayRank==1) {
618        numArray.resize(arrayShape[0]);
619      }
620      if (arrayRank==2) {
621        numArray.resize(arrayShape[0],arrayShape[1]);
622      }
623      if (arrayRank==3) {
624        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
625      }
626      if (arrayRank==4) {
627        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
628      }
629      if (arrayRank==5) {
630        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);
631      }
632    
633      //
634      // loop through each data point in turn, loading the values for that data point
635      // into the numeric array.
636      int dataPoint = 0;
637      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
638        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
639          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
640          if (dataPointRank==0) {
641            numArray[dataPoint]=dataPointView();
642          }
643          if (dataPointRank==1) {
644            for (int i=0; i<dataPointShape[0]; i++) {
645              numArray[make_tuple(dataPoint,i)]=dataPointView(i);
646            }
647          }
648          if (dataPointRank==2) {
649            for (int i=0; i<dataPointShape[0]; i++) {
650              for (int j=0; j<dataPointShape[1]; j++) {
651                numArray[make_tuple(dataPoint,i,j)] = dataPointView(i,j);
652              }
653            }
654          }
655          if (dataPointRank==3) {
656            for (int i=0; i<dataPointShape[0]; i++) {
657              for (int j=0; j<dataPointShape[1]; j++) {
658                for (int k=0; k<dataPointShape[2]; k++) {
659                  numArray[make_tuple(dataPoint,i,j,k)]=dataPointView(i,j,k);
660                }
661              }
662            }
663          }
664          if (dataPointRank==4) {
665            for (int i=0; i<dataPointShape[0]; i++) {
666              for (int j=0; j<dataPointShape[1]; j++) {
667                for (int k=0; k<dataPointShape[2]; k++) {
668                  for (int l=0; l<dataPointShape[3]; l++) {
669                    numArray[make_tuple(dataPoint,i,j,k,l)]=dataPointView(i,j,k,l);
670                  }
671                }
672              }
673            }
674          }
675          dataPoint++;
676        }
677      }
678    
679      //
680      // return the loaded array
681      return numArray;
682    }
683    
684    
685    const
686    boost::python::numeric::array
687    Data:: getValueOfDataPoint(int dataPointNo)
688    {
689      size_t length=0;
690      int i, j, k, l;
691      //
692      // determine the rank and shape of each data point
693      int dataPointRank = getDataPointRank();
694      DataArrayView::ShapeType dataPointShape = getDataPointShape();
695    
696      //
697      // create the numeric array to be returned
698      boost::python::numeric::array numArray(0.0);
699    
700      //
701      // the shape of the returned numeric array will be the same
702      // as that of the data point
703      int arrayRank = dataPointRank;
704      DataArrayView::ShapeType arrayShape = dataPointShape;
705    
706      //
707      // resize the numeric array to the shape just calculated
708      if (arrayRank==0) {
709        numArray.resize(1);
710      }
711      if (arrayRank==1) {
712        numArray.resize(arrayShape[0]);
713      }
714      if (arrayRank==2) {
715        numArray.resize(arrayShape[0],arrayShape[1]);
716      }
717      if (arrayRank==3) {
718        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
719      }
720      if (arrayRank==4) {
721        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
722      }
723    
724      if (getNumDataPointsPerSample()>0) {
725           int sampleNo = dataPointNo/getNumDataPointsPerSample();
726           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
727           //
728           // Check a valid sample number has been supplied
729           if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
730               throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
731           }
732                  
733           //
734           // Check a valid data point number has been supplied
735           if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
736               throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
737           }
738           // TODO: global error handling
739           // create a view of the data if it is stored locally
740           DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
741            
742           switch( dataPointRank ){
743                case 0 :
744                    numArray[0] = dataPointView();
745                    break;
746                case 1 :        
747                    for( i=0; i<dataPointShape[0]; i++ )
748                        numArray[i]=dataPointView(i);
749                    break;
750                case 2 :        
751                    for( i=0; i<dataPointShape[0]; i++ )
752                        for( j=0; j<dataPointShape[1]; j++)
753                            numArray[make_tuple(i,j)]=dataPointView(i,j);
754                    break;
755                case 3 :        
756                    for( i=0; i<dataPointShape[0]; i++ )
757                        for( j=0; j<dataPointShape[1]; j++ )
758                            for( k=0; k<dataPointShape[2]; k++)
759                                numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
760                    break;
761                case 4 :
762                    for( i=0; i<dataPointShape[0]; i++ )
763                        for( j=0; j<dataPointShape[1]; j++ )
764                            for( k=0; k<dataPointShape[2]; k++ )
765                                for( l=0; l<dataPointShape[3]; l++)
766                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
767                    break;
768        }
769      }
770      //
771      // return the array
772      return numArray;
773    
774    }
775    void
776    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array num_array)
777    {
778      if (isProtected()) {
779            throw DataException("Error - attempt to update protected Data object.");
780      }
781      //
782      // check rank
783      if (num_array.getrank()<getDataPointRank())
784          throw DataException("Rank of numarray does not match Data object rank");
785    
786      //
787      // check shape of num_array
788      for (int i=0; i<getDataPointRank(); i++) {
789        if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
790           throw DataException("Shape of numarray does not match Data object rank");
791      }
792      //
793      // make sure data is expanded:
794      if (!isExpanded()) {
795        expand();
796      }
797      if (getNumDataPointsPerSample()>0) {
798           int sampleNo = dataPointNo/getNumDataPointsPerSample();
799           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
800           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
801      } else {
802           m_data->copyToDataPoint(-1, 0,num_array);
803      }
804    }
805    
806    void
807    Data::setValueOfDataPoint(int dataPointNo, const double value)
808    {
809      if (isProtected()) {
810            throw DataException("Error - attempt to update protected Data object.");
811      }
812      //
813      // make sure data is expanded:
814      if (!isExpanded()) {
815        expand();
816      }
817      if (getNumDataPointsPerSample()>0) {
818           int sampleNo = dataPointNo/getNumDataPointsPerSample();
819           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
820           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
821      } else {
822           m_data->copyToDataPoint(-1, 0,value);
823      }
824    }
825    
826    const
827    boost::python::numeric::array
828    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
829    {
830      size_t length=0;
831      int i, j, k, l, pos;
832      //
833      // determine the rank and shape of each data point
834      int dataPointRank = getDataPointRank();
835      DataArrayView::ShapeType dataPointShape = getDataPointShape();
836    
837      //
838      // create the numeric array to be returned
839      boost::python::numeric::array numArray(0.0);
840    
841      //
842      // the shape of the returned numeric array will be the same
843      // as that of the data point
844      int arrayRank = dataPointRank;
845      DataArrayView::ShapeType arrayShape = dataPointShape;
846    
847      //
848      // resize the numeric array to the shape just calculated
849      if (arrayRank==0) {
850        numArray.resize(1);
851      }
852      if (arrayRank==1) {
853        numArray.resize(arrayShape[0]);
854      }
855      if (arrayRank==2) {
856        numArray.resize(arrayShape[0],arrayShape[1]);
857      }
858      if (arrayRank==3) {
859        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);
860      }
861      if (arrayRank==4) {
862        numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
863      }
864    
865      // added for the MPI communication
866      length=1;
867      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
868      double *tmpData = new double[length];
869    
870      //
871      // load the values for the data point into the numeric array.
872    
873        // updated for the MPI case
874        if( get_MPIRank()==procNo ){
875                 if (getNumDataPointsPerSample()>0) {
876                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
877                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
878                    //
879                    // Check a valid sample number has been supplied
880                    if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
881                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
882                    }
883                  
884                    //
885                    // Check a valid data point number has been supplied
886                    if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
887                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
888                    }
889                    // TODO: global error handling
890            // create a view of the data if it is stored locally
891            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
892            
893            // pack the data from the view into tmpData for MPI communication
894            pos=0;
895            switch( dataPointRank ){
896                case 0 :
897                    tmpData[0] = dataPointView();
898                    break;
899                case 1 :        
900                    for( i=0; i<dataPointShape[0]; i++ )
901                        tmpData[i]=dataPointView(i);
902                    break;
903                case 2 :        
904                    for( i=0; i<dataPointShape[0]; i++ )
905                        for( j=0; j<dataPointShape[1]; j++, pos++ )
906                            tmpData[pos]=dataPointView(i,j);
907                    break;
908                case 3 :        
909                    for( i=0; i<dataPointShape[0]; i++ )
910                        for( j=0; j<dataPointShape[1]; j++ )
911                            for( k=0; k<dataPointShape[2]; k++, pos++ )
912                                tmpData[pos]=dataPointView(i,j,k);
913                    break;
914                case 4 :
915                    for( i=0; i<dataPointShape[0]; i++ )
916                        for( j=0; j<dataPointShape[1]; j++ )
917                            for( k=0; k<dataPointShape[2]; k++ )
918                                for( l=0; l<dataPointShape[3]; l++, pos++ )
919                                    tmpData[pos]=dataPointView(i,j,k,l);
920                    break;
921            }
922                }
923        }
924            #ifdef PASO_MPI
925            // broadcast the data to all other processes
926        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
927            #endif
928    
929        // unpack the data
930        switch( dataPointRank ){
931            case 0 :
932                numArray[0]=tmpData[0];
933                break;
934            case 1 :        
935                for( i=0; i<dataPointShape[0]; i++ )
936                    numArray[i]=tmpData[i];
937                break;
938            case 2 :        
939                for( i=0; i<dataPointShape[0]; i++ )
940                    for( j=0; j<dataPointShape[1]; j++ )
941                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
942                break;
943            case 3 :        
944                for( i=0; i<dataPointShape[0]; i++ )
945                    for( j=0; j<dataPointShape[1]; j++ )
946                        for( k=0; k<dataPointShape[2]; k++ )
947                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
948                break;
949            case 4 :
950                for( i=0; i<dataPointShape[0]; i++ )
951                    for( j=0; j<dataPointShape[1]; j++ )
952                        for( k=0; k<dataPointShape[2]; k++ )
953                            for( l=0; l<dataPointShape[3]; l++ )
954                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
955                break;
956        }
957    
958        delete [] tmpData;  
959      //
960      // return the loaded array
961      return numArray;
962    }
963    
964    
965    
966  boost::python::numeric::array  boost::python::numeric::array
967  Data::integrate() const  Data::integrate() const
968  {  {
# Line 450  Data::integrate() const Line 970  Data::integrate() const
970    int rank = getDataPointRank();    int rank = getDataPointRank();
971    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
972    
973    #if defined DOPROF
974      profData->integrate++;
975    #endif
976    
977    //    //
978    // calculate the integral values    // calculate the integral values
979    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 472  Data::integrate() const Line 996  Data::integrate() const
996      }      }
997    }    }
998    if (rank==2) {    if (rank==2) {
999      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
1000      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
1001        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
1002          index = i + shape[0] * j;             index = i + shape[0] * j;
1003          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
1004        }           }
1005      }         }
1006    }    }
1007    if (rank==3) {    if (rank==3) {
1008      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 486  Data::integrate() const Line 1010  Data::integrate() const
1010        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
1011          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1012            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
1013            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
1014          }          }
1015        }        }
1016      }      }
# Line 498  Data::integrate() const Line 1022  Data::integrate() const
1022          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
1023            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
1024              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
1025              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
1026            }            }
1027          }          }
1028        }        }
# Line 511  Data::integrate() const Line 1035  Data::integrate() const
1035  }  }
1036    
1037  Data  Data
1038    Data::erf() const
1039    {
1040    #if defined DOPROF
1041      profData->unary++;
1042    #endif
1043      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1044    }
1045    
1046    Data
1047  Data::sin() const  Data::sin() const
1048  {  {
1049    #if defined DOPROF
1050      profData->unary++;
1051    #endif
1052    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sin);
1053  }  }
1054    
1055  Data  Data
1056  Data::cos() const  Data::cos() const
1057  {  {
1058    #if defined DOPROF
1059      profData->unary++;
1060    #endif
1061    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cos);
1062  }  }
1063    
1064  Data  Data
1065  Data::tan() const  Data::tan() const
1066  {  {
1067    #if defined DOPROF
1068      profData->unary++;
1069    #endif
1070    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tan);
1071  }  }
1072    
1073  Data  Data
1074  Data::log() const  Data::asin() const
1075    {
1076    #if defined DOPROF
1077      profData->unary++;
1078    #endif
1079      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
1080    }
1081    
1082    Data
1083    Data::acos() const
1084    {
1085    #if defined DOPROF
1086      profData->unary++;
1087    #endif
1088      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1089    }
1090    
1091    Data
1092    Data::atan() const
1093    {
1094    #if defined DOPROF
1095      profData->unary++;
1096    #endif
1097      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
1098    }
1099    
1100    Data
1101    Data::sinh() const
1102    {
1103    #if defined DOPROF
1104      profData->unary++;
1105    #endif
1106      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
1107    }
1108    
1109    Data
1110    Data::cosh() const
1111    {
1112    #if defined DOPROF
1113      profData->unary++;
1114    #endif
1115      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
1116    }
1117    
1118    Data
1119    Data::tanh() const
1120    {
1121    #if defined DOPROF
1122      profData->unary++;
1123    #endif
1124      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1125    }
1126    
1127    Data
1128    Data::asinh() const
1129    {
1130    #if defined DOPROF
1131      profData->unary++;
1132    #endif
1133      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1134    }
1135    
1136    Data
1137    Data::acosh() const
1138  {  {
1139    #if defined DOPROF
1140      profData->unary++;
1141    #endif
1142      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1143    }
1144    
1145    Data
1146    Data::atanh() const
1147    {
1148    #if defined DOPROF
1149      profData->unary++;
1150    #endif
1151      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1152    }
1153    
1154    Data
1155    Data::log10() const
1156    {
1157    #if defined DOPROF
1158      profData->unary++;
1159    #endif
1160    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
1161  }  }
1162    
1163  Data  Data
1164  Data::ln() const  Data::log() const
1165  {  {
1166    #if defined DOPROF
1167      profData->unary++;
1168    #endif
1169    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
1170  }  }
1171    
1172  Data  Data
1173  Data::sign() const  Data::sign() const
1174  {  {
1175    #if defined DOPROF
1176      profData->unary++;
1177    #endif
1178    return escript::unaryOp(*this,escript::fsign);    return escript::unaryOp(*this,escript::fsign);
1179  }  }
1180    
1181  Data  Data
1182  Data::abs() const  Data::abs() const
1183  {  {
1184    #if defined DOPROF
1185      profData->unary++;
1186    #endif
1187    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::fabs);
1188  }  }
1189    
1190  Data  Data
1191  Data::neg() const  Data::neg() const
1192  {  {
1193    #if defined DOPROF
1194      profData->unary++;
1195    #endif
1196    return escript::unaryOp(*this,negate<double>());    return escript::unaryOp(*this,negate<double>());
1197  }  }
1198    
1199  Data  Data
1200  Data::pos() const  Data::pos() const
1201  {  {
1202    return (*this);  #if defined DOPROF
1203      profData->unary++;
1204    #endif
1205      Data result;
1206      // perform a deep copy
1207      result.copy(*this);
1208      return result;
1209  }  }
1210    
1211  Data  Data
1212  Data::exp() const  Data::exp() const
1213  {  {
1214    #if defined DOPROF
1215      profData->unary++;
1216    #endif
1217    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::exp);
1218  }  }
1219    
1220  Data  Data
1221  Data::sqrt() const  Data::sqrt() const
1222  {  {
1223    #if defined DOPROF
1224      profData->unary++;
1225    #endif
1226    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sqrt);
1227  }  }
1228    
1229  double  double
1230  Data::Lsup() const  Data::Lsup() const
1231  {  {
1232      double localValue, globalValue;
1233    #if defined DOPROF
1234      profData->reduction1++;
1235    #endif
1236    //    //
1237    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
1238    return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
1239      AbsMax abs_max_func;
1240      localValue = algorithm(abs_max_func,0);
1241    #ifdef PASO_MPI
1242      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1243      return globalValue;
1244    #else
1245      return localValue;
1246    #endif
1247    }
1248    
1249    double
1250    Data::Linf() const
1251    {
1252      double localValue, globalValue;
1253    #if defined DOPROF
1254      profData->reduction1++;
1255    #endif
1256      //
1257      // set the initial absolute minimum value to max double
1258      AbsMin abs_min_func;
1259      localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1260    
1261    #ifdef PASO_MPI
1262      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1263      return globalValue;
1264    #else
1265      return localValue;
1266    #endif
1267  }  }
1268    
1269  double  double
1270  Data::sup() const  Data::sup() const
1271  {  {
1272      double localValue, globalValue;
1273    #if defined DOPROF
1274      profData->reduction1++;
1275    #endif
1276    //    //
1277    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1278    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1279      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1280    #ifdef PASO_MPI
1281      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1282      return globalValue;
1283    #else
1284      return localValue;
1285    #endif
1286  }  }
1287    
1288  double  double
1289  Data::inf() const  Data::inf() const
1290  {  {
1291      double localValue, globalValue;
1292    #if defined DOPROF
1293      profData->reduction1++;
1294    #endif
1295    //    //
1296    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1297    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1298      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1299    #ifdef PASO_MPI
1300      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1301      return globalValue;
1302    #else
1303      return localValue;
1304    #endif
1305  }  }
1306    
1307    /* TODO */
1308    /* global reduction */
1309  Data  Data
1310  Data::maxval() const  Data::maxval() const
1311  {  {
1312    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));  #if defined DOPROF
1313      profData->reduction2++;
1314    #endif
1315      //
1316      // set the initial maximum value to min possible double
1317      FMax fmax_func;
1318      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1319  }  }
1320    
1321  Data  Data
1322  Data::minval() const  Data::minval() const
1323  {  {
1324    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));  #if defined DOPROF
1325      profData->reduction2++;
1326    #endif
1327      //
1328      // set the initial minimum value to max possible double
1329      FMin fmin_func;
1330      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1331  }  }
1332    
1333  Data  Data
1334  Data::length() const  Data::swapaxes(const int axis0, const int axis1) const
1335  {  {
1336    return dp_algorithm(DataAlgorithmAdapter<Length>(0));       int axis0_tmp,axis1_tmp;
1337         #if defined DOPROF
1338         profData->unary++;
1339         #endif
1340         DataArrayView::ShapeType s=getDataPointShape();
1341         DataArrayView::ShapeType ev_shape;
1342         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1343         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1344         int rank=getDataPointRank();
1345         if (rank<2) {
1346            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1347         }
1348         if (axis0<0 || axis0>rank-1) {
1349            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1350         }
1351         if (axis1<0 || axis1>rank-1) {
1352             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1353         }
1354         if (axis0 == axis1) {
1355             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1356         }
1357         if (axis0 > axis1) {
1358             axis0_tmp=axis1;
1359             axis1_tmp=axis0;
1360         } else {
1361             axis0_tmp=axis0;
1362             axis1_tmp=axis1;
1363         }
1364         for (int i=0; i<rank; i++) {
1365           if (i == axis0_tmp) {
1366              ev_shape.push_back(s[axis1_tmp]);
1367           } else if (i == axis1_tmp) {
1368              ev_shape.push_back(s[axis0_tmp]);
1369           } else {
1370              ev_shape.push_back(s[i]);
1371           }
1372         }
1373         Data ev(0.,ev_shape,getFunctionSpace());
1374         ev.typeMatchRight(*this);
1375         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1376         return ev;
1377    
1378    }
1379    
1380    Data
1381    Data::symmetric() const
1382    {
1383         #if defined DOPROF
1384            profData->unary++;
1385         #endif
1386         // check input
1387         DataArrayView::ShapeType s=getDataPointShape();
1388         if (getDataPointRank()==2) {
1389            if(s[0] != s[1])
1390               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1391         }
1392         else if (getDataPointRank()==4) {
1393            if(!(s[0] == s[2] && s[1] == s[3]))
1394               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1395         }
1396         else {
1397            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1398         }
1399         Data ev(0.,getDataPointShape(),getFunctionSpace());
1400         ev.typeMatchRight(*this);
1401         m_data->symmetric(ev.m_data.get());
1402         return ev;
1403    }
1404    
1405    Data
1406    Data::nonsymmetric() const
1407    {
1408         #if defined DOPROF
1409            profData->unary++;
1410         #endif
1411         // check input
1412         DataArrayView::ShapeType s=getDataPointShape();
1413         if (getDataPointRank()==2) {
1414            if(s[0] != s[1])
1415               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1416            DataArrayView::ShapeType ev_shape;
1417            ev_shape.push_back(s[0]);
1418            ev_shape.push_back(s[1]);
1419            Data ev(0.,ev_shape,getFunctionSpace());
1420            ev.typeMatchRight(*this);
1421            m_data->nonsymmetric(ev.m_data.get());
1422            return ev;
1423         }
1424         else if (getDataPointRank()==4) {
1425            if(!(s[0] == s[2] && s[1] == s[3]))
1426               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1427            DataArrayView::ShapeType ev_shape;
1428            ev_shape.push_back(s[0]);
1429            ev_shape.push_back(s[1]);
1430            ev_shape.push_back(s[2]);
1431            ev_shape.push_back(s[3]);
1432            Data ev(0.,ev_shape,getFunctionSpace());
1433            ev.typeMatchRight(*this);
1434            m_data->nonsymmetric(ev.m_data.get());
1435            return ev;
1436         }
1437         else {
1438            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1439         }
1440    }
1441    
1442    Data
1443    Data::trace(int axis_offset) const
1444    {
1445         #if defined DOPROF
1446            profData->unary++;
1447         #endif
1448         DataArrayView::ShapeType s=getDataPointShape();
1449         if (getDataPointRank()==2) {
1450            DataArrayView::ShapeType ev_shape;
1451            Data ev(0.,ev_shape,getFunctionSpace());
1452            ev.typeMatchRight(*this);
1453            m_data->trace(ev.m_data.get(), axis_offset);
1454            return ev;
1455         }
1456         if (getDataPointRank()==3) {
1457            DataArrayView::ShapeType ev_shape;
1458            if (axis_offset==0) {
1459              int s2=s[2];
1460              ev_shape.push_back(s2);
1461            }
1462            else if (axis_offset==1) {
1463              int s0=s[0];
1464              ev_shape.push_back(s0);
1465            }
1466            Data ev(0.,ev_shape,getFunctionSpace());
1467            ev.typeMatchRight(*this);
1468            m_data->trace(ev.m_data.get(), axis_offset);
1469            return ev;
1470         }
1471         if (getDataPointRank()==4) {
1472            DataArrayView::ShapeType ev_shape;
1473            if (axis_offset==0) {
1474              ev_shape.push_back(s[2]);
1475              ev_shape.push_back(s[3]);
1476            }
1477            else if (axis_offset==1) {
1478              ev_shape.push_back(s[0]);
1479              ev_shape.push_back(s[3]);
1480            }
1481        else if (axis_offset==2) {
1482          ev_shape.push_back(s[0]);
1483          ev_shape.push_back(s[1]);
1484        }
1485            Data ev(0.,ev_shape,getFunctionSpace());
1486            ev.typeMatchRight(*this);
1487        m_data->trace(ev.m_data.get(), axis_offset);
1488            return ev;
1489         }
1490         else {
1491            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1492         }
1493    }
1494    
1495    Data
1496    Data::transpose(int axis_offset) const
1497    {
1498         #if defined DOPROF
1499         profData->unary++;
1500         #endif
1501         DataArrayView::ShapeType s=getDataPointShape();
1502         DataArrayView::ShapeType ev_shape;
1503         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1504         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1505         int rank=getDataPointRank();
1506         if (axis_offset<0 || axis_offset>rank) {
1507            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1508         }
1509         for (int i=0; i<rank; i++) {
1510           int index = (axis_offset+i)%rank;
1511           ev_shape.push_back(s[index]); // Append to new shape
1512         }
1513         Data ev(0.,ev_shape,getFunctionSpace());
1514         ev.typeMatchRight(*this);
1515         m_data->transpose(ev.m_data.get(), axis_offset);
1516         return ev;
1517    }
1518    
1519    Data
1520    Data::eigenvalues() const
1521    {
1522         #if defined DOPROF
1523            profData->unary++;
1524         #endif
1525         // check input
1526         DataArrayView::ShapeType s=getDataPointShape();
1527         if (getDataPointRank()!=2)
1528            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1529         if(s[0] != s[1])
1530            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1531         // create return
1532         DataArrayView::ShapeType ev_shape(1,s[0]);
1533         Data ev(0.,ev_shape,getFunctionSpace());
1534         ev.typeMatchRight(*this);
1535         m_data->eigenvalues(ev.m_data.get());
1536         return ev;
1537    }
1538    
1539    const boost::python::tuple
1540    Data::eigenvalues_and_eigenvectors(const double tol) const
1541    {
1542         #if defined DOPROF
1543            profData->unary++;
1544         #endif
1545         DataArrayView::ShapeType s=getDataPointShape();
1546         if (getDataPointRank()!=2)
1547            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1548         if(s[0] != s[1])
1549            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1550         // create return
1551         DataArrayView::ShapeType ev_shape(1,s[0]);
1552         Data ev(0.,ev_shape,getFunctionSpace());
1553         ev.typeMatchRight(*this);
1554         DataArrayView::ShapeType V_shape(2,s[0]);
1555         Data V(0.,V_shape,getFunctionSpace());
1556         V.typeMatchRight(*this);
1557         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1558         return make_tuple(boost::python::object(ev),boost::python::object(V));
1559    }
1560    
1561    const boost::python::tuple
1562    Data::minGlobalDataPoint() const
1563    {
1564      // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1565      // abort (for unknown reasons) if there are openmp directives with it in the
1566      // surrounding function
1567    
1568      int DataPointNo;
1569      int ProcNo;
1570      calc_minGlobalDataPoint(ProcNo,DataPointNo);
1571      return make_tuple(ProcNo,DataPointNo);
1572  }  }
1573    
1574  Data  void
1575  Data::trace() const  Data::calc_minGlobalDataPoint(int& ProcNo,
1576                            int& DataPointNo) const
1577  {  {
1578    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    int i,j;
1579      int lowi=0,lowj=0;
1580      double min=numeric_limits<double>::max();
1581    
1582      Data temp=minval();
1583    
1584      int numSamples=temp.getNumSamples();
1585      int numDPPSample=temp.getNumDataPointsPerSample();
1586    
1587      double next,local_min;
1588      int local_lowi,local_lowj;
1589    
1590      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1591      {
1592        local_min=min;
1593        #pragma omp for private(i,j) schedule(static)
1594        for (i=0; i<numSamples; i++) {
1595          for (j=0; j<numDPPSample; j++) {
1596            next=temp.getDataPoint(i,j)();
1597            if (next<local_min) {
1598              local_min=next;
1599              local_lowi=i;
1600              local_lowj=j;
1601            }
1602          }
1603        }
1604        #pragma omp critical
1605        if (local_min<min) {
1606          min=local_min;
1607          lowi=local_lowi;
1608          lowj=local_lowj;
1609        }
1610      }
1611    
1612    #ifdef PASO_MPI
1613        // determine the processor on which the minimum occurs
1614        next = temp.getDataPoint(lowi,lowj)();
1615        int lowProc = 0;
1616        double *globalMins = new double[get_MPISize()+1];
1617        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1618        
1619        if( get_MPIRank()==0 ){
1620            next = globalMins[lowProc];
1621            for( i=1; i<get_MPISize(); i++ )
1622                if( next>globalMins[i] ){
1623                    lowProc = i;
1624                    next = globalMins[i];
1625                }
1626        }
1627        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1628    
1629        delete [] globalMins;
1630        ProcNo = lowProc;
1631    #else
1632        ProcNo = 0;
1633    #endif
1634      DataPointNo = lowj + lowi * numDPPSample;
1635  }  }
1636    
1637  Data  void
1638  Data::transpose(int axis) const  Data::saveDX(std::string fileName) const
1639  {  {
1640    // not implemented    boost::python::dict args;
1641    throw DataException("Error - Data::transpose not implemented yet.");    args["data"]=boost::python::object(this);
1642    return Data();    getDomain().saveDX(fileName,args);
1643      return;
1644  }  }
1645    
1646  void  void
1647  Data::saveDX(std::string fileName) const  Data::saveVTK(std::string fileName) const
1648  {  {
1649    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1650      args["data"]=boost::python::object(this);
1651      getDomain().saveVTK(fileName,args);
1652    return;    return;
1653  }  }
1654    
1655  Data&  Data&
1656  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1657  {  {
1658      if (isProtected()) {
1659            throw DataException("Error - attempt to update protected Data object.");
1660      }
1661    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1662    return (*this);    return (*this);
1663  }  }
# Line 649  Data::operator+=(const Data& right) Line 1665  Data::operator+=(const Data& right)
1665  Data&  Data&
1666  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1667  {  {
1668    binaryOp(right,plus<double>());    Data tmp(right,getFunctionSpace(),false);
1669      binaryOp(tmp,plus<double>());
1670    return (*this);    return (*this);
1671  }  }
1672    
1673  Data&  Data&
1674  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1675  {  {
1676      if (isProtected()) {
1677            throw DataException("Error - attempt to update protected Data object.");
1678      }
1679    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1680    return (*this);    return (*this);
1681  }  }
# Line 663  Data::operator-=(const Data& right) Line 1683  Data::operator-=(const Data& right)
1683  Data&  Data&
1684  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1685  {  {
1686    binaryOp(right,minus<double>());    Data tmp(right,getFunctionSpace(),false);
1687      binaryOp(tmp,minus<double>());
1688    return (*this);    return (*this);
1689  }  }
1690    
1691  Data&  Data&
1692  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1693  {  {
1694      if (isProtected()) {
1695            throw DataException("Error - attempt to update protected Data object.");
1696      }
1697    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1698    return (*this);    return (*this);
1699  }  }
# Line 677  Data::operator*=(const Data& right) Line 1701  Data::operator*=(const Data& right)
1701  Data&  Data&
1702  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1703  {  {
1704    binaryOp(right,multiplies<double>());    Data tmp(right,getFunctionSpace(),false);
1705      binaryOp(tmp,multiplies<double>());
1706    return (*this);    return (*this);
1707  }  }
1708    
1709  Data&  Data&
1710  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1711  {  {
1712      if (isProtected()) {
1713            throw DataException("Error - attempt to update protected Data object.");
1714      }
1715    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1716    return (*this);    return (*this);
1717  }  }
# Line 691  Data::operator/=(const Data& right) Line 1719  Data::operator/=(const Data& right)
1719  Data&  Data&
1720  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1721  {  {
1722    binaryOp(right,divides<double>());    Data tmp(right,getFunctionSpace(),false);
1723      binaryOp(tmp,divides<double>());
1724    return (*this);    return (*this);
1725  }  }
1726    
1727  Data  Data
1728    Data::rpowO(const boost::python::object& left) const
1729    {
1730      Data left_d(left,*this);
1731      return left_d.powD(*this);
1732    }
1733    
1734    Data
1735  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1736  {  {
1737    Data result;    Data tmp(right,getFunctionSpace(),false);
1738    result.copy(*this);    return powD(tmp);
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1739  }  }
1740    
1741  Data  Data
1742  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1743  {  {
1744    Data result;    Data result;
1745    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1746    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1747         result.binaryOp(*this,escript::rpow);
1748      } else {
1749         result.copy(*this);
1750         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1751      }
1752    return result;    return result;
1753  }  }
1754    
1755    
1756  //  //
1757  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1758  Data  Data
1759  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1760  {  {
1761    Data result;    Data result;
1762    //    //
1763    // perform a deep copy    // perform a deep copy
1764    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1765    result+=right;       result.copy(right);
1766         result+=left;
1767      } else {
1768         result.copy(left);
1769         result+=right;
1770      }
1771    return result;    return result;
1772  }  }
1773    
1774  //  //
1775  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1776  Data  Data
1777  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1778  {  {
1779    Data result;    Data result;
1780    //    //
1781    // perform a deep copy    // perform a deep copy
1782    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1783    result-=right;       result=right.neg();
1784         result+=left;
1785      } else {
1786         result.copy(left);
1787         result-=right;
1788      }
1789    return result;    return result;
1790  }  }
1791    
1792  //  //
1793  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1794  Data  Data
1795  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1796  {  {
1797    Data result;    Data result;
1798    //    //
1799    // perform a deep copy    // perform a deep copy
1800    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1801    result*=right;       result.copy(right);
1802         result*=left;
1803      } else {
1804         result.copy(left);
1805         result*=right;
1806      }
1807    return result;    return result;
1808  }  }
1809    
1810  //  //
1811  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1812  Data  Data
1813  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1814  {  {
1815    Data result;    Data result;
1816    //    //
1817    // perform a deep copy    // perform a deep copy
1818    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1819    result/=right;       result=right.oneOver();
1820         result*=left;
1821      } else {
1822         result.copy(left);
1823         result/=right;
1824      }
1825    return result;    return result;
1826  }  }
1827    
1828  //  //
1829  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1830  Data  Data
1831  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1832  {  {
1833    //    return left+Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result+=right;  
   return result;  
1834  }  }
1835    
1836  //  //
1837  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1838  Data  Data
1839  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1840  {  {
1841    //    return left-Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result-=right;  
   return result;  
1842  }  }
1843    
1844  //  //
1845  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1846  Data  Data
1847  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1848  {  {
1849    //    return left*Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result*=right;  
   return result;  
1850  }  }
1851    
1852  //  //
1853  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1854  Data  Data
1855  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1856  {  {
1857    //    return left/Data(right,left.getFunctionSpace(),false);
   // Convert to DataArray format if possible  
   DataArray temp(right);  
   Data result;  
   //  
   // perform a deep copy  
   result.copy(left);  
   result/=right;  
   return result;  
1858  }  }
1859    
1860  //  //
1861  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1862  Data  Data
1863  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1864  {  {
1865    //    return Data(left,right.getFunctionSpace(),false)+right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result+=right;  
   return result;  
1866  }  }
1867    
1868  //  //
1869  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1870  Data  Data
1871  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1872  {  {
1873    //    return Data(left,right.getFunctionSpace(),false)-right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result-=right;  
   return result;  
1874  }  }
1875    
1876  //  //
1877  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1878  Data  Data
1879  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1880  {  {
1881    //    return Data(left,right.getFunctionSpace(),false)*right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result*=right;  
   return result;  
1882  }  }
1883    
1884  //  //
1885  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1886  Data  Data
1887  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1888  {  {
1889    //    return Data(left,right.getFunctionSpace(),false)/right;
   // Construct the result using the given value and the other parameters  
   // from right  
   Data result(left,right);  
   result/=right;  
   return result;  
1890  }  }
1891    
1892  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1893  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1894  //{  //{
1895  //  /*  //  /*
# Line 927  escript::operator/(const boost::python:: Line 1934  escript::operator/(const boost::python::
1934  //  return ret;  //  return ret;
1935  //}  //}
1936    
1937    /* TODO */
1938    /* global reduction */
1939  Data  Data
1940  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1941  {  {
# Line 941  Data::getItem(const boost::python::objec Line 1950  Data::getItem(const boost::python::objec
1950    return getSlice(slice_region);    return getSlice(slice_region);
1951  }  }
1952    
1953    /* TODO */
1954    /* global reduction */
1955  Data  Data
1956  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1957  {  {
1958    #if defined DOPROF
1959      profData->slicing++;
1960    #endif
1961    return Data(*this,region);    return Data(*this,region);
1962  }  }
1963    
1964    /* TODO */
1965    /* global reduction */
1966  void  void
1967  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1968                 const boost::python::object& value)                 const boost::python::object& value)
# Line 955  Data::setItemO(const boost::python::obje Line 1971  Data::setItemO(const boost::python::obje
1971    setItemD(key,tempData);    setItemD(key,tempData);
1972  }  }
1973    
1974    /* TODO */
1975    /* global reduction */
1976  void  void
1977  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1978                 const Data& value)                 const Data& value)
1979  {  {
1980    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1981    
1982    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1983    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1984      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 971  Data::setItemD(const boost::python::obje Line 1990  Data::setItemD(const boost::python::obje
1990    }    }
1991  }  }
1992    
1993    /* TODO */
1994    /* global reduction */
1995  void  void
1996  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1997                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1998  {  {
1999      if (isProtected()) {
2000            throw DataException("Error - attempt to update protected Data object.");
2001      }
2002    #if defined DOPROF
2003      profData->slicing++;
2004    #endif
2005    Data tempValue(value);    Data tempValue(value);
2006    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
2007    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1009  Data::typeMatchRight(const Data& right) Line 2036  Data::typeMatchRight(const Data& right)
2036    }    }
2037  }  }
2038    
2039    /* TODO */
2040    /* global reduction */
2041  void  void
2042  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
2043                       const boost::python::object& value)                       const boost::python::object& value)
2044  {  {
2045      if (isProtected()) {
2046            throw DataException("Error - attempt to update protected Data object.");
2047      }
2048    //    //
2049    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2050    tag();    tag();
# Line 1030  Data::setTaggedValue(int tagKey, Line 2062  Data::setTaggedValue(int tagKey,
2062    m_data->setTaggedValue(tagKey,valueDataArray.getView());    m_data->setTaggedValue(tagKey,valueDataArray.getView());
2063  }  }
2064    
2065  /*  /* TODO */
2066  Note: this version removed for now. Not needed, and breaks escript.cpp  /* global reduction */
2067  void  void
2068  Data::setTaggedValue(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
2069                       const DataArrayView& value)                              const DataArrayView& value)
2070  {  {
2071      if (isProtected()) {
2072            throw DataException("Error - attempt to update protected Data object.");
2073      }
2074    //    //
2075    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
2076    tag();    tag();
# Line 1048  Data::setTaggedValue(int tagKey, Line 2083  Data::setTaggedValue(int tagKey,
2083    // Call DataAbstract::setTaggedValue    // Call DataAbstract::setTaggedValue
2084    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
2085  }  }
2086  */  
2087    /* TODO */
2088    /* global reduction */
2089    int
2090    Data::getTagNumber(int dpno)
2091    {
2092      return m_data->getTagNumber(dpno);
2093    }
2094    
2095    /* TODO */
2096    /* global reduction */
2097    void
2098    Data::setRefValue(int ref,
2099                      const boost::python::numeric::array& value)
2100    {
2101      if (isProtected()) {
2102            throw DataException("Error - attempt to update protected Data object.");
2103      }
2104      //
2105      // Construct DataArray from boost::python::object input value
2106      DataArray valueDataArray(value);
2107    
2108      //
2109      // Call DataAbstract::setRefValue
2110      m_data->setRefValue(ref,valueDataArray);
2111    }
2112    
2113    /* TODO */
2114    /* global reduction */
2115    void
2116    Data::getRefValue(int ref,
2117                      boost::python::numeric::array& value)
2118    {
2119      //
2120      // Construct DataArray for boost::python::object return value
2121      DataArray valueDataArray(value);
2122    
2123      //
2124      // Load DataArray with values from data-points specified by ref
2125      m_data->getRefValue(ref,valueDataArray);
2126    
2127      //
2128      // Load values from valueDataArray into return numarray
2129    
2130      // extract the shape of the numarray
2131      int rank = value.getrank();
2132      DataArrayView::ShapeType shape;
2133      for (int i=0; i < rank; i++) {
2134        shape.push_back(extract<int>(value.getshape()[i]));
2135      }
2136    
2137      // and load the numarray with the data from the DataArray
2138      DataArrayView valueView = valueDataArray.getView();
2139    
2140      if (rank==0) {
2141          boost::python::numeric::array temp_numArray(valueView());
2142          value = temp_numArray;
2143      }
2144      if (rank==1) {
2145        for (int i=0; i < shape[0]; i++) {
2146          value[i] = valueView(i);
2147        }
2148      }
2149      if (rank==2) {
2150        for (int i=0; i < shape[0]; i++) {
2151          for (int j=0; j < shape[1]; j++) {
2152            value[make_tuple(i,j)] = valueView(i,j);
2153          }
2154        }
2155      }
2156      if (rank==3) {
2157        for (int i=0; i < shape[0]; i++) {
2158          for (int j=0; j < shape[1]; j++) {
2159            for (int k=0; k < shape[2]; k++) {
2160              value[make_tuple(i,j,k)] = valueView(i,j,k);
2161            }
2162          }
2163        }
2164      }
2165      if (rank==4) {
2166        for (int i=0; i < shape[0]; i++) {
2167          for (int j=0; j < shape[1]; j++) {
2168            for (int k=0; k < shape[2]; k++) {
2169              for (int l=0; l < shape[3]; l++) {
2170                value[make_tuple(i,j,k,l)] = valueView(i,j,k,l);
2171              }
2172            }
2173          }
2174        }
2175      }
2176    
2177    }
2178    
2179    void
2180    Data::archiveData(const std::string fileName)
2181    {
2182      cout << "Archiving Data object to: " << fileName << endl;
2183    
2184      //
2185      // Determine type of this Data object
2186      int dataType = -1;
2187    
2188      if (isEmpty()) {
2189        dataType = 0;
2190        cout << "\tdataType: DataEmpty" << endl;
2191      }
2192      if (isConstant()) {
2193        dataType = 1;
2194        cout << "\tdataType: DataConstant" << endl;
2195      }
2196      if (isTagged()) {
2197        dataType = 2;
2198        cout << "\tdataType: DataTagged" << endl;
2199      }
2200      if (isExpanded()) {
2201        dataType = 3;
2202        cout << "\tdataType: DataExpanded" << endl;
2203      }
2204    
2205      if (dataType == -1) {
2206        throw DataException("archiveData Error: undefined dataType");
2207      }
2208    
2209      //
2210      // Collect data items common to all Data types
2211      int noSamples = getNumSamples();
2212      int noDPPSample = getNumDataPointsPerSample();
2213      int functionSpaceType = getFunctionSpace().getTypeCode();
2214      int dataPointRank = getDataPointRank();
2215      int dataPointSize = getDataPointSize();
2216      int dataLength = getLength();
2217      DataArrayView::ShapeType dataPointShape = getDataPointShape();
2218      vector<int> referenceNumbers(noSamples);
2219      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2220        referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);
2221      }
2222      vector<int> tagNumbers(noSamples);
2223      if (isTagged()) {
2224        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2225          tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
2226        }
2227      }
2228    
2229      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2230      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2231      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2232    
2233      //
2234      // Flatten Shape to an array of integers suitable for writing to file
2235      int flatShape[4] = {0,0,0,0};
2236      cout << "\tshape: < ";
2237      for (int dim=0; dim<dataPointRank; dim++) {
2238        flatShape[dim] = dataPointShape[dim];
2239        cout << dataPointShape[dim] << " ";
2240      }
2241      cout << ">" << endl;
2242    
2243      //
2244      // Open archive file
2245      ofstream archiveFile;
2246      archiveFile.open(fileName.data(), ios::out);
2247    
2248      if (!archiveFile.good()) {
2249        throw DataException("archiveData Error: problem opening archive file");
2250      }
2251    
2252      //
2253      // Write common data items to archive file
2254      archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
2255      archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
2256      archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2257      archiveFile.write(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2258      archiveFile.write(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2259      archiveFile.write(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2260      archiveFile.write(reinterpret_cast<char *>(&dataLength),sizeof(int));
2261      for (int dim = 0; dim < 4; dim++) {
2262        archiveFile.write(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2263      }
2264      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2265        archiveFile.write(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2266      }
2267      if (isTagged()) {
2268        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2269          archiveFile.write(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2270        }
2271      }
2272    
2273      if (!archiveFile.good()) {
2274        throw DataException("archiveData Error: problem writing to archive file");
2275      }
2276    
2277      //
2278      // Archive underlying data values for each Data type
2279      int noValues;
2280      switch (dataType) {
2281        case 0:
2282          // DataEmpty
2283          noValues = 0;
2284          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2285          cout << "\tnoValues: " << noValues << endl;
2286          break;
2287        case 1:
2288          // DataConstant
2289          noValues = m_data->getLength();
2290          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2291          cout << "\tnoValues: " << noValues << endl;
2292          if (m_data->archiveData(archiveFile,noValues)) {
2293            throw DataException("archiveData Error: problem writing data to archive file");
2294          }
2295          break;
2296        case 2:
2297          // DataTagged
2298          noValues = m_data->getLength();
2299          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2300          cout << "\tnoValues: " << noValues << endl;
2301          if (m_data->archiveData(archiveFile,noValues)) {
2302            throw DataException("archiveData Error: problem writing data to archive file");
2303          }
2304          break;
2305        case 3:
2306          // DataExpanded
2307          noValues = m_data->getLength();
2308          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
2309          cout << "\tnoValues: " << noValues << endl;
2310          if (m_data->archiveData(archiveFile,noValues)) {
2311            throw DataException("archiveData Error: problem writing data to archive file");
2312          }
2313          break;
2314      }
2315    
2316      if (!archiveFile.good()) {
2317        throw DataException("archiveData Error: problem writing data to archive file");
2318      }
2319    
2320      //
2321      // Close archive file
2322      archiveFile.close();
2323    
2324      if (!archiveFile.good()) {
2325        throw DataException("archiveData Error: problem closing archive file");
2326      }
2327    
2328    }
2329    
2330    void
2331    Data::extractData(const std::string fileName,
2332                      const FunctionSpace& fspace)
2333    {
2334      //
2335      // Can only extract Data to an object which is initially DataEmpty
2336      if (!isEmpty()) {
2337        throw DataException("extractData Error: can only extract to DataEmpty object");
2338      }
2339    
2340      cout << "Extracting Data object from: " << fileName << endl;
2341    
2342      int dataType;
2343      int noSamples;
2344      int noDPPSample;
2345      int functionSpaceType;
2346      int dataPointRank;
2347      int dataPointSize;
2348      int dataLength;
2349      DataArrayView::ShapeType dataPointShape;
2350      int flatShape[4];
2351    
2352      //
2353      // Open the archive file
2354      ifstream archiveFile;
2355      archiveFile.open(fileName.data(), ios::in);
2356    
2357      if (!archiveFile.good()) {
2358        throw DataException("extractData Error: problem opening archive file");
2359      }
2360    
2361      //
2362      // Read common data items from archive file
2363      archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
2364      archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
2365      archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
2366      archiveFile.read(reinterpret_cast<char *>(&functionSpaceType),sizeof(int));
2367      archiveFile.read(reinterpret_cast<char *>(&dataPointRank),sizeof(int));
2368      archiveFile.read(reinterpret_cast<char *>(&dataPointSize),sizeof(int));
2369      archiveFile.read(reinterpret_cast<char *>(&dataLength),sizeof(int));
2370      for (int dim = 0; dim < 4; dim++) {
2371        archiveFile.read(reinterpret_cast<char *>(&flatShape[dim]),sizeof(int));
2372        if (flatShape[dim]>0) {
2373          dataPointShape.push_back(flatShape[dim]);
2374        }
2375      }
2376      vector<int> referenceNumbers(noSamples);
2377      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2378        archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2379      }
2380      vector<int> tagNumbers(noSamples);
2381      if (dataType==2) {
2382        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2383          archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
2384        }
2385      }
2386    
2387      if (!archiveFile.good()) {
2388        throw DataException("extractData Error: problem reading from archive file");
2389      }
2390    
2391      //
2392      // Verify the values just read from the archive file
2393      switch (dataType) {
2394        case 0:
2395          cout << "\tdataType: DataEmpty" << endl;
2396          break;
2397        case 1:
2398          cout << "\tdataType: DataConstant" << endl;
2399          break;
2400        case 2:
2401          cout << "\tdataType: DataTagged" << endl;
2402          break;
2403        case 3:
2404          cout << "\tdataType: DataExpanded" << endl;
2405          break;
2406        default:
2407          throw DataException("extractData Error: undefined dataType read from archive file");
2408          break;
2409      }
2410    
2411      cout << "\tnoSamples: " << noSamples << " noDPPSample: " << noDPPSample << endl;
2412      cout << "\tfunctionSpaceType: " << functionSpaceType << endl;
2413      cout << "\trank: " << dataPointRank << " size: " << dataPointSize << " length: " << dataLength << endl;
2414      cout << "\tshape: < ";
2415      for (int dim = 0; dim < dataPointRank; dim++) {
2416        cout << dataPointShape[dim] << " ";
2417      }
2418      cout << ">" << endl;
2419    
2420      //
2421      // Verify that supplied FunctionSpace object is compatible with this Data object.
2422      if ( (fspace.getTypeCode()!=functionSpaceType) ||
2423           (fspace.getNumSamples()!=noSamples) ||
2424           (fspace.getNumDPPSample()!=noDPPSample)
2425         ) {
2426        throw DataException("extractData Error: incompatible FunctionSpace");
2427      }
2428      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2429        if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {
2430          throw DataException("extractData Error: incompatible FunctionSpace");
2431        }
2432      }
2433      if (dataType==2) {
2434        for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2435          if (tagNumbers[sampleNo] != fspace.getTagFromSampleNo(sampleNo)) {
2436            throw DataException("extractData Error: incompatible FunctionSpace");
2437          }
2438        }
2439      }
2440    
2441      //
2442      // Construct a DataVector to hold underlying data values
2443      DataVector dataVec(dataLength);
2444    
2445      //
2446      // Load this DataVector with the appropriate values
2447      int noValues;
2448      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2449      cout << "\tnoValues: " << noValues << endl;
2450      switch (dataType) {
2451        case 0:
2452          // DataEmpty
2453          if (noValues != 0) {
2454            throw DataException("extractData Error: problem reading data from archive file");
2455          }
2456          break;
2457        case 1:
2458          // DataConstant
2459          if (dataVec.extractData(archiveFile,noValues)) {
2460            throw DataException("extractData Error: problem reading data from archive file");
2461          }
2462          break;
2463        case 2:
2464          // DataTagged
2465          if (dataVec.extractData(archiveFile,noValues)) {
2466            throw DataException("extractData Error: problem reading data from archive file");
2467          }
2468          break;
2469        case 3:
2470          // DataExpanded
2471          if (dataVec.extractData(archiveFile,noValues)) {
2472            throw DataException("extractData Error: problem reading data from archive file");
2473          }
2474          break;
2475      }
2476    
2477      if (!archiveFile.good()) {
2478        throw DataException("extractData Error: problem reading from archive file");
2479      }
2480    
2481      //
2482      // Close archive file
2483      archiveFile.close();
2484    
2485      if (!archiveFile.good()) {
2486        throw DataException("extractData Error: problem closing archive file");
2487      }
2488    
2489      //
2490      // Construct an appropriate Data object
2491      DataAbstract* tempData;
2492      switch (dataType) {
2493        case 0:
2494          // DataEmpty
2495          tempData=new DataEmpty();
2496          break;
2497        case 1:
2498          // DataConstant
2499          tempData=new DataConstant(fspace,dataPointShape,dataVec);
2500          break;
2501        case 2:
2502          // DataTagged
2503          tempData=new DataTagged(fspace,dataPointShape,tagNumbers,dataVec);
2504          break;
2505        case 3:
2506          // DataExpanded
2507          tempData=new DataExpanded(fspace,dataPointShape,dataVec);
2508          break;
2509      }
2510      shared_ptr<DataAbstract> temp_data(tempData);
2511      m_data=temp_data;
2512    }
2513    
2514  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
2515  {  {
2516    o << data.toString();    o << data.toString();
2517    return o;    return o;
2518  }  }
2519    
2520    Data
2521    escript::C_GeneralTensorProduct(Data& arg_0,
2522                         Data& arg_1,
2523                         int axis_offset,
2524                         int transpose)
2525    {
2526      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2527      // SM is the product of the last axis_offset entries in arg_0.getShape().
2528    
2529      #if defined DOPROF
2530        // profData->binary++;
2531      #endif
2532    
2533      // Interpolate if necessary and find an appropriate function space
2534      Data arg_0_Z, arg_1_Z;
2535      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2536        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2537          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2538          arg_1_Z = Data(arg_1);
2539        }
2540        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2541          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2542          arg_0_Z =Data(arg_0);
2543        }
2544        else {
2545          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2546        }
2547      } else {
2548          arg_0_Z = Data(arg_0);
2549          arg_1_Z = Data(arg_1);
2550      }
2551      // Get rank and shape of inputs
2552      int rank0 = arg_0_Z.getDataPointRank();
2553      int rank1 = arg_1_Z.getDataPointRank();
2554      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2555      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2556    
2557      // Prepare for the loops of the product and verify compatibility of shapes
2558      int start0=0, start1=0;
2559      if (transpose == 0)       {}
2560      else if (transpose == 1)  { start0 = axis_offset; }
2561      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2562      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2563    
2564      // Adjust the shapes for transpose
2565      DataArrayView::ShapeType tmpShape0;
2566      DataArrayView::ShapeType tmpShape1;
2567      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2568      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2569    
2570    #if 0
2571      // For debugging: show shape after transpose
2572      char tmp[100];
2573      std::string shapeStr;
2574      shapeStr = "(";
2575      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2576      shapeStr += ")";
2577      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2578      shapeStr = "(";
2579      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2580      shapeStr += ")";
2581      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2582    #endif
2583    
2584      // Prepare for the loops of the product
2585      int SL=1, SM=1, SR=1;
2586      for (int i=0; i<rank0-axis_offset; i++)   {
2587        SL *= tmpShape0[i];
2588      }
2589      for (int i=rank0-axis_offset; i<rank0; i++)   {
2590        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2591          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2592        }
2593        SM *= tmpShape0[i];
2594      }
2595      for (int i=axis_offset; i<rank1; i++)     {
2596        SR *= tmpShape1[i];
2597      }
2598    
2599      // Define the shape of the output
2600      DataArrayView::ShapeType shape2;
2601      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2602      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2603    
2604      // Declare output Data object
2605      Data res;
2606    
2607      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2608        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2609        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2610        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2611        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2612        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2613      }
2614      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2615    
2616        // Prepare the DataConstant input
2617        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2618        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2619    
2620        // Borrow DataTagged input from Data object
2621        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2622        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2623    
2624        // Prepare a DataTagged output 2
2625        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2626        res.tag();
2627        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2628        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2629    
2630        // Prepare offset into DataConstant
2631        int offset_0 = tmp_0->getPointOffset(0,0);
2632        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2633        // Get the views
2634        DataArrayView view_1 = tmp_1->getDefaultValue();
2635        DataArrayView view_2 = tmp_2->getDefaultValue();
2636        // Get the pointers to the actual data
2637        double *ptr_1 = &((view_1.getData())[0]);
2638        double *ptr_2 = &((view_2.getData())[0]);
2639        // Compute an MVP for the default
2640        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2641        // Compute an MVP for each tag
2642        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2643        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2644        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2645          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2646          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2647          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2648          double *ptr_1 = &view_1.getData(0);
2649          double *ptr_2 = &view_2.getData(0);
2650          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2651        }
2652    
2653      }
2654      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2655    
2656        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2657        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2658        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2659        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2660        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2661        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2662        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2663        int sampleNo_1,dataPointNo_1;
2664        int numSamples_1 = arg_1_Z.getNumSamples();
2665        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2666        int offset_0 = tmp_0->getPointOffset(0,0);
2667        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2668        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2669          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2670            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2671            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2672            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2673            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2674            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2675            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2676          }
2677        }
2678    
2679      }
2680      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2681    
2682        // Borrow DataTagged input from Data object
2683        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2684        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2685    
2686        // Prepare the DataConstant input
2687        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2688        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2689    
2690        // Prepare a DataTagged output 2
2691        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2692        res.tag();
2693        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2694        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2695    
2696        // Prepare offset into DataConstant
2697        int offset_1 = tmp_1->getPointOffset(0,0);
2698        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2699        // Get the views
2700        DataArrayView view_0 = tmp_0->getDefaultValue();
2701        DataArrayView view_2 = tmp_2->getDefaultValue();
2702        // Get the pointers to the actual data
2703        double *ptr_0 = &((view_0.getData())[0]);
2704        double *ptr_2 = &((view_2.getData())[0]);
2705        // Compute an MVP for the default
2706        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2707        // Compute an MVP for each tag
2708        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2709        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2710        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2711          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2712          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2713          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2714          double *ptr_0 = &view_0.getData(0);
2715          double *ptr_2 = &view_2.getData(0);
2716          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2717        }
2718    
2719      }
2720      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2721    
2722        // Borrow DataTagged input from Data object
2723        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2724        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2725    
2726        // Borrow DataTagged input from Data object
2727        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2728        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2729    
2730        // Prepare a DataTagged output 2
2731        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2732        res.tag();  // DataTagged output
2733        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2734        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2735    
2736        // Get the views
2737        DataArrayView view_0 = tmp_0->getDefaultValue();
2738        DataArrayView view_1 = tmp_1->getDefaultValue();
2739        DataArrayView view_2 = tmp_2->getDefaultValue();
2740        // Get the pointers to the actual data
2741        double *ptr_0 = &((view_0.getData())[0]);
2742        double *ptr_1 = &((view_1.getData())[0]);
2743        double *ptr_2 = &((view_2.getData())[0]);
2744        // Compute an MVP for the default
2745        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2746        // Merge the tags
2747        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2748        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2749        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2750        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2751          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2752        }
2753        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2754          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2755        }
2756        // Compute an MVP for each tag
2757        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2758        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2759          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2760          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2761          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2762          double *ptr_0 = &view_0.getData(0);
2763          double *ptr_1 = &view_1.getData(0);
2764          double *ptr_2 = &view_2.getData(0);
2765          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2766        }
2767    
2768      }
2769      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2770    
2771        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2772        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2773        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2774        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2775        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2776        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2777        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2778        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2779        int sampleNo_0,dataPointNo_0;
2780        int numSamples_0 = arg_0_Z.getNumSamples();
2781        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2782        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2783        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2784          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2785          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2786          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2787            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2788            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2789            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2790            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2791            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2792          }
2793        }
2794    
2795      }
2796      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2797    
2798        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2799        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2800        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2801        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2802        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2803        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2804        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2805        int sampleNo_0,dataPointNo_0;
2806        int numSamples_0 = arg_0_Z.getNumSamples();
2807        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2808        int offset_1 = tmp_1->getPointOffset(0,0);
2809        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2810        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2811          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2812            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2813            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2814            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2815            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2816            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2817            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2818          }
2819        }
2820    
2821    
2822      }
2823      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2824    
2825        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2826        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2827        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2828        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2829        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2830        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2831        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2832        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2833        int sampleNo_0,dataPointNo_0;
2834        int numSamples_0 = arg_0_Z.getNumSamples();
2835        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2836        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2839          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2840          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2841            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2842            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2843            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2844            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2845            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2846          }
2847        }
2848    
2849      }
2850      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2851    
2852        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2853        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2856        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2858        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2859        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2860        int sampleNo_0,dataPointNo_0;
2861        int numSamples_0 = arg_0_Z.getNumSamples();
2862        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2863        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2864        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2865          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2866            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2867            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2868            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2869            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2870            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2871            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2872            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2873          }
2874        }
2875    
2876      }
2877      else {
2878        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2879      }
2880    
2881      return res;
2882    }
2883    
2884    DataAbstract*
2885    Data::borrowData() const
2886    {
2887      return m_data.get();
2888    }
2889    
2890    /* Member functions specific to the MPI implementation */
2891    
2892    void
2893    Data::print()
2894    {
2895      int i,j;
2896      
2897      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2898      for( i=0; i<getNumSamples(); i++ )
2899      {
2900        printf( "[%6d]", i );
2901        for( j=0; j<getNumDataPointsPerSample(); j++ )
2902          printf( "\t%10.7g", (getSampleData(i))[j] );
2903        printf( "\n" );
2904      }
2905    }
2906    
2907    int
2908    Data::get_MPISize() const
2909    {
2910        int error, size;
2911    #ifdef PASO_MPI
2912        error = MPI_Comm_size( get_MPIComm(), &size );
2913    #else
2914        size = 1;
2915    #endif
2916        return size;
2917    }
2918    
2919    int
2920    Data::get_MPIRank() const
2921    {
2922        int error, rank;
2923    #ifdef PASO_MPI
2924        error = MPI_Comm_rank( get_MPIComm(), &rank );
2925    #else
2926        rank = 0;
2927    #endif
2928        return rank;
2929    }
2930    
2931    MPI_Comm
2932    Data::get_MPIComm() const
2933    {
2934    #ifdef PASO_MPI
2935        return MPI_COMM_WORLD;
2936    #else
2937        return -1;
2938    #endif
2939    }
2940    

Legend:
Removed from v.108  
changed lines
  Added in v.922

  ViewVC Help
Powered by ViewVC 1.1.26