/[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 876 by ksteube, Thu Oct 19 03:50:23 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 550  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 642  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 656  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 666  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 681  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 705  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 732  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 838  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 881  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 899  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 1556  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 1630  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 2149  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 2157  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 2167  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        }        }

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

  ViewVC Help
Powered by ViewVC 1.1.26