/[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 1034 by gross, Wed Mar 14 23:49:20 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    
 void  
 Data::fillFromNumArray(const boost::python::numeric::array num_array)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
   //  
   // check rank  
   if (num_array.getrank()<getDataPointRank())  
       throw DataException("Rank of numarray does not match Data object rank");  
   
   //  
   // check shape of num_array  
   for (int i=0; i<getDataPointRank(); i++) {  
     if (extract<int>(num_array.getshape()[i+1])!=getDataPointShape()[i])  
        throw DataException("Shape of numarray does not match Data object rank");  
   }  
   
   //  
   // make sure data is expanded:  
   if (!isExpanded()) {  
     expand();  
   }  
551    
   //  
   // and copy over  
   m_data->copyAll(num_array);  
 }  
552    
553  const  const
554  boost::python::numeric::array  boost::python::numeric::array
555  Data::convertToNumArray()  Data:: getValueOfDataPoint(int dataPointNo)
556  {  {
557    //    size_t length=0;
558    // determine the total number of data points    int i, j, k, l;
   int numSamples = getNumSamples();  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   int numDataPoints = numSamples * numDataPointsPerSample;  
   
559    //    //
560    // determine the rank and shape of each data point    // determine the rank and shape of each data point
561    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 596  Data::convertToNumArray() Line 566  Data::convertToNumArray()
566    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
567    
568    //    //
569    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
570    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
571    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
572    // 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(numDataPoints);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
573    
574    //    //
575    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
576      if (arrayRank==0) {
577        numArray.resize(1);
578      }
579    if (arrayRank==1) {    if (arrayRank==1) {
580      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
581    }    }
# Line 623  Data::convertToNumArray() Line 588  Data::convertToNumArray()
588    if (arrayRank==4) {    if (arrayRank==4) {
589      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
590    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
591    
592    //    if (getNumDataPointsPerSample()>0) {
593    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
594    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
595    int dataPoint = 0;         //
596    for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {         // Check a valid sample number has been supplied
597      for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {         if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
598        DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
599        if (dataPointRank==0) {         }
600          numArray[dataPoint]=dataPointView();                
601        }         //
602        if (dataPointRank==1) {         // Check a valid data point number has been supplied
603          for (int i=0; i<dataPointShape[0]; i++) {         if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
604            numArray[dataPoint][i]=dataPointView(i);             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
605          }         }
606        }         // TODO: global error handling
607        if (dataPointRank==2) {         // create a view of the data if it is stored locally
608          for (int i=0; i<dataPointShape[0]; i++) {         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
609            for (int j=0; j<dataPointShape[1]; j++) {          
610              numArray[dataPoint][i][j] = dataPointView(i,j);         switch( dataPointRank ){
611            }              case 0 :
612          }                  numArray[0] = dataPointView();
613        }                  break;
614        if (dataPointRank==3) {              case 1 :        
615          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
616            for (int j=0; j<dataPointShape[1]; j++) {                      numArray[i]=dataPointView(i);
617              for (int k=0; k<dataPointShape[2]; k++) {                  break;
618                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);              case 2 :        
619              }                  for( i=0; i<dataPointShape[0]; i++ )
620            }                      for( j=0; j<dataPointShape[1]; j++)
621          }                          numArray[make_tuple(i,j)]=dataPointView(i,j);
622        }                  break;
623        if (dataPointRank==4) {              case 3 :        
624          for (int i=0; i<dataPointShape[0]; i++) {                  for( i=0; i<dataPointShape[0]; i++ )
625            for (int j=0; j<dataPointShape[1]; j++) {                      for( j=0; j<dataPointShape[1]; j++ )
626              for (int k=0; k<dataPointShape[2]; k++) {                          for( k=0; k<dataPointShape[2]; k++)
627                for (int l=0; l<dataPointShape[3]; l++) {                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
628                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  break;
629                }              case 4 :
630              }                  for( i=0; i<dataPointShape[0]; i++ )
631            }                      for( j=0; j<dataPointShape[1]; j++ )
632          }                          for( k=0; k<dataPointShape[2]; k++ )
633        }                              for( l=0; l<dataPointShape[3]; l++)
634        dataPoint++;                                  numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
635      }                  break;
636        }
637    }    }
   
