/[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

revision 2037 by jfenwick, Thu Nov 13 06:17:12 2008 UTC revision 2105 by jfenwick, Fri Nov 28 01:52:12 2008 UTC
# Line 26  Line 26 
26  #include "EscriptParams.h"  #include "EscriptParams.h"
27    
28  extern "C" {  extern "C" {
29  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
30  }  }
31    
32  #include <fstream>  #include <fstream>
33  #include <algorithm>  #include <algorithm>
34  #include <vector>  #include <vector>
35  #include <functional>  #include <functional>
36    #include <sstream>  // so we can throw messages about ranks
37    
38  #include <boost/python/dict.hpp>  #include <boost/python/dict.hpp>
39  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
# Line 130  Data::Data(const Data& inData, Line 131  Data::Data(const Data& inData,
131      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space      if (inData.isConstant()) {  // for a constant function, we just need to use the new function space
132        if (!inData.probeInterpolation(functionspace))        if (!inData.probeInterpolation(functionspace))
133        {           // Even though this is constant, we still need to check whether interpolation is allowed        {           // Even though this is constant, we still need to check whether interpolation is allowed
134      throw FunctionSpaceException("Call to probeInterpolation returned false for DataConstant.");      throw FunctionSpaceException("Cannot interpolate across to the domain of the specified FunctionSpace. (DataConstant)");
135        }        }
136        // if the data is not lazy, this will just be a cast to DataReady        // if the data is not lazy, this will just be a cast to DataReady
137        DataReady_ptr dr=inData.m_data->resolve();        DataReady_ptr dr=inData.m_data->resolve();
# Line 176  Data::Data(const numeric::array& value, Line 177  Data::Data(const numeric::array& value,
177    initialise(value,what,expanded);    initialise(value,what,expanded);
178    m_protected=false;    m_protected=false;
179  }  }
 /*  
 Data::Data(const DataArrayView& value,  
        const FunctionSpace& what,  
            bool expanded)  
 {  
   initialise(value,what,expanded);  
   m_protected=false;  
 }*/  
180    
181  Data::Data(const DataTypes::ValueType& value,  Data::Data(const DataTypes::ValueType& value,
182           const DataTypes::ShapeType& shape,           const DataTypes::ShapeType& shape,
# Line 212  Data::Data(const object& value, Line 205  Data::Data(const object& value,
205    
206    // extract the shape of the numarray    // extract the shape of the numarray
207    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);    DataTypes::ShapeType tempShape=DataTypes::shapeFromNumArray(asNumArray);
208  // /*  for (int i=0; i < asNumArray.getrank(); i++) {    if (DataTypes::getRank(tempShape)==0) {
 //     tempShape.push_back(extract<int>(asNumArray.getshape()[i]));  
 //   }*/  
 //   // get the space for the data vector  
 //   int len = DataTypes::noValues(tempShape);  
 //   DataVector temp_data(len, 0.0, len);  
 // /*  DataArrayView temp_dataView(temp_data, tempShape);  
 //   temp_dataView.copy(asNumArray);*/  
 //   temp_data.copyFromNumArray(asNumArray);  
   
   //  
   // Create DataConstant using the given value and all other parameters  
   // copied from other. If value is a rank 0 object this Data  
   // will assume the point data shape of other.  
   
   if (DataTypes::getRank(tempShape)/*temp_dataView.getRank()*/==0) {  
209    
210    
211      // get the space for the data vector      // get the space for the data vector
212      int len1 = DataTypes::noValues(tempShape);      int len1 = DataTypes::noValues(tempShape);
213      DataVector temp_data(len1, 0.0, len1);      DataVector temp_data(len1, 0.0, len1);
214      temp_data.copyFromNumArray(asNumArray);      temp_data.copyFromNumArray(asNumArray,1);
215    
216      int len = DataTypes::noValues(other.getDataPointShape());      int len = DataTypes::noValues(other.getDataPointShape());
217    
218      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);      DataVector temp2_data(len, temp_data[0]/*temp_dataView()*/, len);
     //DataArrayView temp2_dataView(temp2_data, other.getPointDataView().getShape());  
 //     initialise(temp2_dataView, other.getFunctionSpace(), false);  
   
219      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);      DataConstant* t=new DataConstant(other.getFunctionSpace(),other.getDataPointShape(),temp2_data);
 //     boost::shared_ptr<DataAbstract> sp(t);  
 //     m_data=sp;  
