/[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 921 by gross, Fri Jan 5 00:54:37 2007 UTC
# Line 350  Data::isTagged() const Line 350  Data::isTagged() const
350    return (temp!=0);    return (temp!=0);
351  }  }
352    
 /* TODO */  
 /* global reduction -- the local data being empty does not imply that it is empty on other processers*/  
353  bool  bool
354  Data::isEmpty() const  Data::isEmpty() const
355  {  {
# Line 422  Data::tag() Line 420  Data::tag()
420    }    }
421  }  }
422    
423  void  Data
424  Data::reshapeDataPoint(const DataArrayView::ShapeType& shape)  Data::oneOver() const
425  {  {
426    m_data->reshapeDataPoint(shape);  #if defined DOPROF
427      profData->where++;
428    #endif
429      return escript::unaryOp(*this,bind1st(divides<double>(),1.));
430  }  }
431    
432  Data  Data
# Line 547  Data::getDataPointShape() const Line 548  Data::getDataPointShape() const
548    return getPointDataView().getShape();    return getPointDataView().getShape();
549  }  }
550    
551    
552    
553  void  void
554  Data::fillFromNumArray(const boost::python::numeric::array num_array)  Data::fillFromNumArray(const boost::python::numeric::array num_array)
555  {  {
# Line 639  Data::convertToNumArray() Line 642  Data::convertToNumArray()
642        }        }
643        if (dataPointRank==1) {        if (dataPointRank==1) {
644          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
645            numArray[dataPoint][i]=dataPointView(i);            numArray[make_tuple(dataPoint,i)]=dataPointView(i);
646          }          }
647        }        }
648        if (dataPointRank==2) {        if (dataPointRank==2) {
649          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
650            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
651              numArray[dataPoint][i][j] = dataPointView(i,j);              numArray[make_tuple(dataPoint,i,j)] = dataPointView(i,j);
652            }            }
653          }          }
654        }        }
# Line 653  Data::convertToNumArray() Line 656  Data::convertToNumArray()
656          for (int i=0; i<dataPointShape[0]; i++) {          for (int i=0; i<dataPointShape[0]; i++) {
657            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
658              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
659                numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                numArray[make_tuple(dataPoint,i,j,k)]=dataPointView(i,j,k);
660              }              }
661            }            }
662          }          }
# Line 663  Data::convertToNumArray() Line 666  Data::convertToNumArray()
666            for (int j=0; j<dataPointShape[1]; j++) {            for (int j=0; j<dataPointShape[1]; j++) {
667              for (int k=0; k<dataPointShape[2]; k++) {              for (int k=0; k<dataPointShape[2]; k++) {
668                for (int l=0; l<dataPointShape[3]; l++) {                for (int l=0; l<dataPointShape[3]; l++) {
669                  numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                  numArray[make_tuple(dataPoint,i,j,k,l)]=dataPointView(i,j,k,l);
670                }                }
671              }              }
672            }            }
# Line 678  Data::convertToNumArray() Line 681  Data::convertToNumArray()
681    return numArray;    return numArray;
682  }  }
683    
684  const  
685    const
686  boost::python::numeric::array  boost::python::numeric::array
687  Data::convertToNumArrayFromSampleNo(int sampleNo)  Data:: getValueOfDataPoint(int dataPointNo)
688  {  {
689    //    size_t length=0;
690    // Check a valid sample number has been supplied    int i, j, k, l;
   if (sampleNo >= getNumSamples()) {  
     throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");  
   }  
   
   //  
   // determine the number of data points per sample  
   int numDataPointsPerSample = getNumDataPointsPerSample();  
   
691    //    //
692    // determine the rank and shape of each data point    // determine the rank and shape of each data point
693    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 702  Data::convertToNumArrayFromSampleNo(int Line 698  Data::convertToNumArrayFromSampleNo(int
698    boost::python::numeric::array numArray(0.0);    boost::python::numeric::array numArray(0.0);
699    
700    //    //
701    // the rank of the returned numeric array will be the rank of    // the shape of the returned numeric array will be the same
702    // the data points, plus one. Where the rank of the array is n,    // as that of the data point
703    // the last n-1 dimensions will be equal to the shape of the    int arrayRank = dataPointRank;
704    // data points, whilst the first dimension will be equal to the    DataArrayView::ShapeType arrayShape = dataPointShape;
   // total number of data points. Thus the array will consist of  
   // a serial vector of the data points.  
   int arrayRank = dataPointRank + 1;  
   DataArrayView::ShapeType arrayShape;  
   arrayShape.push_back(numDataPointsPerSample);  
   for (int d=0; d<dataPointRank; d++) {  
      arrayShape.push_back(dataPointShape[d]);  
   }  
705    
706    //    //
707    // resize the numeric array to the shape just calculated    // resize the numeric array to the shape just calculated
708      if (arrayRank==0) {
709        numArray.resize(1);
710      }
711    if (arrayRank==1) {    if (arrayRank==1) {
712      numArray.resize(arrayShape[0]);      numArray.resize(arrayShape[0]);
713    }    }
# Line 729  Data::convertToNumArrayFromSampleNo(int Line 720  Data::convertToNumArrayFromSampleNo(int
720    if (arrayRank==4) {    if (arrayRank==4) {
721      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
722    }    }
   if (arrayRank==5) {  
     numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3],arrayShape[4]);  
   }  
723    
724    //    if (getNumDataPointsPerSample()>0) {
725    // loop through each data point in turn, loading the values for that data point         int sampleNo = dataPointNo/getNumDataPointsPerSample();
726    // into the numeric array.         int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
727    for (int dataPoint = 0; dataPoint < numDataPointsPerSample; dataPoint++) {         //
728      DataArrayView dataPointView = getDataPoint(sampleNo, dataPoint);         // Check a valid sample number has been supplied
729      if (dataPointRank==0) {         if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
730        numArray[dataPoint]=dataPointView();             throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
731      }         }
732      if (dataPointRank==1) {                
733        for (int i=0; i<dataPointShape[0]; i++) {         //
734          numArray[dataPoint][i]=dataPointView(i);         // Check a valid data point number has been supplied
735        }         if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
736      }             throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
737      if (dataPointRank==2) {         }
738        for (int i=0; i<dataPointShape[0]; i++) {         // TODO: global error handling
739          for (int j=0; j<dataPointShape[1]; j++) {         // create a view of the data if it is stored locally
740            numArray[dataPoint][i][j] = dataPointView(i,j);         DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
741          }          
742        }         switch( dataPointRank ){
743      }              case 0 :
744      if (dataPointRank==3) {                  numArray[0] = dataPointView();
745        for (int i=0; i<dataPointShape[0]; i++) {                  break;
746          for (int j=0; j<dataPointShape[1]; j++) {              case 1 :        
747            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
748              numArray[dataPoint][i][j][k]=dataPointView(i,j,k);                      numArray[i]=dataPointView(i);
749            }                  break;
750          }              case 2 :        
751        }                  for( i=0; i<dataPointShape[0]; i++ )
752      }                      for( j=0; j<dataPointShape[1]; j++)
753      if (dataPointRank==4) {                          numArray[make_tuple(i,j)]=dataPointView(i,j);
754        for (int i=0; i<dataPointShape[0]; i++) {                  break;
755          for (int j=0; j<dataPointShape[1]; j++) {              case 3 :        
756            for (int k=0; k<dataPointShape[2]; k++) {                  for( i=0; i<dataPointShape[0]; i++ )
757              for (int l=0; l<dataPointShape[3]; l++) {                      for( j=0; j<dataPointShape[1]; j++ )
758                numArray[dataPoint][i][j][k][l]=dataPointView(i,j,k,l);                          for( k=0; k<dataPointShape[2]; k++)
759              }                              numArray[make_tuple(i,j,k)]=dataPointView(i,j,k);
760            }                  break;
761          }              case 4 :
762        }                  for( i=0; i<dataPointShape[0]; i++ )
763      }                      for( j=0; j<dataPointShape[1]; j++ )
764                            for( k=0; k<dataPointShape[2]; k++ )
765                                for( l=0; l<dataPointShape[3]; l++)
766                                    numArray[make_tuple(i,j,k,l)]=dataPointView(i,j,k,l);
767                    break;
768        }
769    }    }
   
