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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2742 - (hide annotations)
Thu Nov 12 06:03:37 2009 UTC (10 years, 1 month ago) by jfenwick
File size: 12087 byte(s)
Merging changes from the lapack branch.

The inverse() operation has been moved into c++. [No lazy support for this operation yet.]
Optional Lapack support has been added for matrices larger than 3x3. 
service0 is set to use mkl_lapack.



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