/[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 757 by woo409, Mon Jun 26 13:12:56 2006 UTC
# Line 1  Line 1 
1  // $Id$  // $Id$
2  /*  /*
3   ******************************************************************************   ************************************************************
4   *                                                                            *   *          Copyright 2006 by ACcESS MNRF                   *
5   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *   *                                                          *
6   *                                                                            *   *              http://www.access.edu.au                    *
7   * This software is the property of ACcESS. No part of this code              *   *       Primary Business: Queensland, Australia            *
8   * may be copied in any form or by any means without the expressed written    *   *  Licensed under the Open Software License version 3.0    *
9   * consent of ACcESS.  Copying, use or modification of this software          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
10   * by any unauthorised person is illegal unless that person has a software    *   *                                                          *
11   * license agreement with ACcESS.                                             *   ************************************************************
  *                                                                            *  
  ******************************************************************************  
12  */  */
13    
14  #include "escript/Data/DataException.h"  #include "DataExpanded.h"
15  #include "escript/Data/DataExpanded.h"  #include "DataException.h"
16  #include "escript/Data/DataConstant.h"  #include "DataConstant.h"
17  #include "escript/Data/DataTagged.h"  #include "DataTagged.h"
 #include "escript/Data/DataArrayView.h"  
18    
19  #include <boost/python/extract.hpp>  #include <boost/python/extract.hpp>
20    
 #include <iostream>  
   
21  using namespace std;  using namespace std;
22  using namespace boost::python;  using namespace boost::python;
23  using namespace boost;  using namespace boost;
# Line 72  DataExpanded::DataExpanded(const DataTag Line 67  DataExpanded::DataExpanded(const DataTag
67    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
68  {  {
69    //    //
70    // initialise the data for this object    // initialise the data array for this object
71    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
72    //    //
73    // 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 75  DataExpanded::DataExpanded(const DataTag
75    int i,j;    int i,j;
76    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
77    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
78  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
79    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
80      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
81        try {        try {
82          getPointDataView().copy(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j));          getPointDataView().copy(getPointOffset(i,j),
83                                    other.getPointDataView(),
84                                    other.getPointOffset(i,j));
85        }        }
86        catch (std::exception& e) {        catch (std::exception& e) {
87          cout << e.what() << endl;          cout << e.what() << endl;
# Line 109  DataExpanded::DataExpanded(const DataExp Line 106  DataExpanded::DataExpanded(const DataExp
106    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
107    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
108    int i,j;    int i,j;
109  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
110    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
111      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
112        try {        try {
113          getPointDataView().copySlice(getPointOffset(i,j),other.getPointDataView(),other.getPointOffset(i,j),region_loop_range);          getPointDataView().copySlice(getPointOffset(i,j),
114                                         other.getPointDataView(),
115                                         other.getPointOffset(i,j),
116                                         region_loop_range);
117        }        }
118        catch (std::exception& e) {        catch (std::exception& e) {
119          cout << e.what() << endl;          cout << e.what() << endl;
# Line 137  DataExpanded::DataExpanded(const DataArr Line 137  DataExpanded::DataExpanded(const DataArr
137    copy(value);    copy(value);
138  }  }
139    
140    DataExpanded::DataExpanded(const FunctionSpace& what,
141                               const DataArrayView::ShapeType &shape,
142                               const DataArrayView::ValueType &data)
143      : DataAbstract(what)
144    {
145      //
146      // create the view of the data
147      initialise(shape,what.getNumSamples(),what.getNumDPPSample());
148      //
149      // copy the data in the correct format
150      m_data.getData()=data;
151    }
152    
153  DataExpanded::~DataExpanded()  DataExpanded::~DataExpanded()
154  {  {
155  }  }
# Line 160  DataExpanded::reshapeDataPoint(const Dat Line 173  DataExpanded::reshapeDataPoint(const Dat
173    int i,j;    int i,j;
174    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
175    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
176  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
177    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
178      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
179        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
# Line 200  DataExpanded::setSlice(const DataAbstrac Line 213  DataExpanded::setSlice(const DataAbstrac
213    if (getPointDataView().getRank()!=region.size()) {    if (getPointDataView().getRank()!=region.size()) {
214      throw DataException("Error - Invalid slice region.");      throw DataException("Error - Invalid slice region.");
215    }    }
216    if (tempDataExp->getPointDataView().getRank()>0 and !value->getPointDataView().checkShape(shape)) {    if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
217      throw DataException (value->getPointDataView().createShapeErrorMessage(      throw DataException (value->getPointDataView().createShapeErrorMessage(
218          "Error - Couldn't copy slice due to shape mismatch.",shape));          "Error - Couldn't copy slice due to shape mismatch.",shape));
219    }    }
# Line 209  DataExpanded::setSlice(const DataAbstrac Line 222  DataExpanded::setSlice(const DataAbstrac
222    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();    DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
223    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();    DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
224    int i, j;    int i, j;
225  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
226    for (i=0;i<numRows;i++) {    for (i=0;i<numRows;i++) {
227      for (j=0;j<numCols;j++) {      for (j=0;j<numCols;j++) {
228        getPointDataView().copySliceFrom(getPointOffset(i,j),tempDataExp->getPointDataView(),tempDataExp->getPointOffset(i,j),region_loop_range);        getPointDataView().copySliceFrom(getPointOffset(i,j),
229                                           tempDataExp->getPointDataView(),
230                                           tempDataExp->getPointOffset(i,j),
231                                           region_loop_range);
232      }      }
233    }    }
234  }  }
# Line 225  DataExpanded::copy(const DataArrayView& Line 241  DataExpanded::copy(const DataArrayView&
241    int nRows=m_data.getNumRows();    int nRows=m_data.getNumRows();
242    int nCols=m_data.getNumCols();    int nCols=m_data.getNumCols();
243    int i,j;    int i,j;
244  #pragma omp parallel for private(i,j) schedule(static)    #pragma omp parallel for private(i,j) schedule(static)
245    for (i=0;i<nRows;i++) {    for (i=0;i<nRows;i++) {
246      for (j=0;j<nCols;j++) {      for (j=0;j<nCols;j++) {
247        // 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 353  DataExpanded::setRefValue(int ref,
353    
354    for (int n=0; n<numDPPSample; n++) {    for (int n=0; n<numDPPSample; n++) {
355      //      //
356      // Get each data-point in the sample in turn.      // Get *each* data-point in the sample in turn.
357      DataArrayView pointView = getDataPoint(sampleNo, n);      DataArrayView pointView = getDataPoint(sampleNo, n);
358      //      //
359      // 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 386  DataExpanded::getRefValue(int ref,
386    }    }
387    
388    //    //
389    // Get the first data-point associated with this sample number.    // Get the *first* data-point associated with this sample number.
390    DataArrayView pointView = getDataPoint(sampleNo, 0);    DataArrayView pointView = getDataPoint(sampleNo, 0);
391    
392    //    //
# Line 378  DataExpanded::getRefValue(int ref, Line 394  DataExpanded::getRefValue(int ref,
394    value.getView().copy(pointView);    value.getView().copy(pointView);
395  }  }
396    
397    int
398    DataExpanded::archiveData(ofstream& archiveFile,
399                              const DataArrayView::ValueType::size_type noValues) const
400    {
401      return(m_data.archiveData(archiveFile, noValues));
402    }
403    
404    int
405    DataExpanded::extractData(ifstream& archiveFile,
406                              const DataArrayView::ValueType::size_type noValues)
407    {
408      return(m_data.extractData(archiveFile, noValues));
409    }
410    
411    void
412    DataExpanded::copyAll(const boost::python::numeric::array& value) {
413      //
414      // Get the number of samples and data-points per sample.
415      int numSamples = getNumSamples();
416      int numDataPointsPerSample = getNumDPPSample();
417      int dataPointRank = getPointDataView().getRank();
418      ShapeType dataPointShape = getPointDataView().getShape();
419      //
420      // check rank:
421      if (value.getrank()!=dataPointRank+1)
422           throw DataException("Rank of numarray does not match Data object rank");
423      if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
424           throw DataException("leading dimension of numarray is too small");
425      //
426      int dataPoint = 0;
427      for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
428        for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
429          DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
430          if (dataPointRank==0) {
431             dataPointView()=extract<double>(value[dataPoint]);
432          } else if (dataPointRank==1) {
433             for (int i=0; i<dataPointShape[0]; i++) {
434                dataPointView(i)=extract<double>(value[dataPoint][i]);
435             }
436          } else if (dataPointRank==2) {
437             for (int i=0; i<dataPointShape[0]; i++) {
438               for (int j=0; j<dataPointShape[1]; j++) {
439                 dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
440               }
441             }
442           } else if (dataPointRank==3) {
443             for (int i=0; i<dataPointShape[0]; i++) {
444               for (int j=0; j<dataPointShape[1]; j++) {
445                 for (int k=0; k<dataPointShape[2]; k++) {
446                     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
447                 }
448               }
449             }
450           } else if (dataPointRank==4) {
451             for (int i=0; i<dataPointShape[0]; i++) {
452               for (int j=0; j<dataPointShape[1]; j++) {
453                 for (int k=0; k<dataPointShape[2]; k++) {
454                   for (int l=0; l<dataPointShape[3]; l++) {
455                     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
456                   }
457                 }
458               }
459             }
460          }
461          dataPoint++;
462        }
463      }
464    }
465    void
466    DataExpanded::eigenvalues(DataAbstract* ev)
467    {
468      int sampleNo,dataPointNo;
469      int numSamples = getNumSamples();
470      int numDataPointsPerSample = getNumDPPSample();
471      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
472      if (temp_ev==0) {
473        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
474      }
475      DataArrayView& thisView=getPointDataView();
476      DataArrayView& evView=ev->getPointDataView();
477      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
478      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
479        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
480             DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
481                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
482        }
483      }
484    }
485    void
486    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
487    {
488      int numSamples = getNumSamples();
489      int numDataPointsPerSample = getNumDPPSample();
490      int sampleNo,dataPointNo;
491      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492      if (temp_ev==0) {
493        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
494      }
495      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
496      if (temp_V==0) {
497        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
498      }
499      DataArrayView& thisView=getPointDataView();
500      DataArrayView& evView=ev->getPointDataView();
501      DataArrayView& VView=V->getPointDataView();
502      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
503      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
504        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
505             DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
506                                                         evView,ev->getPointOffset(sampleNo,dataPointNo),
507                                                         VView,V->getPointOffset(sampleNo,dataPointNo),
508                                                         tol);
509        }
510      }
511    }
512    
513  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26