638    //    //
639    // return the loaded array    // return the array
640    return numArray;    return numArray;
 }  
641    
642  const  }
643  boost::python::numeric::array  void
644  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data::setValueOfDataPointToPyObject(int dataPointNo, const boost::python::object& py_object)
645  {  {
646    //      // this will throw if the value cannot be represented
647    // Check a valid sample number has been supplied      boost::python::numeric::array num_array(py_object);
648    if (sampleNo >= getNumSamples()) {      setValueOfDataPointToArray(dataPointNo,num_array);
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
649    
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
650    
651    //  }
   // determine the rank and shape of each data point  
   int dataPointRank = getDataPointRank();  
   DataArrayView::ShapeType dataPointShape = getDataPointShape();  
652    
653    void
654    Data::setValueOfDataPointToArray(int dataPointNo, const boost::python::numeric::array& num_array)
655    {
656      if (isProtected()) {
657            throw DataException("Error - attempt to update protected Data object.");
658      }
659    //    //
660    // create the numeric array to be returned    // check rank
661    boost::python::numeric::array numArray(0.0);    if (num_array.getrank()<getDataPointRank())
662          throw DataException("Rank of numarray does not match Data object rank");
663    
664    //    //
665    // the rank of the returned numeric array will be the rank of    // check shape of num_array
666    // the data points, plus one. Where the rank of the array is n,    for (int i=0; i<getDataPointRank(); i++) {
667    // the last n-1 dimensions will be equal to the shape of the      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
668    // data points, whilst the first dimension will be equal to the         throw DataException("Shape of numarray does not match Data object rank");
   // 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]);  
669    }    }
   
670    //    //
671    // resize the numeric array to the shape just calculated    // make sure data is expanded:
672    if (arrayRank==1) {    if (!isExpanded()) {
673      numArray.resize(arrayShape[0]);      expand();
   }  
   if (arrayRank==2) {  
     numArray.resize(arrayShape[0],arrayShape[1]);  
   }  
   if (arrayRank==3) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2]);  
   }  
   if (arrayRank==4) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);  
   }  
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
674    }    }
675      if (getNumDataPointsPerSample()>0) {
676    //         int sampleNo = dataPointNo/getNumDataPointsPerSample();
677    // loop through each data point in turn, loading the values for that data point         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
678    // into the numeric array.         m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
679    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {    } else {
680      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         m_data->copyToDataPoint(-1, 0,num_array);
     if (dataPointRank==0) {  
       numArray[dataPoint]=dataPointView();  
     }  
     if (dataPointRank==1) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         numArray[dataPoint][i]=dataPointView(i);  
       }  
     }  
     if (dataPointRank==2) {  
       for (int i=0; i<dataPointShape[0]; i++) {  
         for (int j=0; j<dataPointShape[1]; j++) {  
           numArray[dataPoint][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[dataPoint][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[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);  
             }  
           }  
         }  
       }  
     }  
681    }    }
   
   //  
   // return the loaded array  
   return numArray;  
682  }  }
683    
684  const  void
685  boost::python::numeric::array  Data::setValueOfDataPoint(int dataPointNo, const double value)
 Data::convertToNumArrayFromDPNo(int procNo,  
                                 int sampleNo,  
                                                                 int dataPointNo)  
   