770    //    //
771    // return the loaded array    // return the array
772    return numArray;    return numArray;
 }  
773    
774  const  }
 boost::python::numeric::array  
 Data::convertToNumArrayFromDPNo(int procNo,  
                                 int sampleNo,  
                                                                 int dataPointNo)  
775    
776    void
777    Data::setValueOfDataPoint(int dataPointNo, const boost::python::numeric::array num_array)
778  {  {
779      size_t length=0;    if (isProtected()) {
780      int i, j, k, l, pos;          throw DataException("Error - attempt to update protected Data object.");
781      }
782      //
783      // check rank
784      if (num_array.getrank()<getDataPointRank())
785          throw DataException("Rank of numarray does not match Data object rank");
786    
787    //    //
788    // Check a valid sample number has been supplied    // check shape of num_array
789    if (sampleNo >= getNumSamples()) {    for (int i=0; i<getDataPointRank(); i++) {
790      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");      if (extract<int>(num_array.getshape()[i])!=getDataPointShape()[i])
791           throw DataException("Shape of numarray does not match Data object rank");
792    }    }
   
793    //    //
794    // Check a valid data point number has been supplied    // make sure data is expanded:
795    if (dataPointNo >= getNumDataPointsPerSample()) {    if (!isExpanded()) {
796      throw DataException("Error - Data::convertToNumArray: invalid dataPointNo.");      expand();
797      }
798      if (getNumDataPointsPerSample()>0) {
799           int sampleNo = dataPointNo/getNumDataPointsPerSample();
800           int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
801           m_data->copyToDataPoint(sampleNo, dataPointNoInSample,num_array);
802      } else {
803           m_data->copyToDataPoint(-1, 0,num_array);
804    }    }
805    }
806    
807    const
808    boost::python::numeric::array
809    Data::getValueOfGlobalDataPoint(int procNo, int dataPointNo)
810    {
811      size_t length=0;
812      int i, j, k, l, pos;
813    //    //
814    // determine the rank and shape of each data point    // determine the rank and shape of each data point
815    int dataPointRank = getDataPointRank();    int dataPointRank = getDataPointRank();
# Line 835  Data::convertToNumArrayFromDPNo(int proc Line 843  Data::convertToNumArrayFromDPNo(int proc
843      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);      numArray.resize(arrayShape[0],arrayShape[1],arrayShape[2],arrayShape[3]);
844    }    }
845    
846      // added for the MPI communication    // added for the MPI communication
847      length=1;    length=1;
848      for( i=0; i<arrayRank; i++ )    for( i=0; i<arrayRank; i++ ) length *= arrayShape[i];
849          length *= arrayShape[i];    double *tmpData = new double[length];
     double *tmpData = new double[length];  
