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