/[escript]/branches/clazy/escriptcore/src/DataLazy.cpp
ViewVC logotype

Diff of /branches/clazy/escriptcore/src/DataLazy.cpp

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

revision 5997 by caltinay, Mon Feb 29 07:24:47 2016 UTC revision 6074 by jfenwick, Fri Mar 18 02:42:48 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
   
17  #include "DataLazy.h"  #include "DataLazy.h"
18  #include "Data.h"  #include "Data.h"
19  #include "DataTypes.h"  #include "DataTypes.h"
20  #include "EscriptParams.h"  #include "EscriptParams.h"
21  #include "FunctionSpace.h"  #include "FunctionSpace.h"
 #include "UnaryFuncs.h"    // for escript::fsign  
22  #include "Utils.h"  #include "Utils.h"
23    #include "DataMaths.h"
24    
25  #ifdef USE_NETCDF  #ifdef USE_NETCDF
26  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 993  DataLazy::collapse() const Line 990  DataLazy::collapse() const
990    m_op=IDENTITY;    m_op=IDENTITY;
991  }  }
992    
   
   
   
   
   
 #define PROC_OP(TYPE,X)                               \  
         for (int j=0;j<onumsteps;++j)\  
         {\  
           for (int i=0;i<numsteps;++i,resultp+=resultStep) \  
           { \  
 LAZYDEBUG(cout << "[left,right]=[" << lroffset << "," << rroffset << "]" << endl;)\  
 LAZYDEBUG(cout << "{left,right}={" << (*left)[lroffset] << "," << (*right)[rroffset] << "}\n";)\  
              tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X); \  
 LAZYDEBUG(cout << " result=      " << resultp[0] << endl;) \  
              lroffset+=leftstep; \  
              rroffset+=rightstep; \  
           }\  
           lroffset+=oleftstep;\  
           rroffset+=orightstep;\  
         }  
   
   
993  // The result will be stored in m_samples  // The result will be stored in m_samples
994  // The return value is a pointer to the DataVector, offset is the offset within the return value  // The return value is a pointer to the DataVector, offset is the offset within the return value
995  const DataTypes::RealVectorType*  const DataTypes::RealVectorType*
# Line 1028  LAZYDEBUG(cout << "Resolve sample " << t Line 1003  LAZYDEBUG(cout << "Resolve sample " << t
1003    }    }
1004    if (m_op==IDENTITY)      if (m_op==IDENTITY)  
1005    {    {
1006      const ValueType& vec=m_id->getVectorRO();      const RealVectorType& vec=m_id->getVectorRO();
1007      roffset=m_id->getPointOffset(sampleNo, 0);      roffset=m_id->getPointOffset(sampleNo, 0);
1008  #ifdef LAZY_STACK_PROF  #ifdef LAZY_STACK_PROF
1009  int x;  int x;
# Line 1085  DataLazy::resolveNodeUnary(int tid, int Line 1060  DataLazy::resolveNodeUnary(int tid, int
1060    const double* left=&((*leftres)[roffset]);    const double* left=&((*leftres)[roffset]);
1061    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1062    double* result=&(m_samples[roffset]);    double* result=&(m_samples[roffset]);
1063      escript::ESFunction operation=SINF;
1064    switch (m_op)    switch (m_op)
1065    {    {
1066      case SIN:        case SIN:
1067          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);      operation=SINF;
1068          break;      break;
1069      case COS:      case COS:
1070          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);          operation=COSF;
1071          break;      break;
1072      case TAN:      case TAN:
1073          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);          operation=TANF;
1074          break;      break;
1075      case ASIN:      case ASIN:
1076          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);          operation=ASINF;
1077          break;      break;
1078      case ACOS:      case ACOS:
1079          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);          operation=ACOSF;
1080          break;      break;
1081      case ATAN:      case ATAN:
1082          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);          operation=ATANF;
1083          break;      break;
1084      case SINH:      case SINH:
1085          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);          operation=SINHF;
1086          break;      break;
1087      case COSH:      case COSH:
1088          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);          operation=COSHF;
1089          break;      break;
1090      case TANH:      case TANH:
1091          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);          operation=TANHF;
1092          break;      break;
1093      case ERF:      case ERF:
1094  #if defined (_WIN32) && !defined(__INTEL_COMPILER)          operation=ERFF;
1095          throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");      break;
 #else  
         tensor_unary_operation(m_samplesize, left, result, ::erf);  
         break;  
 #endif  
