/[escript]/trunk/escript/src/DataConstant.cpp
ViewVC logotype

Contents of /trunk/escript/src/DataConstant.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File size: 12029 byte(s)
Merging dudley and scons updates from branches

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26