/[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

trunk/esys2/escript/src/Data/DataExpanded.cpp revision 117 by jgs, Fri Apr 1 05:48:57 2005 UTC trunk/escript/src/DataExpanded.cpp revision 584 by gross, Thu Mar 9 23:03:38 2006 UTC
# Line 13  Line 13 
13   ******************************************************************************   ******************************************************************************
14  */  */
15    
16  #include "escript/Data/DataException.h"  #include "DataExpanded.h"
17  #include "escript/Data/DataExpanded.h"  #include "DataException.h"
18  #include "escript/Data/DataConstant.h"  #include "DataConstant.h"
19  #include "escript/Data/DataTagged.h"  #include "DataTagged.h"
 #include "escript/Data/DataArrayView.h"  
20    
21  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
22    
 #include <iostream>  
   
23  using namespace std;  using namespace std;
24  using namespace boost::python;  using namespace boost::python;
25  using namespace boost;  using namespace boost;
# Line 72  DataExpanded::DataExpanded(const DataTag Line 69  DataExpanded::DataExpanded(const DataTag
69    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
70  {  {
71    //    //
72    // initialise the data for this object    // initialise the data array for this object
73    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
74    //    //
75    // for each data point in this object, extract and copy the corresponding data    // for each data point in this object, extract and copy the corresponding data
# Line 80  DataExpanded::DataExpanded(const DataTag Line 77  DataExpanded::DataExpanded(const DataTag
77    int i,j;    int i,j;
78    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
79    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
80  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
81    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
82      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
83        try {        try {
84          getPointDataView().copy(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j));          getPointDataView().copy(getPointOffset(i,j),
85                                    other.getPointDataView(),
86                                    other.getPointOffset(i,j));
87        }        }
88        catch (std::exception& e) {        catch (std::exception& e) {
89          cout << e.what() << endl;          cout << e.what() << endl;
# Line 109  DataExpanded::DataExpanded(const DataExp Line 108  DataExpanded::DataExpanded(const DataExp
108    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
109    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
110    int i,j;    int i,j;
111  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
112    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
113      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
114        try {        try {
115          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region_loop_range);          getPointDataView().copySlice(getPointOffset(i,j),
116                                         other.getPointDataView(),
117                                         other.getPointOffset(i,j),
118                                         region_loop_range);
119        }        }
120        catch (std::exception& e) {        catch (std::exception& e) {
121          cout << e.what() << endl;          cout << e.what() << endl;
# Line 137  DataExpanded::DataExpanded(const DataArr Line 139  DataExpanded::DataExpanded(const DataArr
139    copy(value);    copy(value);
140  }  }
141    
142    DataExpanded::DataExpanded(const FunctionSpace& what,
143                               const DataArrayView::ShapeType &shape,
144                               const DataArrayView::ValueType &data)
145      : DataAbstract(what)
146    {
147      //
148      // create the view of the data
149      initialise(shape,what.getNumSamples(),what.getNumDPPSample());
150      //
151      // copy the data in the correct format
152      m_data.getData()=data;
153    }
154    
155  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
156  {  {
157  }  }
# Line 160  DataExpanded::reshapeDataPoint(const Dat Line 175  DataExpanded::reshapeDataPoint(const Dat
175    int i,j;    int i,j;
176    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
177    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
178  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
179    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
180      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
181        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
# Line 209  DataExpanded::setSlice(const DataAbstrac Line 224  DataExpanded::setSlice(const DataAbstrac
224    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
225    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
226    int i, j;    int i, j;
227  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
228    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
229      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
230        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region_loop_range);        getPointDataView().copySliceFrom(getPointOffset(i,j),
231                                           tempDataExp->getPointDataView(),
232                                           tempDataExp->getPointOffset(i,j),
233                                           region_loop_range);
234      }      }
235    }    }
236  }  }
# Line 225  DataExpanded::copy(const DataArrayView& Line 243  DataExpanded::copy(const DataArrayView&
243    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
244    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
245    int i,j;    int i,j;
246  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
247    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
248      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
249        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
# Line 337  DataExpanded::setRefValue(int ref, Line 355  DataExpanded::setRefValue(int ref,
355    
356    for (int n=0; n<numDPPSample; n++) {    for (int n=0; n<numDPPSample; n++) {
357      //      //
358      // Get each data-point in the sample in turn.      // Get *each* data-point in the sample in turn.
359      DataArrayView pointView = getDataPoint(sampleNo, n);      DataArrayView pointView = getDataPoint(sampleNo, n);
360      //      //
361      // Assign the values in the DataArray to this data-point.      // Assign the values in the DataArray to this data-point.
# Line 370  DataExpanded::getRefValue(int ref, Line 388  DataExpanded::getRefValue(int ref,
388    }    }
389    
390    //    //
391    // Get the first data-point associated with this sample number.    // Get the *first* data-point associated with this sample number.
392    DataArrayView pointView = getDataPoint(sampleNo, 0);    DataArrayView pointView = getDataPoint(sampleNo, 0);
393    
394    //    //
# Line 378  DataExpanded::getRefValue(int ref, Line 396  DataExpanded::getRefValue(int ref,
396    value.getView().copy(pointView);    value.getView().copy(pointView);
397  }  }
398    
399    int
400    DataExpanded::archiveData(ofstream& archiveFile,
401                              const DataArrayView::ValueType::size_type noValues) const
402    {
403      return(m_data.archiveData(archiveFile, noValues));
404    }
405    
406    int
407    DataExpanded::extractData(ifstream& archiveFile,
408                              const DataArrayView::ValueType::size_type noValues)
409    {
410      return(m_data.extractData(archiveFile, noValues));
411    }
412    
413    void
414    DataExpanded::copyAll(const boost::python::numeric::array& value) {
415      //
416      // Get the number of samples and data-points per sample.
417      int numSamples = getNumSamples();
418      int numDataPointsPerSample = getNumDPPSample();
419      int dataPointRank = getPointDataView().getRank();
420      ShapeType dataPointShape = getPointDataView().getShape();
421      //
422      // check rank:
423      if (value.getrank()!=dataPointRank+1)
424           throw DataException("Rank of numarray does not match Data object rank");
425      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
426           throw DataException("leading dimension of numarray is too small");
427      //
428      int dataPoint = 0;
429      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
430        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
431          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
432          if (dataPointRank==0) {
433             dataPointView()=extract<double>(value[dataPoint]);
434          } else if (dataPointRank==1) {
435             for (int i=0; i<dataPointShape[0]; i++) {
436                dataPointView(i)=extract<double>(value[dataPoint][i]);
437             }
438          } else if (dataPointRank==2) {
439             for (int i=0; i<dataPointShape[0]; i++) {
440               for (int j=0; j<dataPointShape[1]; j++) {
441                 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
442               }
443             }
444           } else if (dataPointRank==3) {
445             for (int i=0; i<dataPointShape[0]; i++) {
446               for (int j=0; j<dataPointShape[1]; j++) {
447                 for (int k=0; k<dataPointShape[2]; k++) {
448                     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
449                 }
450               }
451             }
452           } else if (dataPointRank==4) {
453             for (int i=0; i<dataPointShape[0]; i++) {
454               for (int j=0; j<dataPointShape[1]; j++) {
455                 for (int k=0; k<dataPointShape[2]; k++) {
456                   for (int l=0; l<dataPointShape[3]; l++) {
457                     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
458                   }
459                 }
460               }
461             }
462          }
463          dataPoint++;
464        }
465      }
466    }
467    void
468    DataExpanded::eigenvalues(DataAbstract* ev)
469    {
470      int sampleNo,dataPointNo;
471      int numSamples = getNumSamples();
472      int numDataPointsPerSample = getNumDPPSample();
473      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
474      if (temp_ev==0) {
475        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
476      }
477      DataArrayView& thisView=getPointDataView();
478      DataArrayView& evView=ev->getPointDataView();
479      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
480      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
481        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
482             DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
483                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
484        }
485      }
486    }
487    void
488    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
489    {
490      int numSamples = getNumSamples();
491      int numDataPointsPerSample = getNumDPPSample();
492      int sampleNo,dataPointNo;
493      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
494      if (temp_ev==0) {
495        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
496      }
497      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
498      if (temp_V==0) {
499        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
500      }
501      DataArrayView& thisView=getPointDataView();
502      DataArrayView& evView=ev->getPointDataView();
503      DataArrayView& VView=V->getPointDataView();
504      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
505      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
506        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
507             DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
508                                                         evView,ev->getPointOffset(sampleNo,dataPointNo),
509                                                         VView,V->getPointOffset(sampleNo,dataPointNo),
510                                                         tol);
511        }
512      }
513    }
514    
515  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.117  
changed lines
  Added in v.584

  ViewVC Help
Powered by ViewVC 1.1.26