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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26