850    
851    //    //
852    // load the values for the data point into the numeric array.    // load the values for the data point into the numeric array.
853    
854      // updated for the MPI case      // updated for the MPI case
855      if( get_MPIRank()==procNo ){      if( get_MPIRank()==procNo ){
856                 if (getNumDataPointsPerSample()>0) {
857                    int sampleNo = dataPointNo/getNumDataPointsPerSample();
858                    int dataPointNoInSample = dataPointNo - sampleNo * getNumDataPointsPerSample();
859                    //
860                    // Check a valid sample number has been supplied
861                    if (sampleNo >= getNumSamples() or sampleNo < 0 ) {
862                      throw DataException("Error - Data::convertToNumArray: invalid sampleNo.");
863                    }
864                  
865                    //
866                    // Check a valid data point number has been supplied
867                    if (dataPointNoInSample >= getNumDataPointsPerSample() or dataPointNoInSample < 0) {
868                      throw DataException("Error - Data::convertToNumArray: invalid dataPointNoInSample.");
869                    }
870                    // TODO: global error handling
871          // create a view of the data if it is stored locally          // create a view of the data if it is stored locally
872          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNoInSample);
873                    
874          // pack the data from the view into tmpData for MPI communication          // pack the data from the view into tmpData for MPI communication
875          pos=0;          pos=0;
# Line 878  Data::convertToNumArrayFromDPNo(int proc Line 900  Data::convertToNumArrayFromDPNo(int proc
900                                  tmpData[pos]=dataPointView(i,j,k,l);                                  tmpData[pos]=dataPointView(i,j,k,l);
901                  break;                  break;
902          }          }
903                }
904      }      }
905  #ifdef PASO_MPI          #ifdef PASO_MPI
906          // broadcast the data to all other processes          // broadcast the data to all other processes
907          MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );      MPI_Bcast( tmpData, length, MPI_DOUBLE, procNo, get_MPIComm() );
908  #endif          #endif
909    
910      // unpack the data      // unpack the data
911      switch( dataPointRank ){      switch( dataPointRank ){
912          case 0 :          case 0 :
913              numArray[i]=tmpData[0];              numArray[0]=tmpData[0];
914              break;              break;
915          case 1 :                  case 1 :        
916              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
# Line 896  Data::convertToNumArrayFromDPNo(int proc Line 919  Data::convertToNumArrayFromDPNo(int proc
919          case 2 :                  case 2 :        
920              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
921                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
922                      tmpData[i+j*dataPointShape[0]];                     numArray[make_tuple(i,j)]=tmpData[i+j*dataPointShape[0]];
923              break;              break;
924          case 3 :                  case 3 :        
925              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
926                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
927                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
928                          tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];                          numArray[make_tuple(i,j,k)]=tmpData[i+dataPointShape[0]*(j*+k*dataPointShape[1])];
929              break;              break;
930          case 4 :          case 4 :
931              for( i=0; i<dataPointShape[0]; i++ )              for( i=0; i<dataPointShape[0]; i++ )
932                  for( j=0; j<dataPointShape[1]; j++ )                  for( j=0; j<dataPointShape[1]; j++ )
933                      for( k=0; k<dataPointShape[2]; k++ )                      for( k=0; k<dataPointShape[2]; k++ )
934                          for( l=0; l<dataPointShape[3]; l++ )                          for( l=0; l<dataPointShape[3]; l++ )
935                              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]))];
936              break;              break;
937      }      }
938    
939      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);  
           }  
         }  
       }  
     }  
   }  
 */  
   
940    //    //
941    // return the loaded array    // return the loaded array
942    return numArray;    return numArray;
943  }  }
944    
945    
946    
947  boost::python::numeric::array  boost::python::numeric::array
948  Data::integrate() const  Data::integrate() const
949  {  {
# Line 1029  Data::integrate() const Line 1016  Data::integrate() const
1016  }  }
1017    
1018  Data  Data
1019    Data::erf() const
1020    {
1021    #if defined DOPROF
1022      profData->unary++;
1023    #endif
1024      return escript::unaryOp(*this,(Data::UnaryDFunPtr)::erf);
1025    }
1026    
1027    Data
1028  Data::sin() const  Data::sin() const
1029  {  {
1030  #if defined DOPROF  #if defined DOPROF
# Line 1316  Data::minval() const Line 1312  Data::minval() const
1312  }  }
1313    
1314  Data  Data
1315  Data::trace() const  Data::swapaxes(const int axis0, const int axis1) const
1316  {  {
1317  #if defined DOPROF       int axis0_tmp,axis1_tmp;
1318    profData->reduction2++;       #if defined DOPROF
1319  #endif       profData->unary++;
1320    Trace trace_func;       #endif
1321    return dp_algorithm(trace_func,0);       DataArrayView::ShapeType s=getDataPointShape();
1322         DataArrayView::ShapeType ev_shape;
1323         // Here's the equivalent of python s_out=s[axis_offset:]+s[:axis_offset]
1324         // which goes thru all shape vector elements starting with axis_offset (at index=rank wrap around to 0)
1325         int rank=getDataPointRank();
1326         if (rank<2) {
1327            throw DataException("Error - Data::swapaxes argument must have at least rank 2.");
1328         }
1329         if (axis0<0 || axis0>rank-1) {
1330            throw DataException("Error - Data::swapaxes: axis0 must be between 0 and rank-1=" + rank-1);
1331         }
1332         if (axis1<0 || axis1>rank-1) {
1333             throw DataException("Error - Data::swapaxes: axis1 must be between 0 and rank-1=" + rank-1);
1334         }
1335         if (axis0 == axis1) {
1336             throw DataException("Error - Data::swapaxes: axis indices must be different.");
1337         }
1338         if (axis0 > axis1) {
1339             axis0_tmp=axis1;
1340             axis1_tmp=axis0;
1341         } else {
1342             axis0_tmp=axis0;
1343             axis1_tmp=axis1;
1344         }
1345         for (int i=0; i<rank; i++) {
1346           if (i == axis0_tmp) {
1347              ev_shape.push_back(s[axis1_tmp]);
1348           } else if (i == axis1_tmp) {
1349              ev_shape.push_back(s[axis0_tmp]);
1350           } else {
1351              ev_shape.push_back(s[i]);
1352           }
1353         }
1354         Data ev(0.,ev_shape,getFunctionSpace());
1355         ev.typeMatchRight(*this);
1356         m_data->swapaxes(ev.m_data.get(), axis0_tmp, axis1_tmp);
1357         return ev;
1358    
1359  }  }
1360    
1361  Data  Data
# Line 1388  Data::nonsymmetric() const Line 1421  Data::nonsymmetric() const
1421  }  }
1422    
1423  Data  Data
1424  Data::matrixtrace(int axis_offset) const  Data::trace(int axis_offset) const
1425  {  {
1426       #if defined DOPROF       #if defined DOPROF
1427          profData->unary++;          profData->unary++;
# Line 1398  Data::matrixtrace(int axis_offset) const Line 1431  Data::matrixtrace(int axis_offset) const
1431          DataArrayView::ShapeType ev_shape;          DataArrayView::ShapeType ev_shape;
1432          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1433          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1434          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1435          return ev;          return ev;
1436       }       }
1437       if (getDataPointRank()==3) {       if (getDataPointRank()==3) {
# Line 1413  Data::matrixtrace(int axis_offset) const Line 1446  Data::matrixtrace(int axis_offset) const
1446          }          }
1447          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1448          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1449          m_data->matrixtrace(ev.m_data.get(), axis_offset);          m_data->trace(ev.m_data.get(), axis_offset);
1450          return ev;          return ev;
1451       }       }
1452       if (getDataPointRank()==4) {       if (getDataPointRank()==4) {
# Line 1432  Data::matrixtrace(int axis_offset) const Line 1465  Data::matrixtrace(int axis_offset) const
1465      }      }
1466          Data ev(0.,ev_shape,getFunctionSpace());          Data ev(0.,ev_shape,getFunctionSpace());
1467          ev.typeMatchRight(*this);          ev.typeMatchRight(*this);
1468      m_data->matrixtrace(ev.m_data.get(), axis_offset);      m_data->trace(ev.m_data.get(), axis_offset);
1469          return ev;          return ev;
1470       }       }
1471       else {       else {
1472          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.");
1473       }       }
1474  }  }
1475    
1476  Data  Data
1477  Data::transpose(int axis_offset) const  Data::transpose(int axis_offset) const
1478  {  {
1479  #if defined DOPROF       #if defined DOPROF
1480       profData->reduction2++;       profData->unary++;
1481  #endif       #endif
1482       DataArrayView::ShapeType s=getDataPointShape();       DataArrayView::ShapeType s=getDataPointShape();
1483       DataArrayView::ShapeType ev_shape;       DataArrayView::ShapeType ev_shape;
1484       // 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 1540  Data::eigenvalues_and_eigenvectors(const
1540  }  }
1541    
1542  const boost::python::tuple  const boost::python::tuple
1543  Data::mindp() const  Data::minGlobalDataPoint() const
1544  {  {
1545    // 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
1546    // 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
1547    // surrounding function    // surrounding function
1548    
   int SampleNo;  
1549    int DataPointNo;    int DataPointNo;
1550      int ProcNo;    int ProcNo;
1551    calc_mindp(ProcNo,SampleNo,DataPointNo);    calc_minGlobalDataPoint(ProcNo,DataPointNo);
1552    return make_tuple(ProcNo,SampleNo,DataPointNo);    return make_tuple(ProcNo,DataPointNo);
1553  }  }
1554    
1555  void  void
1556  Data::calc_mindp(   int& ProcNo,  Data::calc_minGlobalDataPoint(int& ProcNo,
1557                  int& SampleNo,                          int& DataPointNo) const
         int& DataPointNo) const  