686  {  {
687      size_t length=0;    if (isProtected()) {
688      int i, j, k, l, pos;          throw DataException("Error - attempt to update protected Data object.");
   
   //  
   // Check a valid sample number has been supplied  
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
689    }    }
   
690    //    //
691    // Check a valid data point number has been supplied    // make sure data is expanded:
692    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
693      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
694      }
695      if (getNumDataPointsPerSample()>0) {
696           int sampleNo = dataPointNo/getNumDataPointsPerSample();
697           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
698           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,value);
699      } else {
700           m_data->copyToDataPoint(-1, 0,value);
701    }    }
702    }
703    
704    const
705    boost::python::numeric::array
706    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
707    {
708      size_t length=0;
709      int i, j, k, l, pos;
710    //    //
711    // determine the rank and shape of each data point    // determine the rank and shape of each data point
712    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 740  Data::convertToNumArrayFromDPNo(int proc
740      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
741    }    }
742    
743      // added for the MPI communication    // added for the MPI communication
744      length=1;    length=1;
745      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
746          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
747    
748    //    //
749    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
750    
751      // updated for the MPI case      // updated for the MPI case
752      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
753                 if (getNumDataPointsPerSample()>0) {
754                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
755                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
756                    //
757                    // Check a valid sample number has been supplied
758                    if ((sampleNo >= getNumSamples()) || (sampleNo < 0 )) {
759                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
760                    }
761                  
762                    //
763                    // Check a valid data point number has been supplied
764                    if ((dataPointNoInSample >= getNumDataPointsPerSample()) || (dataPointNoInSample < 0)) {
765                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
766                    }
767                    // TODO: global error handling
768          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
769          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
770                    
771          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
772          pos=0;          pos=0;
# Line 878  Data::convertToNumArrayFromDPNo(int proc Line 797  Data::convertToNumArrayFromDPNo(int proc
797                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=dataPointView(i,j,k,l);
798                  break;                  break;
799          }          }
800                }
801      }      }
802  #ifdef PASO_MPI          #ifdef PASO_MPI
803          // broadcast the data to all other processes          // broadcast the data to all other processes
804          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
805  #endif          #endif
806    
807      // unpack the data      // unpack the data
808      switch( dataPointRank ){      switch( dataPointRank ){
809          case 0 :          case 0 :
810              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
811              break;              break;
812          case 1 :                  case 1 :        
813              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
# Line 896  Data::convertToNumArrayFromDPNo(int proc Line 816  Data::convertToNumArrayFromDPNo(int proc
816          case 2 :                  case 2 :        
817              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
818                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
819                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
820              break;              break;
821          case 3 :                  case 3 :        
822              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
823                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
824                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
825                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
826              break;              break;
827          case 4 :          case 4 :
828              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
829                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
830                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
831                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
832                              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]))];
833              break;              break;
834      }      }
835    
836      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);  
           }  
         }  
       }  
     }  
   }  
 */  
   
