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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2682 - (show annotations)
Fri Sep 25 01:35:55 2009 UTC (9 years, 9 months ago) by jfenwick
File size: 11450 byte(s)
Experimenting with Lapack

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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 // #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 // create the view of the data
80 // DataArrayView tempView(m_data,shape);
81 // setPointDataView(tempView);
82 }
83
84 string
85 DataConstant::toString() const
86 {
87 return DataTypes::pointToString(m_data,getShape(),0,"");
88 }
89
90
91 DataAbstract*
92 DataConstant::deepCopy()
93 {
94 return new DataConstant(*this);
95 }
96
97
98 DataTypes::ValueType::size_type
99 DataConstant::getPointOffset(int sampleNo,
100 int dataPointNo) const
101 {
102 // We avoid this check for constant data due to issues manipulating
103 // data with no samples.
104
105 // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
106 // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
107 //
108 // Whatever the coord's always return the same value as this is constant data.
109 return 0;
110 }
111
112 DataTypes::ValueType::size_type
113 DataConstant::getPointOffset(int sampleNo,
114 int dataPointNo)
115 {
116 // We avoid this check for constant data due to issues manipulating
117 // data with no samples.
118
119 // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
120 // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
121 //
122 // Whatever the coord's always return the same value as this is constant data.
123 //
124
125 return 0;
126 }
127
128 DataTypes::ValueType::size_type
129 DataConstant::getLength() const
130 {
131 return m_data.size();
132 }
133
134 DataAbstract*
135 DataConstant::getSlice(const DataTypes::RegionType& region) const
136 {
137 return new DataConstant(*this,region);
138 }
139
140 void
141 DataConstant::setSlice(const DataAbstract* value,
142 const DataTypes::RegionType& region)
143 {
144 const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
145 if (tempDataConst==0) {
146 throw DataException("Programming error - casting to DataConstant.");
147 }
148 CHECK_FOR_EX_WRITE
149 //
150 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
151 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
152 //
153 // check shape:
154 if (getRank()!=region.size()) {
155 throw DataException("Error - Invalid slice region.");
156 }
157 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
158 throw DataException (DataTypes::createShapeErrorMessage(
159 "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
160 }
161 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
162 DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVectorRO(), tempDataConst->getShape(),0,region_loop_range);
163 }
164
165
166
167 void
168 DataConstant::symmetric(DataAbstract* ev)
169 {
170 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
171 if (temp_ev==0) {
172 throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
173 }
174 DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
175 }
176
177 void
178 DataConstant::nonsymmetric(DataAbstract* ev)
179 {
180 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
181 if (temp_ev==0) {
182 throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
183 }
184 DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
185 }
186
187 void
188 DataConstant::trace(DataAbstract* ev, int axis_offset)
189 {
190 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
191 if (temp_ev==0) {
192 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
193 }
194 ValueType& evVec=temp_ev->getVectorRW();
195 const ShapeType& evShape=temp_ev->getShape();
196 DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
197 }
198
199 void
200 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
201 {
202 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
203 if (temp_ev==0) {
204 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
205 }
206 DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0,axis0,axis1);
207 }
208
209 void
210 DataConstant::transpose(DataAbstract* ev, int axis_offset)
211 {
212 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
213 if (temp_ev==0) {
214 throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
215 }
216 DataMaths::transpose(m_data, getShape(),0, temp_ev->getVectorRW(),temp_ev->getShape(),0,axis_offset);
217 }
218
219 void
220 DataConstant::eigenvalues(DataAbstract* ev)
221 {
222 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
223 if (temp_ev==0) {
224 throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
225 }
226 DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0);
227 }
228 void
229 DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
230 {
231 DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
232 if (temp_ev==0) {
233 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
234 }
235 DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
236 if (temp_V==0) {
237 throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
238 }
239 DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVectorRW(), temp_ev->getShape(),0,temp_V->getVectorRW(), temp_V->getShape(),0,tol);
240 }
241
242 void
243 DataConstant::setToZero()
244 {
245 CHECK_FOR_EX_WRITE
246 DataTypes::ValueType::size_type n=m_data.size();
247 for (int i=0; i<n ;++i) m_data[i]=0.;
248 }
249
250 void
251 DataConstant::dump(const std::string fileName) const
252 {
253 #ifdef USE_NETCDF
254 const NcDim* ncdims[DataTypes::maxRank];
255 NcVar* var;
256 int rank = getRank();
257 int type= getFunctionSpace().getTypeCode();
258 int ndims =0;
259 long dims[DataTypes::maxRank];
260 const double* d_ptr=&(m_data[0]);
261 DataTypes::ShapeType shape = getShape();
262 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
263 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
264 #ifdef PASO_MPI
265 MPI_Status status;
266 #endif
267
268 #ifdef PASO_MPI
269 /* Serialize NetCDF I/O */
270 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
271 #endif
272
273 // netCDF error handler
274 NcError err(NcError::verbose_nonfatal);
275 // Create the file.
276 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
277 NcFile dataFile(newFileName, NcFile::Replace);
278 // check if writing was successful
279 if (!dataFile.is_valid())
280 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
281 if (!dataFile.add_att("type_id",0) )
282 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
283 if (!dataFile.add_att("rank",rank) )
284 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
285 if (!dataFile.add_att("function_space_type",type))
286 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
287
288 if (rank == 0) {
289 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
290 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
291 dims[0]=1,
292 ndims=1;
293 } else {
294 ndims=rank;
295 dims[0]=shape[0];
296 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
297 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
298 if ( rank >1 ) {
299 dims[1]=shape[1];
300 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
301 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
302 }
303 if ( rank >2 ) {
304 dims[2]=shape[2];
305 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
306 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
307 }
308 if ( rank >3 ) {
309 dims[3]=shape[3];
310 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
311 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
312 }
313 }
314
315 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
316 throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
317 if (! (var->put(d_ptr,dims)) )
318 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
319 #ifdef PASO_MPI
320 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
321 #endif
322 #else
323 throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
324 #endif
325 }
326
327 // These used to be marked as inline in DataConstant.
328 // But they are marked virtual in DataReady
329 DataTypes::ValueType&
330 DataConstant::getVectorRW()
331 {
332 CHECK_FOR_EX_WRITE
333 return m_data;
334 }
335
336 const DataTypes::ValueType&
337 DataConstant::getVectorRO() const
338 {
339 return m_data;
340 }
341
342 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26