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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (show annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 5 months ago) by caltinay
File size: 12374 byte(s)
Assorted spelling fixes.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26