/[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 122 by jgs, Thu Jun 9 05:38:05 2005 UTC trunk/escript/src/Data.cpp revision 1092 by gross, Fri Apr 13 03:39:49 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 "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    
 ******************************************************************************/  
   
 #include "escript/Data/Data.h"  
   
 #include <iostream>  
27  #include <fstream>  #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/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;
# Line 51  Data::Data() Line 45  Data::Data()
45    DataAbstract* temp=new DataEmpty();    DataAbstract* temp=new DataEmpty();
46    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
47    m_data=temp_data;    m_data=temp_data;
48      m_protected=false;
49  }  }
50    
51  Data::Data(double value,  Data::Data(double value,
# Line 64  Data::Data(double value, Line 59  Data::Data(double value,
59    }    }
60    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
61    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
62      m_protected=false;
63  }  }
64    
65  Data::Data(double value,  Data::Data(double value,
# Line 74  Data::Data(double value, Line 70  Data::Data(double value,
70    DataArray temp(dataPointShape,value);    DataArray temp(dataPointShape,value);
71    pair<int,int> dataShape=what.getDataShape();    pair<int,int> dataShape=what.getDataShape();
72    initialise(temp.getView(),what,expanded);    initialise(temp.getView(),what,expanded);
73      m_protected=false;
74  }  }
75    
76  Data::Data(const Data& inData)  Data::Data(const Data& inData)
77  {  {
78    m_data=inData.m_data;    m_data=inData.m_data;
79      m_protected=inData.isProtected();
80  }  }
81    
82  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 89  Data::Data(const Data& inData, Line 87  Data::Data(const Data& inData,
87    DataAbstract* tmp = inData.m_data->getSlice(region);    DataAbstract* tmp = inData.m_data->getSlice(region);
88    shared_ptr<DataAbstract> temp_data(tmp);    shared_ptr<DataAbstract> temp_data(tmp);
89    m_data=temp_data;    m_data=temp_data;
90      m_protected=false;
91  }  }
92    
93  Data::Data(const Data& inData,  Data::Data(const Data& inData,
# Line 98  Data::Data(const Data& inData, Line 97  Data::Data(const Data& inData,
97      m_data=inData.m_data;      m_data=inData.m_data;
98    } else {    } else {
99      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);      Data tmp(0,inData.getPointDataView().getShape(),functionspace,true);
100      // Note for Lutz, Must use a reference or pointer to a derived object      // Note: Must use a reference or pointer to a derived object
101      // in order to get polymorphic behaviour. Shouldn't really      // in order to get polymorphic behaviour. Shouldn't really
102      // be able to create an instance of AbstractDomain but that was done      // be able to create an instance of AbstractDomain but that was done
103      // as a boost python work around which may no longer be required.      // as a boost:python work around which may no longer be required.
104      const AbstractDomain& inDataDomain=inData.getDomain();      const AbstractDomain& inDataDomain=inData.getDomain();
105      if  (inDataDomain==functionspace.getDomain()) {      if  (inDataDomain==functionspace.getDomain()) {
106        inDataDomain.interpolateOnDomain(tmp,inData);        inDataDomain.interpolateOnDomain(tmp,inData);
# Line 110  Data::Data(const Data& inData, Line 109  Data::Data(const Data& inData,
109      }      }
110      m_data=tmp.m_data;      m_data=tmp.m_data;
111    }    }
112      m_protected=false;
113  }  }
114    
115  Data::Data(const DataTagged::TagListType& tagKeys,  Data::Data(const DataTagged::TagListType& tagKeys,
# Line 121  Data::Data(const DataTagged::TagListType Line 121  Data::Data(const DataTagged::TagListType
121    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);    DataAbstract* temp=new DataTagged(tagKeys,values,defaultValue,what);
122    shared_ptr<DataAbstract> temp_data(temp);    shared_ptr<DataAbstract> temp_data(temp);
123    m_data=temp_data;    m_data=temp_data;
124      m_protected=false;
125    if (expanded) {    if (expanded) {
126      expand();      expand();
127    }    }
# Line 131  Data::Data(const numeric::array& value, Line 132  Data::Data(const numeric::array& value,
132             bool expanded)             bool expanded)
133  {  {
134    initialise(value,what,expanded);    initialise(value,what,expanded);
135      m_protected=false;
136  }  }
137    
138  Data::Data(const DataArrayView& value,  Data::Data(const DataArrayView& value,
# Line 138  Data::Data(const DataArrayView& value, Line 140  Data::Data(const DataArrayView& value,
140             bool expanded)             bool expanded)
141  {  {
142    initialise(value,what,expanded);    initialise(value,what,expanded);
143      m_protected=false;
144  }  }
145    
146  Data::Data(const object& value,  Data::Data(const object& value,
# Line 146  Data::Data(const object& value, Line 149  Data::Data(const object& value,
149  {  {
150    numeric::array asNumArray(value);    numeric::array asNumArray(value);
151    initialise(asNumArray,what,expanded);    initialise(asNumArray,what,expanded);
152      m_protected=false;
153  }  }
154    
155  Data::Data(const object& value,  Data::Data(const object& value,
# Line 166  Data::Data(const object& value, Line 170  Data::Data(const object& value,
170      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
171      initialise(temp.getView(),other.getFunctionSpace(),false);      initialise(temp.getView(),other.getFunctionSpace(),false);
172    }    }
173      m_protected=false;
174    }
175    
176    Data::~Data()
177    {
178    
179  }  }
180    
181  escriptDataC  escriptDataC
# Line 302  Data::isConstant() const Line 312  Data::isConstant() const
312  }  }
313    
314  void  void
315    Data::setProtection()
316    {
317       m_protected=true;
318    }
319    
320    bool
321    Data::isProtected() const
322    {
323       return m_protected;
324    }
325    
326    
327    
328    void
329  Data::expand()  Data::expand()
330  {  {
331    if (isConstant()) {    if (isConstant()) {
# Line 343  Data::tag() Line 367  Data::tag()
367    }    }
368  }  }
369    
370  void  Data
371  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
372  {  {
373    m_data->reshapeDataPoint(shape);    return escript::unaryOp(*this,bind1st(divides<double>(),1.));
374  }  }
375    
376  Data  Data
# Line 374  Data::whereNonPositive() const Line 398  Data::whereNonPositive() const
398  }  }
399    
400  Data  Data
401  Data::whereZero() const  Data::whereZero(double tol) const
402  {  {
403    return escript::unaryOp(*this,bind2nd(equal_to<double>(),0.0));    Data dataAbs=abs();
404      return escript::unaryOp(dataAbs,bind2nd(less_equal<double>(),tol));
405  }  }
406    
407  Data  Data
408  Data::whereNonZero() const  Data::whereNonZero(double tol) const
409  {  {
410    return escript::unaryOp(*this,bind2nd(not_equal_to<double>(),0.0));    Data dataAbs=abs();
411      return escript::unaryOp(dataAbs,bind2nd(greater<double>(),tol));
412  }  }
413    
414  Data  Data
# Line 442  Data::getDataPointShape() const Line 468  Data::getDataPointShape() const
468    return getPointDataView().getShape();    return getPointDataView().getShape();
469  }  }
470    
471  const  
472    
473    const
474  boost::python::numeric::array  boost::python::numeric::array
475  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
476  {  {
477    //    size_t length=0;
478    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
479    //    //
480    // determine the rank and shape of each data point    // determine the rank and shape of each data point
481    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 462  Data::convertToNumArray() Line 486  Data::convertToNumArray()
486    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
487    
488    //    //
489    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
490    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
491    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
492    // data points, whilst the first dimension will be equal to the    DataArrayView::ShapeType arrayShape = dataPointShape;
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPoints);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
493    
494    //    //
495    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
496      if (arrayRank==0) {
497        numArray.resize(1);
498      }
499    if (arrayRank==1) {    if (arrayRank==1) {
500      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
501    }    }
# Line 489  Data::convertToNumArray() Line 508  Data::convertToNumArray()
508    if (arrayRank==4) {    if (arrayRank==4) {
509      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
510    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
511    
512    //    if (getNumDataPointsPerSample()>0) {
513    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
514    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
515    int dataPoint = 0;         //
516    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
517      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
518        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
519        if (dataPointRank==0) {         }
520          numArray[dataPoint]=dataPointView();                
521        }         //
522        if (dataPointRank==1) {         // Check a valid data point number has been supplied
523          for (int i=0; i<dataPointShape[0]; i++) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
524            numArray[dataPoint][i]=dataPointView(i);             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
525          }         }
526        }         // TODO: global error handling
527        if (dataPointRank==2) {         // create a view of the data if it is stored locally
528          for (int i=0; i<dataPointShape[0]; i++) {         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
529            for (int j=0; j<dataPointShape[1]; j++) {          
530              numArray[dataPoint][i][j] = dataPointView(i,j);         switch( dataPointRank ){
531            }              case 0 :
532          }                  numArray[0] = dataPointView();
533        }                  break;
534        if (dataPointRank==3) {              case 1 :        
535          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
536            for (int j=0; j<dataPointShape[1]; j++) {                      numArray[i]=dataPointView(i);
537              for (int k=0; k<dataPointShape[2]; k++) {                  break;
538                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);              case 2 :        
539              }                  for( i=0; i<dataPointShape[0]; i++ )
540            }                      for( j=0; j<dataPointShape[1]; j++)
541          }                          numArray[make_tuple(i,j)]=dataPointView(i,j);
542        }                  break;
543        if (dataPointRank==4) {              case 3 :        
544          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
545            for (int j=0; j<dataPointShape[1]; j++) {                      for( j=0; j<dataPointShape[1]; j++ )
546              for (int k=0; k<dataPointShape[2]; k++) {                          for( k=0; k<dataPointShape[2]; k++)
547                for (int l=0; l<dataPointShape[3]; l++) {                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
548                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  break;
549                }              case 4 :
550              }                  for( i=0; i<dataPointShape[0]; i++ )
551            }                      for( j=0; j<dataPointShape[1]; j++ )
552          }                          for( k=0; k<dataPointShape[2]; k++ )
553        }                              for( l=0; l<dataPointShape[3]; l++)
554        dataPoint++;                                  numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
555      }                  break;
556        }
557    }    }
   
