/[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 2497 by jfenwick, Mon Jun 29 00:04:45 2009 UTC revision 2500 by jfenwick, Tue Jun 30 00:42:38 2009 UTC
# Line 435  opToString(ES_optype op) Line 435  opToString(ES_optype op)
435    return ES_opstrings[op];    return ES_opstrings[op];
436  }  }
437    
438    #ifdef LAZY_NODE_STORAGE
439    void DataLazy::LazyNodeSetup()
440    {
441    #ifdef _OPENMP
442        int numthreads=omp_get_max_threads();
443        m_samples.resize(numthreads*m_samplesize);
444        m_sampleids=new int[numthreads];
445        for (int i=0;i<numthreads;++i)
446        {
447            m_sampleids[i]=-1;  
448        }
449    #else
450        m_samples.resize(m_samplesize);
451        m_sampleids=new int[1];
452        m_sampleids[0]=-1;
453    #endif  // _OPENMP
454    }
455    #endif   // LAZY_NODE_STORAGE
456    
457    
458    // Creates an identity node
459  DataLazy::DataLazy(DataAbstract_ptr p)  DataLazy::DataLazy(DataAbstract_ptr p)
460      : parent(p->getFunctionSpace(),p->getShape())      : parent(p->getFunctionSpace(),p->getShape())
461    #ifdef LAZY_NODE_STORAGE
462        ,m_sampleids(0),
463        m_samples(1)
464    #endif
465  {  {
466     if (p->isLazy())     if (p->isLazy())
467     {     {
# Line 456  LAZYDEBUG(cout << "Wrapping " << dr.get( Line 480  LAZYDEBUG(cout << "Wrapping " << dr.get(
480  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(1)Lazy created with " << m_samplesize << endl;)
481  }  }
482    
   
   
   
483  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)  DataLazy::DataLazy(DataAbstract_ptr left, ES_optype op)
484      : parent(left->getFunctionSpace(),left->getShape()),      : parent(left->getFunctionSpace(),left->getShape()),
485      m_op(op),      m_op(op),
# Line 487  DataLazy::DataLazy(DataAbstract_ptr left Line 508  DataLazy::DataLazy(DataAbstract_ptr left
508     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
509     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
510     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
511    #ifdef LAZY_NODE_STORAGE
512       LazyNodeSetup();
513    #endif
514     SIZELIMIT     SIZELIMIT
515  }  }
516    
# Line 557  LAZYDEBUG(cout << "Right " << right.get( Line 581  LAZYDEBUG(cout << "Right " << right.get(
581     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
582     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
583     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
584    #ifdef LAZY_NODE_STORAGE
585       LazyNodeSetup();
586    #endif
587     SIZELIMIT     SIZELIMIT
588  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(3)Lazy created with " << m_samplesize << endl;)
589  }  }
# Line 624  DataLazy::DataLazy(DataAbstract_ptr left Line 651  DataLazy::DataLazy(DataAbstract_ptr left
651     m_buffsRequired=calcBuffs(m_left, m_right,m_op);     m_buffsRequired=calcBuffs(m_left, m_right,m_op);
652     m_children=m_left->m_children+m_right->m_children+2;     m_children=m_left->m_children+m_right->m_children+2;
653     m_height=max(m_left->m_height,m_right->m_height)+1;     m_height=max(m_left->m_height,m_right->m_height)+1;
654    #ifdef LAZY_NODE_STORAGE
655       LazyNodeSetup();
656    #endif
657     SIZELIMIT     SIZELIMIT
658  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(4)Lazy created with " << m_samplesize << endl;)
659  }  }
# Line 656  DataLazy::DataLazy(DataAbstract_ptr left Line 686  DataLazy::DataLazy(DataAbstract_ptr left
686     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
687     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
688     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
689    #ifdef LAZY_NODE_STORAGE
690       LazyNodeSetup();
691    #endif
692     SIZELIMIT     SIZELIMIT
693  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(5)Lazy created with " << m_samplesize << endl;)
694  }  }
# Line 687  DataLazy::DataLazy(DataAbstract_ptr left Line 720  DataLazy::DataLazy(DataAbstract_ptr left
720     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
721     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
722     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
723    #ifdef LAZY_NODE_STORAGE
724       LazyNodeSetup();
725    #endif
726     SIZELIMIT     SIZELIMIT
727  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(6)Lazy created with " << m_samplesize << endl;)
728  }  }
# Line 719  DataLazy::DataLazy(DataAbstract_ptr left Line 755  DataLazy::DataLazy(DataAbstract_ptr left
755     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());     m_maxsamplesize=max(m_samplesize,m_left->getMaxSampleSize());
756     m_children=m_left->m_children+1;     m_children=m_left->m_children+1;
757     m_height=m_left->m_height+1;     m_height=m_left->m_height+1;
758    #ifdef LAZY_NODE_STORAGE
759       LazyNodeSetup();
760    #endif
761     SIZELIMIT     SIZELIMIT
762  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)  LAZYDEBUG(cout << "(7)Lazy created with " << m_samplesize << endl;)
763  }  }
764    
765  DataLazy::~DataLazy()  DataLazy::~DataLazy()
766  {  {
767    #ifdef LAZY_NODE_SETUP
768       delete[] m_sampleids;
769       delete[] m_samples;
770    #endif
771  }  }
772    
773    
# Line 1543  cout << "\nWritten to: " << resultp << " Line 1586  cout << "\nWritten to: " << resultp << "
1586  }  }
1587    
1588    
1589    #ifdef LAZY_NODE_STORAGE
1590    
1591    // The result will be stored in m_samples
1592    // The return value is a pointer to the DataVector, offset is the offset within the return value
1593    const DataTypes::ValueType*
1594    DataLazy::resolveNodeSample(int tid, int sampleNo, size_t& roffset)
1595    {
1596    ENABLEDEBUG
1597    LAZYDEBUG(cout << "Resolve sample " << toString() << endl;)
1598        // collapse so we have a 'E' node or an IDENTITY for some other type
1599      if (m_readytype!='E' && m_op!=IDENTITY)
1600      {
1601        collapse();
1602      }
1603      if (m_op==IDENTITY)  
1604      {
1605        const ValueType& vec=m_id->getVectorRO();
1606    //     if (m_readytype=='C')
1607    //     {
1608    //  roffset=0;      // all samples read from the same position
1609    //  return &(m_samples);
1610    //     }
1611        roffset=m_id->getPointOffset(sampleNo, 0);
1612        return &(vec);
1613      }
1614      if (m_readytype!='E')
1615      {
1616        throw DataException("Programmer Error - Collapse did not produce an expanded node.");
1617      }
1618      switch (getOpgroup(m_op))
1619      {
1620      case G_UNARY:
1621      case G_UNARY_P: return resolveNodeUnary(tid, sampleNo, roffset);
1622      case G_BINARY: return resolveNodeBinary(tid, sampleNo, roffset);
1623      case G_NP1OUT: return resolveNodeNP1OUT(tid, sampleNo, roffset);
1624      case G_NP1OUT_P: return resolveNodeNP1OUT_P(tid, sampleNo, roffset);
1625      case G_TENSORPROD: return resolveNodeTProd(tid, sampleNo, roffset);
1626      case G_NP1OUT_2P: return resolveNodeNP1OUT_2P(tid, sampleNo, roffset);
1627      default:
1628        throw DataException("Programmer Error - resolveSample does not know how to process "+opToString(m_op)+".");
1629      }
1630    }
1631    
1632    const DataTypes::ValueType*
1633    DataLazy::resolveNodeUnary(int tid, int sampleNo, size_t& roffset)
1634    {
1635        // we assume that any collapsing has been done before we get here
1636        // since we only have one argument we don't need to think about only
1637        // processing single points.
1638        // we will also know we won't get identity nodes
1639      if (m_readytype!='E')
1640      {
1641        throw DataException("Programmer error - resolveUnary should only be called on expanded Data.");
1642      }
1643      if (m_op==IDENTITY)
1644      {
1645        throw DataException("Programmer error - resolveNodeUnary should not be called on identity nodes.");
1646      }
1647      if (m_sampleids[tid]==sampleNo)
1648      {
1649        roffset=tid*m_samplesize;
1650        return &(m_samples);        // sample is already resolved
1651      }
1652      const DataTypes::ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, roffset);
1653      const double* left=&((*leftres)[roffset]);
1654      roffset=m_samplesize*tid;
1655      double* result=&(m_samples[roffset]);
1656      switch (m_op)
1657      {
1658        case SIN:  
1659        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sin);
1660        break;
1661        case COS:
1662        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cos);
1663        break;
1664        case TAN:
1665        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tan);
1666        break;
1667        case ASIN:
1668        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::asin);
1669        break;
1670        case ACOS:
1671        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::acos);
1672        break;
1673        case ATAN:
1674        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::atan);
1675        break;
1676        case SINH:
1677        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sinh);
1678        break;
1679        case COSH:
1680        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::cosh);
1681        break;
1682        case TANH:
1683        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::tanh);
1684        break;
1685        case ERF:
1686    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1687        throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1688    #else
1689        tensor_unary_operation(m_samplesize, left, result, ::erf);
1690        break;
1691    #endif
1692       case ASINH:
1693    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1694        tensor_unary_operation(m_samplesize, left, result, escript::asinh_substitute);
1695    #else
1696        tensor_unary_operation(m_samplesize, left, result, ::asinh);
1697    #endif  
1698        break;
1699       case ACOSH:
1700    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1701        tensor_unary_operation(m_samplesize, left, result, escript::acosh_substitute);
1702    #else
1703        tensor_unary_operation(m_samplesize, left, result, ::acosh);
1704    #endif  
1705        break;
1706       case ATANH:
1707    #if defined (_WIN32) && !defined(__INTEL_COMPILER)
1708        tensor_unary_operation(m_samplesize, left, result, escript::atanh_substitute);
1709    #else
1710        tensor_unary_operation(m_samplesize, left, result, ::atanh);
1711    #endif  
1712        break;
1713        case LOG10:
1714        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log10);
1715        break;
1716        case LOG:
1717        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::log);
1718        break;
1719        case SIGN:
1720        tensor_unary_operation(m_samplesize, left, result, escript::fsign);
1721        break;
1722        case ABS:
1723        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::fabs);
1724        break;
1725        case NEG:
1726        tensor_unary_operation(m_samplesize, left, result, negate<double>());
1727        break;
1728        case POS:
1729        // it doesn't mean anything for delayed.
1730        // it will just trigger a deep copy of the lazy object
1731        throw DataException("Programmer error - POS not supported for lazy data.");
1732        break;
1733        case EXP:
1734        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::exp);
1735        break;
1736        case SQRT:
1737        tensor_unary_operation<double (*)(double)>(m_samplesize, left, result, ::sqrt);
1738        break;
1739        case RECIP:
1740        tensor_unary_operation(m_samplesize, left, result, bind1st(divides<double>(),1.));
1741        break;
1742        case GZ:
1743        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater<double>(),0.0));
1744        break;
1745        case LZ:
1746        tensor_unary_operation(m_samplesize, left, result, bind2nd(less<double>(),0.0));
1747        break;
1748        case GEZ:
1749        tensor_unary_operation(m_samplesize, left, result, bind2nd(greater_equal<double>(),0.0));
1750        break;
1751        case LEZ:
1752        tensor_unary_operation(m_samplesize, left, result, bind2nd(less_equal<double>(),0.0));
1753        break;
1754    // There are actually G_UNARY_P but I don't see a compelling reason to treat them differently
1755        case NEZ:
1756        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsGT(),m_tol));
1757        break;
1758        case EZ:
1759        tensor_unary_operation(m_samplesize, left, result, bind2nd(AbsLTE(),m_tol));
1760        break;
1761    
1762        default:
1763        throw DataException("Programmer error - resolveUnary can not resolve operator "+opToString(m_op)+".");
1764      }
1765      return &(m_samples);
1766    }
1767    
1768    
1769    const DataTypes::ValueType*
1770    DataLazy::resolveNodeNP1OUT(int tid, int sampleNo, size_t& roffset)
1771    {
1772        // we assume that any collapsing has been done before we get here
1773        // since we only have one argument we don't need to think about only
1774        // processing single points.
1775      if (m_readytype!='E')
1776      {
1777        throw DataException("Programmer error - resolveNodeNP1OUT should only be called on expanded Data.");
1778      }
1779      if (m_op==IDENTITY)
1780      {
1781        throw DataException("Programmer error - resolveNodeNP1OUT should not be called on identity nodes.");
1782      }
1783      size_t subroffset;
1784      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1785      roffset=m_samplesize*tid;
1786      size_t loop=0;
1787      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1788      size_t step=getNoValues();
1789      size_t offset=roffset;
1790      switch (m_op)
1791      {
1792        case SYM:
1793        for (loop=0;loop<numsteps;++loop)
1794        {
1795            DataMaths::symmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1796            subroffset+=step;
1797            offset+=step;
1798        }
1799        break;
1800        case NSYM:
1801        for (loop=0;loop<numsteps;++loop)
1802        {
1803            DataMaths::nonsymmetric(*leftres,m_left->getShape(),subroffset, m_samples, getShape(), offset);
1804            subroffset+=step;
1805            offset+=step;
1806        }
1807        break;
1808        default:
1809        throw DataException("Programmer error - resolveNP1OUT can not resolve operator "+opToString(m_op)+".");
1810      }
1811      return &m_samples;
1812    }
1813    
1814    const DataTypes::ValueType*
1815    DataLazy::resolveNodeNP1OUT_P(int tid, int sampleNo, size_t& roffset)
1816    {
1817        // we assume that any collapsing has been done before we get here
1818        // since we only have one argument we don't need to think about only
1819        // processing single points.
1820      if (m_readytype!='E')
1821      {
1822        throw DataException("Programmer error - resolveNodeNP1OUT_P should only be called on expanded Data.");
1823      }
1824      if (m_op==IDENTITY)
1825      {
1826        throw DataException("Programmer error - resolveNodeNP1OUT_P should not be called on identity nodes.");
1827      }
1828      size_t subroffset;
1829      size_t offset;
1830      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1831      roffset=m_samplesize*tid;
1832      offset=roffset;
1833      size_t loop=0;
1834      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1835      size_t outstep=getNoValues();
1836      size_t instep=m_left->getNoValues();
1837      switch (m_op)
1838      {
1839        case TRACE:
1840        for (loop=0;loop<numsteps;++loop)
1841        {
1842                DataMaths::trace(*leftres,m_left->getShape(),subroffset, m_samples ,getShape(),offset,m_axis_offset);
1843            subroffset+=instep;
1844            offset+=outstep;
1845        }
1846        break;
1847        case TRANS:
1848        for (loop=0;loop<numsteps;++loop)
1849        {
1850                DataMaths::transpose(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset,m_axis_offset);
1851            subroffset+=instep;
1852            offset+=outstep;
1853        }
1854        break;
1855        default:
1856        throw DataException("Programmer error - resolveNP1OUTP can not resolve operator "+opToString(m_op)+".");
1857      }
1858      return &m_samples;
1859    }
1860    
1861    
1862    const DataTypes::ValueType*
1863    DataLazy::resolveNodeNP1OUT_2P(int tid, int sampleNo, size_t& roffset)
1864    {
1865      if (m_readytype!='E')
1866      {
1867        throw DataException("Programmer error - resolveNodeNP1OUT_2P should only be called on expanded Data.");
1868      }
1869      if (m_op==IDENTITY)
1870      {
1871        throw DataException("Programmer error - resolveNodeNP1OUT_2P should not be called on identity nodes.");
1872      }
1873      size_t subroffset;
1874      size_t offset;
1875      const ValueType* leftres=m_left->resolveNodeSample(tid, sampleNo, subroffset);
1876      roffset=m_samplesize*tid;
1877      offset=roffset;
1878      size_t loop=0;
1879      size_t numsteps=(m_readytype=='E')?getNumDPPSample():1;
1880      size_t outstep=getNoValues();
1881      size_t instep=m_left->getNoValues();
1882      switch (m_op)
1883      {
1884        case SWAP:
1885        for (loop=0;loop<numsteps;++loop)
1886        {
1887                DataMaths::swapaxes(*leftres,m_left->getShape(),subroffset, m_samples, getShape(),offset, m_axis_offset, m_transpose);
1888            subroffset+=instep;
1889            offset+=outstep;
1890        }
1891        break;
1892        default:
1893        throw DataException("Programmer error - resolveNodeNP1OUT2P can not resolve operator "+opToString(m_op)+".");
1894      }
1895      return &m_samples;
1896    }
1897    
1898    
1899    
1900    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
1901    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
1902    // If both children are expanded, then we can process them in a single operation (we treat
1903    // the whole sample as one big datapoint.
1904    // If one of the children is not expanded, then we need to treat each point in the sample
1905    // individually.
1906    // There is an additional complication when scalar operations are considered.
1907    // For example, 2+Vector.
1908    // In this case each double within the point is treated individually
1909    const DataTypes::ValueType*
1910    DataLazy::resolveNodeBinary(int tid, int sampleNo, size_t& roffset)
1911    {
1912    LAZYDEBUG(cout << "Resolve binary: " << toString() << endl;)
1913    
1914      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
1915        // first work out which of the children are expanded
1916      bool leftExp=(m_left->m_readytype=='E');
1917      bool rightExp=(m_right->m_readytype=='E');
1918      if (!leftExp && !rightExp)
1919      {
1920        throw DataException("Programmer Error - please use collapse if neither argument has type 'E'.");
1921      }
1922      bool leftScalar=(m_left->getRank()==0);
1923      bool rightScalar=(m_right->getRank()==0);
1924      if ((m_left->getRank()!=m_right->getRank()) && (!leftScalar && !rightScalar))
1925      {
1926        throw DataException("resolveBinary - ranks of arguments must match unless one of them is scalar.");
1927      }
1928      size_t leftsize=m_left->getNoValues();
1929      size_t rightsize=m_right->getNoValues();
1930      size_t chunksize=1;           // how many doubles will be processed in one go
1931      int leftstep=0;       // how far should the left offset advance after each step
1932      int rightstep=0;
1933      int numsteps=0;       // total number of steps for the inner loop
1934      int oleftstep=0;  // the o variables refer to the outer loop
1935      int orightstep=0; // The outer loop is only required in cases where there is an extended scalar
1936      int onumsteps=1;
1937      
1938      bool LES=(leftExp && leftScalar); // Left is an expanded scalar
1939      bool RES=(rightExp && rightScalar);
1940      bool LS=(!leftExp && leftScalar); // left is a single scalar
1941      bool RS=(!rightExp && rightScalar);
1942      bool LN=(!leftExp && !leftScalar);    // left is a single non-scalar
1943      bool RN=(!rightExp && !rightScalar);
1944      bool LEN=(leftExp && !leftScalar);    // left is an expanded non-scalar
1945      bool REN=(rightExp && !rightScalar);
1946    
1947      if ((LES && RES) || (LEN && REN)) // both are Expanded scalars or both are expanded non-scalars
1948      {
1949        chunksize=m_left->getNumDPPSample()*leftsize;
1950        leftstep=0;
1951        rightstep=0;
1952        numsteps=1;
1953      }
1954      else if (LES || RES)
1955      {
1956        chunksize=1;
1957        if (LES)        // left is an expanded scalar
1958        {
1959            if (RS)
1960            {
1961               leftstep=1;
1962               rightstep=0;
1963               numsteps=m_left->getNumDPPSample();
1964            }
1965            else        // RN or REN
1966            {
1967               leftstep=0;
1968               oleftstep=1;
1969               rightstep=1;
1970               orightstep=(RN ? -(int)rightsize : 0);
1971               numsteps=rightsize;
1972               onumsteps=m_left->getNumDPPSample();
1973            }
1974        }
1975        else        // right is an expanded scalar
1976        {
1977            if (LS)
1978            {
1979               rightstep=1;
1980               leftstep=0;
1981               numsteps=m_right->getNumDPPSample();
1982            }
1983            else
1984            {
1985               rightstep=0;
1986               orightstep=1;
1987               leftstep=1;
1988               oleftstep=(LN ? -(int)leftsize : 0);
1989               numsteps=leftsize;
1990               onumsteps=m_right->getNumDPPSample();
1991            }
1992        }
1993      }
1994      else  // this leaves (LEN, RS), (LEN, RN) and their transposes
1995      {
1996        if (LEN)    // and Right will be a single value
1997        {
1998            chunksize=rightsize;
1999            leftstep=rightsize;
2000            rightstep=0;
2001            numsteps=m_left->getNumDPPSample();
2002            if (RS)
2003            {
2004               numsteps*=leftsize;
2005            }
2006        }
2007        else    // REN
2008        {
2009            chunksize=leftsize;
2010            rightstep=leftsize;
2011            leftstep=0;
2012            numsteps=m_right->getNumDPPSample();
2013            if (LS)
2014            {
2015               numsteps*=rightsize;
2016            }
2017        }
2018      }
2019    
2020      int resultStep=max(leftstep,rightstep);   // only one (at most) should be !=0
2021        // Get the values of sub-expressions
2022      const ValueType* left=m_left->resolveNodeSample(tid,sampleNo,lroffset);  
2023      const ValueType* right=m_right->resolveNodeSample(tid,sampleNo,rroffset);
2024    LAZYDEBUG(cout << "Post sub calls in " << toString() << endl;)
2025    LAZYDEBUG(cout << "shapes=" << DataTypes::shapeToString(m_left->getShape()) << "," << DataTypes::shapeToString(m_right->getShape()) << endl;)
2026    LAZYDEBUG(cout << "chunksize=" << chunksize << endl << "leftstep=" << leftstep << " rightstep=" << rightstep;)
2027    LAZYDEBUG(cout << " numsteps=" << numsteps << endl << "oleftstep=" << oleftstep << " orightstep=" << orightstep;)
2028    LAZYDEBUG(cout << "onumsteps=" << onumsteps << endl;)
2029    LAZYDEBUG(cout << " DPPS=" << m_left->getNumDPPSample() << "," <<m_right->getNumDPPSample() << endl;)
2030    LAZYDEBUG(cout << "" << LS << RS << LN << RN << LES << RES <<LEN << REN <<   endl;)
2031    
2032    LAZYDEBUG(cout << "Left res["<< lroffset<< "]=" << (*left)[lroffset] << endl;)
2033    LAZYDEBUG(cout << "Right res["<< rroffset<< "]=" << (*right)[rroffset] << endl;)
2034    
2035    
2036      roffset=m_samplesize*tid;
2037      double* resultp=&(m_samples[roffset]);        // results are stored at the vector offset we recieved
2038      switch(m_op)
2039      {
2040        case ADD:
2041            PROC_OP(NO_ARG,plus<double>());
2042        break;
2043        case SUB:
2044        PROC_OP(NO_ARG,minus<double>());
2045        break;
2046        case MUL:
2047        PROC_OP(NO_ARG,multiplies<double>());
2048        break;
2049        case DIV:
2050        PROC_OP(NO_ARG,divides<double>());
2051        break;
2052        case POW:
2053           PROC_OP(double (double,double),::pow);
2054        break;
2055        default:
2056        throw DataException("Programmer error - resolveBinary can not resolve operator "+opToString(m_op)+".");
2057      }
2058    LAZYDEBUG(cout << "Result res[" << roffset<< "]" << m_samples[roffset] << endl;)
2059      return &m_samples;
2060    }
2061    
2062    
2063    // This method assumes that any subexpressions which evaluate to Constant or Tagged Data
2064    // have already been collapsed to IDENTITY. So we must have at least one expanded child.
2065    // unlike the other resolve helpers, we must treat these datapoints separately.
2066    const DataTypes::ValueType*
2067    DataLazy::resolveNodeTProd(int tid, int sampleNo, size_t& roffset)
2068    {
2069    LAZYDEBUG(cout << "Resolve TensorProduct: " << toString() << endl;)
2070    
2071      size_t lroffset=0, rroffset=0;    // offsets in the left and right result vectors
2072        // first work out which of the children are expanded
2073      bool leftExp=(m_left->m_readytype=='E');
2074      bool rightExp=(m_right->m_readytype=='E');
2075      int steps=getNumDPPSample();
2076      int leftStep=(leftExp? m_left->getNoValues() : 0);        // do not have scalars as input to this method
2077      int rightStep=(rightExp?m_right->getNoValues() : 0);
2078    
2079      int resultStep=getNoValues();
2080      roffset=m_samplesize*tid;
2081      size_t offset=roffset;
2082    
2083      const ValueType* left=m_left->resolveNodeSample(tid, sampleNo, lroffset);
2084    
2085      const ValueType* right=m_right->resolveNodeSample(tid, sampleNo, rroffset);
2086    
2087    LAZYDEBUG(cerr << "[Left shape]=" << DataTypes::shapeToString(m_left->getShape()) << "\n[Right shape]=" << DataTypes::shapeToString(m_right->getShape()) << " result=" <<DataTypes::shapeToString(getShape()) <<  endl;
2088    cout << getNoValues() << endl;)
2089    
2090    
2091    LAZYDEBUG(cerr << "Post sub calls: " << toString() << endl;)
2092    LAZYDEBUG(cout << "LeftExp=" << leftExp << " rightExp=" << rightExp << endl;)
2093    LAZYDEBUG(cout << "LeftR=" << m_left->getRank() << " rightExp=" << m_right->getRank() << endl;)
2094    LAZYDEBUG(cout << "LeftSize=" << m_left->getNoValues() << " RightSize=" << m_right->getNoValues() << endl;)
2095    LAZYDEBUG(cout << "m_samplesize=" << m_samplesize << endl;)
2096    LAZYDEBUG(cout << "outputshape=" << DataTypes::shapeToString(getShape()) << endl;)
2097    LAZYDEBUG(cout << "DPPS=" << m_right->getNumDPPSample() <<"."<<endl;)
2098    
2099      double* resultp=&(m_samples[offset]);     // results are stored at the vector offset we recieved
2100      switch(m_op)
2101      {
2102        case PROD:
2103        for (int i=0;i<steps;++i,resultp+=resultStep)
2104        {
2105              const double *ptr_0 = &((*left)[lroffset]);
2106              const double *ptr_1 = &((*right)[rroffset]);
2107    
2108    LAZYDEBUG(cout << DataTypes::pointToString(*left, m_left->getShape(),lroffset,"LEFT") << endl;)
2109    LAZYDEBUG(cout << DataTypes::pointToString(*right,m_right->getShape(),rroffset, "RIGHT") << endl;)
2110    
2111              matrix_matrix_product(m_SL, m_SM, m_SR, ptr_0, ptr_1, resultp, m_transpose);
2112    
2113          lroffset+=leftStep;
2114          rroffset+=rightStep;
2115        }
2116        break;
2117        default:
2118        throw DataException("Programmer error - resolveTProduct can not resolve operator "+opToString(m_op)+".");
2119      }
2120      roffset=offset;
2121      return &m_samples;
2122    }
2123    #endif //LAZY_NODE_STORAGE
2124    
2125  /*  /*
2126    \brief Compute the value of the expression for the given sample.    \brief Compute the value of the expression for the given sample.
# Line 1619  DataLazy::resolveToIdentity() Line 2197  DataLazy::resolveToIdentity()
2197  {  {
2198     if (m_op==IDENTITY)     if (m_op==IDENTITY)
2199      return;      return;
2200    #ifndef LAZY_NODE_STORAGE
2201     DataReady_ptr p=resolveVectorWorker();     DataReady_ptr p=resolveVectorWorker();
2202    #else
2203       DataReady_ptr p=resolveNodeWorker();
2204    #endif
2205     makeIdentity(p);     makeIdentity(p);
2206  }  }
2207    
# Line 1650  DataLazy::resolve() Line 2232  DataLazy::resolve()
2232      return m_id;      return m_id;
2233  }  }
2234    
2235    #ifdef LAZY_NODE_STORAGE
2236    
2237    // This version of resolve uses storage in each node to hold results
2238    DataReady_ptr
2239    DataLazy::resolveNodeWorker()
2240    {
2241      if (m_readytype!='E')     // if the whole sub-expression is Constant or Tagged, then evaluate it normally
2242      {
2243        collapse();
2244      }
2245      if (m_op==IDENTITY)       // So a lazy expression of Constant or Tagged data will be returned here.
2246      {
2247        return m_id;
2248      }
2249        // from this point on we must have m_op!=IDENTITY and m_readytype=='E'
2250      DataExpanded* result=new DataExpanded(getFunctionSpace(),getShape(),  ValueType(getNoValues()));
2251      ValueType& resvec=result->getVectorRW();
2252      DataReady_ptr resptr=DataReady_ptr(result);
2253    
2254      int sample;
2255      int totalsamples=getNumSamples();
2256      const ValueType* res=0;   // Storage for answer
2257    LAZYDEBUG(cout << "Total number of samples=" <<totalsamples << endl;)
2258      #pragma omp parallel for private(sample,res) schedule(static)
2259      for (sample=0;sample<totalsamples;++sample)
2260      {
2261        size_t roffset=0;
2262    #ifdef _OPENMP
2263        res=resolveNodeSample(omp_get_thread_num(),sample,roffset);
2264    #else
2265        res=resolveNodeSample(0,sample,roffset);
2266    #endif
2267    LAZYDEBUG(cout << "Sample #" << sample << endl;)
2268    LAZYDEBUG(cout << "Final res[" << roffset<< "]=" << (*res)[roffset] << (*res)[roffset]<< endl; )
2269        DataVector::size_type outoffset=result->getPointOffset(sample,0);
2270        memcpy(&(resvec[outoffset]),&((*res)[roffset]),m_samplesize*sizeof(DataVector::ElementType));
2271      }
2272      return resptr;
2273    }
2274    
2275    #endif // LAZY_NODE_STORAGE
2276    
2277  // To simplify the memory management, all threads operate on one large vector, rather than one each.  // To simplify the memory management, all threads operate on one large vector, rather than one each.
2278  // Each sample is evaluated independently and copied into the result DataExpanded.  // Each sample is evaluated independently and copied into the result DataExpanded.
2279  DataReady_ptr  DataReady_ptr
# Line 1687  LAZYDEBUG(cout << "Total number of sampl Line 2311  LAZYDEBUG(cout << "Total number of sampl
2311    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)    #pragma omp parallel for private(sample,resoffset,outoffset,res) schedule(static)
2312    for (sample=0;sample<totalsamples;++sample)    for (sample=0;sample<totalsamples;++sample)
2313    {    {
       if (sample==0)  {ENABLEDEBUG}  
   
2314  LAZYDEBUG(cout << "################################# " << sample << endl;)  LAZYDEBUG(cout << "################################# " << sample << endl;)
2315  #ifdef _OPENMP  #ifdef _OPENMP
2316      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);      res=resolveSample(v,threadbuffersize*omp_get_thread_num(),sample,resoffset);
# Line 1706  LAZYDEBUG(cerr << "outoffset=" << outoff Line 2328  LAZYDEBUG(cerr << "outoffset=" << outoff
2328      }      }
2329  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)  LAZYDEBUG(cerr << DataTypes::pointToString(resvec,getShape(),outoffset-m_samplesize+DataTypes::noValues(getShape()),"Final result:") << endl;)
2330  LAZYDEBUG(cerr << "*********************************" << endl;)  LAZYDEBUG(cerr << "*********************************" << endl;)
     DISABLEDEBUG  
2331    }    }
2332    return resptr;    return resptr;
2333  }  }
# Line 1884  DataLazy::setToZero() Line 2505  DataLazy::setToZero()
2505  //   m_buffsRequired=1;  //   m_buffsRequired=1;
2506    
2507    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");    throw DataException("Programmer error - setToZero not supported for DataLazy (DataLazy objects should be read only).");
2508      (int)privdebug;   // to stop the compiler complaining about unused privdebug
2509  }  }
2510    
2511  bool  bool

Legend:
Removed from v.2497  
changed lines
  Added in v.2500

  ViewVC Help
Powered by ViewVC 1.1.26