/[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 824 by ksteube, Tue Aug 29 05:22:50 2006 UTC revision 922 by gross, Fri Jan 5 04:23:05 2007 UTC
# Line 350  Data::isTagged() const Line 350  Data::isTagged() const
350    return (temp!=0);    return (temp!=0);
351  }  }
352    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
353  bool  bool
354  Data::isEmpty() const  Data::isEmpty() const
355  {  {
# Line 422  Data::tag() Line 420  Data::tag()
420    }    }
421  }  }
422    
423  void  Data
424  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
425  {  {
426    m_data->reshapeDataPoint(shape);  #if defined DOPROF
427      profData->where++;
428    #endif
429      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
430  }  }
431    
432  Data  Data
# Line 547  Data::getDataPointShape() const Line 548  Data::getDataPointShape() const
548    return getPointDataView().getShape();    return getPointDataView().getShape();
549  }  }
550    
551    
552    
553  void  void
554  Data::fillFromNumArray(const boost::python::numeric::array num_array)  Data::fillFromNumArray(const boost::python::numeric::array num_array)
555  {  {
# Line 639  Data::convertToNumArray() Line 642  Data::convertToNumArray()
642        }        }
643        if (dataPointRank==1) {        if (dataPointRank==1) {
644          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
645            numArray[dataPoint][i]=dataPointView(i);            numArray[make_tuple(dataPoint,i)]=dataPointView(i);
646          }          }
647        }        }
648        if (dataPointRank==2) {        if (dataPointRank==2) {
649          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
650            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
651              numArray[dataPoint][i][j] = dataPointView(i,j);              numArray[make_tuple(dataPoint,i,j)] = dataPointView(i,j);
652            }            }
653          }          }
654        }        }
# Line 653  Data::convertToNumArray() Line 656  Data::convertToNumArray()
656          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
657            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
658              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
659                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                numArray[make_tuple(dataPoint,i,j,k)]=dataPointView(i,j,k);
660              }              }
661            }            }
662          }          }
# Line 663  Data::convertToNumArray() Line 666  Data::convertToNumArray()
666            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
667              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
668                for (int l=0; l<dataPointShape[3]; l++) {                for (int l=0; l<dataPointShape[3]; l++) {
669                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  numArray[make_tuple(dataPoint,i,j,k,l)]=dataPointView(i,j,k,l);
670                }                }
671              }              }
672            }            }
# Line 678  Data::convertToNumArray() Line 681  Data::convertToNumArray()
681    return numArray;    return numArray;
682  }  }
683    
684  const  
685    const
686  boost::python::numeric::array  boost::python::numeric::array
687  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data:: getValueOfDataPoint(int dataPointNo)
688  {  {
689    //    size_t length=0;
690    // Check a valid sample number has been supplied    int i, j, k, l;
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   
691    //    //
692    // determine the rank and shape of each data point    // determine the rank and shape of each data point
693    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 702  Data::convertToNumArrayFromSampleNo(int Line 698  Data::convertToNumArrayFromSampleNo(int
698    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
699    
700    //    //
701    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
702    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
703    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
704    // data points, whilst the first dimension will be equal to the    DataArrayView::ShapeType arrayShape = dataPointShape;
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
705    
706    //    //
707    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
708      if (arrayRank==0) {
709        numArray.resize(1);
710      }
711    if (arrayRank==1) {    if (arrayRank==1) {
712      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
713    }    }
# Line 729  Data::convertToNumArrayFromSampleNo(int Line 720  Data::convertToNumArrayFromSampleNo(int
720    if (arrayRank==4) {    if (arrayRank==4) {
721      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
722    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
723    
724    //    if (getNumDataPointsPerSample()>0) {
725    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
726    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
727    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {         //
728      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         // Check a valid sample number has been supplied
729      if (dataPointRank==0) {         if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
730        numArray[dataPoint]=dataPointView();             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
731      }         }
732      if (dataPointRank==1) {                
733        for (int i=0; i<dataPointShape[0]; i++) {         //
734          numArray[dataPoint][i]=dataPointView(i);         // Check a valid data point number has been supplied
735        }         if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
736      }             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
737      if (dataPointRank==2) {         }
738        for (int i=0; i<dataPointShape[0]; i++) {         // TODO: global error handling
739          for (int j=0; j<dataPointShape[1]; j++) {         // create a view of the data if it is stored locally
740            numArray[dataPoint][i][j] = dataPointView(i,j);         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
741          }          
742        }         switch( dataPointRank ){
743      }              case 0 :
744      if (dataPointRank==3) {                  numArray[0] = dataPointView();
745        for (int i=0; i<dataPointShape[0]; i++) {                  break;
746          for (int j=0; j<dataPointShape[1]; j++) {              case 1 :        
747            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
748              numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                      numArray[i]=dataPointView(i);
749            }                  break;
750          }              case 2 :        
751        }                  for( i=0; i<dataPointShape[0]; i++ )
752      }                      for( j=0; j<dataPointShape[1]; j++)
753      if (dataPointRank==4) {                          numArray[make_tuple(i,j)]=dataPointView(i,j);
754        for (int i=0; i<dataPointShape[0]; i++) {                  break;
755          for (int j=0; j<dataPointShape[1]; j++) {              case 3 :        
756            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
757              for (int l=0; l<dataPointShape[3]; l++) {                      for( j=0; j<dataPointShape[1]; j++ )
758                numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                          for( k=0; k<dataPointShape[2]; k++)
759              }                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
760            }                  break;
761          }              case 4 :
762        }                  for( i=0; i<dataPointShape[0]; i++ )
763      }                      for( j=0; j<dataPointShape[1]; j++ )
764                            for( k=0; k<dataPointShape[2]; k++ )
765                                for( l=0; l<dataPointShape[3]; l++)
766                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
767                    break;
768        }
769    }    }
   
770    //    //
771    // return the loaded array    // return the array
772    return numArray;    return numArray;
 }  
   
 const  
 boost::python::numeric::array  
 Data::convertToNumArrayFromDPNo(int procNo,  
                                 int sampleNo,  
                                                                 int dataPointNo)  
