/[escript]/trunk/escript/src/DataExpanded.cpp
ViewVC logotype

Diff of /trunk/escript/src/DataExpanded.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 967 by gross, Tue Feb 13 09:40:12 2007 UTC revision 1748 by ksteube, Wed Sep 3 06:10:39 2008 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /* $Id$ */
3   ************************************************************  
4   *          Copyright 2006 by ACcESS MNRF                   *  /*******************************************************
5   *                                                          *   *
6   *              http://www.access.edu.au                    *   *           Copyright 2003-2007 by ACceSS MNRF
7   *       Primary Business: Queensland, Australia            *   *       Copyright 2007 by University of Queensland
8   *  Licensed under the Open Software License version 3.0    *   *
9   *     http://www.opensource.org/licenses/osl-3.0.php       *   *                http://esscc.uq.edu.au
10   *                                                          *   *        Primary Business: Queensland, Australia
11   ************************************************************   *  Licensed under the Open Software License version 3.0
12  */   *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  #include "DataExpanded.h"  #include "DataExpanded.h"
17  #include "DataException.h"  #include "DataException.h"
18  #include "DataConstant.h"  #include "DataConstant.h"
19  #include "DataTagged.h"  #include "DataTagged.h"
20    #ifdef USE_NETCDF
21  #include <netcdfcpp.h>  #include <netcdfcpp.h>
22    #endif
23    #ifdef PASO_MPI
24    #include <mpi.h>
25    #endif
26    
27  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
28    
# Line 156  DataExpanded::~DataExpanded() Line 163  DataExpanded::~DataExpanded()
163  }  }
164    
165  DataAbstract*  DataAbstract*
166  DataExpanded::getSlice(const DataArrayView::RegionType& region) const  DataExpanded::getSlice(const DataArrayView::RegionType& region) const
167  {  {
168    return new DataExpanded(*this,region);    return new DataExpanded(*this,region);
169  }  }
170    
171  void  void
172  DataExpanded::setSlice(const DataAbstract* value,  DataExpanded::setSlice(const DataAbstract* value,
173                         const DataArrayView::RegionType& region)                         const DataArrayView::RegionType& region)
174  {  {
175    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);    const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
176    if (tempDataExp==0) {    if (tempDataExp==0) {
# Line 199  DataExpanded::setSlice(const DataAbstrac Line 206  DataExpanded::setSlice(const DataAbstrac
206  }  }
207    
208  void  void
209  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataArrayView& value)
210  {  {
211    //    //
212    // copy a single value to every data point in this object    // copy a single value to every data point in this object
# Line 209  DataExpanded::copy(const DataArrayView& Line 216  DataExpanded::copy(const DataArrayView&
216    #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
217    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
218      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
219        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
220        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
221        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
222        // if using DOASSERT.        // if using DOASSERT.
223        getPointDataView().copy(m_data.index(i,j),value);        getPointDataView().copy(getPointOffset(i,j),value);
224      }      }
225    }    }
226  }  }
227    
228  void  void
229  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
230  {  {
231    //  
232    // first convert the numarray into a DataArray object    // extract the shape of the numarray
233    DataArray temp(value);    DataArrayView::ShapeType tempShape;
234      for (int i=0; i < value.getrank(); i++) {
235        tempShape.push_back(extract<int>(value.getshape()[i]));
236      }
237    
238      // get the space for the data vector
239      int len = DataArrayView::noValues(tempShape);
240      DataVector temp_data(len, 0.0, len);
241      DataArrayView temp_dataView(temp_data, tempShape);
242      temp_dataView.copy(value);
243    
244    //    //
245    // check the input shape matches this shape    // check the input shape matches this shape
246    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!getPointDataView().checkShape(temp_dataView.getShape())) {
247      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(getPointDataView().createShapeErrorMessage(
248                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
249                          temp.getView().getShape()));                          temp_dataView.getShape()));
250    }    }
251    //    //
252    // now copy over the data    // now copy over the data
253    copy(temp.getView());    copy(temp_dataView);
254  }  }
255    
256  void  void
# Line 260  DataExpanded::toString() const Line 277  DataExpanded::toString() const
277    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
278    for (int i=0;i<m_data.getNumRows();i++) {    for (int i=0;i<m_data.getNumRows();i++) {
279      for (int j=0;j<m_data.getNumCols();j++) {      for (int j=0;j<m_data.getNumCols();j++) {
280        tempView.setOffset(m_data.index(i,j));        tempView.setOffset(getPointOffset(i,j));
281        stringstream suffix;        stringstream suffix;
282        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";        suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
283        temp << tempView.toString(suffix.str());        temp << tempView.toString(suffix.str());
# Line 283  DataArrayView Line 300  DataArrayView
300  DataExpanded::getDataPoint(int sampleNo,  DataExpanded::getDataPoint(int sampleNo,
301                             int dataPointNo)                             int dataPointNo)
302  {  {
303    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));    DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
304    return temp;    return temp;
305  }  }
306    
# Line 307  DataExpanded::extractData(ifstream& arch Line 324  DataExpanded::extractData(ifstream& arch
324    return(m_data.extractData(archiveFile, noValues));    return(m_data.extractData(archiveFile, noValues));
325  }  }
326    
327  void  void
328  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
329    //    //
330    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
# Line 354  DataExpanded::copyToDataPoint(const int Line 371  DataExpanded::copyToDataPoint(const int
371               }               }
372             }             }
373           }           }
374       }       }
375    }    }
376  }  }
377  void  void
378  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {  DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
379    //    //
380    // Get the number of samples and data-points per sample.    // Get the number of samples and data-points per sample.
# Line 408  DataExpanded::copyToDataPoint(const int Line 425  DataExpanded::copyToDataPoint(const int
425               }               }
426             }             }
427           }           }
428       }       }
429    }    }
430  }  }
431  void  void
# Line 616  DataExpanded::eigenvalues_and_eigenvecto Line 633  DataExpanded::eigenvalues_and_eigenvecto
633  }  }
634    
635  void  void
636    DataExpanded::setToZero(){
637      int numSamples = getNumSamples();
638      int numDataPointsPerSample = getNumDPPSample();
639      DataArrayView& thisView=getPointDataView();
640      DataArrayView::ValueType::size_type n = thisView.noValues();
641      double* p;
642      int  sampleNo,dataPointNo, i;
643      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
644      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
645        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
646            p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
647            for (i=0; i<n ;++i) p[i]=0.;
648        }
649      }
650    }
651    
652    /* Append MPI rank to file name if multiple MPI processes */
653    char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
654      /* Make plenty of room for the mpi_rank number and terminating '\0' */
655      char *newFileName = (char *)malloc(strlen(fileName)+20);
656      strncpy(newFileName, fileName, strlen(fileName)+1);
657      if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
658      return(newFileName);
659    }
660    
661    void
662  DataExpanded::dump(const std::string fileName) const  DataExpanded::dump(const std::string fileName) const
663  {  {
664     #ifdef PASO_MPI     #ifdef USE_NETCDF
    throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");  
    #endif  
665     const int ldims=2+DataArrayView::maxRank;     const int ldims=2+DataArrayView::maxRank;
666     const NcDim* ncdims[ldims];     const NcDim* ncdims[ldims];
667     NcVar *var, *ids;     NcVar *var, *ids;
# Line 628  DataExpanded::dump(const std::string fil Line 669  DataExpanded::dump(const std::string fil
669     int type=  getFunctionSpace().getTypeCode();     int type=  getFunctionSpace().getTypeCode();
670     int ndims =0;     int ndims =0;
671     long dims[ldims];     long dims[ldims];
672       const double* d_ptr=&(m_data[0]);
673     DataArrayView::ShapeType shape = getPointDataView().getShape();     DataArrayView::ShapeType shape = getPointDataView().getShape();
674       int mpi_iam=0, mpi_num=1;
675    
676     // netCDF error handler     // netCDF error handler
677     NcError err(NcError::verbose_nonfatal);     NcError err(NcError::verbose_nonfatal);
678     // Create the file.     // Create the file.
679     NcFile dataFile(fileName.c_str(), NcFile::Replace);  #ifdef PASO_MPI
680       MPI_Comm_rank(MPI_COMM_WORLD, &mpi_iam);
681       MPI_Comm_size(MPI_COMM_WORLD, &mpi_num);
682    #endif
683       char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
684       NcFile dataFile(newFileName, NcFile::Replace);
685     // check if writing was successful     // check if writing was successful
686     if (!dataFile.is_valid())     if (!dataFile.is_valid())
687          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");          throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
688     if (!dataFile.add_att("type","expanded") )     if (!dataFile.add_att("type_id",2) )
689          throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
690     if (!dataFile.add_att("rank",rank) )     if (!dataFile.add_att("rank",rank) )
691          throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
# Line 679  DataExpanded::dump(const std::string fil Line 727  DataExpanded::dump(const std::string fil
727    
728     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
729          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");          throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
730     if (! (var->put(&m_data[0],dims)) )     if (! (var->put(d_ptr,dims)) )
731          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");          throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
732       #else
733       throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
734       #endif
735    }
736    
737    void
738    DataExpanded::setTaggedValue(int tagKey,
739                                 const DataArrayView& value)
740    {
741      int numSamples = getNumSamples();
742      int numDataPointsPerSample = getNumDPPSample();
743      int sampleNo,dataPointNo, i;
744      DataArrayView& thisView=getPointDataView();
745      DataArrayView::ValueType::size_type n = thisView.noValues();
746      double* p,*in=&(value.getData()[0]);
747      
748      if (value.noValues() != n) {
749        throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
750      }
751    
752      #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
753      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
754        if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
755            for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
756                p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
757                for (i=0; i<n ;++i) p[i]=in[i];
758            }
759        }
760      }
761    }
762    
763    void
764    DataExpanded::reorderByReferenceIDs(int *reference_ids)
765    {
766      int numSamples = getNumSamples();
767      DataArrayView& thisView=getPointDataView();
768      DataArrayView::ValueType::size_type n = thisView.noValues() * getNumDPPSample();
769      int sampleNo, sampleNo2,i;
770      double* p,*p2;
771      register double rtmp;
772      FunctionSpace fs=getFunctionSpace();
773    
774      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
775         const int id_in=reference_ids[sampleNo];
776         const int id=fs.getReferenceIDOfSample(sampleNo);
777         if (id!=id_in) {
778             bool matched=false;
779             for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
780                  if (id == reference_ids[sampleNo2]) {
781                     p=&(m_data[getPointOffset(sampleNo,0)]);
782                     p2=&(m_data[getPointOffset(sampleNo2,0)]);
783                     for (i=0; i<n ;i++) {
784                             rtmp=p[i];
785                             p[i]=p2[i];
786                             p2[i]=rtmp;
787                     }
788                     reference_ids[sampleNo]=id;
789                     reference_ids[sampleNo2]=id_in;
790                     matched=true;
791                     break;
792                  }
793             }
794             if (! matched) {
795                throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
796             }
797         }
798       }
799  }  }
800    
801    

Legend:
Removed from v.967  
changed lines
  Added in v.1748

  ViewVC Help
Powered by ViewVC 1.1.26