escript/src/DataExpanded.cpp

// $Id$
/*
 ************************************************************
 *          Copyright 2006 by ACcESS MNRF                   *
 *                                                          *
 *              http://www.access.edu.au                    *
 *       Primary Business: Queensland, Australia            *
 *  Licensed under the Open Software License version 3.0    *
 *     http://www.opensource.org/licenses/osl-3.0.php       *
 *                                                          *
 ************************************************************
*/

#include "DataExpanded.h"
#include "DataException.h"
#include "DataConstant.h"
#include "DataTagged.h"

#include <boost/python/extract.hpp>

using namespace std;
using namespace boost::python;
using namespace boost;

namespace escript {

DataExpanded::DataExpanded(const boost::python::numeric::array& value,
                           const FunctionSpace& what)
  : DataAbstract(what)
{
  DataArrayView::ShapeType tempShape;
  //
  // extract the shape of the python numarray
  for (int i=0; i<value.getrank(); i++) {
    tempShape.push_back(extract<int>(value.getshape()[i]));
  }
  //
  // initialise the data array for this object
  initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
  //
  // copy the given value to every data point
  copy(value);
}

DataExpanded::DataExpanded(const DataExpanded& other)
  : DataAbstract(other.getFunctionSpace()),
  m_data(other.m_data)
{
  //
  // create the view for the data
  DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
  setPointDataView(temp);
}

DataExpanded::DataExpanded(const DataConstant& other)
  : DataAbstract(other.getFunctionSpace())
{
  //
  // initialise the data array for this object
  initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
  //
  // DataConstant only has one value, copy this to every data point
  copy(other.getPointDataView());
}

DataExpanded::DataExpanded(const DataTagged& other)
  : DataAbstract(other.getFunctionSpace())
{
  //
  // initialise the data array for this object
  initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
  //
  // for each data point in this object, extract and copy the corresponding data
  // value from the given DataTagged object
  int i,j;
  DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
  DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
  #pragma omp parallel for private(i,j) schedule(static)
  for (i=0;i<numRows;i++) {
    for (j=0;j<numCols;j++) {
      try {
        getPointDataView().copy(getPointOffset(i,j),
                                other.getPointDataView(),
                                other.getPointOffset(i,j));
      }
      catch (std::exception& e) {
        cout << e.what() << endl;
      }
    }
  }
}

DataExpanded::DataExpanded(const DataExpanded& other,
                           const DataArrayView::RegionType& region)
  : DataAbstract(other.getFunctionSpace())
{
  //
  // get the shape of the slice
  DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
  //
  // initialise this Data object to the shape of the slice
  initialise(shape,other.getNumSamples(),other.getNumDPPSample());
  //
  // copy the data
  DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
  DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
  DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
  int i,j;
  #pragma omp parallel for private(i,j) schedule(static)
  for (i=0;i<numRows;i++) {
    for (j=0;j<numCols;j++) {
      try {
        getPointDataView().copySlice(getPointOffset(i,j),
                                     other.getPointDataView(),
                                     other.getPointOffset(i,j),
                                     region_loop_range);
      }
      catch (std::exception& e) {
        cout << e.what() << endl;
      }
    }
  }
}

DataExpanded::DataExpanded(const DataArrayView& value,
                           const FunctionSpace& what)
  : DataAbstract(what)
{
  //
  // get the shape of the given data value
  DataArrayView::ShapeType tempShape=value.getShape();
  //
  // initialise this Data object to the shape of the given data value
  initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
  //
  // copy the given value to every data point
  copy(value);
}

DataExpanded::DataExpanded(const FunctionSpace& what,
                           const DataArrayView::ShapeType &shape,
                           const DataArrayView::ValueType &data)
  : DataAbstract(what)
{
  //
  // create the view of the data
  initialise(shape,what.getNumSamples(),what.getNumDPPSample());
  //
  // copy the data in the correct format
  m_data.getData()=data;
}

DataExpanded::~DataExpanded()
{
}

DataAbstract*
DataExpanded::getSlice(const DataArrayView::RegionType& region) const 
{
  return new DataExpanded(*this,region);
}

void
DataExpanded::setSlice(const DataAbstract* value,
                       const DataArrayView::RegionType& region) 
{
  const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
  if (tempDataExp==0) {
    throw DataException("Programming error - casting to DataExpanded.");
  }
  //
  // get shape of slice
  DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
  DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
  //
  // check shape
  if (getPointDataView().getRank()!=region.size()) {
    throw DataException("Error - Invalid slice region.");
  }
  if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
    throw DataException (value->getPointDataView().createShapeErrorMessage(
        "Error - Couldn't copy slice due to shape mismatch.",shape));
  }
  //
  // copy the data from the slice into this object
  DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
  DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
  int i, j;
  #pragma omp parallel for private(i,j) schedule(static)
  for (i=0;i<numRows;i++) {
    for (j=0;j<numCols;j++) {
      getPointDataView().copySliceFrom(getPointOffset(i,j),
                                       tempDataExp->getPointDataView(),
                                       tempDataExp->getPointOffset(i,j),
                                       region_loop_range);
    }
  }
}