558    //    //
559    // return the loaded array    // return the array
560    return numArray;    return numArray;
 }  
561    
562  const  }
563  boost::python::numeric::array  void
564  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
565  {  {
566    //      // this will throw if the value cannot be represented
567    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
568    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
569    
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
570    
571    //  }
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
572    
573    void
574    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
575    {
576      if (isProtected()) {
577            throw DataException("Error - attempt to update protected Data object.");
578      }
579    //    //
580    // create the numeric array to be returned    // check rank
581    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
582          throw DataException("Rank of numarray does not match Data object rank");
583    
584    //    //
585    // the rank of the returned numeric array will be the rank of    // check shape of num_array
586    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
587    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
588    // data points, whilst the first dimension will be equal to the         throw DataException("Shape of numarray does not match Data object rank");
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
589    }    }
   
590    //    //
591    // resize the numeric array to the shape just calculated    // make sure data is expanded:
592    if (arrayRank==1) {    if (!isExpanded()) {
593      numArray.resize(arrayShape[0]);      expand();
   }  
   if (arrayRank==2) {  
     numArray.resize(arrayShape[0],arrayShape[1]);  
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
   }  
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
594    }    }
595      if (getNumDataPointsPerSample()>0) {
596    //         int sampleNo = dataPointNo/getNumDataPointsPerSample();
597    // loop through each data point in turn, loading the values for that data point         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
598    // into the numeric array.         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
599    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    } else {
600      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         m_data->copyToDataPoint(-1, 0,num_array);
     if (dataPointRank==0) {  
       numArray[dataPoint]=dataPointView();  
     }  
     if (dataPointRank==1) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         numArray[dataPoint][i]=dataPointView(i);  
       }  
     }  
     if (dataPointRank==2) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           numArray[dataPoint][i][j] = dataPointView(i,j);  
         }  
       }  
     }  
     if (dataPointRank==3) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             numArray[dataPoint][i][j][k]=dataPointView(i,j,k);  
           }  
         }  
       }  
     }  
     if (dataPointRank==4) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           for (int k=0; k<dataPointShape[2]; k++) {  
             for (int l=0; l<dataPointShape[3]; l++) {  
               numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
             }  
           }  
         }  
       }  
     }  
601    }    }
   
   //  
   // return the loaded array  
   return numArray;  
602  }  }
603    
604  const  void
605  boost::python::numeric::array  Data::setValueOfDataPoint(int dataPointNo, const double value)
 Data::convertToNumArrayFromDPNo(int sampleNo,  
                                 int dataPointNo)  
606  {  {
607    //    if (isProtected()) {
608    // Check a valid sample number has been supplied          throw DataException("Error - attempt to update protected Data object.");
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
609    }    }
   
