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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (hide annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 4 months ago) by phornby
File size: 29036 byte(s)

Merge in /branches/windows_from_1456_trunk_1620_merged_in branch.

You will find a preserved pre-merge trunk in tags under tags/trunk_at_1625.
That will be useful for diffing & checking on my stupidity.

Here is a list of the conflicts and their resolution at this
point in time.


=================================================================================
(LLWS == looks like white space).

finley/src/Assemble_addToSystemMatrix.c - resolve to branch - unused var. may be wrong.....
finley/src/CPPAdapter/SystemMatrixAdapter.cpp - resolve to branch - LLWS
finley/src/CPPAdapter/MeshAdapter.cpp - resolve to branch - LLWS
paso/src/PCG.c - resolve to branch - unused var fixes.
paso/src/SolverFCT.c - resolve to branch - LLWS
paso/src/FGMRES.c - resolve to branch - LLWS
paso/src/Common.h - resolve to trunk version. It's omp.h's include... not sure it's needed,
but for the sake of saftey.....
paso/src/Functions.c - resolve to branch version, indentation/tab removal and return error
on bad unimplemented Paso_FunctionCall.
paso/src/SolverFCT_solve.c - resolve to branch version, unused vars
paso/src/SparseMatrix_MatrixVector.c - resolve to branch version, unused vars.
escript/src/Utils.cpp - resloved to branch, needs WinSock2.h
escript/src/DataExpanded.cpp - resolved to branch version - LLWS
escript/src/DataFactory.cpp - resolve to branch version
=================================================================================

This currently passes tests on linux (debian), but is not checked on windows or Altix yet.

This checkin is to make a trunk I can check out for windows to do tests on it.

Known outstanding problem is in the operator=() method of exceptions
causing warning messages on the intel compilers.

May the God of doughnuts have mercy on my soul.


