/[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 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 9 months ago) by ksteube
Original Path: trunk/escript/src/DataConstant.cpp
File size: 12331 byte(s)
Copyright updated in all files

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