610    //    //
611    // Check a valid data point number has been supplied    // make sure data is expanded:
612    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
613      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
614      }
615      if (getNumDataPointsPerSample()>0) {
616           int sampleNo = dataPointNo/getNumDataPointsPerSample();
617           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
618           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
619      } else {
620           m_data->copyToDataPoint(-1, 0,value);
621    }    }
622    }
623    
624    const
625    boost::python::numeric::array
626    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
627    {
628      size_t length=0;
629      int i, j, k, l, pos;
630    //    //
631    // determine the rank and shape of each data point    // determine the rank and shape of each data point
632    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 696  Data::convertToNumArrayFromDPNo(int samp Line 660  Data::convertToNumArrayFromDPNo(int samp
660      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
661    }    }
662    
663      // added for the MPI communication
664      length=1;
665      for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
666      double *tmpData = new double[length];
667    
668    //    //
669    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
   DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
670    
671        // updated for the MPI case
672        if( get_MPIRank()==procNo ){
673                 if (getNumDataPointsPerSample()>0) {
674                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
675                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
676                    //
677                    // Check a valid sample number has been supplied
678                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
679                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
680                    }
681                  
682                    //
683                    // Check a valid data point number has been supplied
684                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
685                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
686                    }
687                    // TODO: global error handling
688            // create a view of the data if it is stored locally
689            DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
690            
691            // pack the data from the view into tmpData for MPI communication
692            pos=0;
693            switch( dataPointRank ){
694                case 0 :
695                    tmpData[0] = dataPointView();
696                    break;
697                case 1 :        
698                    for( i=0; i<dataPointShape[0]; i++ )
699                        tmpData[i]=dataPointView(i);
700                    break;
701                case 2 :        
702                    for( i=0; i<dataPointShape[0]; i++ )
703                        for( j=0; j<dataPointShape[1]; j++, pos++ )
704                            tmpData[pos]=dataPointView(i,j);
705                    break;
706                case 3 :        
707                    for( i=0; i<dataPointShape[0]; i++ )
708                        for( j=0; j<dataPointShape[1]; j++ )
709                            for( k=0; k<dataPointShape[2]; k++, pos++ )
710                                tmpData[pos]=dataPointView(i,j,k);
711                    break;
712                case 4 :
713                    for( i=0; i<dataPointShape[0]; i++ )
714                        for( j=0; j<dataPointShape[1]; j++ )
715                            for( k=0; k<dataPointShape[2]; k++ )
716                                for( l=0; l<dataPointShape[3]; l++, pos++ )
717                                    tmpData[pos]=dataPointView(i,j,k,l);
718                    break;
719            }
720                }
721        }
722            #ifdef PASO_MPI
723            // broadcast the data to all other processes
724        MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
725            #endif
726    
727        // unpack the data
728        switch( dataPointRank ){
729            case 0 :
730                numArray[0]=tmpData[0];
731                break;
732            case 1 :        
733                for( i=0; i<dataPointShape[0]; i++ )
734                    numArray[i]=tmpData[i];
735                break;
736            case 2 :        
737                for( i=0; i<dataPointShape[0]; i++ )
738                    for( j=0; j<dataPointShape[1]; j++ )
739                       numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
740                break;
741            case 3 :        
742                for( i=0; i<dataPointShape[0]; i++ )
743                    for( j=0; j<dataPointShape[1]; j++ )
744                        for( k=0; k<dataPointShape[2]; k++ )
745                            numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
746                break;
747            case 4 :
748                for( i=0; i<dataPointShape[0]; i++ )
749                    for( j=0; j<dataPointShape[1]; j++ )
750                        for( k=0; k<dataPointShape[2]; k++ )
751                            for( l=0; l<dataPointShape[3]; l++ )
752                                    numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
753                break;
754        }
755    
756        delete [] tmpData;  
757    //    //
758    // return the loaded array    // return the loaded array
759    return numArray;    return numArray;
760  }  }
761    
762    
763    
764  boost::python::numeric::array  boost::python::numeric::array
765  Data::integrate() const  Data::integrate() const
766  {  {
# Line 747  Data::integrate() const Line 768  Data::integrate() const
768    int rank = getDataPointRank();    int rank = getDataPointRank();
769    DataArrayView::ShapeType shape = getDataPointShape();    DataArrayView::ShapeType shape = getDataPointShape();
770    
771    
772    //    //
773    // calculate the integral values    // calculate the integral values
774    vector<double> integrals(getDataPointSize());    vector<double> integrals(getDataPointSize());
# Line 769  Data::integrate() const Line 791  Data::integrate() const
791      }      }
792    }    }
793    if (rank==2) {    if (rank==2) {
794      bp_array.resize(shape[0],shape[1]);         bp_array.resize(shape[0],shape[1]);
795      for (int i=0; i<shape[0]; i++) {         for (int i=0; i<shape[0]; i++) {
796        for (int j=0; j<shape[1]; j++) {           for (int j=0; j<shape[1]; j++) {
797          index = i + shape[0] * j;             index = i + shape[0] * j;
798          bp_array[i,j] = integrals[index];             bp_array[make_tuple(i,j)] = integrals[index];
799        }           }
800      }         }
801    }    }
802    if (rank==3) {    if (rank==3) {
803      bp_array.resize(shape[0],shape[1],shape[2]);      bp_array.resize(shape[0],shape[1],shape[2]);
# Line 783  Data::integrate() const Line 805  Data::integrate() const
805        for (int j=0; j<shape[1]; j++) {        for (int j=0; j<shape[1]; j++) {
806          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
807            index = i + shape[0] * ( j + shape[1] * k );            index = i + shape[0] * ( j + shape[1] * k );
808            bp_array[i,j,k] = integrals[index];            bp_array[make_tuple(i,j,k)] = integrals[index];
809          }          }
810        }        }
811      }      }
# Line 795  Data::integrate() const Line 817  Data::integrate() const
817          for (int k=0; k<shape[2]; k++) {          for (int k=0; k<shape[2]; k++) {
818            for (int l=0; l<shape[3]; l++) {            for (int l=0; l<shape[3]; l++) {
819              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );              index = i + shape[0] * ( j + shape[1] * ( k + shape[2] * l ) );
820              bp_array[i,j,k,l] = integrals[index];              bp_array[make_tuple(i,j,k,l)] = integrals[index];
821            }            }
822          }          }
823        }        }
# Line 826  Data::tan() const Line 848  Data::tan() const
848  }  }
849    
850  Data  Data
851  Data::log() const  Data::asin() const
852    {
853      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asin);
854    }
855    
856    Data
857    Data::acos() const
858    {
859      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
860    }
861    
862    
863    Data
864    Data::atan() const
865    {
866      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atan);
867    }
868    
869    Data
870    Data::sinh() const
871    {
872      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::sinh);
873    }
874    
875    Data
876    Data::cosh() const
877    {
878      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::cosh);
879    }
880    
881    Data
882    Data::tanh() const
883    {
884      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
885    }
886    
887    
888    Data
889    Data::erf() const
890    {
891    #ifdef _WIN32
892      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
893    #else
894      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
895    #endif
896    }
897    
898    Data
899    Data::asinh() const
900    {
901    #ifdef _WIN32
902      return escript::unaryOp(*this,escript::asinh_substitute);
903    #else
904      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
905    #endif
906    }
907    
908    Data
909    Data::acosh() const
910    {
911    #ifdef _WIN32
912      return escript::unaryOp(*this,escript::acosh_substitute);
913    #else
914      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
915    #endif
916    }
917    
918    Data
919    Data::atanh() const
920    {
921    #ifdef _WIN32
922      return escript::unaryOp(*this,escript::atanh_substitute);
923    #else
924      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
925    #endif
926    }
927    
928    Data
929    Data::log10() const
930  {  {
931    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log10);
932  }  }
933    
934  Data  Data
935  Data::ln() const  Data::log() const
936  {  {
937    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::log);
938  }  }
# Line 858  Data::neg() const Line 958  Data::neg() const
958  Data  Data
959  Data::pos() const  Data::pos() const
960  {  {
961    return (*this);    Data result;
962      // perform a deep copy
963      result.copy(*this);
964      return result;
965  }  }
966    
967  Data  Data
# Line 876  Data::sqrt() const Line 979  Data::sqrt() const
979  double  double
980  Data::Lsup() const  Data::Lsup() const
981  {  {
982      double localValue, globalValue;
983    //    //
984    // set the initial absolute maximum value to zero    // set the initial absolute maximum value to zero
985    return algorithm(DataAlgorithmAdapter<AbsMax>(0));  
986      AbsMax abs_max_func;
987      localValue = algorithm(abs_max_func,0);
988    #ifdef PASO_MPI
989      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
990      return globalValue;
991    #else
992      return localValue;
993    #endif
994  }  }
995    
996  double  double
997  Data::Linf() const  Data::Linf() const
998  {  {
999      double localValue, globalValue;
1000    //    //
1001    // set the initial absolute minimum value to max double    // set the initial absolute minimum value to max double
1002    return algorithm(DataAlgorithmAdapter<AbsMin>(numeric_limits<double>::max()));    AbsMin abs_min_func;
1003      localValue = algorithm(abs_min_func,numeric_limits<double>::max());
1004    
1005    #ifdef PASO_MPI
1006      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1007      return globalValue;
1008    #else
1009      return localValue;
1010    #endif
1011  }  }
1012    
1013  double  double
1014  Data::sup() const  Data::sup() const
1015  {  {
1016      double localValue, globalValue;
1017    //    //
1018    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1019    return algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1020      localValue = algorithm(fmax_func,numeric_limits<double>::max()*-1);
1021    #ifdef PASO_MPI
1022      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD );
1023      return globalValue;
1024    #else
1025      return localValue;
1026    #endif
1027  }  }
1028    
1029  double  double
1030  Data::inf() const  Data::inf() const
1031  {  {
1032      double localValue, globalValue;
1033    //    //
1034    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1035    return algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1036      localValue = algorithm(fmin_func,numeric_limits<double>::max());
1037    #ifdef PASO_MPI
1038      MPI_Allreduce( &localValue, &globalValue, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD );
1039      return globalValue;
1040    #else
1041      return localValue;
1042    #endif
1043  }  }
1044    
1045    /* TODO */
1046    /* global reduction */
1047  Data  Data
1048  Data::maxval() const  Data::maxval() const
1049  {  {
1050    //    //
1051    // set the initial maximum value to min possible double    // set the initial maximum value to min possible double
1052    return dp_algorithm(DataAlgorithmAdapter<FMax>(numeric_limits<double>::max()*-1));    FMax fmax_func;
1053      return dp_algorithm(fmax_func,numeric_limits<double>::max()*-1);
1054  }  }
1055    
1056  Data  Data
# Line 918  Data::minval() const Line 1058  Data::minval() const
1058  {  {
1059    //    //
1060    // set the initial minimum value to max possible double    // set the initial minimum value to max possible double
1061    return dp_algorithm(DataAlgorithmAdapter<FMin>(numeric_limits<double>::max()));    FMin fmin_func;
1062      return dp_algorithm(fmin_func,numeric_limits<double>::max());
1063  }  }
1064    
1065  const boost::python::tuple  Data
1066  Data::mindp() const  Data::swapaxes(const int axis0, const int axis1) const
1067  {  {
1068    Data temp=minval();       int axis0_tmp,axis1_tmp;
1069         DataArrayView::ShapeType s=getDataPointShape();
1070    int numSamples=temp.getNumSamples();       DataArrayView::ShapeType ev_shape;
1071    int numDPPSample=temp.getNumDataPointsPerSample();       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1072         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1073    int i,j,lowi=0,lowj=0;       int rank=getDataPointRank();
1074    double min=numeric_limits<double>::max();       if (rank<2) {
1075            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1076    for (i=0; i<numSamples; i++) {       }
1077      for (j=0; j<numDPPSample; j++) {       if (axis0<0 || axis0>rank-1) {
1078        double next=temp.getDataPoint(i,j)();          throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1079        if (next<min) {       }
1080          min=next;       if (axis1<0 || axis1>rank-1) {
1081          lowi=i;           throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1082          lowj=j;       }
1083        }       if (axis0 == axis1) {
1084      }           throw DataException("Error - Data::swapaxes: axis indices must be different.");
1085    }       }
1086         if (axis0 > axis1) {
1087             axis0_tmp=axis1;
1088             axis1_tmp=axis0;
1089         } else {
1090             axis0_tmp=axis0;
1091             axis1_tmp=axis1;
1092         }
1093         for (int i=0; i<rank; i++) {
1094           if (i == axis0_tmp) {
1095              ev_shape.push_back(s[axis1_tmp]);
1096           } else if (i == axis1_tmp) {
1097              ev_shape.push_back(s[axis0_tmp]);
1098           } else {
1099              ev_shape.push_back(s[i]);
1100           }
1101         }
1102         Data ev(0.,ev_shape,getFunctionSpace());
1103         ev.typeMatchRight(*this);
1104         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1105         return ev;
1106    
1107    }
1108    
1109    Data
1110    Data::symmetric() const
1111    {
1112         // check input
1113         DataArrayView::ShapeType s=getDataPointShape();
1114         if (getDataPointRank()==2) {
1115            if(s[0] != s[1])
1116               throw DataException("Error - Data::symmetric can only be calculated for rank 2 object with equal first and second dimension.");
1117         }
1118         else if (getDataPointRank()==4) {
1119            if(!(s[0] == s[2] && s[1] == s[3]))
1120               throw DataException("Error - Data::symmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1121         }
1122         else {
1123            throw DataException("Error - Data::symmetric can only be calculated for rank 2 or 4 object.");
1124         }
1125         Data ev(0.,getDataPointShape(),getFunctionSpace());
1126         ev.typeMatchRight(*this);
1127         m_data->symmetric(ev.m_data.get());
1128         return ev;
1129    }
1130    
1131    Data
1132    Data::nonsymmetric() const
1133    {
1134         // check input
1135         DataArrayView::ShapeType s=getDataPointShape();
1136         if (getDataPointRank()==2) {
1137            if(s[0] != s[1])
1138               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 object with equal first and second dimension.");
1139            DataArrayView::ShapeType ev_shape;
1140            ev_shape.push_back(s[0]);
1141            ev_shape.push_back(s[1]);
1142            Data ev(0.,ev_shape,getFunctionSpace());
1143            ev.typeMatchRight(*this);
1144            m_data->nonsymmetric(ev.m_data.get());
1145            return ev;
1146         }
1147         else if (getDataPointRank()==4) {
1148            if(!(s[0] == s[2] && s[1] == s[3]))
1149               throw DataException("Error - Data::nonsymmetric can only be calculated for rank 4 object with dim0==dim2 and dim1==dim3.");
1150            DataArrayView::ShapeType ev_shape;
1151            ev_shape.push_back(s[0]);
1152            ev_shape.push_back(s[1]);
1153            ev_shape.push_back(s[2]);
1154            ev_shape.push_back(s[3]);
1155            Data ev(0.,ev_shape,getFunctionSpace());
1156            ev.typeMatchRight(*this);
1157            m_data->nonsymmetric(ev.m_data.get());
1158            return ev;
1159         }
1160         else {
1161            throw DataException("Error - Data::nonsymmetric can only be calculated for rank 2 or 4 object.");
1162         }
1163    }
1164    
1165    Data
1166    Data::trace(int axis_offset) const
1167    {
1168         DataArrayView::ShapeType s=getDataPointShape();
1169         if (getDataPointRank()==2) {
1170            DataArrayView::ShapeType ev_shape;
1171            Data ev(0.,ev_shape,getFunctionSpace());
1172            ev.typeMatchRight(*this);
1173            m_data->trace(ev.m_data.get(), axis_offset);
1174            return ev;
1175         }
1176         if (getDataPointRank()==3) {
1177            DataArrayView::ShapeType ev_shape;
1178            if (axis_offset==0) {
1179              int s2=s[2];
1180              ev_shape.push_back(s2);
1181            }
1182            else if (axis_offset==1) {
1183              int s0=s[0];
1184              ev_shape.push_back(s0);
1185            }
1186            Data ev(0.,ev_shape,getFunctionSpace());
1187            ev.typeMatchRight(*this);
1188            m_data->trace(ev.m_data.get(), axis_offset);
1189            return ev;
1190         }
1191         if (getDataPointRank()==4) {
1192            DataArrayView::ShapeType ev_shape;
1193            if (axis_offset==0) {
1194              ev_shape.push_back(s[2]);
1195              ev_shape.push_back(s[3]);
1196            }
1197            else if (axis_offset==1) {
1198              ev_shape.push_back(s[0]);
1199              ev_shape.push_back(s[3]);
1200            }
1201        else if (axis_offset==2) {
1202          ev_shape.push_back(s[0]);
1203          ev_shape.push_back(s[1]);
1204        }
1205            Data ev(0.,ev_shape,getFunctionSpace());
1206            ev.typeMatchRight(*this);
1207        m_data->trace(ev.m_data.get(), axis_offset);
1208            return ev;
1209         }
1210         else {
1211            throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1212         }
1213    }
1214    
1215    return make_tuple(lowi,lowj);  Data
1216    Data::transpose(int axis_offset) const
1217    {
1218         DataArrayView::ShapeType s=getDataPointShape();
1219         DataArrayView::ShapeType ev_shape;
1220         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1221         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1222         int rank=getDataPointRank();
1223         if (axis_offset<0 || axis_offset>rank) {
1224            throw DataException("Error - Data::transpose must have 0 <= axis_offset <= rank=" + rank);
1225         }
1226         for (int i=0; i<rank; i++) {
1227           int index = (axis_offset+i)%rank;
1228           ev_shape.push_back(s[index]); // Append to new shape
1229         }
1230         Data ev(0.,ev_shape,getFunctionSpace());
1231         ev.typeMatchRight(*this);
1232         m_data->transpose(ev.m_data.get(), axis_offset);
1233         return ev;
1234  }  }
1235    
1236  Data  Data
1237  Data::length() const  Data::eigenvalues() const
1238    {
1239         // check input
1240         DataArrayView::ShapeType s=getDataPointShape();
1241         if (getDataPointRank()!=2)
1242            throw DataException("Error - Data::eigenvalues can only be calculated for rank 2 object.");
1243         if(s[0] != s[1])
1244            throw DataException("Error - Data::eigenvalues can only be calculated for object with equal first and second dimension.");
1245         // create return
1246         DataArrayView::ShapeType ev_shape(1,s[0]);
1247         Data ev(0.,ev_shape,getFunctionSpace());
1248         ev.typeMatchRight(*this);
1249         m_data->eigenvalues(ev.m_data.get());
1250         return ev;
1251    }
1252    
1253    const boost::python::tuple
1254    Data::eigenvalues_and_eigenvectors(const double tol) const
1255  {  {
1256    return dp_algorithm(DataAlgorithmAdapter<Length>(0));       DataArrayView::ShapeType s=getDataPointShape();
1257         if (getDataPointRank()!=2)
1258            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for rank 2 object.");
1259         if(s[0] != s[1])
1260            throw DataException("Error - Data::eigenvalues and eigenvectors can only be calculated for object with equal first and second dimension.");
1261         // create return
1262         DataArrayView::ShapeType ev_shape(1,s[0]);
1263         Data ev(0.,ev_shape,getFunctionSpace());
1264         ev.typeMatchRight(*this);
1265         DataArrayView::ShapeType V_shape(2,s[0]);
1266         Data V(0.,V_shape,getFunctionSpace());
1267         V.typeMatchRight(*this);
1268         m_data->eigenvalues_and_eigenvectors(ev.m_data.get(),V.m_data.get(),tol);
1269         return make_tuple(boost::python::object(ev),boost::python::object(V));
1270  }  }
1271    
1272  Data  const boost::python::tuple
1273  Data::trace() const  Data::minGlobalDataPoint() const
1274  {  {
1275    return dp_algorithm(DataAlgorithmAdapter<Trace>(0));    // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1276      // abort (for unknown reasons) if there are openmp directives with it in the
1277      // surrounding function
1278    
1279      int DataPointNo;
1280      int ProcNo;
1281      calc_minGlobalDataPoint(ProcNo,DataPointNo);
1282      return make_tuple(ProcNo,DataPointNo);
1283  }  }
1284    
1285  Data  void
1286  Data::transpose(int axis) const  Data::calc_minGlobalDataPoint(int& ProcNo,
1287                            int& DataPointNo) const
1288  {  {
1289    // not implemented    int i,j;
1290    throw DataException("Error - Data::transpose not implemented yet.");    int lowi=0,lowj=0;
1291    return Data();    double min=numeric_limits<double>::max();
1292    
1293      Data temp=minval();
1294    
1295      int numSamples=temp.getNumSamples();
1296      int numDPPSample=temp.getNumDataPointsPerSample();
1297    
1298      double next,local_min;
1299      int local_lowi,local_lowj;
1300    
1301      #pragma omp parallel private(next,local_min,local_lowi,local_lowj)
1302      {
1303        local_min=min;
1304        #pragma omp for private(i,j) schedule(static)
1305        for (i=0; i<numSamples; i++) {
1306          for (j=0; j<numDPPSample; j++) {
1307            next=temp.getDataPoint(i,j)();
1308            if (next<local_min) {
1309              local_min=next;
1310              local_lowi=i;
1311              local_lowj=j;
1312            }
1313          }
1314        }
1315        #pragma omp critical
1316        if (local_min<min) {
1317          min=local_min;
1318          lowi=local_lowi;
1319          lowj=local_lowj;
1320        }
1321      }
1322    
1323    #ifdef PASO_MPI
1324        // determine the processor on which the minimum occurs
1325        next = temp.getDataPoint(lowi,lowj)();
1326        int lowProc = 0;
1327        double *globalMins = new double[get_MPISize()+1];
1328        int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1329        
1330        if( get_MPIRank()==0 ){
1331            next = globalMins[lowProc];
1332            for( i=1; i<get_MPISize(); i++ )
1333                if( next>globalMins[i] ){
1334                    lowProc = i;
1335                    next = globalMins[i];
1336                }
1337        }
1338        MPI_Bcast( &lowProc, 1, MPI_DOUBLE, 0, get_MPIComm() );
1339    
1340        delete [] globalMins;
1341        ProcNo = lowProc;
1342    #else
1343        ProcNo = 0;
1344    #endif
1345      DataPointNo = lowj + lowi * numDPPSample;
1346  }  }
1347    
1348  void  void
1349  Data::saveDX(std::string fileName) const  Data::saveDX(std::string fileName) const
1350  {  {
1351    getDomain().saveDX(fileName,*this);    boost::python::dict args;
1352      args["data"]=boost::python::object(this);
1353      getDomain().saveDX(fileName,args);
1354    return;    return;
1355  }  }
1356    
1357  void  void
1358  Data::saveVTK(std::string fileName) const  Data::saveVTK(std::string fileName) const
1359  {  {
1360    getDomain().saveVTK(fileName,*this);    boost::python::dict args;
1361      args["data"]=boost::python::object(this);
1362      getDomain().saveVTK(fileName,args);
1363    return;    return;
1364  }  }
1365    
1366  Data&  Data&
1367  Data::operator+=(const Data& right)  Data::operator+=(const Data& right)
1368  {  {
1369      if (isProtected()) {
1370            throw DataException("Error - attempt to update protected Data object.");
1371      }
1372    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1373    return (*this);    return (*this);
1374  }  }
# Line 990  Data::operator+=(const Data& right) Line 1376  Data::operator+=(const Data& right)
1376  Data&  Data&
1377  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1378  {  {
1379    binaryOp(right,plus<double>());    Data tmp(right,getFunctionSpace(),false);
1380      binaryOp(tmp,plus<double>());
1381    return (*this);    return (*this);
1382  }  }
1383    
1384  Data&  Data&
1385  Data::operator-=(const Data& right)  Data::operator-=(const Data& right)
1386  {  {
1387      if (isProtected()) {
1388            throw DataException("Error - attempt to update protected Data object.");
1389      }
1390    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1391    return (*this);    return (*this);
1392  }  }
# Line 1004  Data::operator-=(const Data& right) Line 1394  Data::operator-=(const Data& right)
1394  Data&  Data&
1395  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1396  {  {
1397    binaryOp(right,minus<double>());    Data tmp(right,getFunctionSpace(),false);
1398      binaryOp(tmp,minus<double>());
1399    return (*this);    return (*this);
1400  }  }
1401    
1402  Data&  Data&
1403  Data::operator*=(const Data& right)  Data::operator*=(const Data& right)
1404  {  {
1405      if (isProtected()) {
1406            throw DataException("Error - attempt to update protected Data object.");
1407      }
1408    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1409    return (*this);    return (*this);
1410  }  }
# Line 1018  Data::operator*=(const Data& right) Line 1412  Data::operator*=(const Data& right)
1412  Data&  Data&
1413  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1414  {  {
1415    binaryOp(right,multiplies<double>());    Data tmp(right,getFunctionSpace(),false);
1416      binaryOp(tmp,multiplies<double>());
1417    return (*this);    return (*this);
1418  }  }
1419    
1420  Data&  Data&
1421  Data::operator/=(const Data& right)  Data::operator/=(const Data& right)
1422  {  {
1423      if (isProtected()) {
1424            throw DataException("Error - attempt to update protected Data object.");
1425      }
1426    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1427    return (*this);    return (*this);
1428  }  }
# Line 1032  Data::operator/=(const Data& right) Line 1430  Data::operator/=(const Data& right)
1430  Data&  Data&
1431  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1432  {  {
1433    binaryOp(right,divides<double>());    Data tmp(right,getFunctionSpace(),false);
1434      binaryOp(tmp,divides<double>());
1435    return (*this);    return (*this);
1436  }  }
1437    
1438  Data  Data
1439    Data::rpowO(const boost::python::object& left) const
1440    {
1441      Data left_d(left,*this);
1442      return left_d.powD(*this);
1443    }
1444    
1445    Data
1446  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1447  {  {
1448    Data result;    Data tmp(right,getFunctionSpace(),false);
1449    result.copy(*this);    return powD(tmp);
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1450  }  }
1451    
1452  Data  Data
1453  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1454  {  {
1455    Data result;    Data result;
1456    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1457    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1458         result.binaryOp(*this,escript::rpow);
1459      } else {
1460         result.copy(*this);
1461         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1462      }
1463    return result;    return result;
1464  }  }
1465    
1466    
1467  //  //
1468  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1469  Data  Data
1470  escript::operator+(const Data& left, const Data& right)  escript::operator+(const Data& left, const Data& right)
1471  {  {
1472    Data result;    Data result;
1473    //    //
1474    // perform a deep copy    // perform a deep copy
1475    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1476    result+=right;       result.copy(right);
1477         result+=left;
1478      } else {
1479         result.copy(left);
1480         result+=right;
1481      }
1482    return result;    return result;
1483  }  }
1484    
1485  //  //
1486  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1487  Data  Data
1488  escript::operator-(const Data& left, const Data& right)  escript::operator-(const Data& left, const Data& right)
1489  {  {
1490    Data result;    Data result;
1491    //    //
1492    // perform a deep copy    // perform a deep copy
1493    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1494    result-=right;       result=right.neg();
1495         result+=left;
1496      } else {
1497         result.copy(left);
1498         result-=right;
1499      }
1500    return result;    return result;
1501  }  }
1502    
1503  //  //
1504  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1505  Data  Data
1506  escript::operator*(const Data& left, const Data& right)  escript::operator*(const Data& left, const Data& right)
1507  {  {
1508    Data result;    Data result;
1509    //    //
1510    // perform a deep copy    // perform a deep copy
1511    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1512    result*=right;       result.copy(right);
1513         result*=left;
1514      } else {
1515         result.copy(left);
1516         result*=right;
1517      }
1518    return result;    return result;
1519  }  }
1520    
1521  //  //
1522  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1523  Data  Data
1524  escript::operator/(const Data& left, const Data& right)  escript::operator/(const Data& left, const Data& right)
1525  {  {
1526    Data result;    Data result;
1527    //    //
1528    // perform a deep copy    // perform a deep copy
1529    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1530    result/=right;       result=right.oneOver();
1531         result*=left;
1532      } else {
1533         result.copy(left);
1534         result/=right;
1535      }
1536    return result;    return result;
1537  }  }
1538    
1539  //  //
1540  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1541  Data  Data
1542  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1543  {  {
1544    //    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;  
1545  }  }
1546    
1547  //  //
1548  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1549  Data  Data
1550  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1551  {  {
1552    //    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;  
1553  }  }
1554    
1555  //  //
1556  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1557  Data  Data
1558  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1559  {  {
1560    //    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;  
1561  }  }
1562    
1563  //  //
1564  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1565  Data  Data
1566  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1567  {  {
1568    //    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;  
1569  }  }
1570    
1571  //  //
1572  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1573  Data  Data
1574  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1575  {  {
1576    //    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;  
1577  }  }
1578    
1579  //  //
1580  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1581  Data  Data
1582  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1583  {  {
1584    //    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;  
1585  }  }
1586    
1587  //  //
1588  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1589  Data  Data
1590  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1591  {  {
1592    //    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;  
1593  }  }
1594    
1595  //  //
1596  // NOTE: It is essential to specify the namepsace this operator belongs to  // NOTE: It is essential to specify the namespace this operator belongs to
1597  Data  Data
1598  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1599  {  {
1600    //    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;  
1601  }  }
1602    
1603  //  //
 // NOTE: It is essential to specify the namepsace this operator belongs to  