220      m_data=DataAbstract_ptr(t);      m_data=DataAbstract_ptr(t);
221    
222    } else {    } else {
223      //      //
224      // Create a DataConstant with the same sample shape as other      // Create a DataConstant with the same sample shape as other
 //     initialise(temp_dataView, other.getFunctionSpace(), false);  
225      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());      DataConstant* t=new DataConstant(asNumArray,other.getFunctionSpace());
226  //     boost::shared_ptr<DataAbstract> sp(t);  //     boost::shared_ptr<DataAbstract> sp(t);
227  //     m_data=sp;  //     m_data=sp;
# Line 435  Data::getShapeTuple() const Line 407  Data::getShapeTuple() const
407  // It can't work out what type the function is based soley on its name.  // It can't work out what type the function is based soley on its name.
408  // There are ways to fix this involving creating function pointer variables for each form  // There are ways to fix this involving creating function pointer variables for each form
409  // but there doesn't seem to be a need given that the methods have the same name from the python point of view  // but there doesn't seem to be a need given that the methods have the same name from the python point of view
410  Data*  Data
411  Data::copySelf()  Data::copySelf()
412  {  {
413     DataAbstract* temp=m_data->deepCopy();     DataAbstract* temp=m_data->deepCopy();
414     return new Data(temp);     return Data(temp);
415  }  }
416    
417  void  void
# Line 1151  Data::integrateWorker() const Line 1123  Data::integrateWorker() const
1123    // calculate the integral values    // calculate the integral values
1124    vector<double> integrals(dataPointSize);    vector<double> integrals(dataPointSize);
1125    vector<double> integrals_local(dataPointSize);    vector<double> integrals_local(dataPointSize);
1126      const AbstractContinuousDomain* dom=dynamic_cast<const AbstractContinuousDomain*>(getDomain().get());
1127      if (dom==0)
1128      {            
1129        throw DataException("Can not integrate over non-continuous domains.");
1130      }
1131  #ifdef PASO_MPI  #ifdef PASO_MPI
1132    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals_local,*this);    dom->setToIntegrals(integrals_local,*this);
1133    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory    // Global sum: use an array instead of a vector because elements of array are guaranteed to be contiguous in memory
1134    double *tmp = new double[dataPointSize];    double *tmp = new double[dataPointSize];
1135    double *tmp_local = new double[dataPointSize];    double *tmp_local = new double[dataPointSize];
# Line 1162  Data::integrateWorker() const Line 1139  Data::integrateWorker() const
1139    delete[] tmp;    delete[] tmp;
1140    delete[] tmp_local;    delete[] tmp_local;
1141  #else  #else
1142    AbstractContinuousDomain::asAbstractContinuousDomain(*getDomain()).setToIntegrals(integrals,*this);    dom->setToIntegrals(integrals,*this);
1143  #endif  #endif
1144    
1145    //    //
# Line 1324  Data::tanh() const Line 1301  Data::tanh() const
1301  Data  Data
1302  Data::erf() const  Data::erf() const
1303  {  {
1304  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1305    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1306  #else  #else
1307    if (isLazy())    if (isLazy())
# Line 1344  Data::asinh() const Line 1321  Data::asinh() const
1321      DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);      DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1322      return Data(c);      return Data(c);
1323    }    }
1324  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1325    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1326  #else  #else
1327    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1359  Data::acosh() const Line 1336  Data::acosh() const
1336      DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);      DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1337      return Data(c);      return Data(c);
1338    }    }
1339  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1340    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1341  #else  #else
1342    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1374  Data::atanh() const Line 1351  Data::atanh() const
1351      DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);      DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1352      return Data(c);      return Data(c);
1353    }    }
1354  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1355    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1356  #else  #else
1357    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1483  Data::Lsup() Line 1460  Data::Lsup()
1460  {  {
1461     if (isLazy())     if (isLazy())
1462     {     {
1463      expand();      resolve();
1464     }     }
1465     return LsupWorker();     return LsupWorker();
1466  }  }
# Line 1503  Data::sup() Line 1480  Data::sup()
1480  {  {
1481     if (isLazy())     if (isLazy())
1482     {     {
1483      expand();      resolve();
1484     }     }
1485     return supWorker();     return supWorker();
1486  }  }
# Line 1523  Data::inf() Line 1500  Data::inf()
1500  {  {
1501     if (isLazy())     if (isLazy())
1502     {     {
1503      expand();      resolve();
1504     }     }
1505     return infWorker();     return infWorker();
1506  }  }
# Line 1615  Data::minval() const Line 1592  Data::minval() const
1592  Data  Data
1593  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1594  {  {
1595         if (isLazy())
1596         {
1597        Data temp(*this);
1598        temp.resolve();
1599        return temp.swapaxes(axis0,axis1);
1600         }
1601       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1602       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1603       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1722  Data::nonsymmetric() const Line 1705  Data::nonsymmetric() const
1705       }       }
1706  }  }
1707    
   
 // Doing a lazy version of this would require some thought.  
 // First it needs a parameter (which DataLazy doesn't support at the moment).  
 // (secondly although it does not apply to trace) we can't handle operations which return  
 // multiple results (like eigenvectors_values) or return values of different shapes to their input  
 // (like eigenvalues).  