void
DataExpanded::copy(const DataArrayView& value) 
{
  //
  // copy a single value to every data point in this object
  int nRows=m_data.getNumRows();
  int nCols=m_data.getNumCols();
  int i,j;
  #pragma omp parallel for private(i,j) schedule(static)
  for (i=0;i<nRows;i++) {
    for (j=0;j<nCols;j++) {
      // NOTE: An exception may be thown from this call if 
      // DOASSERT is on which of course will play
      // havoc with the omp threads. Run single threaded
      // if using DOASSERT. 
      getPointDataView().copy(m_data.index(i,j),value);
    }
  }
}

void
DataExpanded::copy(const boost::python::numeric::array& value) 
{
  //
  // first convert the numarray into a DataArray object
  DataArray temp(value);
  //
  // check the input shape matches this shape
  if (!getPointDataView().checkShape(temp.getView().getShape())) {
    throw DataException(getPointDataView().createShapeErrorMessage(
                        "Error - (DataExpanded) Cannot copy due to shape mismatch.",
                        temp.getView().getShape()));
  }
  //
  // now copy over the data
  copy(temp.getView());
}

void
DataExpanded::initialise(const DataArrayView::ShapeType& shape,
                         int noSamples,
                         int noDataPointsPerSample)
{
  //
  // resize data array to the required size
  m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
  //
  // create the data view of the data array
  DataArrayView temp(m_data.getData(),shape);
  setPointDataView(temp);
}

string
DataExpanded::toString() const
{
  stringstream temp;
  FunctionSpace fs=getFunctionSpace();
  //
  // create a temporary view as the offset will be changed
  DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
  for (int i=0;i<m_data.getNumRows();i++) {
    for (int j=0;j<m_data.getNumCols();j++) {
      tempView.setOffset(m_data.index(i,j));
      stringstream suffix;
      suffix << "( id: " << i << ", ref: " << fs.getReferenceNoFromSampleNo(i) << ", pnt: " << j << ")";
      temp << tempView.toString(suffix.str());
      if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
        temp << endl;
      }
    }
  }
  return temp.str();
}

DataArrayView::ValueType::size_type
DataExpanded::getPointOffset(int sampleNo,
                             int dataPointNo) const
{
  return m_data.index(sampleNo,dataPointNo);
}

DataArrayView
DataExpanded::getDataPoint(int sampleNo,
                           int dataPointNo)
{
  DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
  return temp;
}

DataArrayView::ValueType::size_type
DataExpanded::getLength() const
{
  return m_data.size();
}