1096     case ASINH:     case ASINH:
1097  #if defined (_WIN32) && !defined(__INTEL_COMPILER)          operation=ASINHF;
1098          tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);      break;
 #else  
         tensor_unary_operation(m_samplesize, left, result, ::asinh);  
 #endif    
         break;  
1099     case ACOSH:     case ACOSH:
1100  #if defined (_WIN32) && !defined(__INTEL_COMPILER)          operation=ACOSHF;
1101          tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);      break;
 #else  
         tensor_unary_operation(m_samplesize, left, result, ::acosh);  
 #endif    
         break;  
1102     case ATANH:     case ATANH:
1103  #if defined (_WIN32) && !defined(__INTEL_COMPILER)          operation=ATANHF;
1104          tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);      break;
 #else  
         tensor_unary_operation(m_samplesize, left, result, ::atanh);  
 #endif    
         break;  
1105      case LOG10:      case LOG10:
1106          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);          operation=LOG10F;
1107          break;      break;
1108      case LOG:      case LOG:
1109          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);          operation=LOGF;
1110          break;      break;
1111      case SIGN:      case SIGN:
1112          tensor_unary_operation(m_samplesize, left, result, escript::fsign);          operation=SIGNF;
1113          break;      break;
1114      case ABS:      case ABS:
1115          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);          operation=ABSF;
1116          break;      break;
1117      case NEG:      case NEG:
1118          tensor_unary_operation(m_samplesize, left, result, negate<double>());          operation=NEGF;
1119          break;      break;
1120      case POS:      case POS:
1121          // it doesn't mean anything for delayed.          // it doesn't mean anything for delayed.
1122          // it will just trigger a deep copy of the lazy object          // it will just trigger a deep copy of the lazy object
1123          throw DataException("Programmer error - POS not supported for lazy data.");          throw DataException("Programmer error - POS not supported for lazy data.");
1124          break;          break;
1125      case EXP:      case EXP:
1126          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);          operation=EXPF;
1127          break;      break;
1128      case SQRT:      case SQRT:
1129          tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);          operation=SQRTF;
1130          break;      break;
1131      case RECIP:      case RECIP:
1132          tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));          operation=INVF;
1133          break;      break;
1134      case GZ:      case GZ:
1135          tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));          operation=GTZEROF;
1136          break;      break;
1137      case LZ:      case LZ:
1138          tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));          operation=LTZEROF;
1139          break;      break;
1140      case GEZ:      case GEZ:
1141          tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));          operation=GEZEROF;
1142          break;      break;
1143      case LEZ:      case LEZ:
1144          tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));          operation=LEZEROF;
1145          break;      break;
1146  // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently  // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1147      case NEZ:      case NEZ:
1148          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));          operation=NEQZEROF;
1149          break;      break;
1150      case EZ:      case EZ:
1151          tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));          operation=EQZEROF;
1152          break;      break;
   
