escript/src/DataConstant.cpp


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


#include "Data.h"
#include "DataConstant.h"
#include "DataException.h"
#include "esysUtils/EsysAssert.h"

#include <iostream>
#include <boost/python/extract.hpp>
#ifdef USE_NETCDF
#include <netcdfcpp.h>
#endif
#ifdef PASO_MPI
#include <mpi.h>
#endif

#include <boost/python/extract.hpp>
#include "DataMaths.h"

using namespace std;
using namespace boost::python;

namespace escript {

DataConstant::DataConstant(const boost::python::numeric::array& value,
                           const FunctionSpace& what)
  : DataAbstract(what,DataTypes::shapeFromNumArray(value))
{
//   // extract the shape of the numarray
//   DataTypes::ShapeType tempShape;
//   for (int i=0; i < value.getrank(); i++) {
//     tempShape.push_back(extract<int>(value.getshape()[i]));
//   }

  // get the space for the data vector
//   int len = getNoValues();
//   DataVector temp_data(len, 0.0, len);
//   DataArrayView temp_dataView(temp_data, tempShape);
//   temp_dataView.copy(value);

  m_data.copyFromNumArray(value);
  //

  // copy the data in the correct format
//   m_data=temp_data;
  //
  // create the view of the data
//   DataArrayView tempView(m_data,temp_dataView.getShape());
//   setPointDataView(tempView);
}

// DataConstant::DataConstant(const DataArrayView& value,
//                            const FunctionSpace& what)
//   : DataAbstract(what)
// {
//   //
//   // copy the data in the correct format
//   m_data=value.getData();
//   //
//   // create the view of the data
//   DataArrayView tempView(m_data,value.getShape());
//   setPointDataView(tempView);
// }

DataConstant::DataConstant(const DataConstant& other)
  : DataAbstract(other.getFunctionSpace(),other.getShape())
{  //
  // copy the data in the correct format
  m_data=other.m_data;
  //
//   // create the view of the data
//   DataArrayView tempView(m_data,other.getPointDataView().getShape());
//   setPointDataView(tempView);
}

DataConstant::DataConstant(const DataConstant& other,
                           const DataTypes::RegionType& region)
  : DataAbstract(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
{
  //
  // get the shape of the slice to copy from
//   DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
  //
  // allocate space for this new DataConstant's data
  int len = getNoValues();
  m_data.resize(len,0.,len);
  //
  // create a view of the data with the correct shape
//   DataArrayView tempView(m_data,shape);
  DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
  //
  // load the view with the data from the slice
//   tempView.copySlice(other.getPointDataView(),region_loop_range);
  DataTypes::copySlice(m_data,getShape(),0,other.getVector(),other.getShape(),0,region_loop_range);
//   setPointDataView(tempView);
}

DataConstant::DataConstant(const FunctionSpace& what,
                           const DataTypes::ShapeType &shape,
                           const DataTypes::ValueType &data)
  : DataAbstract(what,shape)
{
  //
  // copy the data in the correct format
  m_data=data;
  //
  // create the view of the data
//   DataArrayView tempView(m_data,shape);
//   setPointDataView(tempView);
}

string
DataConstant::toString() const
{
  return DataTypes::pointToString(m_data,getShape(),0,"");
}


DataAbstract*
DataConstant::deepCopy()
{
  return new DataConstant(*this);
}


DataTypes::ValueType::size_type
DataConstant::getPointOffset(int sampleNo,
                             int dataPointNo) const
{
  EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
              "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
  //
  // Whatever the coord's always return the same value as this is constant data.
  return 0;
}

DataTypes::ValueType::size_type
DataConstant::getLength() const
{
  return m_data.size();
}

// DataArrayView
// DataConstant::getDataPoint(int sampleNo,
//                            int dataPointNo)
// {
//   EsysAssert((validSamplePointNo(dataPointNo) && validSampleNo(sampleNo)),
//              "Invalid index, sampleNo: " << sampleNo << " dataPointNo: " << dataPointNo);
//   //
//   // Whatever the coord's always return the same value as this is constant data.
//   return getPointDataView();
// }

DataAbstract*
DataConstant::getSlice(const DataTypes::RegionType& region) const
{
  return new DataConstant(*this,region);
}

void
DataConstant::setSlice(const DataAbstract* value,
                       const DataTypes::RegionType& region)
{
  const DataConstant* tempDataConst=dynamic_cast<const DataConstant*>(value);
  if (tempDataConst==0) {
    throw DataException("Programming error - casting to DataConstant.");
  }
  //
  DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
  DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
  //
  // check shape:
  if (getRank()!=region.size()) {
    throw DataException("Error - Invalid slice region.");
  }
  if (getRank()>0 && !DataTypes::checkShape(value->getShape(),shape)) {
    throw DataException (DataTypes::createShapeErrorMessage(
                "Error - Couldn't copy slice due to shape mismatch.",shape,value->getShape()));
  }
  //   getPointDataView().copySliceFrom(tempDataConst->getPointDataView(),region_loop_range);
  DataTypes::copySliceFrom(m_data,getShape(),0,tempDataConst->getVector(), tempDataConst->getShape(),0,region_loop_range);
}


void
DataConstant::symmetric(DataAbstract* ev)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::symmetric: casting to DataConstant failed (propably a programming error).");
  }
/*  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();*/
  DataMaths::symmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
}

