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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 11 months ago) by jfenwick
File size: 12073 byte(s)
Don't panic.
Updating copyright stamps

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26