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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1628 - (show annotations)
Fri Jul 11 13:12:46 2008 UTC (11 years, 2 months ago) by phornby
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
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "DataArrayView.h"
17 #include "DataException.h"
18
19 #include <sstream>
20
21 #include <boost/python/extract.hpp>
22
23 using namespace std;
24 using namespace boost::python;
25
26 namespace escript {
27
28 DataArrayView::DataArrayView():
29 m_shape(),
30 m_offset(0),
31 m_noValues(0)
32 {
33 m_data=0;
34 }
35
36 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
63 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
71 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 bool
82 DataArrayView::isEmpty() const
83 {
84 return (m_data==NULL);
85 }
86
87 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
95 EsysAssert((!isEmpty()&&checkShape(tempShape)),
96 createShapeErrorMessage("Error - Couldn't copy due to shape mismatch.",tempShape));
97
98 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 }
129 }
130
131 void
132 DataArrayView::copy(const DataArrayView& other)
133 {
134 copy(m_offset,other);
135 }
136
137 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
152 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
168 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
182 int
183 DataArrayView::getRank() const
184 {
185 return m_shape.size();
186 }
187
188 const
189 DataArrayView::ShapeType&
190 DataArrayView::getShape() const
191 {
192 return m_shape;
193 }
194
195 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
208 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
221 int
222 DataArrayView::noValues() const
223 {
224 return m_noValues;
225 }
226
227 bool
228 DataArrayView::checkShape(const DataArrayView::ShapeType& other) const
229 {
230 return (m_shape==other);
231 }
232
233 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
244 DataArrayView::ValueType::size_type
245 DataArrayView::getOffset() const
246 {
247 return m_offset;
248 }
249
250 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
261 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
272 bool
273 DataArrayView::checkOffset() const
274 {
275 return checkOffset(m_offset);
276 }
277
278 bool
279 DataArrayView::checkOffset(ValueType::size_type offset) const
280 {
281 return (m_data->size() >= (offset+noValues()));
282 }
283
284 DataArrayView::ValueType&
285 DataArrayView::getData() const
286 {
287 EsysAssert(!isEmpty(),"Error - View is empty");
288 return *m_data;
289 }
290
291 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
300 DataArrayView::ShapeType
301 DataArrayView::getResultSliceShape(const RegionType& region)
302 {
303 int dimSize;
304 ShapeType result;
305 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
315 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 } else {
338 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 }
344 }
345 for (i=slice_rank;i<this_rank;i++) {
346 out[i]=std::pair<int,int>(0,getShape()[i]);
347 }
348 return out;
349 }
350
351 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 }
381 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
392 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 }
406 return region_loop_range;
407 }
408
409 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
418 void
419 DataArrayView::copySlice(ValueType::size_type thisOffset,
420 const DataArrayView& other,
421 ValueType::size_type otherOffset,
422 const RegionLoopRangeType& region)
423 {
424
425 //
426 // Make sure views are not empty
427
428 EsysAssert(!isEmpty(),
429 "Error - this view is empty.");
430 EsysAssert(!other.isEmpty(),
431 "Error - other view is empty.");
432
433 //
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
437 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
442 EsysAssert(other.getRank()==region.size(),
443 "Error - slice not same rank as view to be sliced from.");
444
445 EsysAssert(noValues()==noValues(getResultSliceShape(region)),
446 "Error - slice shape not compatible shape for this view.");
447
448 //
449 // copy the values in the specified region of the other view into this view
450
451 // the following loops cannot be parallelised due to the numCopy counter
452 int numCopy=0;
453
454 switch (region.size()) {
455 case 0:
456 /* this case should never be encountered,
457 as python will never pass us an empty region.
458 here for completeness only, allows slicing of a scalar */
459 (*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 }
503 }
504
505 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
514 void
515 DataArrayView::copySliceFrom(ValueType::size_type thisOffset,
516 const DataArrayView& other,
517 ValueType::size_type otherOffset,
518 const RegionLoopRangeType& region)
519 {
520
521 //
522 // Make sure views are not empty
523
524 EsysAssert(!isEmpty(),
525 "Error - this view is empty.");
526 EsysAssert(!other.isEmpty(),
527 "Error - other view is empty.");
528
529 //
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
533 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
538 EsysAssert(getRank()==region.size(),
539 "Error - slice not same rank as this view.");
540
541 EsysAssert(other.getRank()==0 || other.noValues()==noValues(getResultSliceShape(region)),
542 "Error - slice shape not compatible shape for other view.");
543
544 //
545 // copy the values in the other view into the specified region of this view
546
547 // allow for case where other view is a scalar
548 if (other.getRank()==0) {
549
550 // the following loops cannot be parallelised due to the numCopy counter
551 int numCopy=0;
552
553 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
603 } else {
604
605 // the following loops cannot be parallelised due to the numCopy counter
606 int numCopy=0;
607
608 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
658 }
659
660 }
661
662 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 }
671 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 }
681 }
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 }
724 return temp.str();
725 }
726
727 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 }
739 temp << ")";
740 return temp.str();
741 }
742
743 void
744 DataArrayView::matMult(const DataArrayView& left,
745 const DataArrayView& right,
746 DataArrayView& result)
747 {
748
749 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
755 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
764 int outputRank = left.getRank()+right.getRank()-2;
765
766 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
773 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
781 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 }
789
790 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 }
801
802 switch (left.getRank()) {
803
804 case 1:
805 switch (right.getRank()) {
806 case 1:
807 result()=0;
808 for (int i=0;i<left.getShape()[0];i++) {
809 result()+=left(i)*right(i);
810 }
811 break;
812 case 2:
813 for (int i=0;i<result.getShape()[0];i++) {
814 result(i)=0;
815 for (int j=0;j<right.getShape()[0];j++) {
816 result(i)+=left(j)*right(j,i);
817 }
818 }
819 break;
820 default:
821 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
822 throw DataException(temp.str());
823 break;
824 }
825 break;
826
827 case 2:
828 switch (right.getRank()) {
829 case 1:
830 result()=0;
831 for (int i=0;i<left.getShape()[0];i++) {
832 result(i)=0;
833 for (int j=0;j<left.getShape()[1];j++) {
834 result(i)+=left(i,j)*right(i);
835 }
836 }
837 break;
838 case 2:
839 for (int i=0;i<result.getShape()[0];i++) {
840 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 }
847 break;
848 default:
849 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
850 throw DataException(temp.str());
851 break;
852 }
853 break;
854
855 default:
856 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << left.getRank();
857 throw DataException(temp.str());
858 break;
859 }
860
861 }
862
863 DataArrayView::ShapeType
864 DataArrayView::determineResultShape(const DataArrayView& left,
865 const DataArrayView& right)
866 {
867 ShapeType result;
868 for (int i=0; i<(left.getRank()-1); i++) {
869 result.push_back(left.getShape()[i]);
870 }
871 for (int i=1; i<right.getRank(); i++) {
872 result.push_back(right.getShape()[i]);
873 }
874 return result;
875 }
876
877 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 }
888 return result;
889 }
890
891 bool operator!=(const DataArrayView& left, const DataArrayView& right)
892 {
893 return (!operator==(left,right));
894 }
895
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