void
DataExpanded::setRefValue(int ref,
                          const DataArray& value)
{
  //
  // Get the number of samples and data-points per sample.
  int numSamples = getNumSamples();
  int numDPPSample = getNumDPPSample();

  //
  // Determine the sample number which corresponds to this reference number.
  int sampleNo = -1;
  int tempRef = -1;
  for (int n=0; n<numSamples; n++) {
    tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
    if (tempRef == ref) {
      sampleNo = n;
      break;
    }
  }
  if (sampleNo == -1) {
    throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
  }

  for (int n=0; n<numDPPSample; n++) {
    //
    // Get *each* data-point in the sample in turn.
    DataArrayView pointView = getDataPoint(sampleNo, n);
    //
    // Assign the values in the DataArray to this data-point.
    pointView.copy(value.getView());
  }
}

void
DataExpanded::getRefValue(int ref,
                          DataArray& value)
{
  //
  // Get the number of samples and data-points per sample.
  int numSamples = getNumSamples();
  int numDPPSample = getNumDPPSample();

  //
  // Determine the sample number which corresponds to this reference number
  int sampleNo = -1;
  int tempRef = -1;
  for (int n=0; n<numSamples; n++) {
    tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
    if (tempRef == ref) {
      sampleNo = n;
      break;
    }
  }
  if (sampleNo == -1) {
    throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
  }

  //
  // Get the *first* data-point associated with this sample number.
  DataArrayView pointView = getDataPoint(sampleNo, 0);

  //
  // Load the values from this data-point into the DataArray
  value.getView().copy(pointView);
}

int
DataExpanded::archiveData(ofstream& archiveFile,
                          const DataArrayView::ValueType::size_type noValues) const
{
  return(m_data.archiveData(archiveFile, noValues));
}

int
DataExpanded::extractData(ifstream& archiveFile,
                          const DataArrayView::ValueType::size_type noValues)
{
  return(m_data.extractData(archiveFile, noValues));
}

void
DataExpanded::copyAll(const boost::python::numeric::array& value) {
  //
  // Get the number of samples and data-points per sample.
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  int dataPointRank = getPointDataView().getRank();
  ShapeType dataPointShape = getPointDataView().getShape();
  //
  // check rank:
  if (value.getrank()!=dataPointRank+1)
       throw DataException("Rank of numarray does not match Data object rank");
  if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
       throw DataException("leading dimension of numarray is too small");
  //
  int dataPoint = 0;
  for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
      DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
      if (dataPointRank==0) {
         dataPointView()=extract<double>(value[dataPoint]);
      } else if (dataPointRank==1) {
         for (int i=0; i<dataPointShape[0]; i++) {
            dataPointView(i)=extract<double>(value[dataPoint][i]);
         }
      } else if (dataPointRank==2) {
         for (int i=0; i<dataPointShape[0]; i++) {
           for (int j=0; j<dataPointShape[1]; j++) {
             dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
           }
         }
       } else 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++) {
                 dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
             }
           }
         }
       } else 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++) {
                 dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
               }
             }
           }
         }
      }
      dataPoint++;
    }
  }
}
void
DataExpanded::symmetric(DataAbstract* ev)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo));
    }
  }
}
void
DataExpanded::nonsymmetric(DataAbstract* ev)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo));
    }
  }
}
void
DataExpanded::trace(DataAbstract* ev, int axis_offset)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
    }
  }
}

void
DataExpanded::transpose(DataAbstract* ev, int axis_offset)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
    }
  }
}