1 jgs 82
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 jgs 480 #include "DataExpanded.h"
17 jgs 474 #include "DataException.h"
18     #include "DataConstant.h"
19     #include "DataTagged.h"
20 gross 1023 #ifdef USE_NETCDF
21 ksteube 1312 #include <netcdfcpp.h>
22 gross 1023 #endif
23 jgs 82
24     #include <boost/python/extract.hpp>
25    
26     using namespace std;
27     using namespace boost::python;
28     using namespace boost;
29    
30     namespace escript {
31    
32 jgs 102 DataExpanded::DataExpanded(const boost::python::numeric::array& value,
33     const FunctionSpace& what)
34     : DataAbstract(what)
35     {
36     DataArrayView::ShapeType tempShape;
37     //
38 jgs 117 // extract the shape of the python numarray
39     for (int i=0; i<value.getrank(); i++) {
40 jgs 102 tempShape.push_back(extract<int>(value.getshape()[i]));
41 jgs 82 }
42 jgs 117 //
43     // initialise the data array for this object
44 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
45     //
46     // copy the given value to every data point
47     copy(value);
48     }
49 jgs 82
50 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other)
51     : DataAbstract(other.getFunctionSpace()),
52     m_data(other.m_data)
53 jgs 117 {
54 jgs 102 //
55     // create the view for the data
56     DataArrayView temp(m_data.getData(),other.getPointDataView().getShape());
57     setPointDataView(temp);
58     }
59 jgs 82
60 jgs 102 DataExpanded::DataExpanded(const DataConstant& other)
61     : DataAbstract(other.getFunctionSpace())
62 jgs 117 {
63     //
64     // initialise the data array for this object
65 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
66     //
67 jgs 117 // DataConstant only has one value, copy this to every data point
68 jgs 102 copy(other.getPointDataView());
69     }
70 jgs 82
71 jgs 102 DataExpanded::DataExpanded(const DataTagged& other)
72     : DataAbstract(other.getFunctionSpace())
73 jgs 117 {
74     //
75 jgs 119 // initialise the data array for this object
76 jgs 102 initialise(other.getPointDataView().getShape(),other.getNumSamples(),other.getNumDPPSample());
77 jgs 117 //
78     // for each data point in this object, extract and copy the corresponding data
79     // value from the given DataTagged object
80 jgs 102 int i,j;
81     DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
82     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
83 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
84 jgs 117 for (i=0;i<numRows;i++) {
85     for (j=0;j<numCols;j++) {
86 jgs 102 try {
87 jgs 122 getPointDataView().copy(getPointOffset(i,j),
88     other.getPointDataView(),
89     other.getPointOffset(i,j));
90 jgs 82 }
91 jgs 102 catch (std::exception& e) {
92     cout << e.what() << endl;
93     }
94 jgs 82 }
95     }
96 jgs 102 }
97 jgs 82
98 jgs 102 DataExpanded::DataExpanded(const DataExpanded& other,
99     const DataArrayView::RegionType& region)
100     : DataAbstract(other.getFunctionSpace())
101     {
102     //
103     // get the shape of the slice
104     DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
105     //
106     // initialise this Data object to the shape of the slice
107     initialise(shape,other.getNumSamples(),other.getNumDPPSample());
108     //
109 jgs 117 // copy the data
110 jgs 108 DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
111 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
112     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
113     int i,j;
114 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
115 jgs 117 for (i=0;i<numRows;i++) {
116     for (j=0;j<numCols;j++) {
117 jgs 102 try {
118 jgs 122 getPointDataView().copySlice(getPointOffset(i,j),
119     other.getPointDataView(),
120     other.getPointOffset(i,j),
121     region_loop_range);
122 jgs 82 }
123 jgs 102 catch (std::exception& e) {
124     cout << e.what() << endl;
125     }
126 jgs 82 }
127     }
128 jgs 102 }
129 jgs 82
130 jgs 102 DataExpanded::DataExpanded(const DataArrayView& value,
131     const FunctionSpace& what)
132     : DataAbstract(what)
133     {
134 jgs 117 //
135     // get the shape of the given data value
136 jgs 102 DataArrayView::ShapeType tempShape=value.getShape();
137 jgs 117 //
138     // initialise this Data object to the shape of the given data value
139 jgs 102 initialise(tempShape,what.getNumSamples(),what.getNumDPPSample());
140 jgs 117 //
141     // copy the given value to every data point
142 jgs 102 copy(value);
143     }
144 jgs 82
145 jgs 119 DataExpanded::DataExpanded(const FunctionSpace& what,
146     const DataArrayView::ShapeType &shape,
147     const DataArrayView::ValueType &data)
148     : DataAbstract(what)
149     {
150     //
151 jgs 160 // create the view of the data
152     initialise(shape,what.getNumSamples(),what.getNumDPPSample());
153     //
154 jgs 119 // copy the data in the correct format
155     m_data.getData()=data;
156     }
157    
158 jgs 102 DataExpanded::~DataExpanded()
159     {
160     }
161 jgs 82
162 jgs 102 DataAbstract*
163 matt 1319 DataExpanded::getSlice(const DataArrayView::RegionType& region) const
164 jgs 102 {
165     return new DataExpanded(*this,region);
166     }
167    
168     void
169     DataExpanded::setSlice(const DataAbstract* value,
170 matt 1319 const DataArrayView::RegionType& region)
171 jgs 102 {
172     const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
173     if (tempDataExp==0) {
174     throw DataException("Programming error - casting to DataExpanded.");
175     }
176 jgs 108 //
177 jgs 117 // get shape of slice
178 jgs 108 DataArrayView::ShapeType shape(DataArrayView::getResultSliceShape(region));
179     DataArrayView::RegionLoopRangeType region_loop_range=getSliceRegionLoopRange(region);
180     //
181 jgs 117 // check shape
182 jgs 102 if (getPointDataView().getRank()!=region.size()) {
183     throw DataException("Error - Invalid slice region.");
184     }
185 woo409 757 if (tempDataExp->getPointDataView().getRank()>0 && !value->getPointDataView().checkShape(shape)) {
186 jgs 102 throw DataException (value->getPointDataView().createShapeErrorMessage(
187 jgs 108 "Error - Couldn't copy slice due to shape mismatch.",shape));
188 jgs 102 }
189     //
190 jgs 117 // copy the data from the slice into this object
191 jgs 102 DataArrayView::ValueType::size_type numRows=m_data.getNumRows();
192     DataArrayView::ValueType::size_type numCols=m_data.getNumCols();
193     int i, j;
194 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
195 jgs 102 for (i=0;i<numRows;i++) {
196     for (j=0;j<numCols;j++) {
197 jgs 122 getPointDataView().copySliceFrom(getPointOffset(i,j),
198     tempDataExp->getPointDataView(),
199     tempDataExp->getPointOffset(i,j),
200     region_loop_range);
201 jgs 82 }
202     }
203 jgs 102 }
204 jgs 82
205 jgs 102 void
206 matt 1319 DataExpanded::copy(const DataArrayView& value)
207 jgs 102 {
208     //
209 jgs 117 // copy a single value to every data point in this object
210 jgs 102 int nRows=m_data.getNumRows();
211     int nCols=m_data.getNumCols();
212 jgs 117 int i,j;
213 jgs 122 #pragma omp parallel for private(i,j) schedule(static)
214 jgs 117 for (i=0;i<nRows;i++) {
215     for (j=0;j<nCols;j++) {
216 matt 1319 // NOTE: An exception may be thown from this call if
217 jgs 102 // DOASSERT is on which of course will play
218     // havoc with the omp threads. Run single threaded
219 matt 1319 // if using DOASSERT.
220 matt 1323 getPointDataView().copy(getPointOffset(i,j),value);
221 jgs 82 }
222     }
223 jgs 102 }
224 jgs 82
225 jgs 102 void
226 matt 1319 DataExpanded::copy(const boost::python::numeric::array& value)
227 jgs 102 {
228 matt 1319
229     // extract the shape of the numarray
230     DataArrayView::ShapeType tempShape;
231     for (int i=0; i < value.getrank(); i++) {
232     tempShape.push_back(extract<int>(value.getshape()[i]));
233     }
234    
235     // get the space for the data vector
236     int len = DataArrayView::noValues(tempShape);
237     DataVector temp_data(len, 0.0, len);
238     DataArrayView temp_dataView(temp_data, tempShape);
239     temp_dataView.copy(value);
240    
241 jgs 102 //
242 jgs 117 // check the input shape matches this shape
243 matt 1319 if (!getPointDataView().checkShape(temp_dataView.getShape())) {
244 jgs 102 throw DataException(getPointDataView().createShapeErrorMessage(
245     "Error - (DataExpanded) Cannot copy due to shape mismatch.",
246 matt 1319 temp_dataView.getShape()));
247 jgs 82 }
248 jgs 102 //
249 jgs 117 // now copy over the data
250 matt 1319 copy(temp_dataView);
251 jgs 102 }
252 jgs 82
253 jgs 102 void
254     DataExpanded::initialise(const DataArrayView::ShapeType& shape,
255     int noSamples,
256     int noDataPointsPerSample)
257     {
258     //
259 jgs 117 // resize data array to the required size
260 jgs 102 m_data.resize(noSamples,noDataPointsPerSample,DataArrayView::noValues(shape));
261     //
262 jgs 117 // create the data view of the data array
263 jgs 102 DataArrayView temp(m_data.getData(),shape);
264     setPointDataView(temp);
265     }
266 jgs 82
267 jgs 102 string
268     DataExpanded::toString() const
269     {
270     stringstream temp;
271 gross 856 FunctionSpace fs=getFunctionSpace();
272 jgs 102 //
273     // create a temporary view as the offset will be changed
274     DataArrayView tempView(getPointDataView().getData(),getPointDataView().getShape(),getPointDataView().getOffset());
275 jgs 117 for (int i=0;i<m_data.getNumRows();i++) {
276     for (int j=0;j<m_data.getNumCols();j++) {
277 matt 1323 tempView.setOffset(getPointOffset(i,j));
278 jgs 102 stringstream suffix;
279 gross 964 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
280 jgs 102 temp << tempView.toString(suffix.str());
281     if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
282     temp << endl;
283 jgs 82 }
284     }
285     }
286 ksteube 1312 return temp.str();
287 jgs 102 }
288 jgs 82
289 jgs 102 DataArrayView::ValueType::size_type
290     DataExpanded::getPointOffset(int sampleNo,
291     int dataPointNo) const
292     {
293     return m_data.index(sampleNo,dataPointNo);
294     }
295 jgs 82
296 jgs 102 DataArrayView
297     DataExpanded::getDataPoint(int sampleNo,
298     int dataPointNo)
299     {
300 matt 1323 DataArrayView temp(m_data.getData(),getPointDataView().getShape(),getPointOffset(sampleNo,dataPointNo));
301 jgs 102 return temp;
302     }
303 jgs 82
304 jgs 102 DataArrayView::ValueType::size_type
305     DataExpanded::getLength() const
306     {
307 jgs 117 return m_data.size();
308 jgs 102 }
309 jgs 82
310 jgs 123 int
311     DataExpanded::archiveData(ofstream& archiveFile,
312     const DataArrayView::ValueType::size_type noValues) const
313     {
314     return(m_data.archiveData(archiveFile, noValues));
315     }
316    
317     int
318     DataExpanded::extractData(ifstream& archiveFile,
319     const DataArrayView::ValueType::size_type noValues)
320     {
321     return(m_data.extractData(archiveFile, noValues));
322     }
323    
324 matt 1319 void
325 gross 922 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
326     //
327     // Get the number of samples and data-points per sample.
328     int numSamples = getNumSamples();
329     int numDataPointsPerSample = getNumDPPSample();
330     int dataPointRank = getPointDataView().getRank();
331     ShapeType dataPointShape = getPointDataView().getShape();
332     if (numSamples*numDataPointsPerSample > 0) {
333     //TODO: global error handling
334 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
335 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
336     }
337 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
338 gross 922 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
339     }
340     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
341     if (dataPointRank==0) {
342     dataPointView()=value;
343     } else if (dataPointRank==1) {
344     for (int i=0; i<dataPointShape[0]; i++) {
345     dataPointView(i)=value;
346     }
347     } else if (dataPointRank==2) {
348     for (int i=0; i<dataPointShape[0]; i++) {
349     for (int j=0; j<dataPointShape[1]; j++) {
350     dataPointView(i,j)=value;
351     }
352     }
353     } else if (dataPointRank==3) {
354     for (int i=0; i<dataPointShape[0]; i++) {
355     for (int j=0; j<dataPointShape[1]; j++) {
356     for (int k=0; k<dataPointShape[2]; k++) {
357     dataPointView(i,j,k)=value;
358     }
359     }
360     }
361     } else if (dataPointRank==4) {
362     for (int i=0; i<dataPointShape[0]; i++) {
363     for (int j=0; j<dataPointShape[1]; j++) {
364     for (int k=0; k<dataPointShape[2]; k++) {
365     for (int l=0; l<dataPointShape[3]; l++) {
366     dataPointView(i,j,k,l)=value;
367     }
368     }
369     }
370     }
371 matt 1319 }
372 gross 922 }
373     }
374 matt 1319 void
375 gross 921 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const boost::python::numeric::array& value) {
376     //
377     // Get the number of samples and data-points per sample.
378     int numSamples = getNumSamples();
379     int numDataPointsPerSample = getNumDPPSample();
380     int dataPointRank = getPointDataView().getRank();
381     ShapeType dataPointShape = getPointDataView().getShape();
382     //
383     // check rank:
384     if (value.getrank()!=dataPointRank)
385     throw DataException("Rank of numarray does not match Data object rank");
386     if (numSamples*numDataPointsPerSample > 0) {
387     //TODO: global error handling
388 trankine 924 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
389 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
390     }
391 trankine 924 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
392 gross 921 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
393     }
394     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
395     if (dataPointRank==0) {
396     dataPointView()=extract<double>(value[0]);
397     } else if (dataPointRank==1) {
398     for (int i=0; i<dataPointShape[0]; i++) {
399     dataPointView(i)=extract<double>(value[i]);
400     }
401     } else if (dataPointRank==2) {
402     for (int i=0; i<dataPointShape[0]; i++) {
403     for (int j=0; j<dataPointShape[1]; j++) {
404     dataPointView(i,j)=extract<double>(value[i][j]);
405     }
406     }
407     } else if (dataPointRank==3) {
408     for (int i=0; i<dataPointShape[0]; i++) {
409     for (int j=0; j<dataPointShape[1]; j++) {
410     for (int k=0; k<dataPointShape[2]; k++) {
411     dataPointView(i,j,k)=extract<double>(value[i][j][k]);
412     }
413     }
414     }
415     } else if (dataPointRank==4) {
416     for (int i=0; i<dataPointShape[0]; i++) {
417     for (int j=0; j<dataPointShape[1]; j++) {
418     for (int k=0; k<dataPointShape[2]; k++) {
419     for (int l=0; l<dataPointShape[3]; l++) {
420     dataPointView(i,j,k,l)=extract<double>(value[i][j][k][l]);
421     }
422     }
423     }
424     }
425 matt 1319 }
426 gross 921 }
427     }
428 jgs 126 void
429     DataExpanded::copyAll(const boost::python::numeric::array& value) {
430     //
431     // Get the number of samples and data-points per sample.
432     int numSamples = getNumSamples();
433     int numDataPointsPerSample = getNumDPPSample();
434     int dataPointRank = getPointDataView().getRank();
435     ShapeType dataPointShape = getPointDataView().getShape();
436     //
437     // check rank:
438     if (value.getrank()!=dataPointRank+1)
439     throw DataException("Rank of numarray does not match Data object rank");
440     if (value.getshape()[0]!=numSamples*numDataPointsPerSample)
441     throw DataException("leading dimension of numarray is too small");
442     //
443     int dataPoint = 0;
444     for (int sampleNo = 0; sampleNo < numSamples; sampleNo++) {
445     for (int dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
446     DataArrayView dataPointView = getDataPoint(sampleNo, dataPointNo);
447     if (dataPointRank==0) {
448     dataPointView()=extract<double>(value[dataPoint]);
449     } else if (dataPointRank==1) {
450     for (int i=0; i<dataPointShape[0]; i++) {
451     dataPointView(i)=extract<double>(value[dataPoint][i]);
452     }
453     } else if (dataPointRank==2) {
454     for (int i=0; i<dataPointShape[0]; i++) {
455     for (int j=0; j<dataPointShape[1]; j++) {
456     dataPointView(i,j)=extract<double>(value[dataPoint][i][j]);
457     }
458     }
459     } else if (dataPointRank==3) {
460     for (int i=0; i<dataPointShape[0]; i++) {
461     for (int j=0; j<dataPointShape[1]; j++) {
462     for (int k=0; k<dataPointShape[2]; k++) {
463     dataPointView(i,j,k)=extract<double>(value[dataPoint][i][j][k]);
464     }
465     }
466     }
467     } else if (dataPointRank==4) {
468     for (int i=0; i<dataPointShape[0]; i++) {
469     for (int j=0; j<dataPointShape[1]; j++) {
470     for (int k=0; k<dataPointShape[2]; k++) {
471     for (int l=0; l<dataPointShape[3]; l++) {
472     dataPointView(i,j,k,l)=extract<double>(value[dataPoint][i][j][k][l]);
473     }
474     }
475     }
476     }
477     }
478     dataPoint++;
479     }
480     }
481     }
482 gross 580 void
483 ksteube 775 DataExpanded::symmetric(DataAbstract* ev)
484     {
485     int sampleNo,dataPointNo;
486     int numSamples = getNumSamples();
487     int numDataPointsPerSample = getNumDPPSample();
488     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
489     if (temp_ev==0) {
490     throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
491     }
492     DataArrayView& thisView=getPointDataView();
493     DataArrayView& evView=ev->getPointDataView();
494     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
495     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
496     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
497     DataArrayView::symmetric(thisView,getPointOffset(sampleNo,dataPointNo),
498     evView,ev->getPointOffset(sampleNo,dataPointNo));
499     }
500     }
501     }
502     void
503     DataExpanded::nonsymmetric(DataAbstract* ev)
504     {
505     int sampleNo,dataPointNo;
506     int numSamples = getNumSamples();
507     int numDataPointsPerSample = getNumDPPSample();
508     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
509     if (temp_ev==0) {
510     throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
511     }
512     DataArrayView& thisView=getPointDataView();
513     DataArrayView& evView=ev->getPointDataView();
514     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
515     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
516     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
517     DataArrayView::nonsymmetric(thisView,getPointOffset(sampleNo,dataPointNo),
518     evView,ev->getPointOffset(sampleNo,dataPointNo));
519     }
520     }
521     }
522     void
523 gross 800 DataExpanded::trace(DataAbstract* ev, int axis_offset)
524 ksteube 775 {
525     int sampleNo,dataPointNo;
526     int numSamples = getNumSamples();
527     int numDataPointsPerSample = getNumDPPSample();
528     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
529     if (temp_ev==0) {
530 gross 800 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
531 ksteube 775 }
532     DataArrayView& thisView=getPointDataView();
533     DataArrayView& evView=ev->getPointDataView();
534     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
535     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
536     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
537 gross 800 DataArrayView::trace(thisView,getPointOffset(sampleNo,dataPointNo),
538 ksteube 775 evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
539     }
540     }
541     }
542 gross 800
543 ksteube 775 void
544     DataExpanded::transpose(DataAbstract* ev, int axis_offset)
545     {
546     int sampleNo,dataPointNo;
547     int numSamples = getNumSamples();
548     int numDataPointsPerSample = getNumDPPSample();
549     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
550     if (temp_ev==0) {
551     throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
552     }
553     DataArrayView& thisView=getPointDataView();
554     DataArrayView& evView=ev->getPointDataView();
555     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
556     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
557     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
558     DataArrayView::transpose(thisView,getPointOffset(sampleNo,dataPointNo),
559     evView,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
560     }
561     }
562     }
563 gross 800
564 ksteube 775 void
565 gross 804 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
566 gross 800 {
567     int sampleNo,dataPointNo;
568     int numSamples = getNumSamples();
569     int numDataPointsPerSample = getNumDPPSample();
570     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
571     if (temp_ev==0) {
572 gross 804 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
573 gross 800 }
574     DataArrayView& thisView=getPointDataView();
575     DataArrayView& evView=ev->getPointDataView();
576     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
577     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
578     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
579 gross 804 DataArrayView::swapaxes(thisView,getPointOffset(sampleNo,dataPointNo),
580     evView,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
581 gross 800 }
582     }
583     }
584     void
585 gross 580 DataExpanded::eigenvalues(DataAbstract* ev)
586     {
587 gross 584 int sampleNo,dataPointNo;
588 gross 580 int numSamples = getNumSamples();
589     int numDataPointsPerSample = getNumDPPSample();
590     DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
591     if (temp_ev==0) {
592     throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
593     }
594     DataArrayView& thisView=getPointDataView();
595     DataArrayView& evView=ev->getPointDataView();
596     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
597 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
598     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
599 gross 580 DataArrayView::eigenvalues(thisView,getPointOffset(sampleNo,dataPointNo),
600     evView,ev->getPointOffset(sampleNo,dataPointNo));
601     }
602     }
603     }
604     void
605     DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
606     {
607     int numSamples = getNumSamples();
608     int numDataPointsPerSample = getNumDPPSample();
609 gross 584 int sampleNo,dataPointNo;
610 gross 580 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
611     if (temp_ev==0) {
612     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
613     }
614     DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
615     if (temp_V==0) {
616     throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
617     }
618     DataArrayView& thisView=getPointDataView();
619     DataArrayView& evView=ev->getPointDataView();
620     DataArrayView& VView=V->getPointDataView();
621     #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
622 gross 584 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
623     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
624 gross 580 DataArrayView::eigenvalues_and_eigenvectors(thisView,getPointOffset(sampleNo,dataPointNo),
625     evView,ev->getPointOffset(sampleNo,dataPointNo),
626     VView,V->getPointOffset(sampleNo,dataPointNo),
627     tol);
628     }
629     }
630     }
631    
632 gross 950 void
633 gross 1118 DataExpanded::setToZero(){
634     int numSamples = getNumSamples();
635     int numDataPointsPerSample = getNumDPPSample();
636     DataArrayView& thisView=getPointDataView();
637     DataArrayView::ValueType::size_type n = thisView.noValues();
638     double* p;
639     int sampleNo,dataPointNo, i;
640     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
641     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
642     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
643     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
644 phornby 1455 for (i=0; i<n ;++i) p[i]=0.;
645 gross 1118 }
646     }
647     }
648    
649    
650     void
651 gross 950 DataExpanded::dump(const std::string fileName) const
652     {
653 gross 967 #ifdef PASO_MPI
654     throw DataException("Error - DataExpanded:: dump is not implemented for MPI yet.");
655     #endif
656 gross 1023 #ifdef USE_NETCDF
657 gross 967 const int ldims=2+DataArrayView::maxRank;
658     const NcDim* ncdims[ldims];
659     NcVar *var, *ids;
660     int rank = getPointDataView().getRank();
661     int type= getFunctionSpace().getTypeCode();
662     int ndims =0;
663     long dims[ldims];
664 gross 1141 const double* d_ptr=&(m_data[0]);
665 gross 967 DataArrayView::ShapeType shape = getPointDataView().getShape();
666 gross 965
667 gross 967 // netCDF error handler
668     NcError err(NcError::verbose_nonfatal);
669     // Create the file.
670     NcFile dataFile(fileName.c_str(), NcFile::Replace);
671     // check if writing was successful
672     if (!dataFile.is_valid())
673     throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
674 gross 1141 if (!dataFile.add_att("type_id",2) )
675 gross 967 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
676     if (!dataFile.add_att("rank",rank) )
677     throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
678     if (!dataFile.add_att("function_space_type",type))
679     throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
680     ndims=rank+2;
681     if ( rank >0 ) {
682     dims[0]=shape[0];
683     if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
684     throw DataException("Error - DataExpanded:: appending ncdimsion 0 to netCDF file failed.");
685     }
686     if ( rank >1 ) {
687     dims[1]=shape[1];
688     if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
689     throw DataException("Error - DataExpanded:: appending ncdimsion 1 to netCDF file failed.");
690     }
691     if ( rank >2 ) {
692     dims[2]=shape[2];
693     if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
694     throw DataException("Error - DataExpanded:: appending ncdimsion 2 to netCDF file failed.");
695     }
696     if ( rank >3 ) {
697     dims[3]=shape[3];
698     if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
699     throw DataException("Error - DataExpanded:: appending ncdimsion 3 to netCDF file failed.");
700     }
701     dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
702     if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
703     throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
704     dims[rank+1]=getFunctionSpace().getNumSamples();
705     if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
706     throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
707    
708     if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
709     throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
710     const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
711     if (! (ids->put(ids_p,dims[rank+1])) )
712     throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
713    
714     if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
715     throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
716 gross 1141 if (! (var->put(d_ptr,dims)) )
717 gross 967 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
718 gross 1023 #else
719     throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
720     #endif
721 gross 950 }
722 gross 1358
723     void
724     DataExpanded::setTaggedValue(int tagKey,
725     const DataArrayView& value)
726     {
727     int numSamples = getNumSamples();
728     int numDataPointsPerSample = getNumDPPSample();
729     int sampleNo,dataPointNo, i;
730     DataArrayView& thisView=getPointDataView();
731     DataArrayView::ValueType::size_type n = thisView.noValues();
732     double* p,*in=&(value.getData()[0]);
733    
734     if (value.noValues() != n) {
735     throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
736     }
737    
738     #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
739     for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
740     if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
741     for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
742     p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
743 phornby 1455 for (i=0; i<n ;++i) p[i]=in[i];
744 gross 1358 }
745     }
746     }
747     }
748    
749 gross 1487 void
750     DataExpanded::reorderByReferenceIDs(int *reference_ids)
751     {
752     int numSamples = getNumSamples();
753     DataArrayView& thisView=getPointDataView();
754     DataArrayView::ValueType::size_type n = thisView.noValues() * getNumDPPSample();
755     int sampleNo, sampleNo2,i;
756     double* p,*p2;
757     register double rtmp;
758     FunctionSpace fs=getFunctionSpace();
759 gross 1358
760 gross 1487 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
761     const int id_in=reference_ids[sampleNo];
762     const int id=fs.getReferenceIDOfSample(sampleNo);
763     if (id!=id_in) {
764     bool matched=false;
765     for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
766     if (id == reference_ids[sampleNo2]) {
767     p=&(m_data[getPointOffset(sampleNo,0)]);
768     p2=&(m_data[getPointOffset(sampleNo2,0)]);
769 gross 1513 for (i=0; i<n ;i++) {
770 gross 1487 rtmp=p[i];
771     p[i]=p2[i];
772     p2[i]=rtmp;
773     }
774     reference_ids[sampleNo]=id;
775     reference_ids[sampleNo2]=id_in;
776     matched=true;
777     break;
778     }
779     }
780 phornby 1628 if (! matched) {
781 gross 1487 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
782     }
783     }
784     }
785     }
786    
787    
788 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