/[escript]/branches/arrexp_2137_win_merge/escript/src/DataMaths.cpp
ViewVC logotype

Contents of /branches/arrexp_2137_win_merge/escript/src/DataMaths.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2213 - (show annotations)
Wed Jan 14 00:23:39 2009 UTC (10 years, 6 months ago) by jfenwick
File size: 5337 byte(s)
In preparation for merging to trunk

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 #include "DataTypes.h"
16 #include "DataMaths.h"
17 #include <sstream>
18
19 namespace escript
20 {
21 namespace DataMaths
22 {
23
24 void
25 matMult(const DataTypes::ValueType& left,
26 const DataTypes::ShapeType& leftShape,
27 DataTypes::ValueType::size_type leftOffset,
28 const DataTypes::ValueType& right,
29 const DataTypes::ShapeType& rightShape,
30 DataTypes::ValueType::size_type rightOffset,
31 DataTypes::ValueType& result,
32 const DataTypes::ShapeType& resultShape)
33 {
34 using namespace escript::DataTypes;
35 using namespace std;
36
37 int leftRank=getRank(leftShape);
38 int rightRank=getRank(rightShape);
39 int resultRank=getRank(resultShape);
40 if (leftRank==0 || rightRank==0) {
41 stringstream temp;
42 temp << "Error - (matMult) Invalid for rank 0 objects.";
43 throw DataException(temp.str());
44 }
45
46 if (leftShape[leftRank-1] != rightShape[0]) {
47 stringstream temp;
48 temp << "Error - (matMult) Dimension: " << leftRank
49 << ", size: " << leftShape[leftRank-1]
50 << " of LHS and dimension: 1, size: " << rightShape[0]
51 << " of RHS don't match.";
52 throw DataException(temp.str());
53 }
54
55 int outputRank = leftRank+rightRank-2;
56
57 if (outputRank < 0) {
58 stringstream temp;
59 temp << "Error - (matMult) LHS and RHS cannot be multiplied "
60 << "as they have incompatable rank.";
61 throw DataException(temp.str());
62 }
63
64 if (outputRank != resultRank) {
65 stringstream temp;
66 temp << "Error - (matMult) Rank of result array is: "
67 << resultRank
68 << " it must be: " << outputRank;
69 throw DataException(temp.str());
70 }
71
72 for (int i=0; i<(leftRank-1); i++) {
73 if (leftShape[i] != resultShape[i]) {
74 stringstream temp;
75 temp << "Error - (matMult) Dimension: " << i
76 << " of LHS and result array don't match.";
77 throw DataException(temp.str());
78 }
79 }
80
81 for (int i=1; i<rightRank; i++) {
82 if (rightShape[i] != resultShape[i+leftRank-2]) {
83 stringstream temp;
84 temp << "Error - (matMult) Dimension: " << i
85 << ", size: " << rightShape[i]
86 << " of RHS and dimension: " << i+leftRank-1
87 << ", size: " << resultShape[i+leftRank-1]
88 << " of result array don't match.";
89 throw DataException(temp.str());
90 }
91 }
92
93 switch (leftRank) {
94
95 case 1:
96 switch (rightRank) {
97 case 1:
98 result[0]=0;
99 for (int i=0;i<leftShape[0];i++) {
100 result[0]+=left[i+leftOffset]*right[i+rightOffset];
101 }
102 break;
103 case 2:
104 for (int i=0;i<resultShape[0];i++) {
105 result[i]=0;
106 for (int j=0;j<rightShape[0];j++) {
107 result[i]+=left[j+leftOffset]*right[getRelIndex(rightShape,j,i)+rightOffset];
108 }
109 }
110 break;
111 default:
112 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
113 throw DataException(temp.str());
114 break;
115 }
116 break;
117
118 case 2:
119 switch (rightRank) {
120 case 1:
121 result[0]=0;
122 for (int i=0;i<leftShape[0];i++) {
123 result[i]=0;
124 for (int j=0;j<leftShape[1];j++) {
125 result[i]+=left[leftOffset+getRelIndex(leftShape,i,j)]*right[i+rightOffset];
126 }
127 }
128 break;
129 case 2:
130 for (int i=0;i<resultShape[0];i++) {
131 for (int j=0;j<resultShape[1];j++) {
132 result[getRelIndex(resultShape,i,j)]=0;
133 for (int jR=0;jR<rightShape[0];jR++) {
134 result[getRelIndex(resultShape,i,j)]+=left[leftOffset+getRelIndex(leftShape,i,jR)]*right[rightOffset+getRelIndex(rightShape,jR,j)];
135 }
136 }
137 }
138 break;
139 default:
140 stringstream temp; temp << "Error - (matMult) Invalid rank. Programming error.";
141 throw DataException(temp.str());
142 break;
143 }
144 break;
145
146 default:
147 stringstream temp; temp << "Error - (matMult) Not supported for rank: " << leftRank;
148 throw DataException(temp.str());
149 break;
150 }
151
152 }
153
154
155 DataTypes::ShapeType
156 determineResultShape(const DataTypes::ShapeType& left,
157 const DataTypes::ShapeType& right)
158 {
159 DataTypes::ShapeType result;
160 for (int i=0; i<(DataTypes::getRank(left)-1); i++) {
161 result.push_back(left[i]);
162 }
163 for (int i=1; i<DataTypes::getRank(right); i++) {
164 result.push_back(right[i]);
165 }
166 return result;
167 }
168
169
170 } // end namespace
171 } // end namespace
172

  ViewVC Help
Powered by ViewVC 1.1.26