void
DataConstant::nonsymmetric(DataAbstract* ev)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::nonsymmetric: casting to DataConstant failed (propably a programming error).");
  }
/*  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();*/
  DataMaths::nonsymmetric(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
}

void
DataConstant::trace(DataAbstract* ev, int axis_offset)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::trace: casting to DataConstant failed (propably a programming error).");
  }
/*  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();*/
  ValueType& evVec=temp_ev->getVector();
  const ShapeType& evShape=temp_ev->getShape();
  DataMaths::trace(m_data,getShape(),0,evVec,evShape,0,axis_offset);
}

void
DataConstant::swapaxes(DataAbstract* ev, int axis0, int axis1)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::swapaxes: casting to DataConstant failed (propably a programming error).");
  }
//   DataArrayView& thisView=getPointDataView();
//   DataArrayView& evView=ev->getPointDataView();
  DataMaths::swapaxes(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,axis0,axis1);
}

void
DataConstant::transpose(DataAbstract* ev, int axis_offset)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::transpose: casting to DataConstant failed (propably a programming error).");
  }
/*  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();*/
  DataMaths::transpose(m_data, getShape(),0, temp_ev->getVector(),temp_ev->getShape(),0,axis_offset);
}

void
DataConstant::eigenvalues(DataAbstract* ev)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::eigenvalues: casting to DataConstant failed (propably a programming error).");
  }
/*  DataArrayView& thisView=getPointDataView();
  DataArrayView& evView=ev->getPointDataView();*/
  DataMaths::eigenvalues(m_data,getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0);
}
void
DataConstant::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
{
  DataConstant* temp_ev=dynamic_cast<DataConstant*>(ev);
  if (temp_ev==0) {
    throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
  }
  DataConstant* temp_V=dynamic_cast<DataConstant*>(V);
  if (temp_V==0) {
    throw DataException("Error - DataConstant::eigenvalues_and_eigenvectors: casting to DataConstant failed (propably a programming error).");
  }
//   DataArrayView thisView=getPointDataView();
//   DataArrayView evView=ev->getPointDataView();
//   DataArrayView VView=V->getPointDataView();

  DataMaths::eigenvalues_and_eigenvectors(m_data, getShape(),0,temp_ev->getVector(), temp_ev->getShape(),0,temp_V->getVector(), temp_V->getShape(),0,tol);
}

void
DataConstant::setToZero()
{
    DataTypes::ValueType::size_type n=m_data.size();
    for (int i=0; i<n ;++i) m_data[i]=0.;
}

void
DataConstant::dump(const std::string fileName) const
{
   #ifdef USE_NETCDF
   const NcDim* ncdims[DataTypes::maxRank];
   NcVar* var;
   int rank = getRank();
   int type=  getFunctionSpace().getTypeCode();
   int ndims =0;
   long dims[DataTypes::maxRank];
   const double* d_ptr=&(m_data[0]);
   DataTypes::ShapeType shape = getShape();
   int mpi_iam=getFunctionSpace().getDomain().getMPIRank();
   int mpi_num=getFunctionSpace().getDomain().getMPISize();
#ifdef PASO_MPI
   MPI_Status status;
#endif

#ifdef PASO_MPI
   /* Serialize NetCDF I/O */
   if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81802, MPI_COMM_WORLD, &status);
#endif

   // netCDF error handler
   NcError err(NcError::verbose_nonfatal);
   // Create the file.
   char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
   NcFile dataFile(newFileName, NcFile::Replace);
   // check if writing was successful
   if (!dataFile.is_valid())
    throw DataException("Error - DataConstant:: opening of netCDF file for output failed.");
   if (!dataFile.add_att("type_id",0) )
    throw DataException("Error - DataConstant:: appending data type to netCDF file failed.");
   if (!dataFile.add_att("rank",rank) )
    throw DataException("Error - DataConstant:: appending rank attribute to netCDF file failed.");
   if (!dataFile.add_att("function_space_type",type))
    throw DataException("Error - DataConstant:: appending function space attribute to netCDF file failed.");

   if (rank == 0) {
      if( ! (ncdims[0] = dataFile.add_dim("l", 1)) )
        throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
      dims[0]=1,
      ndims=1;
   } else {
       ndims=rank;
       dims[0]=shape[0];
       if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
        throw DataException("Error - DataConstant:: appending ncdimension 0 to netCDF file failed.");
       if ( rank >1 ) {
           dims[1]=shape[1];
           if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
        throw DataException("Error - DataConstant:: appending ncdimension 1 to netCDF file failed.");
       }
       if ( rank >2 ) {
           dims[2]=shape[2];
           if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
        throw DataException("Error - DataConstant:: appending ncdimension 2 to netCDF file failed.");
       }
       if ( rank >3 ) {
           dims[3]=shape[3];
           if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
        throw DataException("Error - DataConstant:: appending ncdimension 3 to netCDF file failed.");
       }
   }

   if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
    throw DataException("Error - DataConstant:: appending variable to netCDF file failed.");
   if (! (var->put(d_ptr,dims)) )
         throw DataException("Error - DataConstant:: copy data to netCDF buffer failed.");
#ifdef PASO_MPI
   if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81802, MPI_COMM_WORLD);
#endif
   #else
   throw DataException("Error - DataConstant:: dump is not configured with netCDF. Please contact your installation manager.");
   #endif
}

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