/[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 1028 by gross, Wed Mar 14 00:15:24 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()) || (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()) || (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()) || (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()) || (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 1073  Data::acos() const Line 1079  Data::acos() const
1079    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acos);
1080  }  }
1081    
1082    
1083  Data  Data
1084  Data::atan() const  Data::atan() const
1085  {  {
# Line 1109  Data::tanh() const Line 1116  Data::tanh() const
1116    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::tanh);
1117  }  }
1118    
1119    
1120    Data
1121    Data::erf() const
1122    {
1123    #if defined DOPROF
1124      profData->unary++;
1125    #endif
1126    #ifdef _WIN32
1127      throw DataException("Error - Data:: erf function is not supported on _WIN32 platforms.");
1128    #else
1129      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1130    #endif
1131    }
1132    
1133  Data  Data
1134  Data::asinh() const  Data::asinh() const
1135  {  {
1136  #if defined DOPROF  #if defined DOPROF
1137    profData->unary++;    profData->unary++;
1138  #endif  #endif
1139    #ifdef _WIN32
1140      return escript::unaryOp(*this,escript::asinh_substitute);
1141    #else
1142    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::asinh);
1143    #endif
1144  }  }
1145    
1146  Data  Data
# Line 1124  Data::acosh() const Line 1149  Data::acosh() const
1149  #if defined DOPROF  #if defined DOPROF
1150    profData->unary++;    profData->unary++;
1151  #endif  #endif
1152    #ifdef _WIN32
1153      return escript::unaryOp(*this,escript::acosh_substitute);
1154    #else
1155    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::acosh);
1156    #endif
1157  }  }
1158    
1159  Data  Data
# Line 1133  Data::atanh() const Line 1162  Data::atanh() const
1162  #if defined DOPROF  #if defined DOPROF
1163    profData->unary++;    profData->unary++;
1164  #endif  #endif
1165    #ifdef _WIN32
1166      return escript::unaryOp(*this,escript::atanh_substitute);
1167    #else
1168    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);    return escript::unaryOp(*this,(Data::UnaryDFunPtr)::atanh);
1169    #endif
1170  }  }
1171    
1172  Data  Data
# Line 1316  Data::minval() const Line 1349  Data::minval() const
1349  }  }
1350    
1351  Data  Data
1352  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1353  {  {
1354  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1355    profData->reduction2++;       #if defined DOPROF
1356  #endif       profData->unary++;
1357    Trace trace_func;       #endif
1358    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1359         DataArrayView::ShapeType ev_shape;
1360         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1361         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1362         int rank=getDataPointRank();
1363         if (rank<2) {
1364            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1365         }
1366         if (axis0<0 || axis0>rank-1) {
1367            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1368         }
1369         if (axis1<0 || axis1>rank-1) {
1370             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1371         }
1372         if (axis0 == axis1) {
1373             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1374         }
1375         if (axis0 > axis1) {
1376             axis0_tmp=axis1;
1377             axis1_tmp=axis0;
1378         } else {
1379             axis0_tmp=axis0;
1380             axis1_tmp=axis1;
1381         }
1382         for (int i=0; i<rank; i++) {
1383           if (i == axis0_tmp) {
1384              ev_shape.push_back(s[axis1_tmp]);
1385           } else if (i == axis1_tmp) {
1386              ev_shape.push_back(s[axis0_tmp]);
1387           } else {
1388              ev_shape.push_back(s[i]);
1389           }
1390         }
1391         Data ev(0.,ev_shape,getFunctionSpace());
1392         ev.typeMatchRight(*this);
1393         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1394         return ev;
1395    
1396  }  }
1397    
1398  Data  Data
# Line 1388  Data::nonsymmetric() const Line 1458  Data::nonsymmetric() const
1458  }  }
1459    
1460  Data  Data
1461  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1462  {  {
1463       #if defined DOPROF       #if defined DOPROF
1464          profData->unary++;          profData->unary++;
# Line 1398  Data::matrixtrace(int axis_offset) const Line 1468  Data::matrixtrace(int axis_offset) const
1468          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1469          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1470          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1471          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1472          return ev;          return ev;
1473       }       }
1474       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1483  Data::matrixtrace(int axis_offset) const
1483          }          }
1484          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1485          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1486          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1487          return ev;          return ev;
1488       }       }
1489       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1502  Data::matrixtrace(int axis_offset) const
1502      }      }
1503          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1504          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1505      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1506          return ev;          return ev;
1507       }       }
1508       else {       else {
1509          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.");
1510       }       }
1511  }  }
1512    
1513  Data  Data
1514  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1515  {  {
1516  #if defined DOPROF       #if defined DOPROF
1517       profData->reduction2++;       profData->unary++;
1518  #endif       #endif
1519       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1520       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1521       // 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 1577  Data::eigenvalues_and_eigenvectors(const
1577  }  }
1578    
1579  const boost::python::tuple  const boost::python::tuple
1580  Data::mindp() const  Data::minGlobalDataPoint() const
1581  {  {
1582    // 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
1583    // 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
1584    // surrounding function    // surrounding function
1585    
   int SampleNo;  
1586    int DataPointNo;    int DataPointNo;
1587      int ProcNo;    int ProcNo;
1588    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1589    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1590  }  }
1591    
1592  void  void
1593  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1594                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1595  {  {
1596    int i,j;    int i,j;
1597    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1649  Data::calc_mindp(  int& ProcNo,
1649  #else  #else
1650      ProcNo = 0;      ProcNo = 0;
1651  #endif  #endif
1652    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1653  }  }
1654    
1655  void  void
# Line 1609  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,plus<double>());    binaryOp(right,plus<double>());
1680    return (*this);    return (*this);
1681  }  }
# Line 1619  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,plus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,plus<double>());  
1688    return (*this);    return (*this);
1689  }  }
1690    
# Line 1635  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,minus<double>());    binaryOp(right,minus<double>());
1698    return (*this);    return (*this);
1699  }  }
# Line 1645  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,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1706    return (*this);    return (*this);
1707  }  }
1708    
# Line 1661  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,multiplies<double>());    binaryOp(right,multiplies<double>());
1716    return (*this);    return (*this);
1717  }  }
# Line 1671  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,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1724    return (*this);    return (*this);
1725  }  }
1726    
# Line 1687  Data::operator/=(const Data& right) Line 1730  Data::operator/=(const Data& right)
1730    if (isProtected()) {    if (isProtected()) {
1731          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1732    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1733    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1734    return (*this);    return (*this);
1735  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1737  Data::operator/=(const Data& right)
1737  Data&  Data&
1738  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1739  {  {
1740    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1741          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1742    return (*this);    return (*this);
1743  }  }
1744    
1745  Data  Data
1746  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1747  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1748    Data left_d(left,*this);    Data left_d(left,*this);
1749    return left_d.powD(*this);    return left_d.powD(*this);
1750  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1752  Data::rpowO(const boost::python::object&
1752  Data  Data
1753  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1754  {  {
1755  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1756    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1757  }  }
1758    
1759  Data  Data
1760  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1761  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1762    Data result;    Data result;
1763    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1764    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1765         result.binaryOp(*this,escript::rpow);
1766      } else {
1767         result.copy(*this);
1768         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1769      }
1770    return result;    return result;
1771  }  }
1772    
# Line 1753  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.copy(right);
1784         result+=left;
1785      } else {
1786         result.copy(left);
1787         result+=right;
1788      }
1789    return result;    return result;
1790  }  }
1791    
# Line 1766  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=right.neg();
1802         result+=left;
1803      } else {
1804         result.copy(left);
1805         result-=right;
1806      }
1807    return result;    return result;
1808  }  }
1809    
# Line 1779  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.copy(right);
1820         result*=left;
1821      } else {
1822         result.copy(left);
1823         result*=right;
1824      }
1825    return result;    return result;
1826  }  }
1827    
# Line 1792  escript::operator/(const Data& left, con Line 1833  escript::operator/(const Data& left, con
1833    Data result;    Data result;
1834    //    //
1835    // perform a deep copy    // perform a deep copy
1836    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1837    result/=right;       result=right.oneOver();
1838         result*=left;
1839      } else {
1840         result.copy(left);
1841         result/=right;
1842      }
1843    return result;    return result;
1844  }  }
1845    
# Line 1802  escript::operator/(const Data& left, con Line 1848  escript::operator/(const Data& left, con
1848  Data  Data
1849  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1850  {  {
1851    //    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;  
1852  }  }
1853    
1854  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1856  escript::operator+(const Data& left, con
1856  Data  Data
1857  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1858  {  {
1859    //    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;  
1860  }  }
1861    
1862  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1864  escript::operator-(const Data& left, con
1864  Data  Data
1865  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1866  {  {
1867    //    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;  
1868  }  }
1869    
1870  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1872  escript::operator*(const Data& left, con
1872  Data  Data
1873  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1874  {  {
1875    //    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;  
1876  }  }
1877    
1878  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1880  escript::operator/(const Data& left, con
1880  Data  Data
1881  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1882  {  {
1883    //    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;  
1884  }  }
1885    
1886  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1888  escript::operator+(const boost::python::
1888  Data  Data
1889  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1890  {  {
1891    //    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;  
1892  }  }
1893    
1894  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1896  escript::operator-(const boost::python::
1896  Data  Data
1897  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1898  {  {
1899    //    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;  
1900  }  }
1901    
1902  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1904  escript::operator*(const boost::python::
1904  Data  Data
1905  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1906  {  {
1907    //    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;  
1908  }  }
1909    
1910  //  //
# Line 2116  Data::getTagNumber(int dpno) Line 2110  Data::getTagNumber(int dpno)
2110    return m_data->getTagNumber(dpno);    return m_data->getTagNumber(dpno);
2111  }  }
2112    
 /* 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);  
           }  
         }  
       }  
     }  
   }  
   
 }  
   
2113  void  void
2114  Data::archiveData(const std::string fileName)  Data::archiveData(const std::string fileName)
2115  {  {
# Line 2241  Data::archiveData(const std::string file Line 2151  Data::archiveData(const std::string file
2151    DataArrayView::ShapeType dataPointShape = getDataPointShape();    DataArrayView::ShapeType dataPointShape = getDataPointShape();
2152    vector<int> referenceNumbers(noSamples);    vector<int> referenceNumbers(noSamples);
2153    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2154      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceNoFromSampleNo(sampleNo);      referenceNumbers[sampleNo] = getFunctionSpace().getReferenceIDOfSample(sampleNo);
2155    }    }
2156    vector<int> tagNumbers(noSamples);    vector<int> tagNumbers(noSamples);
2157    if (isTagged()) {    if (isTagged()) {
# Line 2450  Data::extractData(const std::string file Line 2360  Data::extractData(const std::string file
2360      throw DataException("extractData Error: incompatible FunctionSpace");      throw DataException("extractData Error: incompatible FunctionSpace");
2361    }    }
2362    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {    for (int sampleNo=0; sampleNo<noSamples; sampleNo++) {
2363      if (referenceNumbers[sampleNo] != fspace.getReferenceNoFromSampleNo(sampleNo)) {      if (referenceNumbers[sampleNo] != fspace.getReferenceIDOfSample(sampleNo)) {
2364        throw DataException("extractData Error: incompatible FunctionSpace");        throw DataException("extractData Error: incompatible FunctionSpace");
2365      }      }
2366    }    }
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2451  ostream& escript::operator<<(ostream& o,
2451    return o;    return o;
2452  }  }
2453    
2454    Data
2455    escript::C_GeneralTensorProduct(Data& arg_0,
2456                         Data& arg_1,
2457                         int axis_offset,
2458                         int transpose)
2459    {
2460      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2461      // SM is the product of the last axis_offset entries in arg_0.getShape().
2462    
2463      #if defined DOPROF
2464        // profData->binary++;
2465      #endif
2466    
2467      // Interpolate if necessary and find an appropriate function space
2468      Data arg_0_Z, arg_1_Z;
2469      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2470        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2471          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2472          arg_1_Z = Data(arg_1);
2473        }
2474        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2475          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2476          arg_0_Z =Data(arg_0);
2477        }
2478        else {
2479          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2480        }
2481      } else {
2482          arg_0_Z = Data(arg_0);
2483          arg_1_Z = Data(arg_1);
2484      }
2485      // Get rank and shape of inputs
2486      int rank0 = arg_0_Z.getDataPointRank();
2487      int rank1 = arg_1_Z.getDataPointRank();
2488      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2489      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2490    
2491      // Prepare for the loops of the product and verify compatibility of shapes
2492      int start0=0, start1=0;
2493      if (transpose == 0)       {}
2494      else if (transpose == 1)  { start0 = axis_offset; }
2495      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2496      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2497    
2498      // Adjust the shapes for transpose
2499      DataArrayView::ShapeType tmpShape0;
2500      DataArrayView::ShapeType tmpShape1;
2501      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2502      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2503    
2504    #if 0
2505      // For debugging: show shape after transpose
2506      char tmp[100];
2507      std::string shapeStr;
2508      shapeStr = "(";
2509      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2510      shapeStr += ")";
2511      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2512      shapeStr = "(";
2513      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2514      shapeStr += ")";
2515      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2516    #endif
2517    
2518      // Prepare for the loops of the product
2519      int SL=1, SM=1, SR=1;
2520      for (int i=0; i<rank0-axis_offset; i++)   {
2521        SL *= tmpShape0[i];
2522      }
2523      for (int i=rank0-axis_offset; i<rank0; i++)   {
2524        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2525          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2526        }
2527        SM *= tmpShape0[i];
2528      }
2529      for (int i=axis_offset; i<rank1; i++)     {
2530        SR *= tmpShape1[i];
2531      }
2532    
2533      // Define the shape of the output
2534      DataArrayView::ShapeType shape2;
2535      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2536      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2537    
2538      // Declare output Data object
2539      Data res;
2540    
2541      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2542        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2543        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2544        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2545        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2546        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2547      }
2548      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2549    
2550        // Prepare the DataConstant input
2551        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2552        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2553    
2554        // Borrow DataTagged input from Data object
2555        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2556        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2557    
2558        // Prepare a DataTagged output 2
2559        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2560        res.tag();
2561        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2562        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2563    
2564        // Prepare offset into DataConstant
2565        int offset_0 = tmp_0->getPointOffset(0,0);
2566        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2567        // Get the views
2568        DataArrayView view_1 = tmp_1->getDefaultValue();
2569        DataArrayView view_2 = tmp_2->getDefaultValue();
2570        // Get the pointers to the actual data
2571        double *ptr_1 = &((view_1.getData())[0]);
2572        double *ptr_2 = &((view_2.getData())[0]);
2573        // Compute an MVP for the default
2574        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2575        // Compute an MVP for each tag
2576        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2577        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2578        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2579          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2580          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2581          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2582          double *ptr_1 = &view_1.getData(0);
2583          double *ptr_2 = &view_2.getData(0);
2584          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2585        }
2586    
2587      }
2588      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2589    
2590        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2591        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2592        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2593        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2594        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2595        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2596        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2597        int sampleNo_1,dataPointNo_1;
2598        int numSamples_1 = arg_1_Z.getNumSamples();
2599        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2600        int offset_0 = tmp_0->getPointOffset(0,0);
2601        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2602        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2603          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2604            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2605            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2606            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2607            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2608            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2609            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2610          }
2611        }
2612    
2613      }
2614      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2615    
2616        // Borrow DataTagged input from Data object
2617        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2618        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2619    
2620        // Prepare the DataConstant input
2621        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2622        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2623    
2624        // Prepare a DataTagged output 2
2625        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2626        res.tag();
2627        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2628        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2629    
2630        // Prepare offset into DataConstant
2631        int offset_1 = tmp_1->getPointOffset(0,0);
2632        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2633        // Get the views
2634        DataArrayView view_0 = tmp_0->getDefaultValue();
2635        DataArrayView view_2 = tmp_2->getDefaultValue();
2636        // Get the pointers to the actual data
2637        double *ptr_0 = &((view_0.getData())[0]);
2638        double *ptr_2 = &((view_2.getData())[0]);
2639        // Compute an MVP for the default
2640        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2641        // Compute an MVP for each tag
2642        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2643        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2644        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2645          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2646          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2647          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2648          double *ptr_0 = &view_0.getData(0);
2649          double *ptr_2 = &view_2.getData(0);
2650          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2651        }
2652    
2653      }
2654      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2655    
2656        // Borrow DataTagged input from Data object
2657        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2658        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2659    
2660        // Borrow DataTagged input from Data object
2661        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2662        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2663    
2664        // Prepare a DataTagged output 2
2665        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2666        res.tag();  // DataTagged output
2667        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2668        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2669    
2670        // Get the views
2671        DataArrayView view_0 = tmp_0->getDefaultValue();
2672        DataArrayView view_1 = tmp_1->getDefaultValue();
2673        DataArrayView view_2 = tmp_2->getDefaultValue();
2674        // Get the pointers to the actual data
2675        double *ptr_0 = &((view_0.getData())[0]);
2676        double *ptr_1 = &((view_1.getData())[0]);
2677        double *ptr_2 = &((view_2.getData())[0]);
2678        // Compute an MVP for the default
2679        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2680        // Merge the tags
2681        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2682        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2683        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2684        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2685          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2686        }
2687        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2688          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2689        }
2690        // Compute an MVP for each tag
2691        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2692        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2693          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2694          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2695          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2696          double *ptr_0 = &view_0.getData(0);
2697          double *ptr_1 = &view_1.getData(0);
2698          double *ptr_2 = &view_2.getData(0);
2699          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2700        }
2701    
2702      }
2703      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2704    
2705        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2706        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2707        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2708        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2709        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2710        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2711        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2712        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2713        int sampleNo_0,dataPointNo_0;
2714        int numSamples_0 = arg_0_Z.getNumSamples();
2715        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2716        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2717        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2718          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2719          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2720          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2721            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2722            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2723            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2724            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2725            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2726          }
2727        }
2728    
2729      }
2730      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2731    
2732        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2733        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2734        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2735        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2736        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2737        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2738        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2739        int sampleNo_0,dataPointNo_0;
2740        int numSamples_0 = arg_0_Z.getNumSamples();
2741        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2742        int offset_1 = tmp_1->getPointOffset(0,0);
2743        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2744        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2745          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2746            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2747            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2748            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2749            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2750            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2751            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2752          }
2753        }
2754    
2755    
2756      }
2757      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2758    
2759        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2760        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2761        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2762        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2763        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2764        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2765        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2766        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2767        int sampleNo_0,dataPointNo_0;
2768        int numSamples_0 = arg_0_Z.getNumSamples();
2769        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2770        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2771        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2772          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2773          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2774          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2775            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2776            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2777            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2778            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2779            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2780          }
2781        }
2782    
2783      }
2784      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2785    
2786        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2787        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2788        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2789        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2790        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2791        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2792        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2793        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2794        int sampleNo_0,dataPointNo_0;
2795        int numSamples_0 = arg_0_Z.getNumSamples();
2796        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2797        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2798        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2799          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2800            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2801            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2802            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2803            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2804            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2805            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2806            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2807          }
2808        }
2809    
2810      }
2811      else {
2812        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2813      }
2814    
2815      return res;
2816    }
2817    
2818    DataAbstract*
2819    Data::borrowData() const
2820    {
2821      return m_data.get();
2822    }
2823    
2824  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2825    
2826  void  void

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

  ViewVC Help
Powered by ViewVC 1.1.26