1708  Data  Data
1709  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1710  {  {    
1711       if (isLazy())       if (isLazy())
1712       {       {
1713      Data temp(*this);   // to get around the fact that you can't resolve a const Data      DataLazy* c=new DataLazy(borrowDataPtr(),TRACE,axis_offset);
1714      temp.resolve();      return Data(c);
     return temp.trace(axis_offset);  
1715       }       }
1716       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1717       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1789  Data::transpose(int axis_offset) const Line 1765  Data::transpose(int axis_offset) const
1765  {      {    
1766       if (isLazy())       if (isLazy())
1767       {       {
1768      Data temp(*this);   // to get around the fact that you can't resolve a const Data      DataLazy* c=new DataLazy(borrowDataPtr(),TRANS,axis_offset);
1769      temp.resolve();      return Data(c);
     return temp.transpose(axis_offset);  
1770       }       }
1771       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1772       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1894  Data::calc_minGlobalDataPoint(int& ProcN Line 1869  Data::calc_minGlobalDataPoint(int& ProcN
1869    double next,local_min;    double next,local_min;
1870    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1871    
1872    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1873    {    {
1874      local_min=min;      local_min=min;
1875      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
# Line 1921  Data::calc_minGlobalDataPoint(int& ProcN Line 1896  Data::calc_minGlobalDataPoint(int& ProcN
1896      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPoint(lowi,lowj);
1897      int lowProc = 0;      int lowProc = 0;
1898      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1899      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1900        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1901    
1902      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1903          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 2461  Data::setTaggedValue(int tagKey, Line 2437  Data::setTaggedValue(int tagKey,
2437    }    }
2438    
2439    DataVector temp_data2;    DataVector temp_data2;
2440    temp_data2.copyFromNumArray(asNumArray);    temp_data2.copyFromNumArray(asNumArray,1);
2441    
2442    m_data->setTaggedValue(tagKey,tempShape, temp_data2);    m_data->setTaggedValue(tagKey,tempShape, temp_data2);
2443  }  }
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2488  escript::C_GeneralTensorProduct(Data& ar
2488    // SM is the product of the last axis_offset entries in arg_0.getShape().    // SM is the product of the last axis_offset entries in arg_0.getShape().
2489    
2490    // deal with any lazy data    // deal with any lazy data
2491    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2492    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2493      if (arg_0.isLazy() || arg_1.isLazy())
2494      {
2495        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2496        return Data(c);
2497      }
2498    
2499    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2500    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2571  escript::C_GeneralTensorProduct(Data& ar
2571       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z       for (int i=axis_offset; i<rank1; i++, ++out_index)   { shape2[out_index]=tmpShape1[i]; } // Last part of arg_1_Z
2572    }    }
2573    
2574      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2575      {
2576         ostringstream os;
2577         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2578         throw DataException(os.str());
2579      }
2580    
2581    // Declare output Data object    // Declare output Data object
2582    Data res;    Data res;
2583    
# Line 2923  Data::borrowReadyPtr() const Line 2911  Data::borrowReadyPtr() const
2911  std::string  std::string
2912  Data::toString() const  Data::toString() const
2913  {  {
2914      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2915      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2916        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2917      {      {
2918      stringstream temp;      stringstream temp;
2919      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();      temp << "Summary: inf="<< inf_const() << " sup=" << sup_const() << " data points=" << getNumDataPoints();

Legend:
Removed from v.2037  
changed lines
  Added in v.2105

  ViewVC Help
Powered by ViewVC 1.1.26