1604  //bool escript::operator==(const Data& left, const Data& right)  //bool escript::operator==(const Data& left, const Data& right)
1605  //{  //{
1606  //  /*  //  /*
# Line 1268  escript::operator/(const boost::python:: Line 1645  escript::operator/(const boost::python::
1645  //  return ret;  //  return ret;
1646  //}  //}
1647    
1648    /* TODO */
1649    /* global reduction */
1650  Data  Data
1651  Data::getItem(const boost::python::object& key) const  Data::getItem(const boost::python::object& key) const
1652  {  {
# Line 1282  Data::getItem(const boost::python::objec Line 1661  Data::getItem(const boost::python::objec
1661    return getSlice(slice_region);    return getSlice(slice_region);
1662  }  }
1663    
1664    /* TODO */
1665    /* global reduction */
1666  Data  Data
1667  Data::getSlice(const DataArrayView::RegionType& region) const  Data::getSlice(const DataArrayView::RegionType& region) const
1668  {  {
1669    return Data(*this,region);    return Data(*this,region);
1670  }  }
1671    
1672    /* TODO */
1673    /* global reduction */
1674  void  void
1675  Data::setItemO(const boost::python::object& key,  Data::setItemO(const boost::python::object& key,
1676                 const boost::python::object& value)                 const boost::python::object& value)
# Line 1296  Data::setItemO(const boost::python::obje Line 1679  Data::setItemO(const boost::python::obje
1679    setItemD(key,tempData);    setItemD(key,tempData);
1680  }  }
1681    
1682    /* TODO */
1683    /* global reduction */
1684  void  void
1685  Data::setItemD(const boost::python::object& key,  Data::setItemD(const boost::python::object& key,
1686                 const Data& value)                 const Data& value)
1687  {  {
1688    const DataArrayView& view=getPointDataView();    const DataArrayView& view=getPointDataView();
1689    
1690    DataArrayView::RegionType slice_region=view.getSliceRegion(key);    DataArrayView::RegionType slice_region=view.getSliceRegion(key);
1691    if (slice_region.size()!=view.getRank()) {    if (slice_region.size()!=view.getRank()) {
1692      throw DataException("Error - slice size does not match Data rank.");      throw DataException("Error - slice size does not match Data rank.");
# Line 1312  Data::setItemD(const boost::python::obje Line 1698  Data::setItemD(const boost::python::obje
1698    }    }
1699  }  }
1700    
1701    /* TODO */
1702    /* global reduction */
1703  void  void
1704  Data::setSlice(const Data& value,  Data::setSlice(const Data& value,
1705                 const DataArrayView::RegionType& region)                 const DataArrayView::RegionType& region)
1706  {  {
1707      if (isProtected()) {
1708            throw DataException("Error - attempt to update protected Data object.");
1709      }
1710    Data tempValue(value);    Data tempValue(value);
1711    typeMatchLeft(tempValue);    typeMatchLeft(tempValue);
1712    typeMatchRight(tempValue);    typeMatchRight(tempValue);
# Line 1351  Data::typeMatchRight(const Data& right) Line 1742  Data::typeMatchRight(const Data& right)
1742  }  }
1743    
1744  void  void
1745    Data::setTaggedValueByName(std::string name,
1746                               const boost::python::object& value)
1747    {
1748         if (getFunctionSpace().getDomain().isValidTagName(name)) {
1749            int tagKey=getFunctionSpace().getDomain().getTag(name);
1750            setTaggedValue(tagKey,value);
1751         }
1752    }
1753    void
1754  Data::setTaggedValue(int tagKey,  Data::setTaggedValue(int tagKey,
1755                       const boost::python::object& value)                       const boost::python::object& value)
1756  {  {
1757      if (isProtected()) {
1758            throw DataException("Error - attempt to update protected Data object.");
1759      }
1760    //    //
1761    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1762    tag();    tag();
# Line 1375  void Line 1778  void
1778  Data::setTaggedValueFromCPP(int tagKey,  Data::setTaggedValueFromCPP(int tagKey,
1779                              const DataArrayView& value)                              const DataArrayView& value)
1780  {  {
1781      if (isProtected()) {
1782            throw DataException("Error - attempt to update protected Data object.");
1783      }
1784    //    //
1785    // Ensure underlying data object is of type DataTagged    // Ensure underlying data object is of type DataTagged
1786    tag();    tag();
# Line 1388  Data::setTaggedValueFromCPP(int tagKey, Line 1794  Data::setTaggedValueFromCPP(int tagKey,
1794    m_data->setTaggedValue(tagKey,value);    m_data->setTaggedValue(tagKey,value);
1795  }  }
1796    
1797  void  int
1798  Data::setRefValue(int ref,  Data::getTagNumber(int dpno)
                   const boost::python::numeric::array& value)  
1799  {  {
1800    //    return m_data->getTagNumber(dpno);
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
   
   //  
   // Load DataArray with values from data-points specified by ref  
   m_data->getRefValue(ref,valueDataArray);  
   
   //  
   // Load values from valueDataArray into return numarray  
   
   // extract the shape of the numarray  
   int rank = value.getrank();  
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
   
   if (rank==0) {  
     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");  
   }  
   if (rank==1) {  
     for (int i=0; i < shape[0]; i++) {  
       value[i] = valueView(i);  
     }  
   }  
   if (rank==2) {  
     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");  
   }  
   if (rank==3) {  
     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");  
   }  
   if (rank==4) {  
     throw DataException("Data::getRefValue error: only rank 1 data handled for now.");  
   }  
   
1801  }  }
1802    
1803  void  void
# Line 1471  Data::archiveData(const std::string file Line 1825  Data::archiveData(const std::string file
1825      dataType = 3;      dataType = 3;
1826      cout << "\tdataType: DataExpanded" << endl;      cout << "\tdataType: DataExpanded" << endl;
1827    }    }
1828    
1829    if (dataType == -1) {    if (dataType == -1) {
1830      throw DataException("archiveData Error: undefined dataType");      throw DataException("archiveData Error: undefined dataType");
1831    }    }
# Line 1484  Data::archiveData(const std::string file Line 1839  Data::archiveData(const std::string file
1839    int dataPointSize = getDataPointSize();    int dataPointSize = getDataPointSize();
1840    int dataLength = getLength();    int dataLength = getLength();
1841    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
1842    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
1843    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1844      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
1845    }    }
1846    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
1847    if (isTagged()) {    if (isTagged()) {
1848      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
1849        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);        tagNumbers[sampleNo] = getFunctionSpace().getTagFromSampleNo(sampleNo);
# Line 1510  Data::archiveData(const std::string file Line 1865  Data::archiveData(const std::string file
1865    cout << ">" << endl;    cout << ">" << endl;
1866    
1867    //    //
1868    // Write common data items to archive file    // Open archive file
1869    ofstream archiveFile;    ofstream archiveFile;
1870    archiveFile.open(fileName.data(), ios::out);    archiveFile.open(fileName.data(), ios::out);
1871    
# Line 1518  Data::archiveData(const std::string file Line 1873  Data::archiveData(const std::string file
1873      throw DataException("archiveData Error: problem opening archive file");      throw DataException("archiveData Error: problem opening archive file");
1874    }    }
1875    
1876      //
1877      // Write common data items to archive file
1878    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&dataType),sizeof(int));
1879    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noSamples),sizeof(int));
1880    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.write(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1541  Data::archiveData(const std::string file Line 1898  Data::archiveData(const std::string file
1898      throw DataException("archiveData Error: problem writing to archive file");      throw DataException("archiveData Error: problem writing to archive file");
1899    }    }
1900    
   archiveFile.close();  
   
   if (!archiveFile.good()) {  
     throw DataException("archiveData Error: problem closing archive file");  
   }  
   
1901    //    //
1902    // Collect and archive underlying data values for each Data type    // Archive underlying data values for each Data type
1903      int noValues;
1904    switch (dataType) {    switch (dataType) {
1905      case 0:      case 0:
1906        // DataEmpty        // DataEmpty
1907          noValues = 0;
1908          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1909          cout << "\tnoValues: " << noValues << endl;
1910        break;        break;
1911      case 1:      case 1:
1912        // DataConstant        // DataConstant
1913          noValues = m_data->getLength();
1914          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1915          cout << "\tnoValues: " << noValues << endl;
1916          if (m_data->archiveData(archiveFile,noValues)) {
1917            throw DataException("archiveData Error: problem writing data to archive file");
1918          }
1919        break;        break;
1920      case 2:      case 2:
1921        // DataTagged        // DataTagged
1922          noValues = m_data->getLength();
1923          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1924          cout << "\tnoValues: " << noValues << endl;
1925          if (m_data->archiveData(archiveFile,noValues)) {
1926            throw DataException("archiveData Error: problem writing data to archive file");
1927          }
1928        break;        break;
1929      case 3:      case 3:
1930        // DataExpanded        // DataExpanded
1931          noValues = m_data->getLength();
1932          archiveFile.write(reinterpret_cast<char *>(&noValues),sizeof(int));
1933          cout << "\tnoValues: " << noValues << endl;
1934          if (m_data->archiveData(archiveFile,noValues)) {
1935            throw DataException("archiveData Error: problem writing data to archive file");
1936          }
1937        break;        break;
1938    }    }
1939    
1940      if (!archiveFile.good()) {
1941        throw DataException("archiveData Error: problem writing data to archive file");
1942      }
1943    
1944      //
1945      // Close archive file
1946      archiveFile.close();
1947    
1948      if (!archiveFile.good()) {
1949        throw DataException("archiveData Error: problem closing archive file");
1950      }
1951    
1952  }  }
1953    
1954  void  void
# Line 1589  Data::extractData(const std::string file Line 1974  Data::extractData(const std::string file
1974    int flatShape[4];    int flatShape[4];
1975    
1976    //    //
1977    // Open the archive file and read common data items    // Open the archive file
1978    ifstream archiveFile;    ifstream archiveFile;
1979    archiveFile.open(fileName.data(), ios::in);    archiveFile.open(fileName.data(), ios::in);
1980    
# Line 1597  Data::extractData(const std::string file Line 1982  Data::extractData(const std::string file
1982      throw DataException("extractData Error: problem opening archive file");      throw DataException("extractData Error: problem opening archive file");
1983    }    }
1984    
1985      //
1986      // Read common data items from archive file
1987    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&dataType),sizeof(int));
1988    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noSamples),sizeof(int));
1989    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));    archiveFile.read(reinterpret_cast<char *>(&noDPPSample),sizeof(int));
# Line 1610  Data::extractData(const std::string file Line 1997  Data::extractData(const std::string file
1997        dataPointShape.push_back(flatShape[dim]);        dataPointShape.push_back(flatShape[dim]);
1998      }      }
1999    }    }
2000    int referenceNumbers[noSamples];    vector<int> referenceNumbers(noSamples);
2001    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2002      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));      archiveFile.read(reinterpret_cast<char *>(&referenceNumbers[sampleNo]),sizeof(int));
2003    }    }
2004    int tagNumbers[noSamples];    vector<int> tagNumbers(noSamples);
2005    if (dataType==2) {    if (dataType==2) {
2006      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {      for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2007        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));        archiveFile.read(reinterpret_cast<char *>(&tagNumbers[sampleNo]),sizeof(int));
# Line 1625  Data::extractData(const std::string file Line 2012  Data::extractData(const std::string file
2012      throw DataException("extractData Error: problem reading from archive file");      throw DataException("extractData Error: problem reading from archive file");
2013    }    }
2014    
2015    archiveFile.close();    //
2016      // Verify the values just read from the archive file
   if (!archiveFile.good()) {  
     throw DataException("extractData Error: problem closing archive file");  
   }  
   