void
DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
    }
  }
}
void
DataExpanded::eigenvalues(DataAbstract* ev)
{
  int sampleNo,dataPointNo;
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
                                    evView,ev->getPointOffset(sampleNo,dataPointNo));
    }
  }
}
void
DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
{
  int numSamples = getNumSamples();
  int numDataPointsPerSample = getNumDPPSample();
  int sampleNo,dataPointNo;
  DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
  }
  DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
  if (temp_V==0) {
    throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
  }
  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();
  DataArrayView& VView=V->getPointDataView();
  #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
  for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
    for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
         DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
                                                     evView,ev->getPointOffset(sampleNo,dataPointNo),
                                                     VView,V->getPointOffset(sampleNo,dataPointNo),
                                                     tol);
    }
  }
}

}  // end of namespace
1	jgs	102	// $Id$
2	jgs	82	/*
3	elspeth	615	************************************************************
4			* Copyright 2006 by ACcESS MNRF *
5			* *
6			* http://www.access.edu.au *
7			* Primary Business: Queensland, Australia *
8			* Licensed under the Open Software License version 3.0 *
9			* http://www.opensource.org/licenses/osl-3.0.php *
10			* *
11			************************************************************
12	jgs	82	*/
13
14	jgs	480	#include "DataExpanded.h"
15	jgs	474	#include "DataException.h"
16			#include "DataConstant.h"
17			#include "DataTagged.h"
18	jgs	82
19			#include <boost/python/extract.hpp>
20
21			using namespace std;
22			using namespace boost::python;
23			using namespace boost;
24
25			namespace escript {
26
27	jgs	102	DataExpanded::DataExpanded(const boost::python::numeric::array& value,
28			const FunctionSpace& what)
29			: DataAbstract(what)
30			{
31			DataArrayView::ShapeType tempShape;
32			//
33	jgs	117	// extract the shape of the python numarray
34			for (int i=0; i<value.getrank(); i++) {
35	jgs	102	tempShape.push_back(extract<int>(value.getshape()[i]));
36	jgs	82	}
37	jgs	117	//
38			// initialise the data array for this object
39	jgs	102	initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
40			//
41			// copy the given value to every data point
42			copy(value);
43			}
44	jgs	82
45	jgs	102	DataExpanded::DataExpanded(const DataExpanded& other)
46			: DataAbstract(other.getFunctionSpace()),
47			m_data(other.m_data)
48	jgs	117	{
49	jgs	102	//
50			// create the view for the data
51			DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
52			setPointDataView(temp);
53			}
54	jgs	82
55	jgs	102	DataExpanded::DataExpanded(const DataConstant& other)
56			: DataAbstract(other.getFunctionSpace())
57	jgs	117	{
58			//
59			// initialise the data array for this object
60	jgs	102	initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
61			//
62	jgs	117	// DataConstant only has one value, copy this to every data point
63	jgs	102	copy(other.getPointDataView());
64			}
65	jgs	82
66	jgs	102	DataExpanded::DataExpanded(const DataTagged& other)
67			: DataAbstract(other.getFunctionSpace())
68	jgs	117	{
69			//
70	jgs	119	// initialise the data array for this object
71	jgs	102	initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
72	jgs	117	//
73			// for each data point in this object, extract and copy the corresponding data
74			// value from the given DataTagged object
75	jgs	102	int i,j;
76			DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
77			DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
78	jgs	122	#pragma omp parallel for private(i,j) schedule(static)
79	jgs	117	for (i=0;i<numRows;i++) {
80			for (j=0;j<numCols;j++) {
81	jgs	102	try {
82	jgs	122	getPointDataView().copy(getPointOffset(i,j),
83			other.getPointDataView(),
84			other.getPointOffset(i,j));
85	jgs	82	}
86	jgs	102	catch (std::exception& e) {
87			cout << e.what() << endl;
88			}
89	jgs	82	}
90			}
91	jgs	102	}
92	jgs	82
93	jgs	102	DataExpanded::DataExpanded(const DataExpanded& other,
94			const DataArrayView::RegionType& region)
95			: DataAbstract(other.getFunctionSpace())
96			{
97			//
98			// get the shape of the slice
99			DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
100			//
101			// initialise this Data object to the shape of the slice
102			initialise(shape,other.getNumSamples(),other.getNumDPPSample());
103			//
104	jgs	117	// copy the data
105	jgs	108	DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
106	jgs	102	DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
107			DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
108			int i,j;
109	jgs	122	#pragma omp parallel for private(i,j) schedule(static)
110	jgs	117	for (i=0;i<numRows;i++) {
111			for (j=0;j<numCols;j++) {
112	jgs	102	try {
113	jgs	122	getPointDataView().copySlice(getPointOffset(i,j),
114			other.getPointDataView(),
115			other.getPointOffset(i,j),
116			region_loop_range);
117	jgs	82	}
118	jgs	102	catch (std::exception& e) {
119			cout << e.what() << endl;
120			}
121	jgs	82	}
122			}
123	jgs	102	}
124	jgs	82
125	jgs	102	DataExpanded::DataExpanded(const DataArrayView& value,
126			const FunctionSpace& what)
127			: DataAbstract(what)
128			{
129	jgs	117	//
130			// get the shape of the given data value
131	jgs	102	DataArrayView::ShapeType tempShape=value.getShape();
132	jgs	117	//
133			// initialise this Data object to the shape of the given data value
134	jgs	102	initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
135	jgs	117	//
136			// copy the given value to every data point
137	jgs	102	copy(value);
138			}
139	jgs	82
140	jgs	119	DataExpanded::DataExpanded(const FunctionSpace& what,
141			const DataArrayView::ShapeType &shape,
142			const DataArrayView::ValueType &data)
143			: DataAbstract(what)
144			{
145			//
146	jgs	160	// create the view of the data
147			initialise(shape,what.getNumSamples(),what.getNumDPPSample());
148			//
149	jgs	119	// copy the data in the correct format
150			m_data.getData()=data;
151			}
152
153	jgs	102	DataExpanded::~DataExpanded()
154			{
155			}
156	jgs	82
157	jgs	102	DataAbstract*
158			DataExpanded::getSlice(const DataArrayView::RegionType& region) const
159			{
160			return new DataExpanded(*this,region);
161			}
162
163			void
164			DataExpanded::setSlice(const DataAbstract* value,
165			const DataArrayView::RegionType& region)
166			{
167			const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
168			if (tempDataExp==0) {
169			throw DataException("Programming error - casting to DataExpanded.");
170			}
171	jgs	108	//
172	jgs	117	// get shape of slice
173	jgs	108	DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
174			DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
175			//
176	jgs	117	// check shape
177	jgs	102	if (getPointDataView().getRank()!=region.size()) {
178			throw DataException("Error - Invalid slice region.");
179			}
180	woo409	757	if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
181	jgs	102	throw DataException (value->getPointDataView().createShapeErrorMessage(
182	jgs	108	"Error - Couldn't copy slice due to shape mismatch.",shape));
183	jgs	102	}
184			//
185	jgs	117	// copy the data from the slice into this object
186	jgs	102	DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
187			DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
188			int i, j;
189	jgs	122	#pragma omp parallel for private(i,j) schedule(static)
190	jgs	102	for (i=0;i<numRows;i++) {
191			for (j=0;j<numCols;j++) {
192	jgs	122	getPointDataView().copySliceFrom(getPointOffset(i,j),
193			tempDataExp->getPointDataView(),
194			tempDataExp->getPointOffset(i,j),
195			region_loop_range);
196	jgs	82	}
197			}
198	jgs	102	}
199	jgs	82
200	jgs	102	void
201			DataExpanded::copy(const DataArrayView& value)
202			{
203			//
204	jgs	117	// copy a single value to every data point in this object
205	jgs	102	int nRows=m_data.getNumRows();
206			int nCols=m_data.getNumCols();
207	jgs	117	int i,j;
208	jgs	122	#pragma omp parallel for private(i,j) schedule(static)
209	jgs	117	for (i=0;i<nRows;i++) {
210			for (j=0;j<nCols;j++) {
211	jgs	102	// NOTE: An exception may be thown from this call if
212			// DOASSERT is on which of course will play
213			// havoc with the omp threads. Run single threaded
214			// if using DOASSERT.
215			getPointDataView().copy(m_data.index(i,j),value);
216	jgs	82	}
217			}
218	jgs	102	}
219	jgs	82
220	jgs	102	void
221			DataExpanded::copy(const boost::python::numeric::array& value)
222			{
223			//
224	jgs	117	// first convert the numarray into a DataArray object
225	jgs	102	DataArray temp(value);
226			//
227	jgs	117	// check the input shape matches this shape
228	jgs	102	if (!getPointDataView().checkShape(temp.getView().getShape())) {
229			throw DataException(getPointDataView().createShapeErrorMessage(
230			"Error - (DataExpanded) Cannot copy due to shape mismatch.",
231			temp.getView().getShape()));
232	jgs	82	}
233	jgs	102	//
234	jgs	117	// now copy over the data
235	jgs	102	copy(temp.getView());
236			}
237	jgs	82
238	jgs	102	void
239			DataExpanded::initialise(const DataArrayView::ShapeType& shape,
240			int noSamples,
241			int noDataPointsPerSample)
242			{
243			//
244	jgs	117	// resize data array to the required size
245	jgs	102	m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
246			//
247	jgs	117	// create the data view of the data array
248	jgs	102	DataArrayView temp(m_data.getData(),shape);
249			setPointDataView(temp);
250			}
251	jgs	82
252	jgs	102	string
253			DataExpanded::toString() const
254			{
255			stringstream temp;
256	gross	856	FunctionSpace fs=getFunctionSpace();
257	jgs	102	//
258			// create a temporary view as the offset will be changed
259			DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
260	jgs	117	for (int i=0;i<m_data.getNumRows();i++) {
261			for (int j=0;j<m_data.getNumCols();j++) {
262	jgs	102	tempView.setOffset(m_data.index(i,j));
263			stringstream suffix;
264	gross	856	suffix << "( id: " << i << ", ref: " << fs.getReferenceNoFromSampleNo(i) << ", pnt: " << j << ")";
265	jgs	102	temp << tempView.toString(suffix.str());
266			if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
267			temp << endl;
268	jgs	82	}
269			}
270			}
271	jgs	102	return temp.str();
272			}
273	jgs	82
274	jgs	102	DataArrayView::ValueType::size_type
275			DataExpanded::getPointOffset(int sampleNo,
276			int dataPointNo) const
277			{
278			return m_data.index(sampleNo,dataPointNo);
279			}
280	jgs	82
281	jgs	102	DataArrayView
282			DataExpanded::getDataPoint(int sampleNo,
283			int dataPointNo)
284			{
285			DataArrayView temp(m_data.getData(),getPointDataView().getShape(),m_data.index(sampleNo,dataPointNo));
286			return temp;
287			}
288	jgs	82
289	jgs	102	DataArrayView::ValueType::size_type
290			DataExpanded::getLength() const
291			{
292	jgs	117	return m_data.size();
293	jgs	102	}
294	jgs	82
295	jgs	110	void
296			DataExpanded::setRefValue(int ref,
297			const DataArray& value)
298			{
299			//
300	jgs	113	// Get the number of samples and data-points per sample.
301	jgs	110	int numSamples = getNumSamples();
302			int numDPPSample = getNumDPPSample();
303
304			//
305	jgs	113	// Determine the sample number which corresponds to this reference number.
306	jgs	110	int sampleNo = -1;
307			int tempRef = -1;
308			for (int n=0; n<numSamples; n++) {
309			tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
310			if (tempRef == ref) {
311			sampleNo = n;
312			break;
313			}
314			}
315			if (sampleNo == -1) {
316			throw DataException("DataExpanded::setRefValue error: invalid ref number supplied.");
317			}
318
319	jgs	113	for (int n=0; n<numDPPSample; n++) {
320			//
321	jgs	126	// Get each data-point in the sample in turn.
322	jgs	113	DataArrayView pointView = getDataPoint(sampleNo, n);
323			//
324			// Assign the values in the DataArray to this data-point.
325			pointView.copy(value.getView());
326			}
327	jgs	110	}
328
329			void
330			DataExpanded::getRefValue(int ref,
331			DataArray& value)
332			{
333			//
334	jgs	113	// Get the number of samples and data-points per sample.
335	jgs	110	int numSamples = getNumSamples();
336			int numDPPSample = getNumDPPSample();
337
338			//
339	jgs	113	// Determine the sample number which corresponds to this reference number
340	jgs	110	int sampleNo = -1;
341			int tempRef = -1;
342			for (int n=0; n<numSamples; n++) {
343			tempRef = getFunctionSpace().getReferenceNoFromSampleNo(n);
344			if (tempRef == ref) {
345			sampleNo = n;
346			break;
347			}
348			}
349			if (sampleNo == -1) {
350			throw DataException("DataExpanded::getRefValue error: invalid ref number supplied.");
351			}
352
353	jgs	113	//
354	jgs	126	// Get the first data-point associated with this sample number.
355	jgs	110	DataArrayView pointView = getDataPoint(sampleNo, 0);
356
357			//
358			// Load the values from this data-point into the DataArray
359			value.getView().copy(pointView);
360			}
361
362	jgs	123	int
363			DataExpanded::archiveData(ofstream& archiveFile,
364			const DataArrayView::ValueType::size_type noValues) const
365			{
366			return(m_data.archiveData(archiveFile, noValues));
367			}
368
369			int
370			DataExpanded::extractData(ifstream& archiveFile,
371			const DataArrayView::ValueType::size_type noValues)
372			{
373			return(m_data.extractData(archiveFile, noValues));
374			}
375
376	jgs	126	void
377			DataExpanded::copyAll(const boost::python::numeric::array& value) {
378			//
379			// Get the number of samples and data-points per sample.
380			int numSamples = getNumSamples();
381			int numDataPointsPerSample = getNumDPPSample();
382			int dataPointRank = getPointDataView().getRank();
383			ShapeType dataPointShape = getPointDataView().getShape();
384			//
385			// check rank:
386			if (value.getrank()!=dataPointRank+1)
387			throw DataException("Rank of numarray does not match Data object rank");
388			if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
389			throw DataException("leading dimension of numarray is too small");
390			//
391			int dataPoint = 0;
392			for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
393			for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
394			DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
395			if (dataPointRank==0) {
396			dataPointView()=extract<double>(value[dataPoint]);
397			} else if (dataPointRank==1) {
398			for (int i=0; i<dataPointShape[0]; i++) {
399			dataPointView(i)=extract<double>(value[dataPoint][i]);
400			}
401			} else if (dataPointRank==2) {
402			for (int i=0; i<dataPointShape[0]; i++) {
403			for (int j=0; j<dataPointShape[1]; j++) {
404			dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
405			}
406			}
407			} else if (dataPointRank==3) {
408			for (int i=0; i<dataPointShape[0]; i++) {
409			for (int j=0; j<dataPointShape[1]; j++) {
410			for (int k=0; k<dataPointShape[2]; k++) {
411			dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
412			}
413			}
414			}
415			} else if (dataPointRank==4) {
416			for (int i=0; i<dataPointShape[0]; i++) {
417			for (int j=0; j<dataPointShape[1]; j++) {
418			for (int k=0; k<dataPointShape[2]; k++) {
419			for (int l=0; l<dataPointShape[3]; l++) {
420			dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
421			}
422			}
423			}
424			}
425			}
426			dataPoint++;
427			}
428			}
429			}
430	gross	580	void
431	ksteube	775	DataExpanded::symmetric(DataAbstract* ev)
432			{
433			int sampleNo,dataPointNo;
434			int numSamples = getNumSamples();
435			int numDataPointsPerSample = getNumDPPSample();
436			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
437			if (temp_ev==0) {
438			throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
439			}
440			DataArrayView& thisView=getPointDataView();
441			DataArrayView& evView=ev->getPointDataView();
442			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
443			for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
444			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
445			DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
446			evView,ev->getPointOffset(sampleNo,dataPointNo));
447			}
448			}
449			}
450			void
451			DataExpanded::nonsymmetric(DataAbstract* ev)
452			{
453			int sampleNo,dataPointNo;
454			int numSamples = getNumSamples();
455			int numDataPointsPerSample = getNumDPPSample();
456			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
457			if (temp_ev==0) {
458			throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
459			}
460			DataArrayView& thisView=getPointDataView();
461			DataArrayView& evView=ev->getPointDataView();
462			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
463			for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
464			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
465			DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
466			evView,ev->getPointOffset(sampleNo,dataPointNo));
467			}
468			}
469			}
470			void
471	gross	800	DataExpanded::trace(DataAbstract* ev, int axis_offset)
472	ksteube	775	{
473			int sampleNo,dataPointNo;
474			int numSamples = getNumSamples();
475			int numDataPointsPerSample = getNumDPPSample();
476			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
477			if (temp_ev==0) {
478	gross	800	throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
479	ksteube	775	}
480			DataArrayView& thisView=getPointDataView();
481			DataArrayView& evView=ev->getPointDataView();
482			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
483			for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
484			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
485	gross	800	DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
486	ksteube	775	evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
487			}
488			}
489			}
490	gross	800
491	ksteube	775	void
492			DataExpanded::transpose(DataAbstract* ev, int axis_offset)
493			{
494			int sampleNo,dataPointNo;
495			int numSamples = getNumSamples();
496			int numDataPointsPerSample = getNumDPPSample();
497			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
498			if (temp_ev==0) {
499			throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
500			}
501			DataArrayView& thisView=getPointDataView();
502			DataArrayView& evView=ev->getPointDataView();
503			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
504			for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
505			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
506			DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
507			evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
508			}
509			}
510			}
511	gross	800
512	ksteube	775	void
513	gross	804	DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
514	gross	800	{
515			int sampleNo,dataPointNo;
516			int numSamples = getNumSamples();
517			int numDataPointsPerSample = getNumDPPSample();
518			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
519			if (temp_ev==0) {
520	gross	804	throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
521	gross	800	}
522			DataArrayView& thisView=getPointDataView();
523			DataArrayView& evView=ev->getPointDataView();
524			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
525			for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
526			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
527	gross	804	DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
528			evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
529	gross	800	}
530			}
531			}
532			void
533	gross	580	DataExpanded::eigenvalues(DataAbstract* ev)
534			{
535	gross	584	int sampleNo,dataPointNo;
536	gross	580	int numSamples = getNumSamples();
537			int numDataPointsPerSample = getNumDPPSample();
538			DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
539			if (temp_ev==0) {
540			throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
541			}
542			DataArrayView& thisView=getPointDataView();
543			DataArrayView& evView=ev->getPointDataView();
544			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
545	gross	584	for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
546			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
547	gross	580	DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
548			evView,ev->getPointOffset(sampleNo,dataPointNo));
549			}
550			}
551			}
552			void
553			DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
554			{
555			int numSamples = getNumSamples();
556			int numDataPointsPerSample = getNumDPPSample();
557	gross	584	int sampleNo,dataPointNo;
558	gross	580	DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
559			if (temp_ev==0) {
560			throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
561			}
562			DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
563			if (temp_V==0) {
564			throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
565			}
566			DataArrayView& thisView=getPointDataView();
567			DataArrayView& evView=ev->getPointDataView();
568			DataArrayView& VView=V->getPointDataView();
569			#pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
570	gross	584	for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
571			for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
572	gross	580	DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
573			evView,ev->getPointOffset(sampleNo,dataPointNo),
574			VView,V->getPointOffset(sampleNo,dataPointNo),
575			tol);
576			}
577			}
578			}
579
580	jgs	82	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision