/[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 2086 by jfenwick, Mon Nov 24 02:38:50 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 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
# Line 238  Data::Data(const object& value, Line 216  Data::Data(const object& value,
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 1324  Data::tanh() const Line 1296  Data::tanh() const
1296  Data  Data
1297  Data::erf() const  Data::erf() const
1298  {  {
1299  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1300    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");    throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1301  #else  #else
1302    if (isLazy())    if (isLazy())
# Line 1344  Data::asinh() const Line 1316  Data::asinh() const
1316      DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);      DataLazy* c=new DataLazy(borrowDataPtr(),ASINH);
1317      return Data(c);      return Data(c);
1318    }    }
1319  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1320    return C_TensorUnaryOperation(*this, escript::asinh_substitute);    return C_TensorUnaryOperation(*this, escript::asinh_substitute);
1321  #else  #else
1322    return C_TensorUnaryOperation(*this, ::asinh);    return C_TensorUnaryOperation(*this, ::asinh);
# Line 1359  Data::acosh() const Line 1331  Data::acosh() const
1331      DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);      DataLazy* c=new DataLazy(borrowDataPtr(),ACOSH);
1332      return Data(c);      return Data(c);
1333    }    }
1334  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1335    return C_TensorUnaryOperation(*this, escript::acosh_substitute);    return C_TensorUnaryOperation(*this, escript::acosh_substitute);
1336  #else  #else
1337    return C_TensorUnaryOperation(*this, ::acosh);    return C_TensorUnaryOperation(*this, ::acosh);
# Line 1374  Data::atanh() const Line 1346  Data::atanh() const
1346      DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);      DataLazy* c=new DataLazy(borrowDataPtr(),ATANH);
1347      return Data(c);      return Data(c);
1348    }    }
1349  #ifdef _WIN32  #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1350    return C_TensorUnaryOperation(*this, escript::atanh_substitute);    return C_TensorUnaryOperation(*this, escript::atanh_substitute);
1351  #else  #else
1352    return C_TensorUnaryOperation(*this, ::atanh);    return C_TensorUnaryOperation(*this, ::atanh);
# Line 1483  Data::Lsup() Line 1455  Data::Lsup()
1455  {  {
1456     if (isLazy())     if (isLazy())
1457     {     {
1458      expand();      resolve();
1459     }     }
1460     return LsupWorker();     return LsupWorker();
1461  }  }
# Line 1503  Data::sup() Line 1475  Data::sup()
1475  {  {
1476     if (isLazy())     if (isLazy())
1477     {     {
1478      expand();      resolve();
1479     }     }
1480     return supWorker();     return supWorker();
1481  }  }
# Line 1523  Data::inf() Line 1495  Data::inf()
1495  {  {
1496     if (isLazy())     if (isLazy())
1497     {     {
1498      expand();      resolve();
1499     }     }
1500     return infWorker();     return infWorker();
1501  }  }
# Line 1615  Data::minval() const Line 1587  Data::minval() const
1587  Data  Data
1588  Data::swapaxes(const int axis0, const int axis1) const  Data::swapaxes(const int axis0, const int axis1) const
1589  {  {
1590         if (isLazy())
1591         {
1592        Data temp(*this);
1593        temp.resolve();
1594        return temp.swapaxes(axis0,axis1);
1595         }
1596       int axis0_tmp,axis1_tmp;       int axis0_tmp,axis1_tmp;
1597       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1598       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1722  Data::nonsymmetric() const Line 1700  Data::nonsymmetric() const
1700       }       }
1701  }  }
1702    
   
 // 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).  