1153      default:      default:
1154          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1155    }    }
1156      tensor_unary_array_operation(m_samplesize,
1157                                 left,
1158                                 result,
1159                                 operation,
1160                                 m_tol);  
1161    return &(m_samples);    return &(m_samples);
1162  }  }
1163    
# Line 1265  DataLazy::resolveNodeNP1OUT(int tid, int Line 1229  DataLazy::resolveNodeNP1OUT(int tid, int
1229      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");      throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1230    }    }
1231    size_t subroffset;    size_t subroffset;
1232    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1233    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1234    size_t loop=0;    size_t loop=0;
1235    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;    size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
# Line 1311  DataLazy::resolveNodeNP1OUT_P(int tid, i Line 1275  DataLazy::resolveNodeNP1OUT_P(int tid, i
1275    }    }
1276    size_t subroffset;    size_t subroffset;
1277    size_t offset;    size_t offset;
1278    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1279    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1280    offset=roffset;    offset=roffset;
1281    size_t loop=0;    size_t loop=0;
# Line 1356  DataLazy::resolveNodeNP1OUT_2P(int tid, Line 1320  DataLazy::resolveNodeNP1OUT_2P(int tid,
1320    }    }
1321    size_t subroffset;    size_t subroffset;
1322    size_t offset;    size_t offset;
1323    const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1324    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1325    offset=roffset;    offset=roffset;
1326    size_t loop=0;    size_t loop=0;
# Line 1392  DataLazy::resolveNodeCondEval(int tid, i Line 1356  DataLazy::resolveNodeCondEval(int tid, i
1356    }    }
1357    size_t subroffset;    size_t subroffset;
1358    
1359    const ValueType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);    const RealVectorType* maskres=m_mask->resolveNodeSample(tid, sampleNo, subroffset);
1360    const ValueType* srcres=0;    const RealVectorType* srcres=0;
1361    if ((*maskres)[subroffset]>0)    if ((*maskres)[subroffset]>0)
1362    {    {
1363          srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);          srcres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
# Line 1536  LAZYDEBUG(cout << "Resolve binary: " << Line 1500  LAZYDEBUG(cout << "Resolve binary: " <<
1500    
1501    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0    int resultStep=max(leftstep,rightstep);       // only one (at most) should be !=0
1502          // Get the values of sub-expressions          // Get the values of sub-expressions
1503    const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);          const RealVectorType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);      
1504    const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
1505  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)  LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
1506  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)  LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
1507  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)  LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
# Line 1555  LAZYDEBUG(cout << "Right res["<< rroffse Line 1519  LAZYDEBUG(cout << "Right res["<< rroffse
1519    switch(m_op)    switch(m_op)
1520    {    {
1521      case ADD:      case ADD:
1522          PROC_OP(NO_ARG,plus<double>());          //PROC_OP(NO_ARG,plus<double>());
1523          DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1524                 &(*left)[0],
1525                 &(*right)[0],
1526                 chunksize,
1527                 onumsteps,
1528                 numsteps,
1529                 resultStep,
1530                 leftstep,
1531                 rightstep,
1532                 oleftstep,
1533                 orightstep,
1534                 lroffset,
1535                 rroffset,
1536                 escript::ESFunction::PLUSF);  
1537          break;          break;
1538      case SUB:      case SUB:
1539          PROC_OP(NO_ARG,minus<double>());        DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1540                 &(*left)[0],
1541                 &(*right)[0],
1542                 chunksize,
1543                 onumsteps,
1544                 numsteps,
1545                 resultStep,
1546                 leftstep,
1547                 rightstep,
1548                 oleftstep,
1549                 orightstep,
1550                 lroffset,
1551                 rroffset,
1552                 escript::ESFunction::MINUSF);        
1553            //PROC_OP(NO_ARG,minus<double>());
1554          break;          break;
1555      case MUL:      case MUL:
1556          PROC_OP(NO_ARG,multiplies<double>());          //PROC_OP(NO_ARG,multiplies<double>());
1557          DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1558                 &(*left)[0],
1559                 &(*right)[0],
1560                 chunksize,
1561                 onumsteps,
1562                 numsteps,
1563                 resultStep,
1564                 leftstep,
1565                 rightstep,
1566                 oleftstep,
1567                 orightstep,
1568                 lroffset,
1569                 rroffset,
1570                 escript::ESFunction::MULTIPLIESF);      
1571          break;          break;
1572      case DIV:      case DIV:
1573          PROC_OP(NO_ARG,divides<double>());          //PROC_OP(NO_ARG,divides<double>());
1574          DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1575                 &(*left)[0],
1576                 &(*right)[0],
1577                 chunksize,
1578                 onumsteps,
1579                 numsteps,
1580                 resultStep,
1581                 leftstep,
1582                 rightstep,
1583                 oleftstep,
1584                 orightstep,
1585                 lroffset,
1586                 rroffset,
1587                 escript::ESFunction::DIVIDESF);          
1588          break;          break;
1589      case POW:      case POW:
1590         PROC_OP(double (double,double),::pow);         //PROC_OP(double (double,double),::pow);
1591          DataMaths::binaryOpVectorLazyHelper<real_t, real_t, real_t>(resultp,
1592                 &(*left)[0],
1593                 &(*right)[0],
1594                 chunksize,
1595                 onumsteps,
1596                 numsteps,
1597                 resultStep,
1598                 leftstep,
1599                 rightstep,
1600                 oleftstep,
1601                 orightstep,
1602                 lroffset,
1603                 rroffset,
1604                 escript::ESFunction::POWF);          
1605          break;          break;
1606      default:      default:
1607          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");          throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
# Line 1597  LAZYDEBUG(cout << "Resolve TensorProduct Line 1631  LAZYDEBUG(cout << "Resolve TensorProduct
1631    roffset=m_samplesize*tid;    roffset=m_samplesize*tid;
1632    size_t offset=roffset;    size_t offset=roffset;
1633    
1634    const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);    const RealVectorType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
1635    
1636    const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);    const RealVectorType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
1637    
1638  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;  LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
1639  cout << getNoValues() << endl;)  cout << getNoValues() << endl;)
# Line 1736  DataLazy::resolveGroupWorker(std::vector Line 1770  DataLazy::resolveGroupWorker(std::vector
1770    {             // it is possible that dats[0] is one of the objects which we discarded and    {             // it is possible that dats[0] is one of the objects which we discarded and
1771                  // all the other functionspaces match.                  // all the other functionspaces match.
1772          vector<DataExpanded*> dep;          vector<DataExpanded*> dep;
1773          vector<ValueType*> vecs;          vector<RealVectorType*> vecs;
1774          for (int i=0;i<work.size();++i)          for (int i=0;i<work.size();++i)
1775          {          {
1776                  dep.push_back(new DataExpanded(fs,work[i]->getShape(), ValueType(work[i]->getNoValues())));                  dep.push_back(new DataExpanded(fs,work[i]->getShape(), RealVectorType(work[i]->getNoValues())));
1777                  vecs.push_back(&(dep[i]->getVectorRW()));                  vecs.push_back(&(dep[i]->getVectorRW()));
1778          }          }
1779          int totalsamples=work[0]->getNumSamples();          int totalsamples=work[0]->getNumSamples();
1780          const ValueType* res=0; // Storage for answer          const RealVectorType* res=0; // Storage for answer
1781          int sample;          int sample;
1782          #pragma omp parallel private(sample, res)          #pragma omp parallel private(sample, res)
1783          {          {
# Line 1795  DataLazy::resolveNodeWorker() Line 1829  DataLazy::resolveNodeWorker()
1829      return m_id;      return m_id;
1830    }    }
1831          // from this point on we must have m_op!=IDENTITY and m_readytype=='E'          // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
1832    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));    DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  RealVectorType(getNoValues()));
1833    ValueType& resvec=result->getVectorRW();    RealVectorType& resvec=result->getVectorRW();
1834    DataReady_ptr resptr=DataReady_ptr(result);    DataReady_ptr resptr=DataReady_ptr(result);
1835    
1836    int sample;    int sample;
1837    int totalsamples=getNumSamples();    int totalsamples=getNumSamples();
1838    const ValueType* res=0;       // Storage for answer    const RealVectorType* res=0;       // Storage for answer
1839  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)  LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
1840    #pragma omp parallel private(sample,res)    #pragma omp parallel private(sample,res)
1841    {    {

Legend:
Removed from v.5997  
changed lines
  Added in v.6074

  ViewVC Help
Powered by ViewVC 1.1.26