/[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 790 by bcumming, Wed Jul 26 23:12:34 2006 UTC revision 876 by ksteube, Thu Oct 19 03:50:23 2006 UTC
# Line 422  Data::tag() Line 422  Data::tag()
422    }    }
423  }  }
424    
425  void  Data
426  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
427  {  {
428    m_data->reshapeDataPoint(shape);  #if defined DOPROF
429      profData->where++;
430    #endif
431      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
432  }  }
433    
434  Data  Data
# Line 1029  Data::integrate() const Line 1032  Data::integrate() const
1032  }  }
1033    
1034  Data  Data
1035    Data::erf() const
1036    {
1037    #if defined DOPROF
1038      profData->unary++;
1039    #endif
1040      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1041    }
1042    
1043    Data
1044  Data::sin() const  Data::sin() const
1045  {  {
1046  #if defined DOPROF  #if defined DOPROF
# Line 1316  Data::minval() const Line 1328  Data::minval() const
1328  }  }
1329    
1330  Data  Data
1331  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1332  {  {
1333  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1334    profData->reduction2++;       #if defined DOPROF
1335  #endif       profData->unary++;
1336    Trace trace_func;       #endif
1337    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1338         DataArrayView::ShapeType ev_shape;
1339         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1340         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1341         int rank=getDataPointRank();
1342         if (rank<2) {
1343            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1344         }
1345         if (axis0<0 || axis0>rank-1) {
1346            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1347         }
1348         if (axis1<0 || axis1>rank-1) {
1349             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1350         }
1351         if (axis0 == axis1) {
1352             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1353         }
1354         if (axis0 > axis1) {
1355             axis0_tmp=axis1;
1356             axis1_tmp=axis0;
1357         } else {
1358             axis0_tmp=axis0;
1359             axis1_tmp=axis1;
1360         }
1361         for (int i=0; i<rank; i++) {
1362           if (i == axis0_tmp) {
1363              ev_shape.push_back(s[axis1_tmp]);
1364           } else if (i == axis1_tmp) {
1365              ev_shape.push_back(s[axis0_tmp]);
1366           } else {
1367              ev_shape.push_back(s[i]);
1368           }
1369         }
1370         Data ev(0.,ev_shape,getFunctionSpace());
1371         ev.typeMatchRight(*this);
1372         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1373         return ev;
1374    
1375  }  }
1376    
1377  Data  Data
# Line 1388  Data::nonsymmetric() const Line 1437  Data::nonsymmetric() const
1437  }  }
1438    
1439  Data  Data
1440  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1441  {  {
1442       #if defined DOPROF       #if defined DOPROF
1443          profData->unary++;          profData->unary++;
# Line 1398  Data::matrixtrace(int axis_offset) const Line 1447  Data::matrixtrace(int axis_offset) const
1447          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1448          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1449          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1450          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1451          return ev;          return ev;
1452       }       }
1453       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1462  Data::matrixtrace(int axis_offset) const
1462          }          }
1463          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1464          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1465          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1466          return ev;          return ev;
1467       }       }
1468       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1481  Data::matrixtrace(int axis_offset) const
1481      }      }
1482          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1483          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1484      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1485          return ev;          return ev;
1486       }       }
1487       else {       else {
1488          throw DataException("Error - Data::matrixtrace can only be calculated for rank 2, 3 or 4 object.");          throw DataException("Error - Data::trace can only be calculated for rank 2, 3 or 4 object.");
1489       }       }
1490  }  }
1491    
1492  Data  Data
1493  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1494  {  {
1495  #if defined DOPROF       #if defined DOPROF
1496       profData->reduction2++;       profData->unary++;
1497  #endif       #endif
1498       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1499       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1500       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]       // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
# Line 1609  Data::operator+=(const Data& right) Line 1658  Data::operator+=(const Data& right)
1658    if (isProtected()) {    if (isProtected()) {
1659          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1660    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1661    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1662    return (*this);    return (*this);
1663  }  }
# Line 1619  Data::operator+=(const Data& right) Line 1665  Data::operator+=(const Data& right)
1665  Data&  Data&
1666  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1667  {  {
1668    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1669          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,plus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,plus<double>());  
1670    return (*this);    return (*this);
1671  }  }
1672    
# Line 1635  Data::operator-=(const Data& right) Line 1676  Data::operator-=(const Data& right)
1676    if (isProtected()) {    if (isProtected()) {
1677          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1678    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1679    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1680    return (*this);    return (*this);
1681  }  }
# Line 1645  Data::operator-=(const Data& right) Line 1683  Data::operator-=(const Data& right)
1683  Data&  Data&
1684  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1685  {  {
1686    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1687          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1688    return (*this);    return (*this);
1689  }  }
1690    
# Line 1661  Data::operator*=(const Data& right) Line 1694  Data::operator*=(const Data& right)
1694    if (isProtected()) {    if (isProtected()) {
1695          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1696    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1697    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1698    return (*this);    return (*this);
1699  }  }
# Line 1671  Data::operator*=(const Data& right) Line 1701  Data::operator*=(const Data& right)
1701  Data&  Data&
1702  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1703  {  {
1704    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1705          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1706    return (*this);    return (*this);
1707  }  }
1708    
# Line 1687  Data::operator/=(const Data& right) Line 1712  Data::operator/=(const Data& right)
1712    if (isProtected()) {    if (isProtected()) {
1713          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1714    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1715    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1716    return (*this);    return (*this);
1717  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1719  Data::operator/=(const Data& right)
1719  Data&  Data&
1720  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1721  {  {
1722    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1723          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1724    return (*this);    return (*this);
1725  }  }
1726    
1727  Data  Data
1728  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1729  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1730    Data left_d(left,*this);    Data left_d(left,*this);
1731    return left_d.powD(*this);    return left_d.powD(*this);
1732  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1734  Data::rpowO(const boost::python::object&
1734  Data  Data
1735  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1736  {  {
1737  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1738    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1739  }  }
1740    
1741  Data  Data
1742  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1743  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1744    Data result;    Data result;
1745    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1746    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1747         result.binaryOp(*this,escript::rpow);
1748      } else {
1749         result.copy(*this);
1750         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1751      }
1752    return result;    return result;
1753  }  }
1754    
# Line 1753  escript::operator+(const Data& left, con Line 1761  escript::operator+(const Data& left, con
1761    Data result;    Data result;
1762    //    //
1763    // perform a deep copy    // perform a deep copy
1764    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1765    result+=right;       result.copy(right);
1766         result+=left;
1767      } else {
1768         result.copy(left);
1769         result+=right;
1770      }
1771    return result;    return result;
1772  }  }
1773    
# Line 1766  escript::operator-(const Data& left, con Line 1779  escript::operator-(const Data& left, con
1779    Data result;    Data result;
1780    //    //
1781    // perform a deep copy    // perform a deep copy
1782    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1783    result-=right;       result=right.neg();
1784         result+=left;
1785      } else {
1786         result.copy(left);
1787         result-=right;
1788      }
1789    return result;    return result;
1790  }  }
1791    
# Line 1779  escript::operator*(const Data& left, con Line 1797  escript::operator*(const Data& left, con
1797    Data result;    Data result;
1798    //    //
1799    // perform a deep copy    // perform a deep copy
1800    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1801    result*=right;       result.copy(right);
1802         result*=left;
1803      } else {
1804         result.copy(left);
1805         result*=right;
1806      }
1807    return result;    return result;
1808  }  }
1809    
# Line 1792  escript::operator/(const Data& left, con Line 1815  escript::operator/(const Data& left, con
1815    Data result;    Data result;
1816    //    //
1817    // perform a deep copy    // perform a deep copy
1818    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1819    result/=right;       result=right.oneOver();
1820         result*=left;
1821      } else {
1822         result.copy(left);
1823         result/=right;
1824      }
1825    return result;    return result;
1826  }  }
1827    
# Line 1802  escript::operator/(const Data& left, con Line 1830  escript::operator/(const Data& left, con
1830  Data  Data
1831  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1832  {  {
1833    //    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;  
1834  }  }
1835    
1836  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1838  escript::operator+(const Data& left, con
1838  Data  Data
1839  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1840  {  {
1841    //    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;  
1842  }  }
1843    
1844  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1846  escript::operator-(const Data& left, con
1846  Data  Data
1847  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1848  {  {
1849    //    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;  
1850  }  }
1851    
1852  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1854  escript::operator*(const Data& left, con
1854  Data  Data
1855  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1856  {  {
1857    //    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;  
1858  }  }
1859    
1860  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1862  escript::operator/(const Data& left, con
1862  Data  Data
1863  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1864  {  {
1865    //    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;  
1866  }  }
1867    
1868  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1870  escript::operator+(const boost::python::
1870  Data  Data
1871  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1872  {  {
1873    //    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;  
1874  }  }
1875    
1876  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1878  escript::operator-(const boost::python::
1878  Data  Data
1879  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1880  {  {
1881    //    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;  
1882  }  }
1883    
1884  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1886  escript::operator*(const boost::python::
1886  Data  Data
1887  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1888  {  {
1889    //    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;  
1890  }  }
1891    
1892  //  //
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2517  ostream& escript::operator<<(ostream& o,
2517    return o;    return o;
2518  }  }
2519    
2520    Data
2521    escript::C_GeneralTensorProduct(Data& arg_0,
2522                         Data& arg_1,
2523                         int axis_offset,
2524                         int transpose)
2525    {
2526      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2527      // SM is the product of the last axis_offset entries in arg_0.getShape().
2528    
2529      #if defined DOPROF
2530        // profData->binary++;
2531      #endif
2532    
2533      // Interpolate if necessary and find an appropriate function space
2534      Data arg_0_Z, arg_1_Z;
2535      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2536        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2537          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2538          arg_1_Z = Data(arg_1);
2539        }
2540        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2541          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2542          arg_0_Z =Data(arg_0);
2543        }
2544        else {
2545          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2546        }
2547      } else {
2548          arg_0_Z = Data(arg_0);
2549          arg_1_Z = Data(arg_1);
2550      }
2551      // Get rank and shape of inputs
2552      int rank0 = arg_0_Z.getDataPointRank();
2553      int rank1 = arg_1_Z.getDataPointRank();
2554      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2555      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2556    
2557      // Prepare for the loops of the product and verify compatibility of shapes
2558      int start0=0, start1=0;
2559      if (transpose == 0)       {}
2560      else if (transpose == 1)  { start0 = axis_offset; }
2561      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2562      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2563    
2564      // Adjust the shapes for transpose
2565      DataArrayView::ShapeType tmpShape0;
2566      DataArrayView::ShapeType tmpShape1;
2567      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2568      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2569    
2570    #if 0
2571      // For debugging: show shape after transpose
2572      char tmp[100];
2573      std::string shapeStr;
2574      shapeStr = "(";
2575      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2576      shapeStr += ")";
2577      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2578      shapeStr = "(";
2579      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2580      shapeStr += ")";
2581      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2582    #endif
2583    
2584      // Prepare for the loops of the product
2585      int SL=1, SM=1, SR=1;
2586      for (int i=0; i<rank0-axis_offset; i++)   {
2587        SL *= tmpShape0[i];
2588      }
2589      for (int i=rank0-axis_offset; i<rank0; i++)   {
2590        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2591          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2592        }
2593        SM *= tmpShape0[i];
2594      }
2595      for (int i=axis_offset; i<rank1; i++)     {
2596        SR *= tmpShape1[i];
2597      }
2598    
2599      // Define the shape of the output
2600      DataArrayView::ShapeType shape2;
2601      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2602      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2603    
2604      // Declare output Data object
2605      Data res;
2606    
2607      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2608        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2609        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2610        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2611        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2612        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2613      }
2614      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2615    
2616        // Prepare the DataConstant input
2617        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2618        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2619    
2620        // Borrow DataTagged input from Data object
2621        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2622        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2623    
2624        // Prepare a DataTagged output 2
2625        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2626        res.tag();
2627        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2628        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2629    
2630        // Prepare offset into DataConstant
2631        int offset_0 = tmp_0->getPointOffset(0,0);
2632        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2633        // Get the views
2634        DataArrayView view_1 = tmp_1->getDefaultValue();
2635        DataArrayView view_2 = tmp_2->getDefaultValue();
2636        // Get the pointers to the actual data
2637        double *ptr_1 = &((view_1.getData())[0]);
2638        double *ptr_2 = &((view_2.getData())[0]);
2639        // Compute an MVP for the default
2640        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2641        // Compute an MVP for each tag
2642        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2643        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2644        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2645          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2646          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2647          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2648          double *ptr_1 = &view_1.getData(0);
2649          double *ptr_2 = &view_2.getData(0);
2650          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2651        }
2652    
2653      }
2654      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2655    
2656        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2657        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2658        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2659        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2660        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2661        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2662        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2663        int sampleNo_1,dataPointNo_1;
2664        int numSamples_1 = arg_1_Z.getNumSamples();
2665        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2666        int offset_0 = tmp_0->getPointOffset(0,0);
2667        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2668        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2669          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2670            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2671            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2672            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2673            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2674            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2675            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2676          }
2677        }
2678    
2679      }
2680      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2681    
2682        // Borrow DataTagged input from Data object
2683        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2684        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2685    
2686        // Prepare the DataConstant input
2687        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2688        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2689    
2690        // Prepare a DataTagged output 2
2691        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2692        res.tag();
2693        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2694        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2695    
2696        // Prepare offset into DataConstant
2697        int offset_1 = tmp_1->getPointOffset(0,0);
2698        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2699        // Get the views
2700        DataArrayView view_0 = tmp_0->getDefaultValue();
2701        DataArrayView view_2 = tmp_2->getDefaultValue();
2702        // Get the pointers to the actual data
2703        double *ptr_0 = &((view_0.getData())[0]);
2704        double *ptr_2 = &((view_2.getData())[0]);
2705        // Compute an MVP for the default
2706        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2707        // Compute an MVP for each tag
2708        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2709        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2710        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2711          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2712          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2713          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2714          double *ptr_0 = &view_0.getData(0);
2715          double *ptr_2 = &view_2.getData(0);
2716          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2717        }
2718    
2719      }
2720      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2721    
2722        // Borrow DataTagged input from Data object
2723        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2724        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2725    
2726        // Borrow DataTagged input from Data object
2727        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2728        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2729    
2730        // Prepare a DataTagged output 2
2731        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2732        res.tag();  // DataTagged output
2733        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2734        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2735    
2736        // Get the views
2737        DataArrayView view_0 = tmp_0->getDefaultValue();
2738        DataArrayView view_1 = tmp_1->getDefaultValue();
2739        DataArrayView view_2 = tmp_2->getDefaultValue();
2740        // Get the pointers to the actual data
2741        double *ptr_0 = &((view_0.getData())[0]);
2742        double *ptr_1 = &((view_1.getData())[0]);
2743        double *ptr_2 = &((view_2.getData())[0]);
2744        // Compute an MVP for the default
2745        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2746        // Merge the tags
2747        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2748        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2749        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2750        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2751          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2752        }
2753        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2754          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2755        }
2756        // Compute an MVP for each tag
2757        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2758        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2759          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2760          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2761          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2762          double *ptr_0 = &view_0.getData(0);
2763          double *ptr_1 = &view_1.getData(0);
2764          double *ptr_2 = &view_2.getData(0);
2765          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2766        }
2767    
2768      }
2769      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2770    
2771        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2772        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2773        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2774        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2775        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2776        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2777        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2778        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2779        int sampleNo_0,dataPointNo_0;
2780        int numSamples_0 = arg_0_Z.getNumSamples();
2781        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2782        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2783        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2784          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2785          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2786          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2787            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2788            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2789            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2790            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2791            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2792          }
2793        }
2794    
2795      }
2796      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2797    
2798        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2799        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2800        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2801        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2802        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2803        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2804        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2805        int sampleNo_0,dataPointNo_0;
2806        int numSamples_0 = arg_0_Z.getNumSamples();
2807        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2808        int offset_1 = tmp_1->getPointOffset(0,0);
2809        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2810        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2811          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2812            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2813            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2814            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2815            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2816            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2817            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2818          }
2819        }
2820    
2821    
2822      }
2823      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2824    
2825        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2826        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2827        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2828        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2829        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2830        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2831        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2832        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2833        int sampleNo_0,dataPointNo_0;
2834        int numSamples_0 = arg_0_Z.getNumSamples();
2835        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2836        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2839          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2840          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2841            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2842            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2843            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2844            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2845            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2846          }
2847        }
2848    
2849      }
2850      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2851    
2852        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2853        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2856        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2858        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2859        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2860        int sampleNo_0,dataPointNo_0;
2861        int numSamples_0 = arg_0_Z.getNumSamples();
2862        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2863        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2864        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2865          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2866            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2867            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2868            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2869            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2870            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2871            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2872            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2873          }
2874        }
2875    
2876      }
2877      else {
2878        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2879      }
2880    
2881      return res;
2882    }
2883    
2884    DataAbstract*
2885    Data::borrowData() const
2886    {
2887      return m_data.get();
2888    }
2889    
2890  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2891    
2892  void  void

Legend:
Removed from v.790  
changed lines
  Added in v.876

  ViewVC Help
Powered by ViewVC 1.1.26