/[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 113 by jgs, Mon Feb 28 07:06:33 2005 UTC trunk/escript/src/DataExpanded.cpp revision 775 by ksteube, Mon Jul 10 04:00:08 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 35  DataExpanded::DataExpanded(const boost:: Line 30  DataExpanded::DataExpanded(const boost::
30  {  {
31    DataArrayView::ShapeType tempShape;    DataArrayView::ShapeType tempShape;
32    //    //
33    // extract the shape from the python array    // extract the shape of the python numarray
34    for (int i=0; i<value.getrank(); ++i) {    for (int i=0; i<value.getrank(); i++) {
     //cout << extract<int>(value.getshape()[i]) << endl;  
35      tempShape.push_back(extract<int>(value.getshape()[i]));      tempShape.push_back(extract<int>(value.getshape()[i]));
36    }    }
37      //
38      // initialise the data array for this object
39    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
40    //    //
41    // copy the given value to every data point    // copy the given value to every data point
# Line 49  DataExpanded::DataExpanded(const boost:: Line 45  DataExpanded::DataExpanded(const boost::
45  DataExpanded::DataExpanded(const DataExpanded& other)  DataExpanded::DataExpanded(const DataExpanded& other)
46    : DataAbstract(other.getFunctionSpace()),    : DataAbstract(other.getFunctionSpace()),
47    m_data(other.m_data)    m_data(other.m_data)
48  {      {
49    //    //
50    // create the view for the data    // create the view for the data
51    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());    DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
# Line 58  DataExpanded::DataExpanded(const DataExp Line 54  DataExpanded::DataExpanded(const DataExp
54    
55  DataExpanded::DataExpanded(const DataConstant& other)  DataExpanded::DataExpanded(const DataConstant& other)
56    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
57  {      {
58      //
59      // initialise the data array for this object
60    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());    initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
61    //    //
62    // DataConstant only has one value    // DataConstant only has one value, copy this to every data point
63    copy(other.getPointDataView());    copy(other.getPointDataView());
64  }  }
65    
66  DataExpanded::DataExpanded(const DataTagged& other)  DataExpanded::DataExpanded(const DataTagged& other)
67    : DataAbstract(other.getFunctionSpace())    : DataAbstract(other.getFunctionSpace())
68  {      {
69      //
70      // 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
74      // value from the given DataTagged object
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 96  DataExpanded::DataExpanded(const DataExp Line 101  DataExpanded::DataExpanded(const DataExp
101    // initialise this Data object to the shape of the slice    // initialise this Data object to the shape of the slice
102    initialise(shape,other.getNumSamples(),other.getNumDPPSample());    initialise(shape,other.getNumSamples(),other.getNumDPPSample());
103    //    //
   DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);  
   //  
104    // copy the data    // copy the data
105      DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
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 119  DataExpanded::DataExpanded(const DataArr Line 126  DataExpanded::DataExpanded(const DataArr
126                             const FunctionSpace& what)                             const FunctionSpace& what)
127    : DataAbstract(what)    : DataAbstract(what)
128  {  {
129      //
130      // get the shape of the given data value
131    DataArrayView::ShapeType tempShape=value.getShape();    DataArrayView::ShapeType tempShape=value.getShape();
132      //
133      // initialise this Data object to the shape of the given data value
134    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());    initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
135      //
136      // copy the given value to every data point
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 131  DataExpanded::~DataExpanded() Line 157  DataExpanded::~DataExpanded()
157  void  void
158  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)  DataExpanded::reshapeDataPoint(const DataArrayView::ShapeType& shape)
159  {  {
   //  
   // reshape a rank zero data point  
160    if (getPointDataView().getRank()!=0) {    if (getPointDataView().getRank()!=0) {
161      stringstream temp;      stringstream temp;
162      temp << "Error - Can only reshape Data with data points of rank 0. "      temp << "Error - Can only reshape Data with data points of rank 0. "
# Line 140  DataExpanded::reshapeDataPoint(const Dat Line 164  DataExpanded::reshapeDataPoint(const Dat
164           << getPointDataView().getRank();           << getPointDataView().getRank();
165      throw DataException(temp.str());      throw DataException(temp.str());
166    }    }
167    DataBlocks2D newData(getNumSamples(),getNumDPPSample(), DataArrayView::noValues(shape));    //
168      // create the new DataBlocks2D data array, and a corresponding DataArrayView
169      DataBlocks2D newData(getNumSamples(),getNumDPPSample(),DataArrayView::noValues(shape));
170    DataArrayView newView(newData.getData(),shape);    DataArrayView newView(newData.getData(),shape);
171    //    //
172    // Copy the original data to every value for the new shape    // Copy the original data to every value for the new shape
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++) {
       //  
       // Copy the data into the specified offset  
179        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
180        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
181        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
182        // if using DOASSERT.        // if using DOASSERT.
183        newView.copy(newData.index(i,j),m_data.getData()[m_data.index(i,j)]);        newView.copy(newData.index(i,j),m_data(i,j));
184      }      }
185    }    }
186      // swap the new data array for the original
187    m_data.Swap(newData);    m_data.Swap(newData);
188      // set the corresponding DataArrayView
189    DataArrayView temp(m_data.getData(),shape);    DataArrayView temp(m_data.getData(),shape);
190    setPointDataView(temp);    setPointDataView(temp);
191  }  }
# Line 179  DataExpanded::setSlice(const DataAbstrac Line 205  DataExpanded::setSlice(const DataAbstrac
205      throw DataException("Programming error - casting to DataExpanded.");      throw DataException("Programming error - casting to DataExpanded.");
206    }    }
207    //    //
208      // get shape of slice
209    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));    DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
210    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);    DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
211    //    //
212    // check shape:    // check shape
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    }    }
220    //    //
221    // copy the data    // copy the data from the slice into this object
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 207  void Line 237  void
237  DataExpanded::copy(const DataArrayView& value)  DataExpanded::copy(const DataArrayView& value)
238  {  {
239    //    //
240    // Copy a single value to every data point    // copy a single value to every data point in this object
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++) {
       //  
       // Copy the data into the specified offset  
247        // NOTE: An exception may be thown from this call if        // NOTE: An exception may be thown from this call if
248        // DOASSERT is on which of course will play        // DOASSERT is on which of course will play
249        // havoc with the omp threads. Run single threaded        // havoc with the omp threads. Run single threaded
# Line 229  void Line 257  void
257  DataExpanded::copy(const boost::python::numeric::array& value)  DataExpanded::copy(const boost::python::numeric::array& value)
258  {  {
259    //    //
260    // first convert the numarray into a DataArrayView format    // first convert the numarray into a DataArray object
261    DataArray temp(value);    DataArray temp(value);
262    //    //
263    // check the input shape matches this shape, this will throw an exception    // check the input shape matches this shape
264    if (!getPointDataView().checkShape(temp.getView().getShape())) {    if (!getPointDataView().checkShape(temp.getView().getShape())) {
265      throw DataException(getPointDataView().createShapeErrorMessage(      throw DataException(getPointDataView().createShapeErrorMessage(
266                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",                          "Error - (DataExpanded) Cannot copy due to shape mismatch.",
267                          temp.getView().getShape()));                          temp.getView().getShape()));
268    }    }
269    //    //
270    // now copy over the entire data structure    // now copy over the data
271    copy(temp.getView());    copy(temp.getView());
272  }  }
273    
# Line 249  DataExpanded::initialise(const DataArray Line 277  DataExpanded::initialise(const DataArray
277                           int noDataPointsPerSample)                           int noDataPointsPerSample)
278  {  {
279    //    //
280    // resize data to the required size    // resize data array to the required size
281    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));    m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
282    //    //
283    // create a point data viewer of the data    // create the data view of the data array
284    DataArrayView temp(m_data.getData(),shape);    DataArrayView temp(m_data.getData(),shape);
285    setPointDataView(temp);    setPointDataView(temp);
286  }  }
# Line 264  DataExpanded::toString() const Line 292  DataExpanded::toString() const
292    //    //
293    // create a temporary view as the offset will be changed    // create a temporary view as the offset will be changed
294    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());    DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
295    for (int i=0;i<m_data.getNumRows();++i) {    for (int i=0;i<m_data.getNumRows();i++) {
296      for (int j=0;j<m_data.getNumCols();++j) {      for (int j=0;j<m_data.getNumCols();j++) {
297        tempView.setOffset(m_data.index(i,j));        tempView.setOffset(m_data.index(i,j));
298        stringstream suffix;        stringstream suffix;
299        suffix << "(" << i << "," << j << ")";        suffix << "(" << i << "," << j << ")";
# Line 296  DataExpanded::getDataPoint(int sampleNo, Line 324  DataExpanded::getDataPoint(int sampleNo,
324  DataArrayView::ValueType::size_type  DataArrayView::ValueType::size_type
325  DataExpanded::getLength() const  DataExpanded::getLength() const
326  {  {
327    return m_data.getData().size();    return m_data.size();
328  }  }
329    
330  void  void
# Line 325  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.
360      pointView.copy(value.getView());      pointView.copy(value.getView());
# Line 359  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 367  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::symmetric(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::symmetric: 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::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
481                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
482        }
483      }
484    }
485    void
486    DataExpanded::nonsymmetric(DataAbstract* ev)
487    {
488      int sampleNo,dataPointNo;
489      int numSamples = getNumSamples();
490      int numDataPointsPerSample = getNumDPPSample();
491      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
492      if (temp_ev==0) {
493        throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
494      }
495      DataArrayView& thisView=getPointDataView();
496      DataArrayView& evView=ev->getPointDataView();
497      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
498      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
499        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
500             DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
501                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
502        }
503      }
504    }
505    void
506    DataExpanded::matrixtrace(DataAbstract* ev, int axis_offset)
507    {
508      int sampleNo,dataPointNo;
509      int numSamples = getNumSamples();
510      int numDataPointsPerSample = getNumDPPSample();
511      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
512      if (temp_ev==0) {
513        throw DataException("Error - DataExpanded::matrixtrace: casting to DataExpanded failed (propably a programming error).");
514      }
515      DataArrayView& thisView=getPointDataView();
516      DataArrayView& evView=ev->getPointDataView();
517      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
518      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
519        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
520             DataArrayView::matrixtrace(thisView,getPointOffset(sampleNo,dataPointNo),
521                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
522        }
523      }
524    }
525    void
526    DataExpanded::transpose(DataAbstract* ev, int axis_offset)
527    {
528      int sampleNo,dataPointNo;
529      int numSamples = getNumSamples();
530      int numDataPointsPerSample = getNumDPPSample();
531      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
532      if (temp_ev==0) {
533        throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
534      }
535      DataArrayView& thisView=getPointDataView();
536      DataArrayView& evView=ev->getPointDataView();
537      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
538      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
539        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
540             DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
541                                        evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
542        }
543      }
544    }
545    void
546    DataExpanded::eigenvalues(DataAbstract* ev)
547    {
548      int sampleNo,dataPointNo;
549      int numSamples = getNumSamples();
550      int numDataPointsPerSample = getNumDPPSample();
551      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
552      if (temp_ev==0) {
553        throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
554      }
555      DataArrayView& thisView=getPointDataView();
556      DataArrayView& evView=ev->getPointDataView();
557      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
558      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
559        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
560             DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
561                                        evView,ev->getPointOffset(sampleNo,dataPointNo));
562        }
563      }
564    }
565    void
566    DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
567    {
568      int numSamples = getNumSamples();
569      int numDataPointsPerSample = getNumDPPSample();
570      int sampleNo,dataPointNo;
571      DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
572      if (temp_ev==0) {
573        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
574      }
575      DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
576      if (temp_V==0) {
577        throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
578      }
579      DataArrayView& thisView=getPointDataView();
580      DataArrayView& evView=ev->getPointDataView();
581      DataArrayView& VView=V->getPointDataView();
582      #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
583      for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
584        for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
585             DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
586                                                         evView,ev->getPointOffset(sampleNo,dataPointNo),
587                                                         VView,V->getPointOffset(sampleNo,dataPointNo),
588                                                         tol);
589        }
590      }
591    }
592    
593  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.113  
changed lines
  Added in v.775

  ViewVC Help
Powered by ViewVC 1.1.26