773    
774    }
775    void
776    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array num_array)
777  {  {
778      size_t length=0;    if (isProtected()) {
779      int i, j, k, l, pos;          throw DataException("Error - attempt to update protected Data object.");
780      }
781      //
782      // check rank
783      if (num_array.getrank()<getDataPointRank())
784          throw DataException("Rank of numarray does not match Data object rank");
785    
786    //    //
787    // Check a valid sample number has been supplied    // check shape of num_array
788    if (sampleNo >= getNumSamples()) {    for (int i=0; i<getDataPointRank(); i++) {
789      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
790           throw DataException("Shape of numarray does not match Data object rank");
791      }
792      //
793      // make sure data is expanded:
794      if (!isExpanded()) {
795        expand();
796      }
797      if (getNumDataPointsPerSample()>0) {
798           int sampleNo = dataPointNo/getNumDataPointsPerSample();
799           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
800           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
801      } else {
802           m_data->copyToDataPoint(-1, 0,num_array);
803    }    }
804    }
805    
806    void
807    Data::setValueOfDataPoint(int dataPointNo, const double value)
808    {
809      if (isProtected()) {
810            throw DataException("Error - attempt to update protected Data object.");
811      }
812    //    //
813    // Check a valid data point number has been supplied    // make sure data is expanded:
814    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
815      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
816    }    }
817      if (getNumDataPointsPerSample()>0) {
818           int sampleNo = dataPointNo/getNumDataPointsPerSample();
819           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
820           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
821      } else {
822           m_data->copyToDataPoint(-1, 0,value);
823      }
824    }
825    
826    const
827    boost::python::numeric::array
828    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
829    {
830      size_t length=0;
831      int i, j, k, l, pos;
832    //    //
833    // determine the rank and shape of each data point    // determine the rank and shape of each data point
834    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 862  Data::convertToNumArrayFromDPNo(int proc
862      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
863    }    }
864    
865      // added for the MPI communication    // added for the MPI communication
866      length=1;    length=1;
867      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
868          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
869    
870    //    //
871    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
872    
873      // updated for the MPI case      // updated for the MPI case
874      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
875                 if (getNumDataPointsPerSample()>0) {
876                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
877                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
878                    //
879                    // Check a valid sample number has been supplied
880                    if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
881                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
882                    }
883                  
884                    //
885                    // Check a valid data point number has been supplied
886                    if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
887                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
888                    }
889                    // TODO: global error handling
890          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
891          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
892                    
893          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
894          pos=0;          pos=0;
# Line 878  Data::convertToNumArrayFromDPNo(int proc Line 919  Data::convertToNumArrayFromDPNo(int proc
919                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=dataPointView(i,j,k,l);
920                  break;                  break;
921          }          }
922                }
923      }      }
924  #ifdef PASO_MPI          #ifdef PASO_MPI
925          // broadcast the data to all other processes          // broadcast the data to all other processes
926          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
927  #endif          #endif
928    
929      // unpack the data      // unpack the data
930      switch( dataPointRank ){      switch( dataPointRank ){
931          case 0 :          case 0 :
932              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
933              break;              break;
934          case 1 :                  case 1 :        
935              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
# Line 896  Data::convertToNumArrayFromDPNo(int proc Line 938  Data::convertToNumArrayFromDPNo(int proc
938          case 2 :                  case 2 :        
939              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
940                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
941                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
942              break;              break;
943          case 3 :                  case 3 :        
944              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
945                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
946                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
947                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
948              break;              break;
949          case 4 :          case 4 :
950              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
951                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
952                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
953                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
954                              tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];                                  numArray[make_tuple(i,j,k,l)]=tmpData[i+dataPointShape[0]*(j*+dataPointShape[1]*(k+l*dataPointShape[2]))];
955              break;              break;
956      }      }
957    
958      delete [] tmpData;        delete [] tmpData;  
 /*  
   if (dataPointRank==0) {  
     numArray[0]=dataPointView();  
   }  
   if (dataPointRank==1) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       numArray[i]=dataPointView(i);  
     }  
   }  
   if (dataPointRank==2) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         numArray[i][j] = dataPointView(i,j);  
       }  
     }  
   }  
   if (dataPointRank==3) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           numArray[i][j][k]=dataPointView(i,j,k);  
         }  
       }  
     }  
   }  
   if (dataPointRank==4) {  
     for (int i=0; i<dataPointShape[0]; i++) {  
       for (int j=0; j<dataPointShape[1]; j++) {  
         for (int k=0; k<dataPointShape[2]; k++) {  
           for (int l=0; l<dataPointShape[3]; l++) {  
             numArray[i][j][k][l]=dataPointView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
 */  
   
