/[escript]/branches/subworld2/escriptcore/src/DataConstant.cpp
ViewVC logotype

Contents of /branches/subworld2/escriptcore/src/DataConstant.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5504 - (show annotations)
Wed Mar 4 22:58:13 2015 UTC (4 years, 1 month ago) by jfenwick
File size: 12877 byte(s)
Again with a more up to date copy


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26