837    //    //
838    // return the loaded array    // return the loaded array
839    return numArray;    return numArray;
840  }  }
841    
842    
843    
844  boost::python::numeric::array  boost::python::numeric::array
845  Data::integrate() const  Data::integrate() const
846  {  {
# Line 1073  Data::acos() const Line 957  Data::acos() const
957    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
958  }  }
959    
960    
961  Data  Data
962  Data::atan() const  Data::atan() const
963  {  {
# Line 1109  Data::tanh() const Line 994  Data::tanh() const
994    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
995  }  }
996    
997    
998    Data
999    Data::erf() const
1000    {
1001    #if defined DOPROF
1002      profData->unary++;
1003    #endif
1004    #ifdef _WIN32
1005      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1006    #else
1007      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1008    #endif
1009    }
1010    
1011  Data  Data
1012  Data::asinh() const  Data::asinh() const
1013  {  {
1014  #if defined DOPROF  #if defined DOPROF
1015    profData->unary++;    profData->unary++;
1016  #endif  #endif
1017    #ifdef _WIN32
1018      return escript::unaryOp(*this,escript::asinh_substitute);
1019    #else
1020    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1021    #endif
1022  }  }
1023    
1024  Data  Data
# Line 1124  Data::acosh() const Line 1027  Data::acosh() const
1027  #if defined DOPROF  #if defined DOPROF
1028    profData->unary++;    profData->unary++;
1029  #endif  #endif
1030    #ifdef _WIN32
1031      return escript::unaryOp(*this,escript::acosh_substitute);
1032    #else
1033    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1034    #endif
1035  }  }
1036    
1037  Data  Data
# Line 1133  Data::atanh() const Line 1040  Data::atanh() const
1040  #if defined DOPROF  #if defined DOPROF
1041    profData->unary++;    profData->unary++;
1042  #endif  #endif
1043    #ifdef _WIN32
1044      return escript::unaryOp(*this,escript::atanh_substitute);
1045    #else
1046    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1047    #endif
1048  }  }
1049    
1050  Data  Data
# Line 1316  Data::minval() const Line 1227  Data::minval() const
1227  }  }
1228    
1229  Data  Data
1230  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1231  {  {
1232  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1233    profData->reduction2++;       #if defined DOPROF
1234  #endif       profData->unary++;
1235    Trace trace_func;       #endif
1236    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1237         DataArrayView::ShapeType ev_shape;
1238         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1239         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1240         int rank=getDataPointRank();
1241         if (rank<2) {
1242            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1243         }
1244         if (axis0<0 || axis0>rank-1) {
1245            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1246         }
1247         if (axis1<0 || axis1>rank-1) {
1248             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1249         }
1250         if (axis0 == axis1) {
1251             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1252         }
1253         if (axis0 > axis1) {
1254             axis0_tmp=axis1;
1255             axis1_tmp=axis0;
1256         } else {
1257             axis0_tmp=axis0;
1258             axis1_tmp=axis1;
1259         }
1260         for (int i=0; i<rank; i++) {
1261           if (i == axis0_tmp) {
1262              ev_shape.push_back(s[axis1_tmp]);
1263           } else if (i == axis1_tmp) {
1264              ev_shape.push_back(s[axis0_tmp]);
1265           } else {
1266              ev_shape.push_back(s[i]);
1267           }
1268         }
1269         Data ev(0.,ev_shape,getFunctionSpace());
1270         ev.typeMatchRight(*this);
1271         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1272         return ev;
1273    
1274  }  }
1275    
1276  Data  Data
# Line 1388  Data::nonsymmetric() const Line 1336  Data::nonsymmetric() const
1336  }  }
1337    
1338  Data  Data
1339  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1340  {  {
1341       #if defined DOPROF       #if defined DOPROF
1342          profData->unary++;          profData->unary++;
# Line 1398  Data::matrixtrace(int axis_offset) const Line 1346  Data::matrixtrace(int axis_offset) const
1346          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1347          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1348          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1349          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1350          return ev;          return ev;
1351       }       }
1352       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1361  Data::matrixtrace(int axis_offset) const
1361          }          }
1362          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1363          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1364          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1365          return ev;          return ev;
1366       }       }
1367       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1380  Data::matrixtrace(int axis_offset) const
1380      }      }
1381          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1382          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1383      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1384          return ev;          return ev;
1385       }       }
1386       else {       else {
1387          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.");
1388       }       }
1389  }  }
1390    
1391  Data  Data
1392  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1393  {  {
1394  #if defined DOPROF       #if defined DOPROF
1395       profData->reduction2++;       profData->unary++;
1396  #endif       #endif
1397       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1398       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1399       // 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 1507  Data::eigenvalues_and_eigenvectors(const Line 1455  Data::eigenvalues_and_eigenvectors(const
1455  }  }
1456    
1457  const boost::python::tuple  const boost::python::tuple
1458  Data::mindp() const  Data::minGlobalDataPoint() const
1459  {  {
1460    // 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
1461    // 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
1462    // surrounding function    // surrounding function
1463    
   int SampleNo;  
1464    int DataPointNo;    int DataPointNo;
1465      int ProcNo;    int ProcNo;
1466    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1467    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1468  }  }
1469    
1470  void  void
1471  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1472                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1473  {  {
1474    int i,j;    int i,j;
1475    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1527  Data::calc_mindp(  int& ProcNo,
1527  #else  #else
1528      ProcNo = 0;      ProcNo = 0;
1529  #endif  #endif
1530    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1531  }  }
1532    
1533  void  void
# Line 1609  Data::operator+=(const Data& right) Line 1554  Data::operator+=(const Data& right)
1554    if (isProtected()) {    if (isProtected()) {
1555          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1556    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1557    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1558    return (*this);    return (*this);
1559  }  }
# Line 1619  Data::operator+=(const Data& right) Line 1561  Data::operator+=(const Data& right)
1561  Data&  Data&
1562  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1563  {  {
1564    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1565          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,plus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,plus<double>());  
1566    return (*this);    return (*this);
1567  }  }
1568    
# Line 1635  Data::operator-=(const Data& right) Line 1572  Data::operator-=(const Data& right)
1572    if (isProtected()) {    if (isProtected()) {
1573          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1574    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1575    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1576    return (*this);    return (*this);
1577  }  }
# Line 1645  Data::operator-=(const Data& right) Line 1579  Data::operator-=(const Data& right)
1579  Data&  Data&
1580  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1581  {  {
1582    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1583          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1584    return (*this);    return (*this);
1585  }  }
1586    
# Line 1661  Data::operator*=(const Data& right) Line 1590  Data::operator*=(const Data& right)
1590    if (isProtected()) {    if (isProtected()) {
1591          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1592    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1593    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1594    return (*this);    return (*this);
1595  }  }
# Line 1671  Data::operator*=(const Data& right) Line 1597  Data::operator*=(const Data& right)
1597  Data&  Data&
1598  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1599  {  {
1600    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1601          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1602    return (*this);    return (*this);
1603  }  }
1604    
# Line 1687  Data::operator/=(const Data& right) Line 1608  Data::operator/=(const Data& right)
1608    if (isProtected()) {    if (isProtected()) {
1609          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1610    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1611    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1612    return (*this);    return (*this);
1613  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1615  Data::operator/=(const Data& right)
1615  Data&  Data&
1616  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1617  {  {
1618    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1619          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1620    return (*this);    return (*this);
1621  }  }
1622    
1623  Data  Data
1624  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1625  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1626    Data left_d(left,*this);    Data left_d(left,*this);
1627    return left_d.powD(*this);    return left_d.powD(*this);
1628  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1630  Data::rpowO(const boost::python::object&
1630  Data  Data
1631  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1632  {  {
1633  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1634    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1635  }  }
1636    
1637  Data  Data
1638  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1639  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1640    Data result;    Data result;
1641    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1642    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1643         result.binaryOp(*this,escript::rpow);
1644      } else {
1645         result.copy(*this);
1646         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1647      }
1648    return result;    return result;
1649  }  }
1650    
# Line 1753  escript::operator+(const Data& left, con Line 1657  escript::operator+(const Data& left, con
1657    Data result;    Data result;
1658    //    //
1659    // perform a deep copy    // perform a deep copy
1660    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1661    result+=right;       result.copy(right);
1662         result+=left;
1663      } else {
1664         result.copy(left);
1665         result+=right;
1666      }
1667    return result;    return result;
1668  }  }
1669    
# Line 1766  escript::operator-(const Data& left, con Line 1675  escript::operator-(const Data& left, con
1675    Data result;    Data result;
1676    //    //
1677    // perform a deep copy    // perform a deep copy
1678    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1679    result-=right;       result=right.neg();
1680         result+=left;
1681      } else {
1682         result.copy(left);
1683         result-=right;
1684      }
1685    return result;    return result;
1686  }  }
1687    
# Line 1779  escript::operator*(const Data& left, con Line 1693  escript::operator*(const Data& left, con
1693    Data result;    Data result;
1694    //    //
1695    // perform a deep copy    // perform a deep copy
1696    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1697    result*=right;       result.copy(right);
1698         result*=left;
1699      } else {
1700         result.copy(left);
1701         result*=right;
1702      }
1703    return result;    return result;
1704  }  }
1705    
# Line 1792  escript::operator/(const Data& left, con Line 1711  escript::operator/(const Data& left, con
1711    Data result;    Data result;
1712    //    //
1713    // perform a deep copy    // perform a deep copy
1714    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1715    result/=right;       result=right.oneOver();
1716         result*=left;
1717      } else {
1718         result.copy(left);
1719         result/=right;
1720      }
1721    return result;    return result;
1722  }  }
1723    
# Line 1802  escript::operator/(const Data& left, con Line 1726  escript::operator/(const Data& left, con
1726  Data  Data
1727  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1728  {  {
1729    //    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;  
1730  }  }
1731    
1732  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1734  escript::operator+(const Data& left, con
1734  Data  Data
1735  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1736  {  {
1737    //    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;  
1738  }  }
1739    
1740  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1742  escript::operator-(const Data& left, con
1742  Data  Data
1743  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1744  {  {
1745    //    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;  
1746  }  }
1747    
1748  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1750  escript::operator*(const Data& left, con
1750  Data  Data
1751  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1752  {  {
1753    //    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;  
1754  }  }
1755    
1756  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1758  escript::operator/(const Data& left, con
1758  Data  Data
1759  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1760  {  {
1761    //    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;  
1762  }  }
1763    
1764  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1766  escript::operator+(const boost::python::
1766  Data  Data
1767  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1768  {  {
1769    //    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;  
1770  }  }
1771    
1772  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1774  escript::operator-(const boost::python::
1774  Data  Data
1775  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1776  {  {
1777    //    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;  
1778  }  }
1779    
1780  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1782  escript::operator*(const boost::python::
1782  Data  Data
1783  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1784  {  {
1785    //    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;  
1786  }  }
1787    
1788  //  //
# Line 2116  Data::getTagNumber(int dpno) Line 1988  Data::getTagNumber(int dpno)
1988    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
1989  }  }
1990    
 /* TODO */  
 /* global reduction */  
 void  
 Data::setRefValue(int ref,  
                   const boost::python::numeric::array& value)  
 {  
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
   //  
   // Construct DataArray from boost::python::object input value  
   DataArray valueDataArray(value);  
   
   //  
   // Call DataAbstract::setRefValue  
   m_data->setRefValue(ref,valueDataArray);  
 }  
   
 /* TODO */  
 /* global reduction */  
 void  
 Data::getRefValue(int ref,  
                   boost::python::numeric::array& value)  
 {  
   //  
   // Construct DataArray for boost::python::object return value  
   DataArray valueDataArray(value);  
   
   //  
   // Load DataArray with values from data-points specified by ref  
   m_data->getRefValue(ref,valueDataArray);  
   
   //  
   // Load values from valueDataArray into return numarray  
   
   // extract the shape of the numarray  
   int rank = value.getrank();  
   DataArrayView::ShapeType shape;  
   for (int i=0; i < rank; i++) {  
     shape.push_back(extract<int>(value.getshape()[i]));  
   }  
   
   // and load the numarray with the data from the DataArray  
   DataArrayView valueView = valueDataArray.getView();  
   
   if (rank==0) {  
       boost::python::numeric::array temp_numArray(valueView());  
       value = temp_numArray;  
   }  
   if (rank==1) {  
     for (int i=0; i < shape[0]; i++) {  
       value[i] = valueView(i);  
     }  
   }  
   if (rank==2) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         value[i][j] = valueView(i,j);  
       }  
     }  
   }  
   if (rank==3) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           value[i][j][k] = valueView(i,j,k);  
         }  
       }  
     }  
   }  
   if (rank==4) {  
     for (int i=0; i < shape[0]; i++) {  
       for (int j=0; j < shape[1]; j++) {  
         for (int k=0; k < shape[2]; k++) {  
           for (int l=0; l < shape[3]; l++) {  
             value[i][j][k][l] = valueView(i,j,k,l);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
1991  void  void
1992  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
1993  {  {
# Line 2241  Data::archiveData(const std::string file Line 2029  Data::archiveData(const std::string file
2029    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2030    vector<int> referenceNumbers(noSamples);    vector<int> referenceNumbers(noSamples);
2031    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2032      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
2033    }    }
2034    vector<int> tagNumbers(noSamples);    vector<int> tagNumbers(noSamples);
2035    if (isTagged()) {    if (isTagged()) {
# Line 2450  Data::extractData(const std::string file Line 2238  Data::extractData(const std::string file
2238      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2239    }    }
2240    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2241      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2242        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2243      }      }
2244    }    }
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2329  ostream& escript::operator<<(ostream& o,
2329    return o;    return o;
2330  }  }
2331    
2332    Data
2333    escript::C_GeneralTensorProduct(Data& arg_0,
2334                         Data& arg_1,
2335                         int axis_offset,
2336                         int transpose)
2337    {
2338      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2339      // SM is the product of the last axis_offset entries in arg_0.getShape().
2340    
2341      #if defined DOPROF
2342        // profData->binary++;
2343      #endif
2344    
2345      // Interpolate if necessary and find an appropriate function space
2346      Data arg_0_Z, arg_1_Z;
2347      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2348        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2349          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2350          arg_1_Z = Data(arg_1);
2351        }
2352        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2353          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2354          arg_0_Z =Data(arg_0);
2355        }
2356        else {
2357          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2358        }
2359      } else {
2360          arg_0_Z = Data(arg_0);
2361          arg_1_Z = Data(arg_1);
2362      }
2363      // Get rank and shape of inputs
2364      int rank0 = arg_0_Z.getDataPointRank();
2365      int rank1 = arg_1_Z.getDataPointRank();
2366      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2367      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2368    
2369      // Prepare for the loops of the product and verify compatibility of shapes
2370      int start0=0, start1=0;
2371      if (transpose == 0)       {}
2372      else if (transpose == 1)  { start0 = axis_offset; }
2373      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2374      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2375    
2376      // Adjust the shapes for transpose
2377      DataArrayView::ShapeType tmpShape0;
2378      DataArrayView::ShapeType tmpShape1;
2379      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2380      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2381    
2382    #if 0
2383      // For debugging: show shape after transpose
2384      char tmp[100];
2385      std::string shapeStr;
2386      shapeStr = "(";
2387      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2388      shapeStr += ")";
2389      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2390      shapeStr = "(";
2391      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2392      shapeStr += ")";
2393      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2394    #endif
2395    
2396      // Prepare for the loops of the product
2397      int SL=1, SM=1, SR=1;
2398      for (int i=0; i<rank0-axis_offset; i++)   {
2399        SL *= tmpShape0[i];
2400      }
2401      for (int i=rank0-axis_offset; i<rank0; i++)   {
2402        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2403          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2404        }
2405        SM *= tmpShape0[i];
2406      }
2407      for (int i=axis_offset; i<rank1; i++)     {
2408        SR *= tmpShape1[i];
2409      }
2410    
2411      // Define the shape of the output
2412      DataArrayView::ShapeType shape2;
2413      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2414      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2415    
2416      // Declare output Data object
2417      Data res;
2418    
2419      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2420        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2421        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2422        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2423        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2424        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2425      }
2426      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2427    
2428        // Prepare the DataConstant input
2429        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2430        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2431    
2432        // Borrow DataTagged input from Data object
2433        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2434        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2435    
2436        // Prepare a DataTagged output 2
2437        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2438        res.tag();
2439        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2440        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2441    
2442        // Prepare offset into DataConstant
2443        int offset_0 = tmp_0->getPointOffset(0,0);
2444        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2445        // Get the views
2446        DataArrayView view_1 = tmp_1->getDefaultValue();
2447        DataArrayView view_2 = tmp_2->getDefaultValue();
2448        // Get the pointers to the actual data
2449        double *ptr_1 = &((view_1.getData())[0]);
2450        double *ptr_2 = &((view_2.getData())[0]);
2451        // Compute an MVP for the default
2452        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2453        // Compute an MVP for each tag
2454        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2455        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2456        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2457          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2458          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2459          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2460          double *ptr_1 = &view_1.getData(0);
2461          double *ptr_2 = &view_2.getData(0);
2462          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2463        }
2464    
2465      }
2466      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2467    
2468        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2469        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2470        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2471        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2472        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2473        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2474        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2475        int sampleNo_1,dataPointNo_1;
2476        int numSamples_1 = arg_1_Z.getNumSamples();
2477        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2478        int offset_0 = tmp_0->getPointOffset(0,0);
2479        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2480        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2481          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2482            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2483            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2484            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2485            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2486            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2487            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2488          }
2489        }
2490    
2491      }
2492      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2493    
2494        // Borrow DataTagged input from Data object
2495        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2496        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2497    
2498        // Prepare the DataConstant input
2499        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2500        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2501    
2502        // Prepare a DataTagged output 2
2503        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2504        res.tag();
2505        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2506        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2507    
2508        // Prepare offset into DataConstant
2509        int offset_1 = tmp_1->getPointOffset(0,0);
2510        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2511        // Get the views
2512        DataArrayView view_0 = tmp_0->getDefaultValue();
2513        DataArrayView view_2 = tmp_2->getDefaultValue();
2514        // Get the pointers to the actual data
2515        double *ptr_0 = &((view_0.getData())[0]);
2516        double *ptr_2 = &((view_2.getData())[0]);
2517        // Compute an MVP for the default
2518        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2519        // Compute an MVP for each tag
2520        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2521        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2522        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2523          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2524          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2525          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2526          double *ptr_0 = &view_0.getData(0);
2527          double *ptr_2 = &view_2.getData(0);
2528          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2529        }
2530    
2531      }
2532      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2533    
2534        // Borrow DataTagged input from Data object
2535        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2536        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2537    
2538        // Borrow DataTagged input from Data object
2539        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2540        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2541    
2542        // Prepare a DataTagged output 2
2543        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2544        res.tag();  // DataTagged output
2545        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2546        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2547    
2548        // Get the views
2549        DataArrayView view_0 = tmp_0->getDefaultValue();
2550        DataArrayView view_1 = tmp_1->getDefaultValue();
2551        DataArrayView view_2 = tmp_2->getDefaultValue();
2552        // Get the pointers to the actual data
2553        double *ptr_0 = &((view_0.getData())[0]);
2554        double *ptr_1 = &((view_1.getData())[0]);
2555        double *ptr_2 = &((view_2.getData())[0]);
2556        // Compute an MVP for the default
2557        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2558        // Merge the tags
2559        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2560        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2561        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2562        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2563          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2564        }
2565        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2566          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2567        }
2568        // Compute an MVP for each tag
2569        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2570        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2571          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2572          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2573          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2574          double *ptr_0 = &view_0.getData(0);
2575          double *ptr_1 = &view_1.getData(0);
2576          double *ptr_2 = &view_2.getData(0);
2577          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2578        }
2579    
2580      }
2581      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2582    
2583        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2584        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2585        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2586        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2587        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2588        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2589        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2590        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2591        int sampleNo_0,dataPointNo_0;
2592        int numSamples_0 = arg_0_Z.getNumSamples();
2593        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2594        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2595        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2596          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2597          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2598          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2599            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2600            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2601            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2602            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2603            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2604          }
2605        }
2606    
2607      }
2608      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2609    
2610        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2611        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2612        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2613        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2614        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2615        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2616        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2617        int sampleNo_0,dataPointNo_0;
2618        int numSamples_0 = arg_0_Z.getNumSamples();
2619        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2620        int offset_1 = tmp_1->getPointOffset(0,0);
2621        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2622        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2623          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2624            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2625            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2626            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2627            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2628            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2629            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2630          }
2631        }
2632    
2633    
2634      }
2635      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2636    
2637        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2638        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2639        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2640        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2641        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2642        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2643        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2644        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2645        int sampleNo_0,dataPointNo_0;
2646        int numSamples_0 = arg_0_Z.getNumSamples();
2647        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2648        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2649        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2650          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2651          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2652          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2653            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2654            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2655            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2656            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2657            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2658          }
2659        }
2660    
2661      }
2662      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2663    
2664        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2665        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2666        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2667        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2668        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2669        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2670        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2671        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2672        int sampleNo_0,dataPointNo_0;
2673        int numSamples_0 = arg_0_Z.getNumSamples();
2674        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2675        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2676        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2677          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2678            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2679            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2680            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2681            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2682            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2683            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2684            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2685          }
2686        }
2687    
2688      }
2689      else {
2690        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2691      }
2692    
2693      return res;
2694    }
2695    
2696    DataAbstract*
2697    Data::borrowData() const
2698    {
2699      return m_data.get();
2700    }
2701    
2702  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2703    
2704  void  void

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

  ViewVC Help
Powered by ViewVC 1.1.26