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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3468 - (show annotations)
Tue Feb 22 06:38:57 2011 UTC (8 years, 8 months ago) by jfenwick
File size: 29102 byte(s)
Made constructors more efficient.
Only has a very minor affect on performance
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 "DataExpanded.h"
17 #include "DataException.h"
18 #include "DataConstant.h"
19 #include "DataTagged.h"
20
21 #include "esysUtils/Esys_MPI.h"
22
23 #ifdef USE_NETCDF
24 #include <netcdfcpp.h>
25 #endif
26
27 #include <boost/python/extract.hpp>
28 #include "DataMaths.h"
29
30 using namespace std;
31 using namespace boost::python;
32 using namespace boost;
33 using namespace escript::DataTypes;
34
35
36 // #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {throw DataException("Attempt to modify shared object");}
37
38 #define CHECK_FOR_EX_WRITE if (!checkNoSharing()) {std::ostringstream ss; ss << " Attempt to modify shared object. line " << __LINE__ << " of " << __FILE__; *((int*)0)=17;throw DataException(ss.str());}
39
40 namespace escript {
41
42 DataExpanded::DataExpanded(const WrappedArray& value,
43 const FunctionSpace& what)
44 : parent(what,value.getShape())
45 {
46 //
47 // initialise the data array for this object
48 initialise(what.getNumSamples(),what.getNumDPPSample());
49 //
50 // copy the given value to every data point
51 copy(value);
52 }
53
54 DataExpanded::DataExpanded(const DataExpanded& other)
55 : parent(other.getFunctionSpace(), other.getShape()),
56 m_data(other.m_data)
57 {
58 }
59
60 DataExpanded::DataExpanded(const DataConstant& other)
61 : parent(other.getFunctionSpace(), other.getShape())
62 {
63 //
64 // initialise the data array for this object
65 initialise(other.getNumSamples(),other.getNumDPPSample());
66 //
67 // DataConstant only has one value, copy this to every data point
68 copy(other);
69 }
70
71 DataExpanded::DataExpanded(const DataTagged& other)
72 : parent(other.getFunctionSpace(), other.getShape())
73 {
74 //
75 // initialise the data array for this object
76 initialise(other.getNumSamples(),other.getNumDPPSample());
77 //
78 // for each data point in this object, extract and copy the corresponding data
79 // value from the given DataTagged object
80 int i,j;
81 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
82 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
83 #pragma omp parallel for private(i,j) schedule(static)
84 for (i=0;i<numRows;i++) {
85 for (j=0;j<numCols;j++) {
86 try {
87 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(),
88 other.getVectorRO(),
89 other.getPointOffset(i,j));
90 }
91 catch (std::exception& e) {
92 cout << e.what() << endl;
93 }
94 }
95 }
96 }
97
98 DataExpanded::DataExpanded(const DataExpanded& other,
99 const DataTypes::RegionType& region)
100 : parent(other.getFunctionSpace(),DataTypes::getResultSliceShape(region))
101 {
102 //
103 // get the shape of the slice
104 // DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
105 //
106 // initialise this Data object to the shape of the slice
107 initialise(other.getNumSamples(),other.getNumDPPSample());
108 //
109 // copy the data
110 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
111 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
112 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
113 int i,j;
114 #pragma omp parallel for private(i,j) schedule(static)
115 for (i=0;i<numRows;i++) {
116 for (j=0;j<numCols;j++) {
117 try {
118 DataTypes::copySlice(getVectorRW(),getShape(),getPointOffset(i,j),
119 other.getVectorRO(),
120 other.getShape(),
121 other.getPointOffset(i,j),
122 region_loop_range);
123 }
124 catch (std::exception& e) {
125 cout << e.what() << endl;
126 }
127 }
128 }
129 }
130
131 DataExpanded::DataExpanded(const FunctionSpace& what,
132 const DataTypes::ShapeType &shape,
133 const DataTypes::ValueType &data)
134 : parent(what,shape)
135 {
136 EsysAssert(data.size()%getNoValues()==0,
137 "DataExpanded Constructor - size of supplied data is not a multiple of shape size.");
138
139 if (data.size()==getNoValues())
140 {
141 ValueType& vec=m_data.getData();
142 //
143 // create the view of the data
144 initialise(what.getNumSamples(),what.getNumDPPSample());
145 // now we copy this value to all elements
146 for (int i=0;i<getLength();)
147 {
148 for (unsigned int j=0;j<getNoValues();++j,++i)
149 {
150 vec[i]=data[j];
151 }
152 }
153 }
154 else
155 {
156 //
157 // copy the data in the correct format
158 m_data.getData()=data;
159 }
160
161
162 }
163
164 DataExpanded::DataExpanded(const FunctionSpace& what,
165 const DataTypes::ShapeType &shape,
166 const double v)
167 : parent(what,shape)
168 {
169 ValueType& vec=m_data.getData();
170 //
171 // create the view of the data
172 initialise(what.getNumSamples(),what.getNumDPPSample());
173 // now we copy this value to all elements
174 const int L=getLength();
175 int i;
176 #pragma omp parallel for schedule(static) private(i)
177 for (i=0;i<L;++i)
178 {
179 vec[i]=v;
180 }
181 }
182
183
184
185 DataExpanded::~DataExpanded()
186 {
187 }
188
189 DataAbstract*
190 DataExpanded::deepCopy()
191 {
192 return new DataExpanded(*this);
193 }
194
195
196 DataAbstract*
197 DataExpanded::getSlice(const DataTypes::RegionType& region) const
198 {
199 return new DataExpanded(*this,region);
200 }
201
202 void
203 DataExpanded::setSlice(const DataAbstract* value,
204 const DataTypes::RegionType& region)
205 {
206 const DataExpanded* tempDataExp=dynamic_cast<const DataExpanded*>(value);
207 if (tempDataExp==0) {
208 throw DataException("Programming error - casting to DataExpanded.");
209 }
210 CHECK_FOR_EX_WRITE
211 //
212 // get shape of slice
213 DataTypes::ShapeType shape(DataTypes::getResultSliceShape(region));
214 DataTypes::RegionLoopRangeType region_loop_range=DataTypes::getSliceRegionLoopRange(region);
215 //
216 // check shape
217 if (getRank()!=region.size()) {
218 throw DataException("Error - Invalid slice region.");
219 }
220 if (tempDataExp->getRank()>0 && !DataTypes::checkShape(value->getShape(), shape)) {
221 throw DataException (DataTypes::createShapeErrorMessage(
222 "Error - Couldn't copy slice due to shape mismatch.",shape, value->getShape()));
223 }
224 //
225 // copy the data from the slice into this object
226 DataTypes::ValueType::size_type numRows=m_data.getNumRows();
227 DataTypes::ValueType::size_type numCols=m_data.getNumCols();
228 int i, j;
229 ValueType& vec=getVectorRW();
230 const ShapeType& mshape=getShape();
231 const ValueType& tVec=tempDataExp->getVectorRO();
232 const ShapeType& tShape=tempDataExp->getShape();
233 #pragma omp parallel for private(i,j) schedule(static)
234 for (i=0;i<numRows;i++) {
235 for (j=0;j<numCols;j++) {
236 DataTypes::copySliceFrom(vec,mshape,getPointOffset(i,j),
237 tVec,
238 tShape,
239 tempDataExp->getPointOffset(i,j),
240 region_loop_range);
241
242 }
243 }
244 }
245
246 void
247 DataExpanded::copy(const DataConstant& value)
248 {
249 EsysAssert((checkShape(getShape(), value.getShape())),
250 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",value.getShape(),getShape()));
251
252 //
253 // copy a single value to every data point in this object
254 int nRows=m_data.getNumRows();
255 int nCols=m_data.getNumCols();
256 int i,j;
257 #pragma omp parallel for private(i,j) schedule(static)
258 for (i=0;i<nRows;i++) {
259 for (j=0;j<nCols;j++) {
260 DataTypes::copyPoint(getVectorRW(), getPointOffset(i,j), getNoValues(), value.getVectorRO(), 0);
261 }
262 }
263 }
264
265 void
266 DataExpanded::copy(const WrappedArray& value)
267 {
268 // check the input shape matches this shape
269 if (!DataTypes::checkShape(getShape(),value.getShape())) {
270 throw DataException(DataTypes::createShapeErrorMessage(
271 "Error - (DataExpanded) Cannot copy due to shape mismatch.",
272 value.getShape(),getShape()));
273 }
274 getVectorRW().copyFromArray(value, getNumDPPSample()*getNumSamples());
275 }
276
277
278 void
279 DataExpanded::initialise(int noSamples,
280 int noDataPointsPerSample)
281 {
282 if (noSamples==0) //retain the default empty object
283 {
284 return;
285 }
286 //
287 // resize data array to the required size
288 m_data.resize(noSamples,noDataPointsPerSample,getNoValues());
289 }
290
291 bool
292 DataExpanded::hasNaN() const
293 {
294 const ValueType& v=m_data.getData();
295 for (ValueType::size_type i=0;i<v.size();++i)
296 {
297 if (nancheck(v[i]))
298 {
299 return true;
300 }
301 }
302 return false;
303 }
304
305
306 string
307 DataExpanded::toString() const
308 {
309 stringstream temp;
310 FunctionSpace fs=getFunctionSpace();
311
312 int offset=0;
313 for (int i=0;i<m_data.getNumRows();i++) {
314 for (int j=0;j<m_data.getNumCols();j++) {
315 offset=getPointOffset(i,j);
316 stringstream suffix;
317 suffix << "( id: " << i << ", ref: " << fs.getReferenceIDOfSample(i) << ", pnt: " << j << ")";
318 temp << DataTypes::pointToString(getVectorRO(),getShape(),offset,suffix.str());
319 if (!(i==(m_data.getNumRows()-1) && j==(m_data.getNumCols()-1))) {
320 temp << endl;
321 }
322 }
323 }
324 string result=temp.str();
325 if (result.empty())
326 {
327 return "(data contains no samples)\n";
328 }
329 return temp.str();
330 }
331
332 DataTypes::ValueType::size_type
333 DataExpanded::getPointOffset(int sampleNo,
334 int dataPointNo) const
335 {
336 return m_data.index(sampleNo,dataPointNo);
337 }
338
339 DataTypes::ValueType::size_type
340 DataExpanded::getPointOffset(int sampleNo,
341 int dataPointNo)
342 {
343 return m_data.index(sampleNo,dataPointNo);
344 }
345
346 DataTypes::ValueType::size_type
347 DataExpanded::getLength() const
348 {
349 return m_data.size();
350 }
351
352
353
354 void
355 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const double value) {
356 CHECK_FOR_EX_WRITE
357 //
358 // Get the number of samples and data-points per sample.
359 int numSamples = getNumSamples();
360 int numDataPointsPerSample = getNumDPPSample();
361 int dataPointRank = getRank();
362 ShapeType dataPointShape = getShape();
363 if (numSamples*numDataPointsPerSample > 0) {
364 //TODO: global error handling
365 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
366 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
367 }
368 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
369 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
370 }
371 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
372 ValueType& vec=getVectorRW();
373 if (dataPointRank==0) {
374 vec[offset]=value;
375 } else if (dataPointRank==1) {
376 for (int i=0; i<dataPointShape[0]; i++) {
377 vec[offset+i]=value;
378 }
379 } else if (dataPointRank==2) {
380 for (int i=0; i<dataPointShape[0]; i++) {
381 for (int j=0; j<dataPointShape[1]; j++) {
382 vec[offset+getRelIndex(dataPointShape,i,j)]=value;
383 }
384 }
385 } else if (dataPointRank==3) {
386 for (int i=0; i<dataPointShape[0]; i++) {
387 for (int j=0; j<dataPointShape[1]; j++) {
388 for (int k=0; k<dataPointShape[2]; k++) {
389 vec[offset+getRelIndex(dataPointShape,i,j,k)]=value;
390 }
391 }
392 }
393 } else if (dataPointRank==4) {
394 for (int i=0; i<dataPointShape[0]; i++) {
395 for (int j=0; j<dataPointShape[1]; j++) {
396 for (int k=0; k<dataPointShape[2]; k++) {
397 for (int l=0; l<dataPointShape[3]; l++) {
398 vec[offset+getRelIndex(dataPointShape,i,j,k,l)]=value;
399 }
400 }
401 }
402 }
403 }
404 }
405 }
406
407 void
408 DataExpanded::copyToDataPoint(const int sampleNo, const int dataPointNo, const WrappedArray& value) {
409 CHECK_FOR_EX_WRITE
410 //
411 // Get the number of samples and data-points per sample.
412 int numSamples = getNumSamples();
413 int numDataPointsPerSample = getNumDPPSample();
414 //
415 // check rank:
416 if (value.getRank()!=getRank())
417 throw DataException("Rank of value does not match Data object rank");
418 if (numSamples*numDataPointsPerSample > 0) {
419 //TODO: global error handling
420 if ((sampleNo >= numSamples) || (sampleNo < 0 )) {
421 throw DataException("Error - DataExpanded::copyDataPoint invalid sampleNo.");
422 }
423 if ((dataPointNo >= numDataPointsPerSample) || (dataPointNo < 0)) {
424 throw DataException("Error - DataExpanded::copyDataPoint invalid dataPointNoInSample.");
425 }
426 ValueType::size_type offset = getPointOffset(sampleNo, dataPointNo);
427 ValueType& vec=getVectorRW();
428 vec.copyFromArrayToOffset(value,offset,1);
429 }
430 }
431
432 void
433 DataExpanded::symmetric(DataAbstract* ev)
434 {
435 int sampleNo,dataPointNo;
436 int numSamples = getNumSamples();
437 int numDataPointsPerSample = getNumDPPSample();
438 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
439 if (temp_ev==0) {
440 throw DataException("Error - DataExpanded::symmetric: casting to DataExpanded failed (propably a programming error).");
441 }
442 const ValueType& vec=getVectorRO();
443 const ShapeType& shape=getShape();
444 ValueType& evVec=temp_ev->getVectorRW();
445 const ShapeType& evShape=temp_ev->getShape();
446 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
447 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
448 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
449 DataMaths::symmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
450 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
451 }
452 }
453 }
454 void
455 DataExpanded::nonsymmetric(DataAbstract* ev)
456 {
457 int sampleNo,dataPointNo;
458 int numSamples = getNumSamples();
459 int numDataPointsPerSample = getNumDPPSample();
460 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
461 if (temp_ev==0) {
462 throw DataException("Error - DataExpanded::nonsymmetric: casting to DataExpanded failed (propably a programming error).");
463 }
464 const ValueType& vec=getVectorRO();
465 const ShapeType& shape=getShape();
466 ValueType& evVec=temp_ev->getVectorRW();
467 const ShapeType& evShape=temp_ev->getShape();
468 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
469 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
470 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
471 DataMaths::nonsymmetric(vec,shape,getPointOffset(sampleNo,dataPointNo),
472 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
473 }
474 }
475 }
476 void
477 DataExpanded::trace(DataAbstract* ev, int axis_offset)
478 {
479 int sampleNo,dataPointNo;
480 int numSamples = getNumSamples();
481 int numDataPointsPerSample = getNumDPPSample();
482 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
483 if (temp_ev==0) {
484 throw DataException("Error - DataExpanded::trace: casting to DataExpanded failed (propably a programming error).");
485 }
486 const ValueType& vec=getVectorRO();
487 const ShapeType& shape=getShape();
488 ValueType& evVec=temp_ev->getVectorRW();
489 const ShapeType& evShape=temp_ev->getShape();
490 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
491 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
492 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
493 DataMaths::trace(vec,shape,getPointOffset(sampleNo,dataPointNo),
494 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
495 }
496 }
497 }
498
499 void
500 DataExpanded::transpose(DataAbstract* ev, int axis_offset)
501 {
502 int sampleNo,dataPointNo;
503 int numSamples = getNumSamples();
504 int numDataPointsPerSample = getNumDPPSample();
505 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
506 if (temp_ev==0) {
507 throw DataException("Error - DataExpanded::transpose: casting to DataExpanded failed (propably a programming error).");
508 }
509 const ValueType& vec=getVectorRO();
510 const ShapeType& shape=getShape();
511 ValueType& evVec=temp_ev->getVectorRW();
512 const ShapeType& evShape=temp_ev->getShape();
513 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
514 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
515 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
516 DataMaths::transpose(vec,shape,getPointOffset(sampleNo,dataPointNo),
517 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis_offset);
518 }
519 }
520 }
521
522 void
523 DataExpanded::swapaxes(DataAbstract* ev, int axis0, int axis1)
524 {
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 throw DataException("Error - DataExpanded::swapaxes: casting to DataExpanded failed (propably a programming error).");
531 }
532 const ValueType& vec=getVectorRO();
533 const ShapeType& shape=getShape();
534 ValueType& evVec=temp_ev->getVectorRW();
535 const ShapeType& evShape=temp_ev->getShape();
536 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
537 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
538 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
539 DataMaths::swapaxes(vec,shape,getPointOffset(sampleNo,dataPointNo),
540 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),axis0,axis1);
541 }
542 }
543 }
544 void
545 DataExpanded::eigenvalues(DataAbstract* ev)
546 {
547 int sampleNo,dataPointNo;
548 int numSamples = getNumSamples();
549 int numDataPointsPerSample = getNumDPPSample();
550 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
551 if (temp_ev==0) {
552 throw DataException("Error - DataExpanded::eigenvalues: casting to DataExpanded failed (propably a programming error).");
553 }
554 const ValueType& vec=getVectorRO();
555 const ShapeType& shape=getShape();
556 ValueType& evVec=temp_ev->getVectorRW();
557 const ShapeType& evShape=temp_ev->getShape();
558 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
559 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
560 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
561 DataMaths::eigenvalues(vec,shape,getPointOffset(sampleNo,dataPointNo),
562 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo));
563 }
564 }
565 }
566 void
567 DataExpanded::eigenvalues_and_eigenvectors(DataAbstract* ev,DataAbstract* V,const double tol)
568 {
569 int numSamples = getNumSamples();
570 int numDataPointsPerSample = getNumDPPSample();
571 int sampleNo,dataPointNo;
572 DataExpanded* temp_ev=dynamic_cast<DataExpanded*>(ev);
573 if (temp_ev==0) {
574 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
575 }
576 DataExpanded* temp_V=dynamic_cast<DataExpanded*>(V);
577 if (temp_V==0) {
578 throw DataException("Error - DataExpanded::eigenvalues_and_eigenvectors: casting to DataExpanded failed (propably a programming error).");
579 }
580 const ValueType& vec=getVectorRO();
581 const ShapeType& shape=getShape();
582 ValueType& evVec=temp_ev->getVectorRW();
583 const ShapeType& evShape=temp_ev->getShape();
584 ValueType& VVec=temp_V->getVectorRW();
585 const ShapeType& VShape=temp_V->getShape();
586 #pragma omp parallel for private(sampleNo,dataPointNo) schedule(static)
587 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
588 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
589 DataMaths::eigenvalues_and_eigenvectors(vec,shape,getPointOffset(sampleNo,dataPointNo),
590 evVec,evShape,ev->getPointOffset(sampleNo,dataPointNo),
591 VVec, VShape,V->getPointOffset(sampleNo,dataPointNo),tol);
592 }
593 }
594 }
595
596
597 int
598 DataExpanded::matrixInverse(DataAbstract* out) const
599 {
600 DataExpanded* temp=dynamic_cast<DataExpanded*>(out);
601 if (temp==0)
602 {
603 throw DataException("Error - DataExpanded::matrixInverse: casting to DataExpanded failed (propably a programming error).");
604 }
605
606 if (getRank()!=2)
607 {
608 throw DataException("Error - DataExpanded::matrixInverse: input must be rank 2.");
609
610 }
611 int sampleNo;
612 const int numdpps=getNumDPPSample();
613 const int numSamples = getNumSamples();
614 const ValueType& vec=m_data.getData();
615 int errcode=0;
616 #pragma omp parallel private(sampleNo)
617 {
618 int errorcode=0;
619 LapackInverseHelper h(getShape()[0]);
620 #pragma omp for schedule(static)
621 for (sampleNo = 0; sampleNo < numSamples; sampleNo++)
622 {
623 // not sure I like all those virtual calls to getPointOffset
624 DataTypes::ValueType::size_type offset=getPointOffset(sampleNo,0);
625 int res=DataMaths::matrix_inverse(vec, getShape(), offset, temp->getVectorRW(), temp->getShape(), offset, numdpps, h);
626 if (res>errorcode)
627 {
628 errorcode=res;
629 #pragma omp critical
630 {
631 errcode=errorcode; // I'm not especially concerned which error gets reported as long as one is
632 }
633 }
634 }
635 }
636 return errcode;
637 if (errcode)
638 {
639 DataMaths::matrixInverseError(errcode); // throws exceptions
640 }
641 }
642
643 void
644 DataExpanded::setToZero(){
645 // TODO: Surely there is a more efficient way to do this????
646 // Why is there no memset here? Parallel issues?
647 // A: This ensures that memory is touched by the correct thread.
648 CHECK_FOR_EX_WRITE
649 int numSamples = getNumSamples();
650 int numDataPointsPerSample = getNumDPPSample();
651 DataTypes::ValueType::size_type n = getNoValues();
652 double* p;
653 int sampleNo,dataPointNo, i;
654 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
655 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
656 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
657 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
658 for (i=0; i<n ;++i) p[i]=0.;
659 }
660 }
661 }
662
663 /* Append MPI rank to file name if multiple MPI processes */
664 char *Escript_MPI_appendRankToFileName(const char *fileName, int mpi_size, int mpi_rank) {
665 /* Make plenty of room for the mpi_rank number and terminating '\0' */
666 char *newFileName = (char *)malloc(strlen(fileName)+20);
667 strncpy(newFileName, fileName, strlen(fileName)+1);
668 if (mpi_size>1) sprintf(newFileName+strlen(newFileName), ".%04d", mpi_rank);
669 return(newFileName);
670 }
671
672 void
673 DataExpanded::dump(const std::string fileName) const
674 {
675 #ifdef USE_NETCDF
676 const int ldims=2+DataTypes::maxRank;
677 const NcDim* ncdims[ldims];
678 NcVar *var, *ids;
679 int rank = getRank();
680 int type= getFunctionSpace().getTypeCode();
681 int ndims =0;
682 long dims[ldims];
683 const double* d_ptr=&(m_data[0]);
684 const DataTypes::ShapeType& shape = getShape();
685 int mpi_iam=getFunctionSpace().getDomain()->getMPIRank();
686 int mpi_num=getFunctionSpace().getDomain()->getMPISize();
687 #ifdef ESYS_MPI
688 MPI_Status status;
689 #endif
690
691 #ifdef ESYS_MPI
692 /* Serialize NetCDF I/O */
693 if (mpi_iam>0) MPI_Recv(&ndims, 0, MPI_INT, mpi_iam-1, 81801, MPI_COMM_WORLD, &status);
694 #endif
695
696 // netCDF error handler
697 NcError err(NcError::verbose_nonfatal);
698 // Create the file.
699 char *newFileName = Escript_MPI_appendRankToFileName(fileName.c_str(), mpi_num, mpi_iam);
700 NcFile dataFile(newFileName, NcFile::Replace);
701 // check if writing was successful
702 if (!dataFile.is_valid())
703 throw DataException("Error - DataExpanded:: opening of netCDF file for output failed.");
704 if (!dataFile.add_att("type_id",2) )
705 throw DataException("Error - DataExpanded:: appending data type to netCDF file failed.");
706 if (!dataFile.add_att("rank",rank) )
707 throw DataException("Error - DataExpanded:: appending rank attribute to netCDF file failed.");
708 if (!dataFile.add_att("function_space_type",type))
709 throw DataException("Error - DataExpanded:: appending function space attribute to netCDF file failed.");
710 ndims=rank+2;
711 if ( rank >0 ) {
712 dims[0]=shape[0];
713 if (! (ncdims[0] = dataFile.add_dim("d0",shape[0])) )
714 throw DataException("Error - DataExpanded:: appending ncdimension 0 to netCDF file failed.");
715 }
716 if ( rank >1 ) {
717 dims[1]=shape[1];
718 if (! (ncdims[1] = dataFile.add_dim("d1",shape[1])) )
719 throw DataException("Error - DataExpanded:: appending ncdimension 1 to netCDF file failed.");
720 }
721 if ( rank >2 ) {
722 dims[2]=shape[2];
723 if (! (ncdims[2] = dataFile.add_dim("d2", shape[2])) )
724 throw DataException("Error - DataExpanded:: appending ncdimension 2 to netCDF file failed.");
725 }
726 if ( rank >3 ) {
727 dims[3]=shape[3];
728 if (! (ncdims[3] = dataFile.add_dim("d3", shape[3])) )
729 throw DataException("Error - DataExpanded:: appending ncdimension 3 to netCDF file failed.");
730 }
731 dims[rank]=getFunctionSpace().getNumDataPointsPerSample();
732 if (! (ncdims[rank] = dataFile.add_dim("num_data_points_per_sample", dims[rank])) )
733 throw DataException("Error - DataExpanded:: appending num_data_points_per_sample to netCDF file failed.");
734 dims[rank+1]=getFunctionSpace().getNumSamples();
735 if (! (ncdims[rank+1] = dataFile.add_dim("num_samples", dims[rank+1])) )
736 throw DataException("Error - DataExpanded:: appending num_sample to netCDF file failed.");
737
738 if (getFunctionSpace().getNumSamples()>0)
739 {
740
741 if (! ( ids = dataFile.add_var("id", ncInt, ncdims[rank+1])) )
742 throw DataException("Error - DataExpanded:: appending reference id to netCDF file failed.");
743 const int* ids_p=getFunctionSpace().borrowSampleReferenceIDs();
744 if (! (ids->put(ids_p,dims[rank+1])) )
745 throw DataException("Error - DataExpanded:: copy reference id to netCDF buffer failed.");
746 if (! ( var = dataFile.add_var("data", ncDouble, ndims, ncdims)) )
747 throw DataException("Error - DataExpanded:: appending variable to netCDF file failed.");
748 if (! (var->put(d_ptr,dims)) )
749 throw DataException("Error - DataExpanded:: copy data to netCDF buffer failed.");
750 }
751 #ifdef ESYS_MPI
752 if (mpi_iam<mpi_num-1) MPI_Send(&ndims, 0, MPI_INT, mpi_iam+1, 81801, MPI_COMM_WORLD);
753 #endif
754 #else
755 throw DataException("Error - DataExpanded:: dump is not configured with netCDF. Please contact your installation manager.");
756 #endif
757 }
758
759 void
760 DataExpanded::setTaggedValue(int tagKey,
761 const DataTypes::ShapeType& pointshape,
762 const DataTypes::ValueType& value,
763 int dataOffset)
764 {
765 CHECK_FOR_EX_WRITE
766 int numSamples = getNumSamples();
767 int numDataPointsPerSample = getNumDPPSample();
768 int sampleNo,dataPointNo, i;
769 DataTypes::ValueType::size_type n = getNoValues();
770 double* p;
771 const double* in=&value[0+dataOffset];
772
773 if (value.size() != n) {
774 throw DataException("Error - DataExpanded::setTaggedValue: number of input values does not match number of values per data points.");
775 }
776
777 #pragma omp parallel for private(sampleNo,dataPointNo,p,i) schedule(static)
778 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
779 if (getFunctionSpace().getTagFromSampleNo(sampleNo) == tagKey ) {
780 for (dataPointNo = 0; dataPointNo < numDataPointsPerSample; dataPointNo++) {
781 p=&(m_data[getPointOffset(sampleNo,dataPointNo)]);
782 for (i=0; i<n ;++i) p[i]=in[i];
783 }
784 }
785 }
786 }
787
788
789 void
790 DataExpanded::reorderByReferenceIDs(int *reference_ids)
791 {
792 CHECK_FOR_EX_WRITE
793 int numSamples = getNumSamples();
794 DataTypes::ValueType::size_type n = getNoValues() * getNumDPPSample();
795 int sampleNo, sampleNo2,i;
796 double* p,*p2;
797 register double rtmp;
798 FunctionSpace fs=getFunctionSpace();
799
800 for (sampleNo = 0; sampleNo < numSamples; sampleNo++) {
801 const int id_in=reference_ids[sampleNo];
802 const int id=fs.getReferenceIDOfSample(sampleNo);
803 if (id!=id_in) {
804 bool matched=false;
805 for (sampleNo2 = sampleNo+1; sampleNo2 < numSamples; sampleNo2++) {
806 if (id == reference_ids[sampleNo2]) {
807 p=&(m_data[getPointOffset(sampleNo,0)]);
808 p2=&(m_data[getPointOffset(sampleNo2,0)]);
809 for (i=0; i<n ;i++) {
810 rtmp=p[i];
811 p[i]=p2[i];
812 p2[i]=rtmp;
813 }
814 reference_ids[sampleNo]=id;
815 reference_ids[sampleNo2]=id_in;
816 matched=true;
817 break;
818 }
819 }
820 if (! matched) {
821 throw DataException("Error - DataExpanded::reorderByReferenceIDs: unable to reorder sample data by reference ids");
822 }
823 }
824 }
825 }
826
827 DataTypes::ValueType&
828 DataExpanded::getVectorRW()
829 {
830 CHECK_FOR_EX_WRITE
831 return m_data.getData();
832 }
833
834 const DataTypes::ValueType&
835 DataExpanded::getVectorRO() const
836 {
837 return m_data.getData();
838 }
839
840 void DataExpanded::randomFill(double seed)
841 {
842 CHECK_FOR_EX_WRITE
843 if (seed==0)
844 {
845 time_t s=time(0);
846 seed=s;
847 }
848 srand(seed);
849 DataVector& dv=getVectorRW();
850 for (long i=0;i<dv.size();++i)
851 {
852 dv[i]=(double)rand()/RAND_MAX;
853 }
854 }
855
856 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26