1558  {  {
1559    int i,j;    int i,j;
1560    int lowi=0,lowj=0;    int lowi=0,lowj=0;
# Line 1581  Data::calc_mindp(  int& ProcNo, Line 1612  Data::calc_mindp(  int& ProcNo,
1612  #else  #else
1613      ProcNo = 0;      ProcNo = 0;
1614  #endif  #endif
1615    SampleNo = lowi;    DataPointNo = lowj + lowi * numDPPSample;
   DataPointNo = lowj;  
1616  }  }
1617    
1618  void  void
# Line 1609  Data::operator+=(const Data& right) Line 1639  Data::operator+=(const Data& right)
1639    if (isProtected()) {    if (isProtected()) {
1640          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1641    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1642    binaryOp(right,plus<double>());    binaryOp(right,plus<double>());
1643    return (*this);    return (*this);
1644  }  }
# Line 1619  Data::operator+=(const Data& right) Line 1646  Data::operator+=(const Data& right)
1646  Data&  Data&
1647  Data::operator+=(const boost::python::object& right)  Data::operator+=(const boost::python::object& right)
1648  {  {
1649    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1650          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,plus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,plus<double>());  
1651    return (*this);    return (*this);
1652  }  }
1653    
# Line 1635  Data::operator-=(const Data& right) Line 1657  Data::operator-=(const Data& right)
1657    if (isProtected()) {    if (isProtected()) {
1658          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1659    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1660    binaryOp(right,minus<double>());    binaryOp(right,minus<double>());
1661    return (*this);    return (*this);
1662  }  }
# Line 1645  Data::operator-=(const Data& right) Line 1664  Data::operator-=(const Data& right)
1664  Data&  Data&
1665  Data::operator-=(const boost::python::object& right)  Data::operator-=(const boost::python::object& right)
1666  {  {
1667    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1668          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,minus<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,minus<double>());  
1669    return (*this);    return (*this);
1670  }  }
1671    
# Line 1661  Data::operator*=(const Data& right) Line 1675  Data::operator*=(const Data& right)
1675    if (isProtected()) {    if (isProtected()) {
1676          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1677    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1678    binaryOp(right,multiplies<double>());    binaryOp(right,multiplies<double>());
1679    return (*this);    return (*this);
1680  }  }
# Line 1671  Data::operator*=(const Data& right) Line 1682  Data::operator*=(const Data& right)
1682  Data&  Data&
1683  Data::operator*=(const boost::python::object& right)  Data::operator*=(const boost::python::object& right)
1684  {  {
1685    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1686          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,multiplies<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,multiplies<double>());  
1687    return (*this);    return (*this);
1688  }  }
1689    
# Line 1687  Data::operator/=(const Data& right) Line 1693  Data::operator/=(const Data& right)
1693    if (isProtected()) {    if (isProtected()) {
1694          throw DataException("Error - attempt to update protected Data object.");          throw DataException("Error - attempt to update protected Data object.");
1695    }    }
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1696    binaryOp(right,divides<double>());    binaryOp(right,divides<double>());
1697    return (*this);    return (*this);
1698  }  }
# Line 1697  Data::operator/=(const Data& right) Line 1700  Data::operator/=(const Data& right)
1700  Data&  Data&
1701  Data::operator/=(const boost::python::object& right)  Data::operator/=(const boost::python::object& right)
1702  {  {
1703    if (isProtected()) {    Data tmp(right,getFunctionSpace(),false);
1704          throw DataException("Error - attempt to update protected Data object.");    binaryOp(tmp,divides<double>());
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
   binaryOp(right,divides<double>());  
1705    return (*this);    return (*this);
1706  }  }
1707    
1708  Data  Data
1709  Data::rpowO(const boost::python::object& left) const  Data::rpowO(const boost::python::object& left) const
1710  {  {
   if (isProtected()) {  
         throw DataException("Error - attempt to update protected Data object.");  
   }  
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1711    Data left_d(left,*this);    Data left_d(left,*this);
1712    return left_d.powD(*this);    return left_d.powD(*this);
1713  }  }
# Line 1723  Data::rpowO(const boost::python::object& Line 1715  Data::rpowO(const boost::python::object&
1715  Data  Data
1716  Data::powO(const boost::python::object& right) const  Data::powO(const boost::python::object& right) const
1717  {  {
1718  #if defined DOPROF    Data tmp(right,getFunctionSpace(),false);
1719    profData->binary++;    return powD(tmp);
 #endif  
   Data result;  
   result.copy(*this);  
   result.binaryOp(right,(Data::BinaryDFunPtr)::pow);  
   return result;  
1720  }  }
1721    
1722  Data  Data
1723  Data::powD(const Data& right) const  Data::powD(const Data& right) const
1724  {  {
 #if defined DOPROF  
   profData->binary++;  
 #endif  
1725    Data result;    Data result;
1726    result.copy(*this);    if (getDataPointRank()<right.getDataPointRank()) {
1727    result.binaryOp(right,(Data::BinaryDFunPtr)::pow);       result.copy(right);
1728         result.binaryOp(*this,escript::rpow);
1729      } else {
1730         result.copy(*this);
1731         result.binaryOp(right,(Data::BinaryDFunPtr)::pow);
1732      }
1733    return result;    return result;
1734  }  }
1735    
# Line 1753  escript::operator+(const Data& left, con Line 1742  escript::operator+(const Data& left, con
1742    Data result;    Data result;
1743    //    //
1744    // perform a deep copy    // perform a deep copy
1745    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1746    result+=right;       result.copy(right);
1747         result+=left;
1748      } else {
1749         result.copy(left);
1750         result+=right;
1751      }
1752    return result;    return result;
1753  }  }
1754    
# Line 1766  escript::operator-(const Data& left, con Line 1760  escript::operator-(const Data& left, con
1760    Data result;    Data result;
1761    //    //
1762    // perform a deep copy    // perform a deep copy
1763    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1764    result-=right;       result=right.neg();
1765         result+=left;
1766      } else {
1767         result.copy(left);
1768         result-=right;
1769      }
1770    return result;    return result;
1771  }  }
1772    
# Line 1779  escript::operator*(const Data& left, con Line 1778  escript::operator*(const Data& left, con
1778    Data result;    Data result;
1779    //    //
1780    // perform a deep copy    // perform a deep copy
1781    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1782    result*=right;       result.copy(right);
1783         result*=left;
1784      } else {
1785         result.copy(left);
1786         result*=right;
1787      }
1788    return result;    return result;
1789  }  }
1790    
# Line 1792  escript::operator/(const Data& left, con Line 1796  escript::operator/(const Data& left, con
1796    Data result;    Data result;
1797    //    //
1798    // perform a deep copy    // perform a deep copy
1799    result.copy(left);    if (left.getDataPointRank()<right.getDataPointRank()) {
1800    result/=right;       result=right.oneOver();
1801         result*=left;
1802      } else {
1803         result.copy(left);
1804         result/=right;
1805      }
1806    return result;    return result;
1807  }  }
1808    
# Line 1802  escript::operator/(const Data& left, con Line 1811  escript::operator/(const Data& left, con
1811  Data  Data
1812  escript::operator+(const Data& left, const boost::python::object& right)  escript::operator+(const Data& left, const boost::python::object& right)
1813  {  {
1814    //    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;  
1815  }  }
1816    
1817  //  //
# Line 1818  escript::operator+(const Data& left, con Line 1819  escript::operator+(const Data& left, con
1819  Data  Data
1820  escript::operator-(const Data& left, const boost::python::object& right)  escript::operator-(const Data& left, const boost::python::object& right)
1821  {  {
1822    //    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;  
1823  }  }
1824    
1825  //  //
# Line 1834  escript::operator-(const Data& left, con Line 1827  escript::operator-(const Data& left, con
1827  Data  Data
1828  escript::operator*(const Data& left, const boost::python::object& right)  escript::operator*(const Data& left, const boost::python::object& right)
1829  {  {
1830    //    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;  
1831  }  }
1832    
1833  //  //
# Line 1850  escript::operator*(const Data& left, con Line 1835  escript::operator*(const Data& left, con
1835  Data  Data
1836  escript::operator/(const Data& left, const boost::python::object& right)  escript::operator/(const Data& left, const boost::python::object& right)
1837  {  {
1838    //    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;  
1839  }  }
1840    
1841  //  //
# Line 1866  escript::operator/(const Data& left, con Line 1843  escript::operator/(const Data& left, con
1843  Data  Data
1844  escript::operator+(const boost::python::object& left, const Data& right)  escript::operator+(const boost::python::object& left, const Data& right)
1845  {  {
1846    //    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;  
1847  }  }
1848    
1849  //  //
# Line 1879  escript::operator+(const boost::python:: Line 1851  escript::operator+(const boost::python::
1851  Data  Data
1852  escript::operator-(const boost::python::object& left, const Data& right)  escript::operator-(const boost::python::object& left, const Data& right)
1853  {  {
1854    //    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;  
1855  }  }
1856    
1857  //  //
# Line 1892  escript::operator-(const boost::python:: Line 1859  escript::operator-(const boost::python::
1859  Data  Data
1860  escript::operator*(const boost::python::object& left, const Data& right)  escript::operator*(const boost::python::object& left, const Data& right)
1861  {  {
1862    //    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;  
1863  }  }
1864    
1865  //  //
# Line 1905  escript::operator*(const boost::python:: Line 1867  escript::operator*(const boost::python::
1867  Data  Data
1868  escript::operator/(const boost::python::object& left, const Data& right)  escript::operator/(const boost::python::object& left, const Data& right)
1869  {  {
1870    //    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;  
1871  }  }
1872    
1873  //  //
# Line 2173  Data::getRefValue(int ref, Line 2130  Data::getRefValue(int ref,
2130    if (rank==2) {    if (rank==2) {
2131      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
2132        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2133          value[i][j] = valueView(i,j);          value[make_tuple(i,j)] = valueView(i,j);
2134        }        }
2135      }      }
2136    }    }
# Line 2181  Data::getRefValue(int ref, Line 2138  Data::getRefValue(int ref,
2138      for (int i=0; i < shape[0]; i++) {      for (int i=0; i < shape[0]; i++) {
2139        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2140          for (int k=0; k < shape[2]; k++) {          for (int k=0; k < shape[2]; k++) {
2141            value[i][j][k] = valueView(i,j,k);            value[make_tuple(i,j,k)] = valueView(i,j,k);
2142          }          }
2143        }        }
2144      }      }
# Line 2191  Data::getRefValue(int ref, Line 2148  Data::getRefValue(int ref,
2148        for (int j=0; j < shape[1]; j++) {        for (int j=0; j < shape[1]; j++) {
2149          for (int k=0; k < shape[2]; k++) {          for (int k=0; k < shape[2]; k++) {
2150            for (int l=0; l < shape[3]; l++) {            for (int l=0; l < shape[3]; l++) {
2151              value[i][j][k][l] = valueView(i,j,k,l);              value[make_tuple(i,j,k,l)] = valueView(i,j,k,l);
2152            }            }
2153          }          }
2154        }        }
# Line 2541  ostream& escript::operator<<(ostream& o, Line 2498  ostream& escript::operator<<(ostream& o,
2498    return o;    return o;
2499  }  }
2500    
2501    Data
2502    escript::C_GeneralTensorProduct(Data& arg_0,
2503                         Data& arg_1,
2504                         int axis_offset,
2505                         int transpose)
2506    {
2507      // General tensor product: res(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
2508      // SM is the product of the last axis_offset entries in arg_0.getShape().
2509    
2510      #if defined DOPROF
2511        // profData->binary++;
2512      #endif
2513    
2514      // Interpolate if necessary and find an appropriate function space
2515      Data arg_0_Z, arg_1_Z;
2516      if (arg_0.getFunctionSpace()!=arg_1.getFunctionSpace()) {
2517        if (arg_0.probeInterpolation(arg_1.getFunctionSpace())) {
2518          arg_0_Z = arg_0.interpolate(arg_1.getFunctionSpace());
2519          arg_1_Z = Data(arg_1);
2520        }
2521        else if (arg_1.probeInterpolation(arg_0.getFunctionSpace())) {
2522          arg_1_Z=arg_1.interpolate(arg_0.getFunctionSpace());
2523          arg_0_Z =Data(arg_0);
2524        }
2525        else {
2526          throw DataException("Error - C_GeneralTensorProduct: arguments have incompatible function spaces.");
2527        }
2528      } else {
2529          arg_0_Z = Data(arg_0);
2530          arg_1_Z = Data(arg_1);
2531      }
2532      // Get rank and shape of inputs
2533      int rank0 = arg_0_Z.getDataPointRank();
2534      int rank1 = arg_1_Z.getDataPointRank();
2535      DataArrayView::ShapeType shape0 = arg_0_Z.getDataPointShape();
2536      DataArrayView::ShapeType shape1 = arg_1_Z.getDataPointShape();
2537    
2538      // Prepare for the loops of the product and verify compatibility of shapes
2539      int start0=0, start1=0;
2540      if (transpose == 0)       {}
2541      else if (transpose == 1)  { start0 = axis_offset; }
2542      else if (transpose == 2)  { start1 = rank1-axis_offset; }
2543      else              { throw DataException("C_GeneralTensorProduct: Error - transpose should be 0, 1 or 2"); }
2544    
2545      // Adjust the shapes for transpose
2546      DataArrayView::ShapeType tmpShape0;
2547      DataArrayView::ShapeType tmpShape1;
2548      for (int i=0; i<rank0; i++)   { tmpShape0.push_back( shape0[(i+start0)%rank0] ); }
2549      for (int i=0; i<rank1; i++)   { tmpShape1.push_back( shape1[(i+start1)%rank1] ); }
2550    
2551    #if 0
2552      // For debugging: show shape after transpose
2553      char tmp[100];
2554      std::string shapeStr;
2555      shapeStr = "(";
2556      for (int i=0; i<rank0; i++)   { sprintf(tmp, "%d,", tmpShape0[i]); shapeStr += tmp; }
2557      shapeStr += ")";
2558      cout << "C_GeneralTensorProduct: Shape of arg0 is " << shapeStr << endl;
2559      shapeStr = "(";
2560      for (int i=0; i<rank1; i++)   { sprintf(tmp, "%d,", tmpShape1[i]); shapeStr += tmp; }
2561      shapeStr += ")";
2562      cout << "C_GeneralTensorProduct: Shape of arg1 is " << shapeStr << endl;
2563    #endif
2564    
2565      // Prepare for the loops of the product
2566      int SL=1, SM=1, SR=1;
2567      for (int i=0; i<rank0-axis_offset; i++)   {
2568        SL *= tmpShape0[i];
2569      }
2570      for (int i=rank0-axis_offset; i<rank0; i++)   {
2571        if (tmpShape0[i] != tmpShape1[i-(rank0-axis_offset)]) {
2572          throw DataException("C_GeneralTensorProduct: Error - incompatible shapes");
2573        }
2574        SM *= tmpShape0[i];
2575      }
2576      for (int i=axis_offset; i<rank1; i++)     {
2577        SR *= tmpShape1[i];
2578      }
2579    
2580      // Define the shape of the output
2581      DataArrayView::ShapeType shape2;
2582      for (int i=0; i<rank0-axis_offset; i++) { shape2.push_back(tmpShape0[i]); } // First part of arg_0_Z
2583      for (int i=axis_offset; i<rank1; i++)   { shape2.push_back(tmpShape1[i]); } // Last part of arg_1_Z
2584    
2585      // Declare output Data object
2586      Data res;
2587    
2588      if      (arg_0_Z.isConstant()   && arg_1_Z.isConstant()) {
2589        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataConstant output
2590        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[0]);
2591        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[0]);
2592        double *ptr_2 = &((res.getPointDataView().getData())[0]);
2593        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2594      }
2595      else if (arg_0_Z.isConstant()   && arg_1_Z.isTagged()) {
2596    
2597        // Prepare the DataConstant input
2598        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2599        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2600    
2601        // Borrow DataTagged input from Data object
2602        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2603        if (tmp_1==0) { throw DataException("GTP_1 Programming error - casting to DataTagged."); }
2604    
2605        // Prepare a DataTagged output 2
2606        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());    // DataTagged output
2607        res.tag();
2608        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2609        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2610    
2611        // Prepare offset into DataConstant
2612        int offset_0 = tmp_0->getPointOffset(0,0);
2613        double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2614        // Get the views
2615        DataArrayView view_1 = tmp_1->getDefaultValue();
2616        DataArrayView view_2 = tmp_2->getDefaultValue();
2617        // Get the pointers to the actual data
2618        double *ptr_1 = &((view_1.getData())[0]);
2619        double *ptr_2 = &((view_2.getData())[0]);
2620        // Compute an MVP for the default
2621        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2622        // Compute an MVP for each tag
2623        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2624        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2625        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2626          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2627          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2628          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2629          double *ptr_1 = &view_1.getData(0);
2630          double *ptr_2 = &view_2.getData(0);
2631          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2632        }
2633    
2634      }
2635      else if (arg_0_Z.isConstant()   && arg_1_Z.isExpanded()) {
2636    
2637        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2638        DataConstant* tmp_0=dynamic_cast<DataConstant*>(arg_0_Z.borrowData());
2639        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2640        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2641        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2642        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2643        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2644        int sampleNo_1,dataPointNo_1;
2645        int numSamples_1 = arg_1_Z.getNumSamples();
2646        int numDataPointsPerSample_1 = arg_1_Z.getNumDataPointsPerSample();
2647        int offset_0 = tmp_0->getPointOffset(0,0);
2648        #pragma omp parallel for private(sampleNo_1,dataPointNo_1) schedule(static)
2649        for (sampleNo_1 = 0; sampleNo_1 < numSamples_1; sampleNo_1++) {
2650          for (dataPointNo_1 = 0; dataPointNo_1 < numDataPointsPerSample_1; dataPointNo_1++) {
2651            int offset_1 = tmp_1->getPointOffset(sampleNo_1,dataPointNo_1);
2652            int offset_2 = tmp_2->getPointOffset(sampleNo_1,dataPointNo_1);
2653            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2654            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2655            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2656            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2657          }
2658        }
2659    
2660      }
2661      else if (arg_0_Z.isTagged()     && arg_1_Z.isConstant()) {
2662    
2663        // Borrow DataTagged input from Data object
2664        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2665        if (tmp_0==0) { throw DataException("GTP_0 Programming error - casting to DataTagged."); }
2666    
2667        // Prepare the DataConstant input
2668        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2669        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2670    
2671        // Prepare a DataTagged output 2
2672        res = Data(0.0, shape2, arg_0_Z.getFunctionSpace());    // DataTagged output
2673        res.tag();
2674        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2675        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2676    
2677        // Prepare offset into DataConstant
2678        int offset_1 = tmp_1->getPointOffset(0,0);
2679        double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2680        // Get the views
2681        DataArrayView view_0 = tmp_0->getDefaultValue();
2682        DataArrayView view_2 = tmp_2->getDefaultValue();
2683        // Get the pointers to the actual data
2684        double *ptr_0 = &((view_0.getData())[0]);
2685        double *ptr_2 = &((view_2.getData())[0]);
2686        // Compute an MVP for the default
2687        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2688        // Compute an MVP for each tag
2689        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2690        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2691        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2692          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2693          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2694          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2695          double *ptr_0 = &view_0.getData(0);
2696          double *ptr_2 = &view_2.getData(0);
2697          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2698        }
2699    
2700      }
2701      else if (arg_0_Z.isTagged()     && arg_1_Z.isTagged()) {
2702    
2703        // Borrow DataTagged input from Data object
2704        DataTagged* tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2705        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2706    
2707        // Borrow DataTagged input from Data object
2708        DataTagged* tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2709        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2710    
2711        // Prepare a DataTagged output 2
2712        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace());
2713        res.tag();  // DataTagged output
2714        DataTagged* tmp_2=dynamic_cast<DataTagged*>(res.borrowData());
2715        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2716    
2717        // Get the views
2718        DataArrayView view_0 = tmp_0->getDefaultValue();
2719        DataArrayView view_1 = tmp_1->getDefaultValue();
2720        DataArrayView view_2 = tmp_2->getDefaultValue();
2721        // Get the pointers to the actual data
2722        double *ptr_0 = &((view_0.getData())[0]);
2723        double *ptr_1 = &((view_1.getData())[0]);
2724        double *ptr_2 = &((view_2.getData())[0]);
2725        // Compute an MVP for the default
2726        matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2727        // Merge the tags
2728        DataTagged::DataMapType::const_iterator i; // i->first is a tag, i->second is an offset into memory
2729        const DataTagged::DataMapType& lookup_0=tmp_0->getTagLookup();
2730        const DataTagged::DataMapType& lookup_1=tmp_1->getTagLookup();
2731        for (i=lookup_0.begin();i!=lookup_0.end();i++) {
2732          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue()); // use tmp_2 to get correct shape
2733        }
2734        for (i=lookup_1.begin();i!=lookup_1.end();i++) {
2735          tmp_2->addTaggedValue(i->first,tmp_2->getDefaultValue());
2736        }
2737        // Compute an MVP for each tag
2738        const DataTagged::DataMapType& lookup_2=tmp_2->getTagLookup();
2739        for (i=lookup_2.begin();i!=lookup_2.end();i++) {
2740          DataArrayView view_0 = tmp_0->getDataPointByTag(i->first);
2741          DataArrayView view_1 = tmp_1->getDataPointByTag(i->first);
2742          DataArrayView view_2 = tmp_2->getDataPointByTag(i->first);
2743          double *ptr_0 = &view_0.getData(0);
2744          double *ptr_1 = &view_1.getData(0);
2745          double *ptr_2 = &view_2.getData(0);
2746          matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2747        }
2748    
2749      }
2750      else if (arg_0_Z.isTagged()     && arg_1_Z.isExpanded()) {
2751    
2752        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2753        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2754        DataTagged*   tmp_0=dynamic_cast<DataTagged*>(arg_0_Z.borrowData());
2755        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2756        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2757        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2758        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2759        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2760        int sampleNo_0,dataPointNo_0;
2761        int numSamples_0 = arg_0_Z.getNumSamples();
2762        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2763        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2764        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2765          int offset_0 = tmp_0->getPointOffset(sampleNo_0,0); // They're all the same, so just use #0
2766          double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2767          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2768            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2769            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2770            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2771            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2772            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2773          }
2774        }
2775    
2776      }
2777      else if (arg_0_Z.isExpanded()   && arg_1_Z.isConstant()) {
2778    
2779        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2780        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2781        DataConstant* tmp_1=dynamic_cast<DataConstant*>(arg_1_Z.borrowData());
2782        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2783        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2784        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataConstant."); }
2785        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2786        int sampleNo_0,dataPointNo_0;
2787        int numSamples_0 = arg_0_Z.getNumSamples();
2788        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2789        int offset_1 = tmp_1->getPointOffset(0,0);
2790        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2791        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2792          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2793            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2794            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2795            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2796            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2797            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2798            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2799          }
2800        }
2801    
2802    
2803      }
2804      else if (arg_0_Z.isExpanded()   && arg_1_Z.isTagged()) {
2805    
2806        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2807        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2808        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2809        DataTagged*   tmp_1=dynamic_cast<DataTagged*>(arg_1_Z.borrowData());
2810        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2811        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2812        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataTagged."); }
2813        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2814        int sampleNo_0,dataPointNo_0;
2815        int numSamples_0 = arg_0_Z.getNumSamples();
2816        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2817        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2818        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2819          int offset_1 = tmp_1->getPointOffset(sampleNo_0,0);
2820          double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2821          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2822            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2823            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2824            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2825            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2826            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2827          }
2828        }
2829    
2830      }
2831      else if (arg_0_Z.isExpanded()   && arg_1_Z.isExpanded()) {
2832    
2833        // After finding a common function space above the two inputs have the same numSamples and num DPPS
2834        res = Data(0.0, shape2, arg_1_Z.getFunctionSpace(),true); // DataExpanded output
2835        DataExpanded* tmp_0=dynamic_cast<DataExpanded*>(arg_0_Z.borrowData());
2836        DataExpanded* tmp_1=dynamic_cast<DataExpanded*>(arg_1_Z.borrowData());
2837        DataExpanded* tmp_2=dynamic_cast<DataExpanded*>(res.borrowData());
2838        if (tmp_0==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2839        if (tmp_1==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2840        if (tmp_2==0) { throw DataException("GTP Programming error - casting to DataExpanded."); }
2841        int sampleNo_0,dataPointNo_0;
2842        int numSamples_0 = arg_0_Z.getNumSamples();
2843        int numDataPointsPerSample_0 = arg_0_Z.getNumDataPointsPerSample();
2844        #pragma omp parallel for private(sampleNo_0,dataPointNo_0) schedule(static)
2845        for (sampleNo_0 = 0; sampleNo_0 < numSamples_0; sampleNo_0++) {
2846          for (dataPointNo_0 = 0; dataPointNo_0 < numDataPointsPerSample_0; dataPointNo_0++) {
2847            int offset_0 = tmp_0->getPointOffset(sampleNo_0,dataPointNo_0);
2848            int offset_1 = tmp_1->getPointOffset(sampleNo_0,dataPointNo_0);
2849            int offset_2 = tmp_2->getPointOffset(sampleNo_0,dataPointNo_0);
2850            double *ptr_0 = &((arg_0_Z.getPointDataView().getData())[offset_0]);
2851            double *ptr_1 = &((arg_1_Z.getPointDataView().getData())[offset_1]);
2852            double *ptr_2 = &((res.getPointDataView().getData())[offset_2]);
2853            matrix_matrix_product(SL, SM, SR, ptr_0, ptr_1, ptr_2, transpose);
2854          }
2855        }
2856    
2857      }
2858      else {
2859        throw DataException("Error - C_GeneralTensorProduct: unknown combination of inputs");
2860      }
2861    
2862      return res;
2863    }
2864    
2865    DataAbstract*
2866    Data::borrowData() const
2867    {
2868      return m_data.get();
2869    }
2870    
2871  /* Member functions specific to the MPI implementation */  /* Member functions specific to the MPI implementation */
2872    
2873  void  void

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

  ViewVC Help
Powered by ViewVC 1.1.26