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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2105 - (hide annotations)
Fri Nov 28 01:52:12 2008 UTC (10 years, 7 months ago) by jfenwick
Original Path: trunk/escript/src/DataConstant.cpp
File size: 13154 byte(s)
Data::copySelf() now returns an object instead of a pointer.
Fixed a bug in copyFromArray relating to expanded data.
1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4     * Copyright (c) 2003-2008 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 ksteube 1312
14 ksteube 1811
15 ksteube 1827 #include "Data.h"
16 jgs 474 #include "DataConstant.h"
17     #include "DataException.h"
18 robwdcock 682 #include "esysUtils/EsysAssert.h"
19 jgs 82
20     #include <iostream>
21     #include <boost/python/extract.hpp>
22 gross 1023 #ifdef USE_NETCDF
23 ksteube 1312 #include <netcdfcpp.h>
24 gross 1023 #endif
25 ksteube 1827 #ifdef PASO_MPI
26     #include <mpi.h>
27     #endif
28 jgs 82
29 matt 1319 #include <boost/python/extract.hpp>
30 jfenwick 1796 #include "DataMaths.h"
31 matt 1319
32 jgs 82 using namespace std;
33 matt 1319 using namespace boost::python;
34 jgs 82
35     namespace escript {
36    
37 jgs 102 DataConstant::DataConstant(const boost::python::numeric::array& value,
38     const FunctionSpace& what)
39 jfenwick 2005 : parent(what,DataTypes::shapeFromNumArray(value))
40 jgs 102 {
41 jfenwick 1796 // // extract the shape of the numarray
42     // DataTypes::ShapeType tempShape;
43     // for (int i=0; i < value.getrank(); i++) {
44     // tempShape.push_back(extract<int>(value.getshape()[i]));
45     // }
46 matt 1319
47     // get the space for the data vector
48 jfenwick 1796 // int len = getNoValues();
49     // DataVector temp_data(len, 0.0, len);
50     // DataArrayView temp_dataView(temp_data, tempShape);
51     // temp_dataView.copy(value);
52 matt 1319
53 jfenwick 2105 m_data.copyFromNumArray(value,1);
54 jgs 102 //
55 jgs 82
56 jgs 102 // copy the data in the correct format
57 jfenwick 1796 // m_data=temp_data;
58 jgs 119 //
59     // create the view of the data
60 jfenwick 1796 // DataArrayView tempView(m_data,temp_dataView.getShape());
61     // setPointDataView(tempView);
62 jgs 102 }
63 jgs 82
64 jfenwick 1796 // DataConstant::DataConstant(const DataArrayView& value,
65     // const FunctionSpace& what)
66     // : DataAbstract(what)
67     // {
68     // //
69     // // copy the data in the correct format
70     // m_data=value.getData();
71     // //
72     // // create the view of the data
73     // DataArrayView tempView(m_data,value.getShape());
74     // setPointDataView(tempView);
75     // }
76    
77 jgs 102 DataConstant::DataConstant(const DataConstant& other)
78 jfenwick 2005 : parent(other.getFunctionSpace(),other.getShape())
79 matt 1319 { //
80 jgs 102 // copy the data in the correct format
81     m_data=other.m_data;
82 jgs 119 //
83 jfenwick 1796 // // create the view of the data
84     // DataArrayView tempView(m_data,other.getPointDataView().getShape());
85     // setPointDataView(tempView);
86 jgs 102 }
87 jgs 82
88 jgs 102 DataConstant::DataConstant(const DataConstant& other,
89 jfenwick 1796 const DataTypes::RegionType& region)
90 jfenwick 2005 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
91 jgs 102 {
92     //
93     // get the shape of the slice to copy from
94 jfenwick 1796 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
95 jgs 102 //
96     // allocate space for this new DataConstant's data
97 jfenwick 1796 int len = getNoValues();
98 jgs 151 m_data.resize(len,0.,len);
99 jgs 102 //
100     // create a view of the data with the correct shape
101 jfenwick 1796 // DataArrayView tempView(m_data,shape);
102     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
103 jgs 119 //
104     // load the view with the data from the slice
105 jfenwick 1796 // tempView.copySlice(other.getPointDataView(),region_loop_range);
106     DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
107     // setPointDataView(tempView);
108 jgs 119 }
109    
110     DataConstant::DataConstant(const FunctionSpace& what,
111 jfenwick 1796 const DataTypes::ShapeType &shape,
112     const DataTypes::ValueType &data)
113 jfenwick 2005 : parent(what,shape)
114 jgs 119 {
115 jgs 102 //
116 jgs 119 // copy the data in the correct format
117     m_data=data;
118     //
119     // create the view of the data
120 jfenwick 1796 // DataArrayView tempView(m_data,shape);
121     // setPointDataView(tempView);
122 jgs 102 }
123 jgs 82
124 jgs 102 string
125     DataConstant::toString() const
126     {
127 jfenwick 1796 return DataTypes::pointToString(m_data,getShape(),0,"");
128 jgs 102 }
129 jgs 82
130 jfenwick 1799
131     DataAbstract*
132     DataConstant::deepCopy()
133     {
134     return new DataConstant(*this);
135     }
136    
137    
138 jfenwick 1796 DataTypes::ValueType::size_type
139 jgs 102 DataConstant::getPointOffset(int sampleNo,
140     int dataPointNo) const
141     {
142     EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
143     "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
144 jgs 122 //
145     // Whatever the coord's always return the same value as this is constant data.
146 jgs 102 return 0;
147     }
148 jgs 82
149 jfenwick 1796 DataTypes::ValueType::size_type
150 jfenwick 2005 DataConstant::getPointOffset(int sampleNo,
151     int dataPointNo)
152     {
153     EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
154     "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
155     //
156     // Whatever the coord's always return the same value as this is constant data.
157     return 0;
158     }
159    
160     DataTypes::ValueType::size_type
161 jgs 102 DataConstant::getLength() const
162     {
163     return m_data.size();
164     }
165 jgs 82
166 jfenwick 1796 // DataArrayView
167     // DataConstant::getDataPoint(int sampleNo,
168     // int dataPointNo)
169     // {
170     // EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
171     // "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
172     // //
173     // // Whatever the coord's always return the same value as this is constant data.
174     // return getPointDataView();
175     // }
176 matt 1319
177 jgs 102 DataAbstract*
178 jfenwick 1796 DataConstant::getSlice(const DataTypes::RegionType& region) const
179 jgs 102 {
180     return new DataConstant(*this,region);
181     }
182 jgs 82
183 jgs 102 void
184     DataConstant::setSlice(const DataAbstract* value,
185 jfenwick 1796 const DataTypes::RegionType& region)
186 jgs 102 {
187     const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
188     if (tempDataConst==0) {
189     throw DataException("Programming error - casting to DataConstant.");
190 jgs 82 }
191 matt 1319 //
192 jfenwick 1796 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
193     DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
194 jgs 108 //
195     // check shape:
196 jfenwick 1796 if (getRank()!=region.size()) {
197 jgs 108 throw DataException("Error - Invalid slice region.");
198     }
199 jfenwick 1796 if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
200     throw DataException (DataTypes::createShapeErrorMessage(
201     "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
202 jgs 108 }
203 jfenwick 1796 // getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
204     DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
205 jgs 102 }
206 jgs 82
207 jgs 123
208    
209 gross 580 void
210 ksteube 775 DataConstant::symmetric(DataAbstract* ev)
211     {
212     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
213     if (temp_ev==0) {
214     throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
215     }
216 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
217     DataArrayView& evView=ev->getPointDataView();*/
218     DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
219 ksteube 775 }
220    
221     void
222     DataConstant::nonsymmetric(DataAbstract* ev)
223     {
224     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
225     if (temp_ev==0) {
226     throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
227     }
228 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
229     DataArrayView& evView=ev->getPointDataView();*/
230     DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
231 ksteube 775 }
232    
233     void
234 gross 800 DataConstant::trace(DataAbstract* ev, int axis_offset)
235 ksteube 775 {
236     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
237     if (temp_ev==0) {
238 gross 800 throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
239 ksteube 775 }
240 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
241     DataArrayView& evView=ev->getPointDataView();*/
242     ValueType& evVec=temp_ev->getVector();
243     const ShapeType& evShape=temp_ev->getShape();
244     DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
245 ksteube 775 }
246    
247     void
248 gross 804 DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
249 gross 800 {
250     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
251     if (temp_ev==0) {
252 gross 804 throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
253 gross 800 }
254 jfenwick 1796 // DataArrayView& thisView=getPointDataView();
255     // DataArrayView& evView=ev->getPointDataView();
256     DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
257 gross 800 }
258    
259     void
260 ksteube 775 DataConstant::transpose(DataAbstract* ev, int axis_offset)
261     {
262     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
263     if (temp_ev==0) {
264     throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
265     }
266 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
267     DataArrayView& evView=ev->getPointDataView();*/
268     DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
269 ksteube 775 }
270    
271     void
272 gross 580 DataConstant::eigenvalues(DataAbstract* ev)
273     {
274     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
275     if (temp_ev==0) {
276     throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
277     }
278 jfenwick 1796 /* DataArrayView& thisView=getPointDataView();
279     DataArrayView& evView=ev->getPointDataView();*/
280     DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
281 gross 580 }
282     void
283     DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
284     {
285     DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
286     if (temp_ev==0) {
287     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
288     }
289     DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
290     if (temp_V==0) {
291     throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
292     }
293 jfenwick 1796 // DataArrayView thisView=getPointDataView();
294     // DataArrayView evView=ev->getPointDataView();
295     // DataArrayView VView=V->getPointDataView();
296 gross 580
297 jfenwick 1796 DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
298 gross 580 }
299    
300 gross 950 void
301 gross 1118 DataConstant::setToZero()
302     {
303 jfenwick 1796 DataTypes::ValueType::size_type n=m_data.size();
304 gross 1118 for (int i=0; i<n ;++i) m_data[i]=0.;
305     }
306    
307     void
308 gross 950 DataConstant::dump(const std::string fileName) const
309     {
310 gross 1149 #ifdef USE_NETCDF
311 jfenwick 1796 const NcDim* ncdims[DataTypes::maxRank];
312 gross 950 NcVar* var;
313 jfenwick 1796 int rank = getRank();
314 gross 950 int type= getFunctionSpace().getTypeCode();
315     int ndims =0;
316 jfenwick 1796 long dims[DataTypes::maxRank];
317 gross 1141 const double* d_ptr=&(m_data[0]);
318 jfenwick 1796 DataTypes::ShapeType shape = getShape();
319 jfenwick 1872 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
320     int mpi_num=getFunctionSpace().getDomain()->getMPISize();
321 ksteube 1827 #ifdef PASO_MPI
322     MPI_Status status;
323     #endif
324 matt 1319
325 ksteube 1827 #ifdef PASO_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 gross 950 // netCDF error handler
331     NcError err(NcError::verbose_nonfatal);
332     // Create the file.
333 ksteube 1827 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
334     NcFile dataFile(newFileName, NcFile::Replace);
335 gross 950 // check if writing was successful
336 matt 1319 if (!dataFile.is_valid())
337 gross 950 throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
338 matt 1319 if (!dataFile.add_att("type_id",0) )
339 gross 950 throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
340 matt 1319 if (!dataFile.add_att("rank",rank) )
341 gross 950 throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
342 matt 1319 if (!dataFile.add_att("function_space_type",type))
343 gross 950 throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");
344 gross 580
345 gross 950 if (rank == 0) {
346 matt 1319 if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
347 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
348 gross 950 dims[0]=1,
349     ndims=1;
350     } else {
351     ndims=rank;
352     dims[0]=shape[0];
353 matt 1319 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
354 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
355 gross 950 if ( rank >1 ) {
356     dims[1]=shape[1];
357 matt 1319 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
358 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
359 gross 950 }
360     if ( rank >2 ) {
361     dims[2]=shape[2];
362 matt 1319 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
363 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
364 gross 950 }
365     if ( rank >3 ) {
366     dims[3]=shape[3];
367 matt 1319 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
368 ksteube 1800 throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
369 gross 950 }
370     }
371 matt 1319
372 gross 950 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
373     throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
374 gross 1141 if (! (var->put(d_ptr,dims)) )
375 matt 1319 throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
376 ksteube 1827 #ifdef PASO_MPI
377     if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
378     #endif
379 gross 1023 #else
380     throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
381     #endif
382 gross 950 }
383    
384 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26