2017    switch (dataType) {    switch (dataType) {
2018      case 0:      case 0:
2019        cout << "\tdataType: DataEmpty" << endl;        cout << "\tdataType: DataEmpty" << endl;
# Line 1667  Data::extractData(const std::string file Line 2050  Data::extractData(const std::string file
2050      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2051    }    }
2052    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2053      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2054        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2055      }      }
2056    }    }
# Line 1685  Data::extractData(const std::string file Line 2068  Data::extractData(const std::string file
2068    
2069    //    //
2070    // Load this DataVector with the appropriate values    // Load this DataVector with the appropriate values
2071      int noValues;
2072      archiveFile.read(reinterpret_cast<char *>(&noValues),sizeof(int));
2073      cout << "\tnoValues: " << noValues << endl;
2074    switch (dataType) {    switch (dataType) {
2075      case 0:      case 0:
2076        // DataEmpty        // DataEmpty
2077          if (noValues != 0) {
2078            throw DataException("extractData Error: problem reading data from archive file");
2079          }
2080        break;        break;
2081      case 1:      case 1:
2082        // DataConstant        // DataConstant
2083          if (dataVec.extractData(archiveFile,noValues)) {
2084            throw DataException("extractData Error: problem reading data from archive file");
2085          }
2086        break;        break;
2087      case 2:      case 2:
2088        // DataTagged        // DataTagged
2089          if (dataVec.extractData(archiveFile,noValues)) {
2090            throw DataException("extractData Error: problem reading data from archive file");
2091          }
2092        break;        break;
2093      case 3:      case 3:
2094        // DataExpanded        // DataExpanded
2095          if (dataVec.extractData(archiveFile,noValues)) {
2096            throw DataException("extractData Error: problem reading data from archive file");
2097          }
2098        break;        break;
2099    }    }
2100    
2101      if (!archiveFile.good()) {
2102        throw DataException("extractData Error: problem reading from archive file");
2103      }
2104    
2105      //
2106      // Close archive file
2107      archiveFile.close();
2108    
2109      if (!archiveFile.good()) {
2110        throw DataException("extractData Error: problem closing archive file");
2111      }
2112    
2113    //    //
2114    // Construct an appropriate Data object    // Construct an appropriate Data object
2115    DataAbstract* tempData;    DataAbstract* tempData;
# Line 1723  Data::extractData(const std::string file Line 2133  Data::extractData(const std::string file
2133    }    }
2134    shared_ptr<DataAbstract> temp_data(tempData);    shared_ptr<DataAbstract> temp_data(tempData);
2135    m_data=temp_data;    m_data=temp_data;
   
