/[escript]/branches/arrexp_trunk2098/escript/src/DataConstant.cpp
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/src/DataConstant.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2119 - (show annotations)
Tue Dec 2 06:06:04 2008 UTC (11 years, 10 months ago) by jfenwick
File size: 13931 byte(s)
Branch commit.
Threading Wrapped Array through the code.
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 : parent(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 WrappedArray& value,
65 const FunctionSpace& what)
66 : parent(what,value.getShape())
67 {
68 // // extract the shape of the numarray
69 // DataTypes::ShapeType tempShape;
70 // for (int i=0; i < value.getrank(); i++) {
71 // tempShape.push_back(extract<int>(value.getshape()[i]));
72 // }
73
74 // get the space for the data vector
75 // int len = getNoValues();
76 // DataVector temp_data(len, 0.0, len);
77 // DataArrayView temp_dataView(temp_data, tempShape);
78 // temp_dataView.copy(value);
79
80 m_data.copyFromArray(value);
81 //
82
83 // copy the data in the correct format
84 // m_data=temp_data;
85 //
86 // create the view of the data
87 // DataArrayView tempView(m_data,temp_dataView.getShape());
88 // setPointDataView(tempView);
89 }
90
91
92 // DataConstant::DataConstant(const DataArrayView& value,
93 // const FunctionSpace& what)
94 // : DataAbstract(what)
95 // {
96 // //
97 // // copy the data in the correct format
98 // m_data=value.getData();
99 // //
100 // // create the view of the data
101 // DataArrayView tempView(m_data,value.getShape());
102 // setPointDataView(tempView);
103 // }
104
105 DataConstant::DataConstant(const DataConstant& other)
106 : parent(other.getFunctionSpace(),other.getShape())
107 { //
108 // copy the data in the correct format
109 m_data=other.m_data;
110 //
111 // // create the view of the data
112 // DataArrayView tempView(m_data,other.getPointDataView().getShape());
113 // setPointDataView(tempView);
114 }
115
116 DataConstant::DataConstant(const DataConstant& other,
117 const DataTypes::RegionType& region)
118 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
119 {
120 //
121 // get the shape of the slice to copy from
122 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
123 //
124 // allocate space for this new DataConstant's data
125 int len = getNoValues();
126 m_data.resize(len,0.,len);
127 //
128 // create a view of the data with the correct shape
129 // DataArrayView tempView(m_data,shape);
130 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
131 //
132 // load the view with the data from the slice
133 // tempView.copySlice(other.getPointDataView(),region_loop_range);
134 DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
135 // setPointDataView(tempView);
136 }
137
138 DataConstant::DataConstant(const FunctionSpace& what,
139 const DataTypes::ShapeType &shape,
140 const DataTypes::ValueType &data)
141 : parent(what,shape)
142 {
143 //
144 // copy the data in the correct format
145 m_data=data;
146 //
147 // create the view of the data
148 // DataArrayView tempView(m_data,shape);
149 // setPointDataView(tempView);
150 }
151
152 string
153 DataConstant::toString() const
154 {
155 return DataTypes::pointToString(m_data,getShape(),0,"");
156 }
157
158
159 DataAbstract*
160 DataConstant::deepCopy()
161 {
162 return new DataConstant(*this);
163 }
164
165
166 DataTypes::ValueType::size_type
167 DataConstant::getPointOffset(int sampleNo,
168 int dataPointNo) const
169 {
170 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
171 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
172 //
173 // Whatever the coord's always return the same value as this is constant data.
174 return 0;
175 }
176
177 DataTypes::ValueType::size_type
178 DataConstant::getPointOffset(int sampleNo,
179 int dataPointNo)
180 {
181 EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
182 "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
183 //
184 // Whatever the coord's always return the same value as this is constant data.
185 return 0;
186 }
187
188 DataTypes::ValueType::size_type
189 DataConstant::getLength() const
190 {
191 return m_data.size();
192 }
193
194 // DataArrayView
195 // DataConstant::getDataPoint(int sampleNo,
196 // int dataPointNo)
197 // {
198 // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
199 // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
200 // //
201 // // Whatever the coord's always return the same value as this is constant data.
202 // return getPointDataView();
203 // }
204
205 DataAbstract*
206 DataConstant::getSlice(const DataTypes::RegionType& region) const
207 {
208 return new DataConstant(*this,region);
209 }
210
211 void
212 DataConstant::setSlice(const DataAbstract* value,
213 const DataTypes::RegionType& region)
214 {
215 const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
216 if (tempDataConst==0) {
217 throw DataException("Programming error - casting to DataConstant.");
218 }
219 //
220 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
221 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
222 //
223 // check shape:
224 if (getRank()!=region.size()) {
225 throw DataException("Error - Invalid slice region.");
226 }
227 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
228 throw DataException (DataTypes::createShapeErrorMessage(
229 "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
230 }
231 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
232 DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
233 }
234
235
236
237 void
238 DataConstant::symmetric(DataAbstract* ev)
239 {
240 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
241 if (temp_ev==0) {
242 throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
243 }
244 /* DataArrayView& thisView=getPointDataView();
245 DataArrayView& evView=ev->getPointDataView();*/
246 DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
247 }
248
249 void
250 DataConstant::nonsymmetric(DataAbstract* ev)
251 {
252 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
253 if (temp_ev==0) {
254 throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
255 }
256 /* DataArrayView& thisView=getPointDataView();
257 DataArrayView& evView=ev->getPointDataView();*/
258 DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
259 }
260
261 void
262 DataConstant::trace(DataAbstract* ev, int axis_offset)
263 {
264 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
265 if (temp_ev==0) {
266 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
267 }
268 /* DataArrayView& thisView=getPointDataView();
269 DataArrayView& evView=ev->getPointDataView();*/
270 ValueType& evVec=temp_ev->getVector();
271 const ShapeType& evShape=temp_ev->getShape();
272 DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
273 }
274
275 void
276 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
277 {
278 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
279 if (temp_ev==0) {
280 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
281 }
282 // DataArrayView& thisView=getPointDataView();
283 // DataArrayView& evView=ev->getPointDataView();
284 DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
285 }
286
287 void
288 DataConstant::transpose(DataAbstract* ev, int axis_offset)
289 {
290 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
291 if (temp_ev==0) {
292 throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
293 }
294 /* DataArrayView& thisView=getPointDataView();
295 DataArrayView& evView=ev->getPointDataView();*/
296 DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
297 }
298
299 void
300 DataConstant::eigenvalues(DataAbstract* ev)
301 {
302 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
303 if (temp_ev==0) {
304 throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
305 }
306 /* DataArrayView& thisView=getPointDataView();
307 DataArrayView& evView=ev->getPointDataView();*/
308 DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
309 }
310 void
311 DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
312 {
313 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
314 if (temp_ev==0) {
315 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
316 }
317 DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
318 if (temp_V==0) {
319 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
320 }
321 // DataArrayView thisView=getPointDataView();
322 // DataArrayView evView=ev->getPointDataView();
323 // DataArrayView VView=V->getPointDataView();
324
325 DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
326 }
327
328 void
329 DataConstant::setToZero()
330 {
331 DataTypes::ValueType::size_type n=m_data.size();
332 for (int i=0; i<n ;++i) m_data[i]=0.;
333 }
334
335 void
336 DataConstant::dump(const std::string fileName) const
337 {
338 #ifdef USE_NETCDF
339 const NcDim* ncdims[DataTypes::maxRank];
340 NcVar* var;
341 int rank = getRank();
342 int type= getFunctionSpace().getTypeCode();
343 int ndims =0;
344 long dims[DataTypes::maxRank];
345 const double* d_ptr=&(m_data[0]);
346 DataTypes::ShapeType shape = getShape();
347 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
348 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
349 #ifdef PASO_MPI
350 MPI_Status status;
351 #endif
352
353 #ifdef PASO_MPI
354 /* Serialize NetCDF I/O */
355 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
356 #endif
357
358 // netCDF error handler
359 NcError err(NcError::verbose_nonfatal);
360 // Create the file.
361 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
362 NcFile dataFile(newFileName, NcFile::Replace);
363 // check if writing was successful
364 if (!dataFile.is_valid())
365 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
366 if (!dataFile.add_att("type_id",0) )
367 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
368 if (!dataFile.add_att("rank",rank) )
369 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
370 if (!dataFile.add_att("function_space_type",type))
371 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
372
373 if (rank == 0) {
374 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
375 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
376 dims[0]=1,
377 ndims=1;
378 } else {
379 ndims=rank;
380 dims[0]=shape[0];
381 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
382 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
383 if ( rank >1 ) {
384 dims[1]=shape[1];
385 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
386 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
387 }
388 if ( rank >2 ) {
389 dims[2]=shape[2];
390 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
391 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
392 }
393 if ( rank >3 ) {
394 dims[3]=shape[3];
395 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
396 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
397 }
398 }
399
400 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
401 throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
402 if (! (var->put(d_ptr,dims)) )
403 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
404 #ifdef PASO_MPI
405 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
406 #endif
407 #else
408 throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
409 #endif
410 }
411
412 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26