/[escript]/branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp
ViewVC logotype

Annotation of /branches/arrayview_from_1695_trunk/escript/src/DataArrayView.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (hide annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 6 months ago) by phornby
Original Path: trunk/escript/src/DataArrayView.cpp
File size: 29584 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 jgs 82
4 ksteube 1312 /*******************************************************
5 phornby 1455 *
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 ksteube 1312
16 jgs 478 #include "DataArrayView.h"
17     #include "DataException.h"
18    
19     #include <sstream>
20    
21     #include <boost/python/extract.hpp>
22    
23 jgs 82 using namespace std;
24     using namespace boost::python;
25    
26     namespace escript {
27    
28 phornby 1455 DataArrayView::DataArrayView():
29     m_shape(),
30     m_offset(0),
31     m_noValues(0)
32     {
33     m_data=0;
34     }
35 jgs 82
36 phornby 1455 DataArrayView::DataArrayView(ValueType& data,
37     const ShapeType& viewShape,
38     int offset):
39     m_data(&data),
40     m_shape(viewShape),
41     m_offset(offset),
42     m_noValues(noValues(viewShape))
43     {
44     //
45     // check the shape rank and size and throw an exception if a
46     // shape incompatible with the given data array is specified
47     if (viewShape.size() > m_maxRank) {
48     stringstream temp;
49     temp << "Error- Couldn't construct DataArrayView, invalid rank."
50     << "Input data rank: " << viewShape.size() << " "
51     << "Maximum rank allowed for DataArrayView: " << m_maxRank;
52     throw DataException(temp.str());
53     }
54     if (m_data->size() < (m_offset+noValues())) {
55     stringstream temp;
56     temp << "Error- Couldn't construct DataArrayView, insufficient data values. "
57     << "Shape requires: " << noValues()+m_offset << " values."
58     << "Data values available: " << m_data->size();
59     throw DataException(temp.str());
60     }
61     }
62 jgs 82
63 phornby 1455 DataArrayView::DataArrayView(const DataArrayView& other):
64     m_data(other.m_data),
65     m_offset(other.m_offset),
66     m_shape(other.m_shape),
67     m_noValues(other.m_noValues)
68     {
69     }
70 jgs 82
71 phornby 1628 DataArrayView &
72     DataArrayView::operator=(const DataArrayView& other)
73     {
74     m_data = other.m_data;
75     m_offset = other.m_offset;
76     m_shape = other.m_shape;
77     m_noValues = other.m_noValues;
78     return *this;
79     }
80    
81 phornby 1455 bool
82     DataArrayView::isEmpty() const
83     {
84 phornby 1628 return (m_data==NULL);
85 phornby 1455 }
86 jgs 82
87 phornby 1455 void
88     DataArrayView::copy(const boost::python::numeric::array& value)
89     {
90     ShapeType tempShape;
91     for (int i=0; i<value.getrank(); i++) {
92     tempShape.push_back(extract<int>(value.getshape()[i]));
93     }
94 jgs 108
95 phornby 1455 EsysAssert((!isEmpty()&&checkShape(tempShape)),
96     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
97 jgs 82
98 phornby 1455 if (value.getrank()==0) {
99     (*this)()=extract<double>(value[value.getshape()]);
100     } else if (value.getrank()==1) {
101     for (ValueType::size_type i=0;i<tempShape[0];i++) {
102     (*this)(i)=extract<double>(value[i]);
103     }
104     } else if (value.getrank()==2) {
105     for (ValueType::size_type i=0;i<tempShape[0];i++) {
106     for (ValueType::size_type j=0;j<tempShape[1];j++) {
107     (*this)(i,j)=extract<double>(value[i][j]);
108     }
109     }
110     } else if (value.getrank()==3) {
111     for (ValueType::size_type i=0;i<tempShape[0];i++) {
112     for (ValueType::size_type j=0;j<tempShape[1];j++) {
113     for (ValueType::size_type k=0;k<tempShape[2];k++) {
114     (*this)(i,j,k)=extract<double>(value[i][j][k]);
115     }
116     }
117     }
118     } else if (value.getrank()==4) {
119     for (ValueType::size_type i=0;i<tempShape[0];i++) {
120     for (ValueType::size_type j=0;j<tempShape[1];j++) {
121     for (ValueType::size_type k=0;k<tempShape[2];k++) {
122     for (ValueType::size_type l=0;l<tempShape[3];l++) {
123     (*this)(i,j,k,l)=extract<double>(value[i][j][k][l]);
124     }
125     }
126     }
127     }
128 jgs 82 }
129 phornby 1455 }
130 jgs 102
131 phornby 1455 void
132     DataArrayView::copy(const DataArrayView& other)
133     {
134     copy(m_offset,other);
135     }
136 jgs 82
137 phornby 1455 void
138     DataArrayView::copy(ValueType::size_type offset,
139     const DataArrayView& other)
140     {
141     EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(offset)),
142     "Error - Couldn't copy due to insufficient storage.");
143     EsysAssert((checkShape(other.getShape())),
144     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
145     if (checkOffset(offset)) {
146     memcpy(&(*m_data)[offset],&(*other.m_data)[other.m_offset],sizeof(double)*noValues());
147     } else {
148     throw DataException("Error - invalid offset specified.");
149     }
150     }
151 jgs 82
152 phornby 1455 void
153     DataArrayView::copy(ValueType::size_type thisOffset,
154     const DataArrayView& other,
155     ValueType::size_type otherOffset)
156     {
157     EsysAssert((!isEmpty()&&!other.isEmpty()&&checkOffset(thisOffset)&&other.checkOffset(otherOffset)),
158     "Error - Couldn't copy due to insufficient storage.");
159     EsysAssert((checkShape(other.getShape())),
160     createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",other.getShape()));
161     if (checkOffset(thisOffset)&&other.checkOffset(otherOffset)) {
162     memcpy(&(*m_data)[thisOffset],&(*other.m_data)[otherOffset],sizeof(double)*noValues());
163     } else {
164     throw DataException("Error - invalid offset specified.");
165     }
166     }
167 jgs 82
168 phornby 1455 void
169     DataArrayView::copy(ValueType::size_type offset,
170     ValueType::value_type value)
171     {
172     EsysAssert((!isEmpty()&&checkOffset(offset)),
173     "Error - Couldn't copy due to insufficient storage.");
174     if (checkOffset(offset)) {
175     vector<double> temp(noValues(),value);
176     memcpy(&(*m_data)[offset],&temp[0],sizeof(double)*noValues());
177     } else {
178     throw DataException("Error - invalid offset specified.");
179     }
180     }
181 jgs 82
182 phornby 1455 int
183     DataArrayView::getRank() const
184     {
185     return m_shape.size();
186     }
187 jgs 82
188 phornby 1455 const
189     DataArrayView::ShapeType&
190     DataArrayView::getShape() const
191     {
192     return m_shape;
193     }
194 jgs 82
195 phornby 1455 int
196     DataArrayView::noValues(const ShapeType& shape)
197     {
198     ShapeType::const_iterator i;
199     //
200     // An empty shape vector means rank 0 which contains 1 value
201     int noValues=1;
202     for (i=shape.begin();i!=shape.end();i++) {
203     noValues*=(*i);
204     }
205     return noValues;
206     }
207 jgs 108
208 phornby 1455 int
209     DataArrayView::noValues(const RegionLoopRangeType& region)
210     {
211     //
212     // An empty region vector means rank 0 which contains 1 value
213     int noValues=1;
214     unsigned int i;
215     for (i=0;i<region.size();i++) {
216     noValues*=region[i].second-region[i].first;
217     }
218     return noValues;
219     }
220 jgs 102
221 phornby 1455 int
222     DataArrayView::noValues() const
223     {
224     return m_noValues;
225     }
226 jgs 82
227 phornby 1455 bool
228     DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
229     {
230     return (m_shape==other);
231     }
232 jgs 82
233 phornby 1455 string
234     DataArrayView::createShapeErrorMessage(const string& messagePrefix,
235     const DataArrayView::ShapeType& other) const
236     {
237     stringstream temp;
238     temp << messagePrefix
239     << " This shape: " << shapeToString(m_shape)
240     << " Other shape: " << shapeToString(other);
241     return temp.str();
242     }
243 jgs 108
244 phornby 1455 DataArrayView::ValueType::size_type
245     DataArrayView::getOffset() const
246     {
247     return m_offset;
248     }
249 jgs 108
250 phornby 1455 void
251     DataArrayView::setOffset(ValueType::size_type offset)
252     {
253     EsysAssert((checkOffset(offset)), "Error - Invalid offset.");
254     if (checkOffset(offset)) {
255     m_offset=offset;
256     } else {
257     throw DataException("Error - invalid offset specified.");
258     }
259     }
260 jgs 82
261 phornby 1455 void
262     DataArrayView::incrOffset()
263     {
264     EsysAssert((checkOffset(m_offset+noValues())), "Error - Cannot increment offset.");
265     if (checkOffset(m_offset+noValues())) {
266     m_offset=m_offset+noValues();
267     } else {
268     throw DataException("Error - Cannot increment offset.");
269     }
270     }
271 jgs 108
272 phornby 1455 bool
273     DataArrayView::checkOffset() const
274     {
275     return checkOffset(m_offset);
276     }
277 jgs 108
278 phornby 1455 bool
279     DataArrayView::checkOffset(ValueType::size_type offset) const
280     {
281     return (m_data->size() >= (offset+noValues()));
282     }
283 jgs 108
284 phornby 1455 DataArrayView::ValueType&
285     DataArrayView::getData() const
286     {
287     EsysAssert(!isEmpty(),"Error - View is empty");
288     return *m_data;
289     }
290 jgs 108
291 phornby 1455 DataArrayView::ValueType::reference
292     DataArrayView::getData(ValueType::size_type i) const
293     {
294     //
295     // don't do any index checking so allow one past the end of the
296     // vector providing the equivalent of end()
297     return (*m_data)[m_offset+i];
298     }
299 jgs 108
300 phornby 1455 DataArrayView::ShapeType
301     DataArrayView::getResultSliceShape(const RegionType& region)
302     {
303     int dimSize;
304 phornby 1628 ShapeType result;
305 phornby 1455 RegionType::const_iterator i;
306     for (i=region.begin();i!=region.end();i++) {
307     dimSize=((i->second)-(i->first));
308     if (dimSize!=0) {
309     result.push_back(dimSize);
310     }
311     }
312     return result;
313     }
314 jgs 82
315 phornby 1455 DataArrayView::RegionType
316     DataArrayView::getSliceRegion(const boost::python::object& key) const
317     {
318     int slice_rank, i;
319     int this_rank=getRank();
320     RegionType out(this_rank);
321     /* allow for case where key is singular eg: [1], this implies we
322     want to generate a rank-1 dimension object, as opposed to eg: [1,2]
323     which implies we want to take a rank dimensional object with one
324     dimension of size 1 */
325     extract<tuple> key_tuple(key);
326     if (key_tuple.check()) {
327     slice_rank=extract<int> (key.attr("__len__")());
328     /* ensure slice is correctly dimensioned */
329     if (slice_rank>this_rank) {
330     throw DataException("Error - rank of slices does not match rank of slicee");
331     } else {
332     /* calculate values for slice range */
333     for (i=0;i<slice_rank;i++) {
334     out[i]=getSliceRange(key[i],getShape()[i]);
335     }
336     }
337 jgs 108 } else {
338 phornby 1455 slice_rank=1;
339     if (slice_rank>this_rank) {
340     throw DataException("Error - rank of slices does not match rank of slicee");
341     } else {
342     out[0]=getSliceRange(key,getShape()[0]);
343 jgs 108 }
344     }
345 phornby 1455 for (i=slice_rank;i<this_rank;i++) {
346     out[i]=std::pair<int,int>(0,getShape()[i]);
347 jgs 108 }
348 phornby 1455 return out;
349     }
350 jgs 102
351 phornby 1455 std::pair<int,int>
352     getSliceRange(const boost::python::object& key,
353     const int shape)
354     {
355     /* default slice range is range of entire shape dimension */
356     int s0=0, s1=shape;;
357     extract<int> slice_int(key);
358     if (slice_int.check()) {
359     /* if the key is a single int set start=key and end=key */
360     /* in this case, we want to return a rank-1 dimension object from
361     this object, taken from a single index value for one of this
362     object's dimensions */
363     s0=slice_int();
364     s1=s0;
365     } else {
366     /* if key is a pair extract begin and end values */
367     extract<int> step(key.attr("step"));
368     if (step.check() && step()!=1) {
369     throw DataException("Error - Data does not support increments in slicing ");
370     } else {
371     extract<int> start(key.attr("start"));
372     if (start.check()) {
373     s0=start();
374     }
375     extract<int> stop(key.attr("stop"));
376     if (stop.check()) {
377     s1=stop();
378     }
379     }
380 jgs 102 }
381 phornby 1455 if (s0 < 0)
382     throw DataException("Error - slice index out of range.");
383     if (s0 == s1 && s1 >= shape)
384     throw DataException("Error - slice index out of range.");
385     if (s0 != s1 && s1>shape)
386     throw DataException("Error - slice index out of range.");
387     if (s0 > s1)
388     throw DataException("Error - lower index must less or equal upper index.");
389     return std::pair<int,int>(s0,s1);
390     }
391 jgs 102
392 phornby 1455 DataArrayView::RegionLoopRangeType
393     getSliceRegionLoopRange(const DataArrayView::RegionType& region)
394     {
395     DataArrayView::RegionLoopRangeType region_loop_range(region.size());
396     unsigned int i;
397     for (i=0;i<region.size();i++) {
398     if (region[i].first==region[i].second) {
399     region_loop_range[i].first=region[i].first;
400     region_loop_range[i].second=region[i].second+1;
401     } else {
402     region_loop_range[i].first=region[i].first;
403     region_loop_range[i].second=region[i].second;
404     }
405 jgs 108 }
406 phornby 1455 return region_loop_range;
407     }
408 jgs 108
409 phornby 1455 void
410     DataArrayView::copySlice(const DataArrayView& other,
411     const RegionLoopRangeType& region)
412     {
413     //
414     // version of copySlice that uses the default offsets
415     copySlice(m_offset,other,other.m_offset,region);
416     }
417 jgs 82
418 phornby 1455 void
419     DataArrayView::copySlice(ValueType::size_type thisOffset,
420     const DataArrayView& other,
421     ValueType::size_type otherOffset,
422     const RegionLoopRangeType& region)
423     {
424 jgs 102
425 phornby 1455 //
426     // Make sure views are not empty
427 jgs 102
428 phornby 1455 EsysAssert(!isEmpty(),
429     "Error - this view is empty.");
430     EsysAssert(!other.isEmpty(),
431     "Error - other view is empty.");
432 jgs 108
433 phornby 1455 //
434     // Check the view to be sliced from is compatible with the region to be sliced,
435     // and that the region to be sliced is compatible with this view:
436 jgs 102
437 phornby 1455 EsysAssert(checkOffset(thisOffset),
438     "Error - offset incompatible with this view.");
439     EsysAssert(otherOffset+noValues()<=other.getData().size(),
440     "Error - offset incompatible with other view.");
441 jgs 102
442 phornby 1455 EsysAssert(other.getRank()==region.size(),
443     "Error - slice not same rank as view to be sliced from.");
444 jgs 102
445 phornby 1455 EsysAssert(noValues()==noValues(getResultSliceShape(region)),
446     "Error - slice shape not compatible shape for this view.");
447 jgs 102
448 phornby 1455 //
449     // copy the values in the specified region of the other view into this view
450 jgs 102
451 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
452     int numCopy=0;
453 jgs 108
454 phornby 1455 switch (region.size()) {
455     case 0:
456     /* this case should never be encountered,
457 jgs 108 as python will never pass us an empty region.
458     here for completeness only, allows slicing of a scalar */
459 phornby 1455 (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex()];
460     numCopy++;
461     break;
462     case 1:
463     for (int i=region[0].first;i<region[0].second;i++) {
464     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i)];
465     numCopy++;
466     }
467     break;
468     case 2:
469     for (int j=region[1].first;j<region[1].second;j++) {
470     for (int i=region[0].first;i<region[0].second;i++) {
471     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j)];
472     numCopy++;
473     }
474     }
475     break;
476     case 3:
477     for (int k=region[2].first;k<region[2].second;k++) {
478     for (int j=region[1].first;j<region[1].second;j++) {
479     for (int i=region[0].first;i<region[0].second;i++) {
480     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k)];
481     numCopy++;
482     }
483     }
484     }
485     break;
486     case 4:
487     for (int l=region[3].first;l<region[3].second;l++) {
488     for (int k=region[2].first;k<region[2].second;k++) {
489     for (int j=region[1].first;j<region[1].second;j++) {
490     for (int i=region[0].first;i<region[0].second;i++) {
491     (*m_data)[thisOffset+numCopy]=(*other.m_data)[otherOffset+other.relIndex(i,j,k,l)];
492     numCopy++;
493     }
494     }
495     }
496     }
497     break;
498     default:
499     stringstream mess;
500     mess << "Error - (copySlice) Invalid slice region rank: " << region.size();
501     throw DataException(mess.str());
502 jgs 82 }
503 phornby 1455 }
504 jgs 102
505 phornby 1455 void
506     DataArrayView::copySliceFrom(const DataArrayView& other,
507     const RegionLoopRangeType& region_loop_range)
508     {
509     //
510     // version of copySliceFrom that uses the default offsets
511     copySliceFrom(m_offset,other,other.m_offset,region_loop_range);
512     }
513 jgs 82
514 phornby 1455 void
515     DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
516     const DataArrayView& other,
517     ValueType::size_type otherOffset,
518     const RegionLoopRangeType& region)
519     {
520 jgs 102
521 phornby 1455 //
522     // Make sure views are not empty
523 jgs 102
524 phornby 1455 EsysAssert(!isEmpty(),
525     "Error - this view is empty.");
526     EsysAssert(!other.isEmpty(),
527     "Error - other view is empty.");
528 jgs 108
529 phornby 1455 //
530     // Check this view is compatible with the region to be sliced,
531     // and that the region to be sliced is compatible with the other view:
532 jgs 102
533 phornby 1455 EsysAssert(other.checkOffset(otherOffset),
534     "Error - offset incompatible with other view.");
535     EsysAssert(thisOffset+other.noValues()<=getData().size(),
536     "Error - offset incompatible with this view.");
537 jgs 102
538 phornby 1455 EsysAssert(getRank()==region.size(),
539     "Error - slice not same rank as this view.");
540 jgs 102
541 phornby 1455 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
542     "Error - slice shape not compatible shape for other view.");
543 jgs 102
544 phornby 1455 //
545     // copy the values in the other view into the specified region of this view
546 jgs 102
547 phornby 1455 // allow for case where other view is a scalar
548     if (other.getRank()==0) {
549 jgs 108
550 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
551     int numCopy=0;
552 jgs 108
553 phornby 1455 switch (region.size()) {
554     case 0:
555     /* this case should never be encountered,
556     as python will never pass us an empty region.
557     here for completeness only, allows slicing of a scalar */
558     (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset];
559     numCopy++;
560     break;
561     case 1:
562     for (int i=region[0].first;i<region[0].second;i++) {
563     (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset];
564     numCopy++;
565     }
566     break;
567     case 2:
568     for (int j=region[1].first;j<region[1].second;j++) {
569     for (int i=region[0].first;i<region[0].second;i++) {
570     (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset];
571     numCopy++;
572     }
573     }
574     break;
575     case 3:
576     for (int k=region[2].first;k<region[2].second;k++) {
577     for (int j=region[1].first;j<region[1].second;j++) {
578     for (int i=region[0].first;i<region[0].second;i++) {
579     (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset];
580     numCopy++;
581     }
582     }
583     }
584     break;
585     case 4:
586     for (int l=region[3].first;l<region[3].second;l++) {
587     for (int k=region[2].first;k<region[2].second;k++) {
588     for (int j=region[1].first;j<region[1].second;j++) {
589     for (int i=region[0].first;i<region[0].second;i++) {
590     (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset];
591     numCopy++;
592     }
593     }
594     }
595     }
596     break;
597     default:
598     stringstream mess;
599     mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
600     throw DataException(mess.str());
601     }
602 jgs 108
603 phornby 1455 } else {
604 jgs 108
605 phornby 1455 // the following loops cannot be parallelised due to the numCopy counter
606     int numCopy=0;
607 jgs 108
608 phornby 1455 switch (region.size()) {
609     case 0:
610     /* this case should never be encountered,
611     as python will never pass us an empty region.
612     here for completeness only, allows slicing of a scalar */
613     (*m_data)[thisOffset+relIndex()]=(*other.m_data)[otherOffset+numCopy];
614     numCopy++;
615     break;
616     case 1:
617     for (int i=region[0].first;i<region[0].second;i++) {
618     (*m_data)[thisOffset+relIndex(i)]=(*other.m_data)[otherOffset+numCopy];
619     numCopy++;
620     }
621     break;
622     case 2:
623     for (int j=region[1].first;j<region[1].second;j++) {
624     for (int i=region[0].first;i<region[0].second;i++) {
625     (*m_data)[thisOffset+relIndex(i,j)]=(*other.m_data)[otherOffset+numCopy];
626     numCopy++;
627     }
628     }
629     break;
630     case 3:
631     for (int k=region[2].first;k<region[2].second;k++) {
632     for (int j=region[1].first;j<region[1].second;j++) {
633     for (int i=region[0].first;i<region[0].second;i++) {
634     (*m_data)[thisOffset+relIndex(i,j,k)]=(*other.m_data)[otherOffset+numCopy];
635     numCopy++;
636     }
637     }
638     }
639     break;
640     case 4:
641     for (int l=region[3].first;l<region[3].second;l++) {
642     for (int k=region[2].first;k<region[2].second;k++) {
643     for (int j=region[1].first;j<region[1].second;j++) {
644     for (int i=region[0].first;i<region[0].second;i++) {
645     (*m_data)[thisOffset+relIndex(i,j,k,l)]=(*other.m_data)[otherOffset+numCopy];
646     numCopy++;
647     }
648     }
649     }
650     }
651     break;
652     default:
653     stringstream mess;
654     mess << "Error - (copySliceFrom) Invalid slice region rank: " << region.size();
655     throw DataException(mess.str());
656     }
657 jgs 108
658 phornby 1455 }
659 jgs 108
660 phornby 1455 }
661 jgs 102
662 phornby 1455 string
663     DataArrayView::toString(const string& suffix) const
664     {
665     EsysAssert(!isEmpty(),"Error - View is empty.");
666     stringstream temp;
667     string finalSuffix=suffix;
668     if (suffix.length() > 0) {
669     finalSuffix+=" ";
670 jgs 82 }
671 phornby 1455 switch (getRank()) {
672     case 0:
673     temp << finalSuffix << (*this)();
674     break;
675     case 1:
676     for (int i=0;i<getShape()[0];i++) {
677     temp << finalSuffix << "(" << i << ") " << (*this)(i);
678     if (i!=(getShape()[0]-1)) {
679     temp << endl;
680 jgs 102 }
681 phornby 1455 }
682     break;
683     case 2:
684     for (int i=0;i<getShape()[0];i++) {
685     for (int j=0;j<getShape()[1];j++) {
686     temp << finalSuffix << "(" << i << "," << j << ") " << (*this)(i,j);
687     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1))) {
688     temp << endl;
689     }
690     }
691     }
692     break;
693     case 3:
694     for (int i=0;i<getShape()[0];i++) {
695     for (int j=0;j<getShape()[1];j++) {
696     for (int k=0;k<getShape()[2];k++) {
697     temp << finalSuffix << "(" << i << "," << j << "," << k << ") " << (*this)(i,j,k);
698     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1))) {
699     temp << endl;
700     }
701     }
702     }
703     }
704     break;
705     case 4:
706     for (int i=0;i<getShape()[0];i++) {
707     for (int j=0;j<getShape()[1];j++) {
708     for (int k=0;k<getShape()[2];k++) {
709     for (int l=0;l<getShape()[3];l++) {
710     temp << finalSuffix << "(" << i << "," << j << "," << k << "," << l << ") " << (*this)(i,j,k,l);
711     if (!(i==(getShape()[0]-1) && j==(getShape()[1]-1) && k==(getShape()[2]-1) && l==(getShape()[3]-1))) {
712     temp << endl;
713     }
714     }
715     }
716     }
717     }
718     break;
719     default:
720     stringstream mess;
721     mess << "Error - (toString) Invalid rank: " << getRank();
722     throw DataException(mess.str());
723 jgs 102 }
724 phornby 1455 return temp.str();
725     }
726 jgs 82
727 phornby 1455 string
728     DataArrayView::shapeToString(const DataArrayView::ShapeType& shape)
729     {
730     stringstream temp;
731     temp << "(";
732     unsigned int i;
733     for (i=0;i<shape.size();i++) {
734     temp << shape[i];
735     if (i < shape.size()-1) {
736     temp << ",";
737     }
738 jgs 82 }
739 phornby 1455 temp << ")";
740     return temp.str();
741     }
742 jgs 82
743 phornby 1455 void
744     DataArrayView::matMult(const DataArrayView& left,
745     const DataArrayView& right,
746     DataArrayView& result)
747     {
748 jgs 113
749 phornby 1455 if (left.getRank()==0 || right.getRank()==0) {
750     stringstream temp;
751     temp << "Error - (matMult) Invalid for rank 0 objects.";
752     throw DataException(temp.str());
753     }
754 jgs 82
755 phornby 1455 if (left.getShape()[left.getRank()-1] != right.getShape()[0]) {
756     stringstream temp;
757     temp << "Error - (matMult) Dimension: " << left.getRank()
758     << ", size: " << left.getShape()[left.getRank()-1]
759     << " of LHS and dimension: 1, size: " << right.getShape()[0]
760     << " of RHS don't match.";
761     throw DataException(temp.str());
762     }
763 jgs 113
764 phornby 1455 int outputRank = left.getRank()+right.getRank()-2;
765 jgs 126
766 phornby 1455 if (outputRank < 0) {
767     stringstream temp;
768     temp << "Error - (matMult) LHS and RHS cannot be multiplied "
769     << "as they have incompatable rank.";
770     throw DataException(temp.str());
771     }
772 jgs 82
773 phornby 1455 if (outputRank != result.getRank()) {
774     stringstream temp;
775     temp << "Error - (matMult) Rank of result array is: "
776     << result.getRank()
777     << " it must be: " << outputRank;
778     throw DataException(temp.str());
779     }
780 jgs 82
781 phornby 1455 for (int i=0; i<(left.getRank()-1); i++) {
782     if (left.getShape()[i] != result.getShape()[i]) {
783     stringstream temp;
784     temp << "Error - (matMult) Dimension: " << i
785     << " of LHS and result array don't match.";
786     throw DataException(temp.str());
787     }
788 jgs 82 }
789    
790 phornby 1455 for (int i=1; i<right.getRank(); i++) {
791     if (right.getShape()[i] != result.getShape()[i+left.getRank()-2]) {
792     stringstream temp;
793     temp << "Error - (matMult) Dimension: " << i
794     << ", size: " << right.getShape()[i]
795     << " of RHS and dimension: " << i+left.getRank()-1
796     << ", size: " << result.getShape()[i+left.getRank()-1]
797     << " of result array don't match.";
798     throw DataException(temp.str());
799     }
800 jgs 82 }
801    
802 phornby 1455 switch (left.getRank()) {
803 jgs 113
804 jgs 82 case 1:
805 phornby 1455 switch (right.getRank()) {
806     case 1:
807 jgs 113 result()=0;
808     for (int i=0;i<left.getShape()[0];i++) {
809 phornby 1455 result()+=left(i)*right(i);
810 jgs 113 }
811     break;
812 phornby 1455 case 2:
813 jgs 113 for (int i=0;i<result.getShape()[0];i++) {
814 phornby 1455 result(i)=0;
815     for (int j=0;j<right.getShape()[0];j++) {
816     result(i)+=left(j)*right(j,i);
817     }
818 jgs 113 }
819     break;
820 phornby 1455 default:
821 jgs 113 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
822     throw DataException(temp.str());
823     break;
824 phornby 1455 }
825     break;
826 jgs 113
827 jgs 82 case 2:
828 phornby 1455 switch (right.getRank()) {
829     case 1:
830 jgs 113 result()=0;
831 phornby 1455 for (int i=0;i<left.getShape()[0];i++) {
832 jgs 113 result(i)=0;
833     for (int j=0;j<left.getShape()[1];j++) {
834 phornby 1455 result(i)+=left(i,j)*right(i);
835 jgs 113 }
836     }
837     break;
838 phornby 1455 case 2:
839 jgs 113 for (int i=0;i<result.getShape()[0];i++) {
840 phornby 1455 for (int j=0;j<result.getShape()[1];j++) {
841     result(i,j)=0;
842     for (int jR=0;jR<right.getShape()[0];jR++) {
843     result(i,j)+=left(i,jR)*right(jR,j);
844     }
845     }
846 jgs 113 }
847     break;
848 phornby 1455 default:
849 jgs 113 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
850     throw DataException(temp.str());
851     break;
852 phornby 1455 }
853     break;
854 jgs 113
855 jgs 82 default:
856 phornby 1455 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
857     throw DataException(temp.str());
858     break;
859     }
860 jgs 113
861 phornby 1455 }
862 jgs 82
863 phornby 1455 DataArrayView::ShapeType
864     DataArrayView::determineResultShape(const DataArrayView& left,
865     const DataArrayView& right)
866     {
867 phornby 1628 ShapeType result;
868 phornby 1455 for (int i=0; i<(left.getRank()-1); i++) {
869 phornby 1628 result.push_back(left.getShape()[i]);
870 phornby 1455 }
871     for (int i=1; i<right.getRank(); i++) {
872 phornby 1628 result.push_back(right.getShape()[i]);
873 phornby 1455 }
874 phornby 1628 return result;
875 phornby 1455 }
876 jgs 113
877 phornby 1455 bool operator==(const DataArrayView& left, const DataArrayView& right)
878     {
879     //
880     // The views are considered equal if they have the same shape and each value
881     // of the view is the same.
882     bool result=(left.m_shape==right.m_shape);
883     if (result) {
884     for (int i=0;i<left.noValues();i++) {
885     result=result && (left.getData(i)==right.getData(i));
886     }
887 jgs 82 }
888 phornby 1455 return result;
889     }
890 jgs 82
891 phornby 1455 bool operator!=(const DataArrayView& left, const DataArrayView& right)
892     {
893     return (!operator==(left,right));
894     }
895 jgs 82
896     } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26