2136  }  }
2137    
2138  ostream& escript::operator<<(ostream& o, const Data& data)  ostream& escript::operator<<(ostream& o, const Data& data)
# Line 1731  ostream& escript::operator<<(ostream& o, Line 2140  ostream& escript::operator<<(ostream& o,
2140    o << data.toString();    o << data.toString();
2141    return o;    return o;
2142  }  }
2143    
2144    Data
2145    escript::C_GeneralTensorProduct(Data& arg_0,
2146                         Data& arg_1,
2147                         int axis_offset,
2148                         int transpose)
2149    {
2150      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2151      // SM is the product of the last axis_offset entries in arg_0.getShape().
2152    
2153    
2154      // Interpolate if necessary and find an appropriate function space
2155      Data arg_0_Z, arg_1_Z;
2156      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2157        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2158          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2159          arg_1_Z = Data(arg_1);
2160        }
2161        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2162          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2163          arg_0_Z =Data(arg_0);
2164        }
2165        else {
2166          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2167        }
2168      } else {
2169          arg_0_Z = Data(arg_0);
2170          arg_1_Z = Data(arg_1);
2171      }
2172      // Get rank and shape of inputs
2173      int rank0 = arg_0_Z.getDataPointRank();
2174      int rank1 = arg_1_Z.getDataPointRank();
2175      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2176      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2177    
2178      // Prepare for the loops of the product and verify compatibility of shapes
2179      int start0=0, start1=0;
2180      if (transpose == 0)       {}
2181      else if (transpose == 1)  { start0 = axis_offset; }
2182      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2183      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2184    
2185      // Adjust the shapes for transpose
2186      DataArrayView::ShapeType tmpShape0;
2187      DataArrayView::ShapeType tmpShape1;
2188      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2189      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2190    
2191    #if 0
2192      // For debugging: show shape after transpose
2193      char tmp[100];
2194      std::string shapeStr;
2195      shapeStr = "(";
2196      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2197      shapeStr += ")";
2198      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2199      shapeStr = "(";
2200      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2201      shapeStr += ")";
2202      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2203    #endif
2204    
2205      // Prepare for the loops of the product
2206      int SL=1, SM=1, SR=1;
2207      for (int i=0; i<rank0-axis_offset; i++)   {
2208        SL *= tmpShape0[i];
2209      }
2210      for (int i=rank0-axis_offset; i<rank0; i++)   {
2211        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2212          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2213        }
2214        SM *= tmpShape0[i];
2215      }
2216      for (int i=axis_offset; i<rank1; i++)     {
2217        SR *= tmpShape1[i];
2218      }
2219    
2220      // Define the shape of the output
2221      DataArrayView::ShapeType shape2;
2222      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2223      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2224    
2225      // Declare output Data object
2226      Data res;
2227    
2228      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2229        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2230        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2231        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2232        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2233        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2234      }
2235      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2236    
2237        // Prepare the DataConstant input
2238        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2239        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2240    
2241        // Borrow DataTagged input from Data object
2242        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2243        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2244    
2245        // Prepare a DataTagged output 2
2246        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2247        res.tag();
2248        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2249        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2250    
2251        // Prepare offset into DataConstant
2252        int offset_0 = tmp_0->getPointOffset(0,0);
2253        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2254        // Get the views
2255        DataArrayView view_1 = tmp_1->getDefaultValue();
2256        DataArrayView view_2 = tmp_2->getDefaultValue();
2257        // Get the pointers to the actual data
2258        double *ptr_1 = &((view_1.getData())[0]);
2259        double *ptr_2 = &((view_2.getData())[0]);
2260        // Compute an MVP for the default
2261        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2262        // Compute an MVP for each tag
2263        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2264        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2265        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2266          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2267          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2268          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2269          double *ptr_1 = &view_1.getData(0);
2270          double *ptr_2 = &view_2.getData(0);
2271          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2272        }
2273    
2274      }
2275      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2276    
2277        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2278        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2279        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2280        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2281        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2282        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2283        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2284        int sampleNo_1,dataPointNo_1;
2285        int numSamples_1 = arg_1_Z.getNumSamples();
2286        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2287        int offset_0 = tmp_0->getPointOffset(0,0);
2288        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2289        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2290          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2291            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2292            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2293            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2294            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2295            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2296            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2297          }
2298        }
2299    
2300      }
2301      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2302    
2303        // Borrow DataTagged input from Data object
2304        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2305        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2306    
2307        // Prepare the DataConstant input
2308        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2309        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2310    
2311        // Prepare a DataTagged output 2
2312        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2313        res.tag();
2314        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2315        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2316    
2317        // Prepare offset into DataConstant
2318        int offset_1 = tmp_1->getPointOffset(0,0);
2319        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2320        // Get the views
2321        DataArrayView view_0 = tmp_0->getDefaultValue();
2322        DataArrayView view_2 = tmp_2->getDefaultValue();
2323        // Get the pointers to the actual data
2324        double *ptr_0 = &((view_0.getData())[0]);
2325        double *ptr_2 = &((view_2.getData())[0]);
2326        // Compute an MVP for the default
2327        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2328        // Compute an MVP for each tag
2329        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2330        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2331        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2332          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2333          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2334          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2335          double *ptr_0 = &view_0.getData(0);
2336          double *ptr_2 = &view_2.getData(0);
2337          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2338        }
2339    
2340      }
2341      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2342    
2343        // Borrow DataTagged input from Data object
2344        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2345        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2346    
2347        // Borrow DataTagged input from Data object
2348        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2349        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2350    
2351        // Prepare a DataTagged output 2
2352        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2353        res.tag();  // DataTagged output
2354        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2355        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2356    
2357        // Get the views
2358        DataArrayView view_0 = tmp_0->getDefaultValue();
2359        DataArrayView view_1 = tmp_1->getDefaultValue();
2360        DataArrayView view_2 = tmp_2->getDefaultValue();
2361        // Get the pointers to the actual data
2362        double *ptr_0 = &((view_0.getData())[0]);
2363        double *ptr_1 = &((view_1.getData())[0]);
2364        double *ptr_2 = &((view_2.getData())[0]);
2365        // Compute an MVP for the default
2366        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2367        // Merge the tags
2368        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2369        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2370        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2371        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2372          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2373        }
2374        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2375          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2376        }
2377        // Compute an MVP for each tag
2378        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2379        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2380          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2381          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2382          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2383          double *ptr_0 = &view_0.getData(0);
2384          double *ptr_1 = &view_1.getData(0);
2385          double *ptr_2 = &view_2.getData(0);
2386          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2387        }
2388    
2389      }
2390      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2391    
2392        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2393        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2394        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2395        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2396        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2397        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2398        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2399        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2400        int sampleNo_0,dataPointNo_0;
2401        int numSamples_0 = arg_0_Z.getNumSamples();
2402        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2403        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2404        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2405          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2406          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2407          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2408            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2409            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2410            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2411            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2412            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2413          }
2414        }
2415    
2416      }
2417      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2418    
2419        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2420        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2421        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2422        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2423        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2424        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2425        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2426        int sampleNo_0,dataPointNo_0;
2427        int numSamples_0 = arg_0_Z.getNumSamples();
2428        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2429        int offset_1 = tmp_1->getPointOffset(0,0);
2430        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2431        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2432          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2433            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2434            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2435            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2436            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2437            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2438            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2439          }
2440        }
2441    
2442    
2443      }
2444      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2445    
2446        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2447        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2448        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2449        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2450        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2451        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2452        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2453        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2454        int sampleNo_0,dataPointNo_0;
2455        int numSamples_0 = arg_0_Z.getNumSamples();
2456        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2457        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2458        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2459          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2460          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2461          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2462            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2463            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2464            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2465            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2466            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2467          }
2468        }
2469    
2470      }
2471      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2472    
2473        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2474        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2475        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2476        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2477        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2478        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2479        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2480        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2481        int sampleNo_0,dataPointNo_0;
2482        int numSamples_0 = arg_0_Z.getNumSamples();
2483        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2484        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2485        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2486          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2487            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2488            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2489            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2490            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2491            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2492            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2493            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2494          }
2495        }
2496    
2497      }
2498      else {
2499        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2500      }
2501    
2502      return res;
2503    }
2504    
2505    DataAbstract*
2506    Data::borrowData() const
2507    {
2508      return m_data.get();
2509    }
2510    
2511    /* Member functions specific to the MPI implementation */
2512    
2513    void
2514    Data::print()
2515    {
2516      int i,j;
2517      
2518      printf( "Data is %dX%d\n", getNumSamples(), getNumDataPointsPerSample() );
2519      for( i=0; i<getNumSamples(); i++ )
2520      {
2521        printf( "[%6d]", i );
2522        for( j=0; j<getNumDataPointsPerSample(); j++ )
2523          printf( "\t%10.7g", (getSampleData(i))[j] );
2524        printf( "\n" );
2525      }
2526    }
2527    
2528    int
2529    Data::get_MPISize() const
2530    {
2531        int error, size;
2532    #ifdef PASO_MPI
2533        error = MPI_Comm_size( get_MPIComm(), &size );
2534    #else
2535        size = 1;
2536    #endif
2537        return size;
2538    }
2539    
2540    int
2541    Data::get_MPIRank() const
2542    {
2543        int error, rank;
2544    #ifdef PASO_MPI
2545        error = MPI_Comm_rank( get_MPIComm(), &rank );
2546    #else
2547        rank = 0;
2548    #endif
2549        return rank;
2550    }
2551    
2552    MPI_Comm
2553    Data::get_MPIComm() const
2554    {
2555    #ifdef PASO_MPI
2556        return MPI_COMM_WORLD;
2557    #else
2558        return -1;
2559    #endif
2560    }
2561    

Legend:
Removed from v.122  
changed lines
  Added in v.1092

  ViewVC Help
Powered by ViewVC 1.1.26