1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2008 by University of Queensland |
5 |
* Earth Systems Science Computational Center (ESSCC) |
6 |
* http://www.uq.edu.au/esscc |
7 |
* |
8 |
* Primary Business: Queensland, Australia |
9 |
* Licensed under the Open Software License version 3.0 |
10 |
* http://www.opensource.org/licenses/osl-3.0.php |
11 |
* |
12 |
*******************************************************/ |
13 |
|
14 |
|
15 |
#include "Data.h" |
16 |
#include "DataConstant.h" |
17 |
#include "DataException.h" |
18 |
#include "esysUtils/EsysAssert.h" |
19 |
|
20 |
#include <iostream> |
21 |
#include <boost/python/extract.hpp> |
22 |
#ifdef USE_NETCDF |
23 |
#include <netcdfcpp.h> |
24 |
#endif |
25 |
#ifdef PASO_MPI |
26 |
#include <mpi.h> |
27 |
#endif |
28 |
|
29 |
#include <boost/python/extract.hpp> |
30 |
#include "DataMaths.h" |
31 |
|
32 |
using namespace std; |
33 |
using namespace boost::python; |
34 |
|
35 |
namespace escript { |
36 |
|
37 |
DataConstant::DataConstant(const boost::python::numeric::array& value, |
38 |
const FunctionSpace& what) |
39 |
: DataAbstract(what,DataTypes::shapeFromNumArray(value)) |
40 |
{ |
41 |
// // extract the shape of the numarray |
42 |
// DataTypes::ShapeType tempShape; |
43 |
// for (int i=0; i < value.getrank(); i++) { |
44 |
// tempShape.push_back(extract<int>(value.getshape()[i])); |
45 |
// } |
46 |
|
47 |
// get the space for the data vector |
48 |
// int len = getNoValues(); |
49 |
// DataVector temp_data(len, 0.0, len); |
50 |
// DataArrayView temp_dataView(temp_data, tempShape); |
51 |
// temp_dataView.copy(value); |
52 |
|
53 |
m_data.copyFromNumArray(value); |
54 |
// |
55 |
|
56 |
// copy the data in the correct format |
57 |
// m_data=temp_data; |
58 |
// |
59 |
// create the view of the data |
60 |
// DataArrayView tempView(m_data,temp_dataView.getShape()); |
61 |
// setPointDataView(tempView); |
62 |
} |
63 |
|
64 |
// DataConstant::DataConstant(const DataArrayView& value, |
65 |
// const FunctionSpace& what) |
66 |
// : DataAbstract(what) |
67 |
// { |
68 |
// // |
69 |
// // copy the data in the correct format |
70 |
// m_data=value.getData(); |
71 |
// // |
72 |
// // create the view of the data |
73 |
// DataArrayView tempView(m_data,value.getShape()); |
74 |
// setPointDataView(tempView); |
75 |
// } |
76 |
|
77 |
DataConstant::DataConstant(const DataConstant& other) |
78 |
: DataAbstract(other.getFunctionSpace(),other.getShape()) |
79 |
{ // |
80 |
// copy the data in the correct format |
81 |
m_data=other.m_data; |
82 |
// |
83 |
// // create the view of the data |
84 |
// DataArrayView tempView(m_data,other.getPointDataView().getShape()); |
85 |
// setPointDataView(tempView); |
86 |
} |
87 |
|
88 |
DataConstant::DataConstant(const DataConstant& other, |
89 |
const DataTypes::RegionType& region) |
90 |
: DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region)) |
91 |
{ |
92 |
// |
93 |
// get the shape of the slice to copy from |
94 |
// DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region)); |
95 |
// |
96 |
// allocate space for this new DataConstant's data |
97 |
int len = getNoValues(); |
98 |
m_data.resize(len,0.,len); |
99 |
// |
100 |
// create a view of the data with the correct shape |
101 |
// DataArrayView tempView(m_data,shape); |
102 |
DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region); |
103 |
// |
104 |
// load the view with the data from the slice |
105 |
// tempView.copySlice(other.getPointDataView(),region_loop_range); |
106 |
DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range); |
107 |
// setPointDataView(tempView); |
108 |
} |
109 |
|
110 |
DataConstant::DataConstant(const FunctionSpace& what, |
111 |
const DataTypes::ShapeType &shape, |
112 |
const DataTypes::ValueType &data) |
113 |
: DataAbstract(what,shape) |
114 |
{ |
115 |
// |
116 |
// copy the data in the correct format |
117 |
m_data=data; |
118 |
// |
119 |
// create the view of the data |
120 |
// DataArrayView tempView(m_data,shape); |
121 |
// setPointDataView(tempView); |
122 |
} |
123 |
|
124 |
string |
125 |
DataConstant::toString() const |
126 |
{ |
127 |
return DataTypes::pointToString(m_data,getShape(),0,""); |
128 |
} |
129 |
|
130 |
|
131 |
DataAbstract* |
132 |
DataConstant::deepCopy() |
133 |
{ |
134 |
return new DataConstant(*this); |
135 |
} |
136 |
|
137 |
|
138 |
DataTypes::ValueType::size_type |
139 |
DataConstant::getPointOffset(int sampleNo, |
140 |
int dataPointNo) const |
141 |
{ |
142 |
EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)), |
143 |
"Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo); |
144 |
// |
145 |
// Whatever the coord's always return the same value as this is constant data. |
146 |
return 0; |
147 |
} |
148 |
|
149 |
DataTypes::ValueType::size_type |
150 |
DataConstant::getLength() const |
151 |
{ |
152 |
return m_data.size(); |
153 |
} |
154 |
|
155 |
// DataArrayView |
156 |
// DataConstant::getDataPoint(int sampleNo, |
157 |
// int dataPointNo) |
158 |
// { |
159 |
// EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)), |
160 |
// "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo); |
161 |
// // |
162 |
// // Whatever the coord's always return the same value as this is constant data. |
163 |
// return getPointDataView(); |
164 |
// } |
165 |
|
166 |
DataAbstract* |
167 |
DataConstant::getSlice(const DataTypes::RegionType& region) const |
168 |
{ |
169 |
return new DataConstant(*this,region); |
170 |
} |
171 |
|
172 |
void |
173 |
DataConstant::setSlice(const DataAbstract* value, |
174 |
const DataTypes::RegionType& region) |
175 |
{ |
176 |
const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value); |
177 |
if (tempDataConst==0) { |
178 |
throw DataException("Programming error - casting to DataConstant."); |
179 |
} |
180 |
// |
181 |
DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region)); |
182 |
DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region); |
183 |
// |
184 |
// check shape: |
185 |
if (getRank()!=region.size()) { |
186 |
throw DataException("Error - Invalid slice region."); |
187 |
} |
188 |
if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) { |
189 |
throw DataException (DataTypes::createShapeErrorMessage( |
190 |
"Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape())); |
191 |
} |
192 |
// getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range); |
193 |
DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range); |
194 |
} |
195 |
|
196 |
|
197 |
|
198 |
void |
199 |
DataConstant::symmetric(DataAbstract* ev) |
200 |
{ |
201 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
202 |
if (temp_ev==0) { |
203 |
throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error)."); |
204 |
} |
205 |
/* DataArrayView& thisView=getPointDataView(); |
206 |
DataArrayView& evView=ev->getPointDataView();*/ |
207 |
DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0); |
208 |
} |
209 |
|
210 |
void |
211 |
DataConstant::nonsymmetric(DataAbstract* ev) |
212 |
{ |
213 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
214 |
if (temp_ev==0) { |
215 |
throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error)."); |
216 |
} |
217 |
/* DataArrayView& thisView=getPointDataView(); |
218 |
DataArrayView& evView=ev->getPointDataView();*/ |
219 |
DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0); |
220 |
} |
221 |
|
222 |
void |
223 |
DataConstant::trace(DataAbstract* ev, int axis_offset) |
224 |
{ |
225 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
226 |
if (temp_ev==0) { |
227 |
throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error)."); |
228 |
} |
229 |
/* DataArrayView& thisView=getPointDataView(); |
230 |
DataArrayView& evView=ev->getPointDataView();*/ |
231 |
ValueType& evVec=temp_ev->getVector(); |
232 |
const ShapeType& evShape=temp_ev->getShape(); |
233 |
DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset); |
234 |
} |
235 |
|
236 |
void |
237 |
DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1) |
238 |
{ |
239 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
240 |
if (temp_ev==0) { |
241 |
throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error)."); |
242 |
} |
243 |
// DataArrayView& thisView=getPointDataView(); |
244 |
// DataArrayView& evView=ev->getPointDataView(); |
245 |
DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1); |
246 |
} |
247 |
|
248 |
void |
249 |
DataConstant::transpose(DataAbstract* ev, int axis_offset) |
250 |
{ |
251 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
252 |
if (temp_ev==0) { |
253 |
throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error)."); |
254 |
} |
255 |
/* DataArrayView& thisView=getPointDataView(); |
256 |
DataArrayView& evView=ev->getPointDataView();*/ |
257 |
DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset); |
258 |
} |
259 |
|
260 |
void |
261 |
DataConstant::eigenvalues(DataAbstract* ev) |
262 |
{ |
263 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
264 |
if (temp_ev==0) { |
265 |
throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error)."); |
266 |
} |
267 |
/* DataArrayView& thisView=getPointDataView(); |
268 |
DataArrayView& evView=ev->getPointDataView();*/ |
269 |
DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0); |
270 |
} |
271 |
void |
272 |
DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol) |
273 |
{ |
274 |
DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev); |
275 |
if (temp_ev==0) { |
276 |
throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error)."); |
277 |
} |
278 |
DataConstant* temp_V=dynamic_cast<DataConstant*>(V); |
279 |
if (temp_V==0) { |
280 |
throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error)."); |
281 |
} |
282 |
// DataArrayView thisView=getPointDataView(); |
283 |
// DataArrayView evView=ev->getPointDataView(); |
284 |
// DataArrayView VView=V->getPointDataView(); |
285 |
|
286 |
DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol); |
287 |
} |
288 |
|
289 |
void |
290 |
DataConstant::setToZero() |
291 |
{ |
292 |
DataTypes::ValueType::size_type n=m_data.size(); |
293 |
for (int i=0; i<n ;++i) m_data[i]=0.; |
294 |
} |
295 |
|
296 |
void |
297 |
DataConstant::dump(const std::string fileName) const |
298 |
{ |
299 |
#ifdef USE_NETCDF |
300 |
const NcDim* ncdims[DataTypes::maxRank]; |
301 |
NcVar* var; |
302 |
int rank = getRank(); |
303 |
int type= getFunctionSpace().getTypeCode(); |
304 |
int ndims =0; |
305 |
long dims[DataTypes::maxRank]; |
306 |
const double* d_ptr=&(m_data[0]); |
307 |
DataTypes::ShapeType shape = getShape(); |
308 |
int mpi_iam=getFunctionSpace().getDomain().getMPIRank(); |
309 |
int mpi_num=getFunctionSpace().getDomain().getMPISize(); |
310 |
#ifdef PASO_MPI |
311 |
MPI_Status status; |
312 |
#endif |
313 |
|
314 |
#ifdef PASO_MPI |
315 |
/* Serialize NetCDF I/O */ |
316 |
if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status); |
317 |
#endif |
318 |
|
319 |
// netCDF error handler |
320 |
NcError err(NcError::verbose_nonfatal); |
321 |
// Create the file. |
322 |
char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam); |
323 |
NcFile dataFile(newFileName, NcFile::Replace); |
324 |
// check if writing was successful |
325 |
if (!dataFile.is_valid()) |
326 |
throw DataException("Error - DataConstant:: opening of netCDF file for output failed."); |
327 |
if (!dataFile.add_att("type_id",0) ) |
328 |
throw DataException("Error - DataConstant:: appending data type to netCDF file failed."); |
329 |
if (!dataFile.add_att("rank",rank) ) |
330 |
throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed."); |
331 |
if (!dataFile.add_att("function_space_type",type)) |
332 |
throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed."); |
333 |
|
334 |
if (rank == 0) { |
335 |
if( ! (ncdims[0] = dataFile.add_dim("l", 1)) ) |
336 |
throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed."); |
337 |
dims[0]=1, |
338 |
ndims=1; |
339 |
} else { |
340 |
ndims=rank; |
341 |
dims[0]=shape[0]; |
342 |
if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) ) |
343 |
throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed."); |
344 |
if ( rank >1 ) { |
345 |
dims[1]=shape[1]; |
346 |
if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) ) |
347 |
throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed."); |
348 |
} |
349 |
if ( rank >2 ) { |
350 |
dims[2]=shape[2]; |
351 |
if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) ) |
352 |
throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed."); |
353 |
} |
354 |
if ( rank >3 ) { |
355 |
dims[3]=shape[3]; |
356 |
if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) ) |
357 |
throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed."); |
358 |
} |
359 |
} |
360 |
|
361 |
if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) ) |
362 |
throw DataException("Error - DataConstant:: appending variable to netCDF file failed."); |
363 |
if (! (var->put(d_ptr,dims)) ) |
364 |
throw DataException("Error - DataConstant:: copy data to netCDF buffer failed."); |
365 |
#ifdef PASO_MPI |
366 |
if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD); |
367 |
#endif |
368 |
#else |
369 |
throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager."); |
370 |
#endif |
371 |
} |
372 |
|
373 |
} // end of namespace |