1703  Data  Data
1704  Data::trace(int axis_offset) const  Data::trace(int axis_offset) const
1705  {  {    
1706       if (isLazy())       if (isLazy())
1707       {       {
1708      Data temp(*this);   // to get around the fact that you can't resolve a const Data      DataLazy* c=new DataLazy(borrowDataPtr(),TRACE,axis_offset);
1709      temp.resolve();      return Data(c);
     return temp.trace(axis_offset);  
1710       }       }
1711       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1712       if (getDataPointRank()==2) {       if (getDataPointRank()==2) {
# Line 1789  Data::transpose(int axis_offset) const Line 1760  Data::transpose(int axis_offset) const
1760  {      {    
1761       if (isLazy())       if (isLazy())
1762       {       {
1763      Data temp(*this);   // to get around the fact that you can't resolve a const Data      DataLazy* c=new DataLazy(borrowDataPtr(),TRANS,axis_offset);
1764      temp.resolve();      return Data(c);
     return temp.transpose(axis_offset);  
1765       }       }
1766       DataTypes::ShapeType s=getDataPointShape();       DataTypes::ShapeType s=getDataPointShape();
1767       DataTypes::ShapeType ev_shape;       DataTypes::ShapeType ev_shape;
# Line 1894  Data::calc_minGlobalDataPoint(int& ProcN Line 1864  Data::calc_minGlobalDataPoint(int& ProcN
1864    double next,local_min;    double next,local_min;
1865    int local_lowi=0,local_lowj=0;        int local_lowi=0,local_lowj=0;    
1866    
1867    #pragma omp parallel private(next,local_min,local_lowi,local_lowj)    #pragma omp parallel firstprivate(local_lowi,local_lowj) private(next,local_min)
1868    {    {
1869      local_min=min;      local_min=min;
1870      #pragma omp for private(i,j) schedule(static)      #pragma omp for private(i,j) schedule(static)
# Line 1921  Data::calc_minGlobalDataPoint(int& ProcN Line 1891  Data::calc_minGlobalDataPoint(int& ProcN
1891      next = temp.getDataPoint(lowi,lowj);      next = temp.getDataPoint(lowi,lowj);
1892      int lowProc = 0;      int lowProc = 0;
1893      double *globalMins = new double[get_MPISize()+1];      double *globalMins = new double[get_MPISize()+1];
1894      int error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );      int error;
1895        error = MPI_Gather ( &next, 1, MPI_DOUBLE, globalMins, 1, MPI_DOUBLE, 0, get_MPIComm() );
1896    
1897      if( get_MPIRank()==0 ){      if( get_MPIRank()==0 ){
1898          next = globalMins[lowProc];          next = globalMins[lowProc];
# Line 2512  escript::C_GeneralTensorProduct(Data& ar Line 2483  escript::C_GeneralTensorProduct(Data& ar
2483    // 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().
2484    
2485    // deal with any lazy data    // deal with any lazy data
2486    if (arg_0.isLazy()) {arg_0.resolve();}  //   if (arg_0.isLazy()) {arg_0.resolve();}
2487    if (arg_1.isLazy()) {arg_1.resolve();}  //   if (arg_1.isLazy()) {arg_1.resolve();}
2488      if (arg_0.isLazy() || arg_1.isLazy())
2489      {
2490        DataLazy* c=new DataLazy(arg_0.borrowDataPtr(), arg_1.borrowDataPtr(), PROD, axis_offset,transpose);
2491        return Data(c);
2492      }
2493    
2494    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2495    Data arg_0_Z, arg_1_Z;    Data arg_0_Z, arg_1_Z;
# Line 2590  escript::C_GeneralTensorProduct(Data& ar Line 2566  escript::C_GeneralTensorProduct(Data& ar
2566       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
2567    }    }
2568    
2569      if (shape2.size()>ESCRIPT_MAX_DATA_RANK)
2570      {
2571         ostringstream os;
2572         os << "C_GeneralTensorProduct: Error - Attempt to create a rank " << shape2.size() << " object. The maximum rank is " << ESCRIPT_MAX_DATA_RANK << ".";
2573         throw DataException(os.str());
2574      }
2575    
2576    // Declare output Data object    // Declare output Data object
2577    Data res;    Data res;
2578    
# Line 2923  Data::borrowReadyPtr() const Line 2906  Data::borrowReadyPtr() const
2906  std::string  std::string
2907  Data::toString() const  Data::toString() const
2908  {  {
2909      if (!m_data->isEmpty() &&      if (!m_data->isEmpty() &&
2910      getNumDataPoints()*getDataPointSize()>escriptParams.getInt("TOO_MANY_LINES"))      !m_data->isLazy() &&
2911        getLength()>escriptParams.getInt("TOO_MANY_LINES"))
2912      {      {
2913      stringstream temp;      stringstream temp;
2914      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.2086

  ViewVC Help
Powered by ViewVC 1.1.26