959    //    //
960    // return the loaded array    // return the loaded array
961    return numArray;    return numArray;
962  }  }
963    
964    
965    
966  boost::python::numeric::array  boost::python::numeric::array
967  Data::integrate() const  Data::integrate() const
968  {  {
# Line 1029  Data::integrate() const Line 1035  Data::integrate() const
1035  }  }
1036    
1037  Data  Data
1038    Data::erf() const
1039    {
1040    #if defined DOPROF
1041      profData->unary++;
1042    #endif
1043      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1044    }
1045    
1046    Data
1047  Data::sin() const  Data::sin() const
1048  {  {
1049  #if defined DOPROF  #if defined DOPROF
# Line 1544  Data::eigenvalues_and_eigenvectors(const Line 1559  Data::eigenvalues_and_eigenvectors(const
1559  }  }
1560    
1561  const boost::python::tuple  const boost::python::tuple
1562  Data::mindp() const  Data::minGlobalDataPoint() const
1563  {  {
1564    // NB: calc_mindp had to be split off from mindp as boost::make_tuple causes an    // NB: calc_minGlobalDataPoint( had to be split off from minGlobalDataPoint( as boost::make_tuple causes an
1565    // abort (for unknown reasons) if there are openmp directives with it in the    // abort (for unknown reasons) if there are openmp directives with it in the
1566    // surrounding function    // surrounding function
1567    
   int SampleNo;  
1568    int DataPointNo;    int DataPointNo;
1569      int ProcNo;    int ProcNo;
1570    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1571    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1572  }  }
1573    
1574  void  void
1575  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1576                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1577  {  {
1578    int i,j;    int i,j;
1579    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1618  Data::calc_mindp(  int& ProcNo, Line 1631  Data::calc_mindp(  int& ProcNo,
1631  #else  #else
1632      ProcNo = 0;      ProcNo = 0;
1633  #endif  #endif
1634    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1635  }  }
1636    
1637  void  void
# Line 1646  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 1656  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 1672  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 1682  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 1698  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 1708  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 1724  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 1734  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 1760  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 1790  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 1803  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 1816  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 1829  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 1839  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 1855  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 1871  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 1887  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 1903  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 1916  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 1929  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 1942  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 2210  Data::getRefValue(int ref, Line 2149  Data::getRefValue(int ref,
2149    if (rank==2) {    if (rank==2) {
2150      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
2151        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2152          value[i][j] = valueView(i,j);          value[make_tuple(i,j)] = valueView(i,j);
2153        }        }
2154      }      }
2155    }    }
# Line 2218  Data::getRefValue(int ref, Line 2157  Data::getRefValue(int ref,
2157      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
2158        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2159          for (int k=0; k < shape[2]; k++) {          for (int k=0; k < shape[2]; k++) {
2160            value[i][j][k] = valueView(i,j,k);            value[make_tuple(i,j,k)] = valueView(i,j,k);
2161          }          }
2162        }        }
2163      }      }
# Line 2228  Data::getRefValue(int ref, Line 2167  Data::getRefValue(int ref,
2167        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2168          for (int k=0; k < shape[2]; k++) {          for (int k=0; k < shape[2]; k++) {
2169            for (int l=0; l < shape[3]; l++) {            for (int l=0; l < shape[3]; l++) {
2170              value[i][j][k][l] = valueView(i,j,k,l);              value[make_tuple(i,j,k,l)] = valueView(i,j,k,l);
2171            }            }
2172          }          }
2173        }        }
# Line 2584  escript::C_GeneralTensorProduct(Data& ar Line 2523  escript::C_GeneralTensorProduct(Data& ar
2523                       int axis_offset,                       int axis_offset,
2524                       int transpose)                       int transpose)
2525  {  {
2526    // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)    // 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().    // SM is the product of the last axis_offset entries in arg_0.getShape().
2528    
2529    #if defined DOPROF    #if defined DOPROF
# Line 2592  escript::C_GeneralTensorProduct(Data& ar Line 2531  escript::C_GeneralTensorProduct(Data& ar
2531    #endif    #endif
2532    
2533    // Interpolate if necessary and find an appropriate function space    // Interpolate if necessary and find an appropriate function space
2534    Data arg_Z;    Data arg_0_Z, arg_1_Z;
2535    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {    if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2536      if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {      if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2537        arg_Z = arg_0.interpolate(arg_1.getFunctionSpace());        arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2538        // arg_0=arg_Z;        arg_1_Z = Data(arg_1);
2539      }      }
2540      else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {      else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2541        arg_1=arg_1.interpolate(arg_0.getFunctionSpace());        arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2542          arg_0_Z =Data(arg_0);
2543      }      }
2544      else {      else {
2545        throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");        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    // Get rank and shape of inputs
2552    int rank0 = arg_Z.getDataPointRank();    int rank0 = arg_0_Z.getDataPointRank();
2553    int rank1 = arg_1.getDataPointRank();    int rank1 = arg_1_Z.getDataPointRank();
2554    DataArrayView::ShapeType shape0 = arg_Z.getDataPointShape();    DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2555    DataArrayView::ShapeType shape1 = arg_1.getDataPointShape();    DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2556    
2557    // Prepare for the loops of the product and verify compatibility of shapes    // Prepare for the loops of the product and verify compatibility of shapes
2558    int start0=0, start1=0;    int start0=0, start1=0;
# Line 2655  escript::C_GeneralTensorProduct(Data& ar Line 2598  escript::C_GeneralTensorProduct(Data& ar
2598    
2599    // Define the shape of the output    // Define the shape of the output
2600    DataArrayView::ShapeType shape2;    DataArrayView::ShapeType shape2;
2601    for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_Z    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    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    // Declare output Data object
2605    Data arg_2;    Data res;
2606    
2607    if      (arg_Z.isConstant()   && arg_1.isConstant()) {    if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2608      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataConstant output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2609      double *ptr_0 = &((arg_Z.getPointDataView().getData())[0]);      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2610      double *ptr_1 = &((arg_1.getPointDataView().getData())[0]);      double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2611      double *ptr_2 = &((arg_2.getPointDataView().getData())[0]);      double *ptr_2 = &((res.getPointDataView().getData())[0]);
2612      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);      matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2613    }    }
2614    else if (arg_Z.isConstant()   && arg_1.isTagged()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2615    
2616      // Prepare the DataConstant input      // Prepare the DataConstant input
2617      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_Z.borrowData());      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2618      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2619    
2620      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2621      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2622      if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }      if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2623    
2624      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2625      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());    // DataTagged output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2626      arg_2.tag();      res.tag();
2627      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2628      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2629    
2630      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2631      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2632      double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);      double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2633      // Get the views      // Get the views
2634      DataArrayView view_1 = tmp_1->getDefaultValue();      DataArrayView view_1 = tmp_1->getDefaultValue();
2635      DataArrayView view_2 = tmp_2->getDefaultValue();      DataArrayView view_2 = tmp_2->getDefaultValue();
# Line 2708  escript::C_GeneralTensorProduct(Data& ar Line 2651  escript::C_GeneralTensorProduct(Data& ar
2651      }      }
2652    
2653    }    }
2654    else if (arg_Z.isConstant()   && arg_1.isExpanded()) {    else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2655    
2656      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2657      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_Z.borrowData());      DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2658      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2659      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2660      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2663      int sampleNo_1,dataPointNo_1;      int sampleNo_1,dataPointNo_1;
2664      int numSamples_1 = arg_1.getNumSamples();      int numSamples_1 = arg_1_Z.getNumSamples();
2665      int numDataPointsPerSample_1 = arg_1.getNumDataPointsPerSample();      int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2666      int offset_0 = tmp_0->getPointOffset(0,0);      int offset_0 = tmp_0->getPointOffset(0,0);
2667      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)      #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2668      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {      for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2669        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {        for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2670          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2671          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);          int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2672          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2673          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2674          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2675          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2676        }        }
2677      }      }
2678    
2679    }    }
2680    else if (arg_Z.isTagged()     && arg_1.isConstant()) {    else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2681    
2682      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2683      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2684      if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }      if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2685    
2686      // Prepare the DataConstant input      // Prepare the DataConstant input
2687      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2688      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2689    
2690      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2691      arg_2 = Data(0.0, shape2, arg_Z.getFunctionSpace());    // DataTagged output      res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2692      arg_2.tag();      res.tag();
2693      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2694      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2695    
2696      // Prepare offset into DataConstant      // Prepare offset into DataConstant
2697      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2698      double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);      double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2699      // Get the views      // Get the views
2700      DataArrayView view_0 = tmp_0->getDefaultValue();      DataArrayView view_0 = tmp_0->getDefaultValue();
2701      DataArrayView view_2 = tmp_2->getDefaultValue();      DataArrayView view_2 = tmp_2->getDefaultValue();
# Line 2774  escript::C_GeneralTensorProduct(Data& ar Line 2717  escript::C_GeneralTensorProduct(Data& ar
2717      }      }
2718    
2719    }    }
2720    else if (arg_Z.isTagged()     && arg_1.isTagged()) {    else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2721    
2722      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2723      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2724      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2725    
2726      // Borrow DataTagged input from Data object      // Borrow DataTagged input from Data object
2727      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2728      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2729    
2730      // Prepare a DataTagged output 2      // Prepare a DataTagged output 2
2731      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace());      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2732      arg_2.tag();    // DataTagged output      res.tag();  // DataTagged output
2733      DataTagged* tmp_2=dynamic_cast<DataTagged*>(arg_2.borrowData());      DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2734      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2735    
2736      // Get the views      // Get the views
# Line 2823  escript::C_GeneralTensorProduct(Data& ar Line 2766  escript::C_GeneralTensorProduct(Data& ar
2766      }      }
2767    
2768    }    }
2769    else if (arg_Z.isTagged()     && arg_1.isExpanded()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2772      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2773      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_Z.borrowData());      DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2774      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2775      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2776      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2779      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2780      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2781      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2782      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2783      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      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        int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2785        double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2786        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2787          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2788          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2789          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2790          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2791          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2792        }        }
2793      }      }
2794    
2795    }    }
2796    else if (arg_Z.isExpanded()   && arg_1.isConstant()) {    else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2797    
2798      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2799      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2800      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1.borrowData());      DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2801      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2802      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2805      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2806      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2807      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2808      int offset_1 = tmp_1->getPointOffset(0,0);      int offset_1 = tmp_1->getPointOffset(0,0);
2809      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2810      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2811        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2812          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2813          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2814          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2815          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2816          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2817          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2818        }        }
2819      }      }
2820    
2821    
2822    }    }
2823    else if (arg_Z.isExpanded()   && arg_1.isTagged()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2826      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2827      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2828      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1.borrowData());      DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2829      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2830      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2833      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2834      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2835      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2836      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2837      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2838        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);        int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2839        double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2840        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2841          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2842          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2843          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2844          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2845          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2846        }        }
2847      }      }
2848    
2849    }    }
2850    else if (arg_Z.isExpanded()   && arg_1.isExpanded()) {    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      // After finding a common function space above the two inputs have the same numSamples and num DPPS
2853      arg_2 = Data(0.0, shape2, arg_1.getFunctionSpace(),true); // DataExpanded output      res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2854      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_Z.borrowData());      DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2855      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1.borrowData());      DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2856      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(arg_2.borrowData());      DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2857      if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }      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."); }      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."); }      if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2860      int sampleNo_0,dataPointNo_0;      int sampleNo_0,dataPointNo_0;
2861      int numSamples_0 = arg_Z.getNumSamples();      int numSamples_0 = arg_0_Z.getNumSamples();
2862      int numDataPointsPerSample_0 = arg_Z.getNumDataPointsPerSample();      int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2863      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)      #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2864      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {      for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2865        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {        for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2866          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2867          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2868          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);          int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2869          double *ptr_0 = &((arg_Z.getPointDataView().getData())[offset_0]);          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2870          double *ptr_1 = &((arg_1.getPointDataView().getData())[offset_1]);          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2871          double *ptr_2 = &((arg_2.getPointDataView().getData())[offset_2]);          double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2872          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2873        }        }
2874      }      }
# Line 2935  escript::C_GeneralTensorProduct(Data& ar Line 2878  escript::C_GeneralTensorProduct(Data& ar
2878      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");      throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2879    }    }
2880    
2881    return arg_2;    return res;
2882  }  }
2883    
2884  DataAbstract*  DataAbstract*

Legend:
Removed from v.824  
changed lines
  Added in v.922

  ViewVC Help
Powered by ViewVC 1.1.26