|
@@ -0,0 +1,406 @@
|
|
1
|
+/*!
|
|
2
|
+ * \file src/utils/Matrix.cpp
|
|
3
|
+ * \brief 3D Matrix
|
|
4
|
+ *
|
|
5
|
+ * \author Mongoose
|
|
6
|
+ */
|
|
7
|
+
|
|
8
|
+#include <stdio.h>
|
|
9
|
+#include <math.h>
|
|
10
|
+
|
|
11
|
+#include "utils/Matrix.h"
|
|
12
|
+
|
|
13
|
+Matrix::Matrix() {
|
|
14
|
+ setIdentity();
|
|
15
|
+}
|
|
16
|
+
|
|
17
|
+Matrix::Matrix(matrix_t m) {
|
|
18
|
+ setMatrix(m);
|
|
19
|
+}
|
|
20
|
+
|
|
21
|
+Matrix::Matrix(Quaternion &q) {
|
|
22
|
+ matrix_t m;
|
|
23
|
+ q.getMatrix(m);
|
|
24
|
+ setMatrix(m);
|
|
25
|
+}
|
|
26
|
+
|
|
27
|
+bool Matrix::getInvert(matrix_t out) {
|
|
28
|
+ matrix_t m;
|
|
29
|
+
|
|
30
|
+#ifdef COLUMN_ORDER
|
|
31
|
+ getMatrix(m);
|
|
32
|
+#else
|
|
33
|
+ getTransposeMatrix(m);
|
|
34
|
+#endif
|
|
35
|
+
|
|
36
|
+ /* Mongoose: This code was from a Jeff Lander tutorial which was based
|
|
37
|
+ on MESA GL's InvertMatrix */
|
|
38
|
+
|
|
39
|
+ /* NB. OpenGL Matrices are COLUMN major. */
|
|
40
|
+#define SWAP_ROWS(a, b) { float *_tmp = a; (a)=(b); (b)=_tmp; }
|
|
41
|
+#define MAT(m,r,c) (m)[(c)*4+(r)]
|
|
42
|
+
|
|
43
|
+ float wtmp[4][8];
|
|
44
|
+ float m0, m1, m2, m3, s;
|
|
45
|
+ float *r0, *r1, *r2, *r3;
|
|
46
|
+
|
|
47
|
+ r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
|
|
48
|
+
|
|
49
|
+ r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
|
|
50
|
+ r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
|
|
51
|
+ r0[4] = 1.0f, r0[5] = r0[6] = r0[7] = 0.0f,
|
|
52
|
+
|
|
53
|
+ r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
|
|
54
|
+ r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
|
|
55
|
+ r1[5] = 1.0f, r1[4] = r1[6] = r1[7] = 0.0f,
|
|
56
|
+
|
|
57
|
+ r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
|
|
58
|
+ r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
|
|
59
|
+ r2[6] = 1.0f, r2[4] = r2[5] = r2[7] = 0.0f,
|
|
60
|
+
|
|
61
|
+ r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
|
|
62
|
+ r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
|
|
63
|
+ r3[7] = 1.0f, r3[4] = r3[5] = r3[6] = 0.0f;
|
|
64
|
+
|
|
65
|
+ /* choose pivot - or die */
|
|
66
|
+ if (fabs(r3[0])>fabs(r2[0])) SWAP_ROWS(r3, r2);
|
|
67
|
+ if (fabs(r2[0])>fabs(r1[0])) SWAP_ROWS(r2, r1);
|
|
68
|
+ if (fabs(r1[0])>fabs(r0[0])) SWAP_ROWS(r1, r0);
|
|
69
|
+ if (0.0f == r0[0]) return false;
|
|
70
|
+
|
|
71
|
+ /* eliminate first variable */
|
|
72
|
+ m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
|
|
73
|
+ s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
|
|
74
|
+ s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
|
|
75
|
+ s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
|
|
76
|
+ s = r0[4];
|
|
77
|
+ if (s != 0.0f) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
|
|
78
|
+ s = r0[5];
|
|
79
|
+ if (s != 0.0f) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
|
|
80
|
+ s = r0[6];
|
|
81
|
+ if (s != 0.0f) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
|
|
82
|
+ s = r0[7];
|
|
83
|
+ if (s != 0.0f) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }
|
|
84
|
+
|
|
85
|
+ /* choose pivot - or die */
|
|
86
|
+ if (fabs(r3[1])>fabs(r2[1])) SWAP_ROWS(r3, r2);
|
|
87
|
+ if (fabs(r2[1])>fabs(r1[1])) SWAP_ROWS(r2, r1);
|
|
88
|
+ if (0.0f == r1[1]) return false;
|
|
89
|
+
|
|
90
|
+ /* eliminate second variable */
|
|
91
|
+ m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
|
|
92
|
+ r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
|
|
93
|
+ r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
|
|
94
|
+ s = r1[4]; if (0.0f != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
|
|
95
|
+ s = r1[5]; if (0.0f != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
|
|
96
|
+ s = r1[6]; if (0.0f != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
|
|
97
|
+ s = r1[7]; if (0.0f != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }
|
|
98
|
+
|
|
99
|
+ /* choose pivot - or die */
|
|
100
|
+ if (fabs(r3[2])>fabs(r2[2])) SWAP_ROWS(r3, r2);
|
|
101
|
+ if (0.0f == r2[2]) return false;
|
|
102
|
+
|
|
103
|
+ /* eliminate third variable */
|
|
104
|
+ m3 = r3[2]/r2[2];
|
|
105
|
+ r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
|
|
106
|
+ r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
|
|
107
|
+ r3[7] -= m3 * r2[7];
|
|
108
|
+
|
|
109
|
+ /* last check */
|
|
110
|
+ if (0.0f == r3[3]) return false;
|
|
111
|
+
|
|
112
|
+ s = 1.0f/r3[3]; /* now back substitute row 3 */
|
|
113
|
+ r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;
|
|
114
|
+
|
|
115
|
+ m2 = r2[3]; /* now back substitute row 2 */
|
|
116
|
+ s = 1.0f/r2[2];
|
|
117
|
+ r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
|
|
118
|
+ r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
|
|
119
|
+ m1 = r1[3];
|
|
120
|
+ r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
|
|
121
|
+ r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
|
|
122
|
+ m0 = r0[3];
|
|
123
|
+ r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
|
|
124
|
+ r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
|
|
125
|
+
|
|
126
|
+ m1 = r1[2]; /* now back substitute row 1 */
|
|
127
|
+ s = 1.0f/r1[1];
|
|
128
|
+ r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
|
|
129
|
+ r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
|
|
130
|
+ m0 = r0[2];
|
|
131
|
+ r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
|
|
132
|
+ r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
|
|
133
|
+
|
|
134
|
+ m0 = r0[1]; /* now back substitute row 0 */
|
|
135
|
+ s = 1.0f/r0[0];
|
|
136
|
+ r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
|
|
137
|
+ r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
|
|
138
|
+
|
|
139
|
+ MAT(out,0,0) = r0[4];
|
|
140
|
+ MAT(out,0,1) = r0[5], MAT(out,0,2) = r0[6];
|
|
141
|
+ MAT(out,0,3) = r0[7], MAT(out,1,0) = r1[4];
|
|
142
|
+ MAT(out,1,1) = r1[5], MAT(out,1,2) = r1[6];
|
|
143
|
+ MAT(out,1,3) = r1[7], MAT(out,2,0) = r2[4];
|
|
144
|
+ MAT(out,2,1) = r2[5], MAT(out,2,2) = r2[6];
|
|
145
|
+ MAT(out,2,3) = r2[7], MAT(out,3,0) = r3[4];
|
|
146
|
+ MAT(out,3,1) = r3[5], MAT(out,3,2) = r3[6];
|
|
147
|
+ MAT(out,3,3) = r3[7];
|
|
148
|
+
|
|
149
|
+ return true;
|
|
150
|
+#undef MAT
|
|
151
|
+#undef SWAP_ROWS
|
|
152
|
+}
|
|
153
|
+
|
|
154
|
+void Matrix::getMatrix(matrix_t mat) {
|
|
155
|
+ copy(mMatrix, mat);
|
|
156
|
+}
|
|
157
|
+
|
|
158
|
+void Matrix::getTransposeMatrix(matrix_t m) {
|
|
159
|
+ m[ 0]= mMatrix[0]; m[ 1]= mMatrix[4]; m[ 2]= mMatrix[ 8]; m[ 3]=mMatrix[12];
|
|
160
|
+ m[ 4]= mMatrix[1]; m[ 5]= mMatrix[5]; m[ 6]= mMatrix[ 9]; m[ 7]=mMatrix[13];
|
|
161
|
+ m[ 8]= mMatrix[2]; m[ 9]= mMatrix[6]; m[10]= mMatrix[10]; m[11]=mMatrix[14];
|
|
162
|
+ m[12]= mMatrix[3]; m[13]= mMatrix[7]; m[14]= mMatrix[11]; m[15]=mMatrix[15];
|
|
163
|
+}
|
|
164
|
+
|
|
165
|
+Matrix Matrix::multiply(const Matrix &a, const Matrix &b) {
|
|
166
|
+ Matrix c;
|
|
167
|
+ multiply(a.mMatrix, b.mMatrix, c.mMatrix);
|
|
168
|
+ return c;
|
|
169
|
+}
|
|
170
|
+
|
|
171
|
+Matrix Matrix::operator *(const Matrix &a) {
|
|
172
|
+ return multiply(a, *this);
|
|
173
|
+}
|
|
174
|
+
|
|
175
|
+Vector3d Matrix::operator *(Vector3d v) {
|
|
176
|
+ vec_t x = v.mVec[0], y = v.mVec[1], z = v.mVec[2];
|
|
177
|
+
|
|
178
|
+#ifdef COLUMN_ORDER
|
|
179
|
+ return Vector3d(mMatrix[0]*x + mMatrix[4]*y + mMatrix[ 8]*z + mMatrix[12],
|
|
180
|
+ mMatrix[1]*x + mMatrix[5]*y + mMatrix[ 9]*z + mMatrix[13],
|
|
181
|
+ mMatrix[2]*x + mMatrix[6]*y + mMatrix[10]*z + mMatrix[14]);
|
|
182
|
+#else
|
|
183
|
+ return Vector3d(mMatrix[0]*x + mMatrix[1]*y + mMatrix[ 2]*z + mMatrix[ 3],
|
|
184
|
+ mMatrix[4]*x + mMatrix[5]*y + mMatrix[ 6]*z + mMatrix[ 7],
|
|
185
|
+ mMatrix[8]*x + mMatrix[9]*y + mMatrix[10]*z + mMatrix[11]);
|
|
186
|
+#endif
|
|
187
|
+}
|
|
188
|
+
|
|
189
|
+void Matrix::multiply3v(vec3_t v, vec3_t result) {
|
|
190
|
+ vec_t x = v[0], y = v[1], z = v[2];
|
|
191
|
+
|
|
192
|
+ result[0] = mMatrix[0]*x + mMatrix[1]*y + mMatrix[ 2]*z + mMatrix[ 3];
|
|
193
|
+ result[1] = mMatrix[4]*x + mMatrix[5]*y + mMatrix[ 6]*z + mMatrix[ 7];
|
|
194
|
+ result[2] = mMatrix[8]*x + mMatrix[9]*y + mMatrix[10]*z + mMatrix[11];
|
|
195
|
+}
|
|
196
|
+
|
|
197
|
+void Matrix::multiply4v(vec4_t v, vec4_t result) {
|
|
198
|
+ vec_t x = v[0], y = v[1], z = v[2], w = v[3];
|
|
199
|
+
|
|
200
|
+ result[0] = mMatrix[ 0]*x + mMatrix[ 1]*y + mMatrix[ 2]*z + mMatrix[ 3]*w;
|
|
201
|
+ result[1] = mMatrix[ 4]*x + mMatrix[ 5]*y + mMatrix[ 6]*z + mMatrix[ 7]*w;
|
|
202
|
+ result[2] = mMatrix[ 8]*x + mMatrix[ 9]*y + mMatrix[10]*z + mMatrix[11]*w;
|
|
203
|
+ result[3] = mMatrix[12]*x + mMatrix[13]*y + mMatrix[14]*z + mMatrix[15]*w;
|
|
204
|
+}
|
|
205
|
+
|
|
206
|
+void Matrix::print() {
|
|
207
|
+ printf("{\n%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n%f %f %f %f\n}\n",
|
|
208
|
+#ifdef COLUMN_ORDER
|
|
209
|
+ mMatrix[0], mMatrix[4], mMatrix[ 8], mMatrix[12],
|
|
210
|
+ mMatrix[1], mMatrix[5], mMatrix[ 9], mMatrix[13],
|
|
211
|
+ mMatrix[2], mMatrix[6], mMatrix[10], mMatrix[14],
|
|
212
|
+ mMatrix[3], mMatrix[7], mMatrix[11], mMatrix[15]);
|
|
213
|
+#else
|
|
214
|
+ mMatrix[ 0], mMatrix[ 1], mMatrix[ 2], mMatrix[ 3],
|
|
215
|
+ mMatrix[ 4], mMatrix[ 5], mMatrix[ 6], mMatrix[ 7],
|
|
216
|
+ mMatrix[ 8], mMatrix[ 9], mMatrix[10], mMatrix[11],
|
|
217
|
+ mMatrix[12], mMatrix[13], mMatrix[14], mMatrix[15]);
|
|
218
|
+#endif
|
|
219
|
+}
|
|
220
|
+
|
|
221
|
+bool Matrix::isIdentity() {
|
|
222
|
+ // Hhhmm... floating point using direct comparisons
|
|
223
|
+ /*
|
|
224
|
+ if (mMatrix[ 0] == 1 && mMatrix[ 1] == 0 && mMatrix[ 2] == 0 &&
|
|
225
|
+ mMatrix[ 3] == 0 && mMatrix[ 4] == 0 && mMatrix[ 5] == 1 &&
|
|
226
|
+ mMatrix[ 6] == 0 && mMatrix[ 7] == 0 && mMatrix[ 8] == 0 &&
|
|
227
|
+ mMatrix[ 9] == 0 && mMatrix[10] == 1 && mMatrix[11] == 0 &&
|
|
228
|
+ mMatrix[12] == 0 && mMatrix[13] == 0 && mMatrix[14] == 0 &&
|
|
229
|
+ mMatrix[15] == 1)
|
|
230
|
+ return true;
|
|
231
|
+ */
|
|
232
|
+ if (equalEpsilon(mMatrix[ 0], 1.0) && equalEpsilon(mMatrix[ 1], 0.0) && equalEpsilon(mMatrix[ 2], 0.0) &&
|
|
233
|
+ equalEpsilon(mMatrix[ 3], 0.0) && equalEpsilon(mMatrix[ 4], 0.0) && equalEpsilon(mMatrix[ 5], 1.0) &&
|
|
234
|
+ equalEpsilon(mMatrix[ 6], 0.0) && equalEpsilon(mMatrix[ 7], 0.0) && equalEpsilon(mMatrix[ 8], 0.0) &&
|
|
235
|
+ equalEpsilon(mMatrix[ 9], 0.0) && equalEpsilon(mMatrix[10], 1.0) && equalEpsilon(mMatrix[11], 0.0) &&
|
|
236
|
+ equalEpsilon(mMatrix[12], 0.0) && equalEpsilon(mMatrix[13], 0.0) && equalEpsilon(mMatrix[14], 0.0) &&
|
|
237
|
+ equalEpsilon(mMatrix[15], 1.0))
|
|
238
|
+ return true;
|
|
239
|
+
|
|
240
|
+ return false;
|
|
241
|
+}
|
|
242
|
+
|
|
243
|
+void Matrix::setMatrix(matrix_t mat) {
|
|
244
|
+ copy(mat, mMatrix);
|
|
245
|
+}
|
|
246
|
+
|
|
247
|
+void Matrix::setIdentity() {
|
|
248
|
+ mMatrix[ 0] = 1; mMatrix[ 1] = 0; mMatrix[ 2] = 0; mMatrix[ 3] = 0;
|
|
249
|
+ mMatrix[ 4] = 0; mMatrix[ 5] = 1; mMatrix[ 6] = 0; mMatrix[ 7] = 0;
|
|
250
|
+ mMatrix[ 8] = 0; mMatrix[ 9] = 0; mMatrix[10] = 1; mMatrix[11] = 0;
|
|
251
|
+ mMatrix[12] = 0; mMatrix[13] = 0; mMatrix[14] = 0; mMatrix[15] = 1;
|
|
252
|
+}
|
|
253
|
+
|
|
254
|
+void Matrix::scale(const vec_t *xyz) {
|
|
255
|
+ scale(xyz[0], xyz[1], xyz[2]);
|
|
256
|
+}
|
|
257
|
+
|
|
258
|
+void Matrix::scale(vec_t sx, vec_t sy, vec_t sz) {
|
|
259
|
+ matrix_t smatrix;
|
|
260
|
+ matrix_t tmp;
|
|
261
|
+
|
|
262
|
+ smatrix[ 0] = sx; smatrix[ 1] = 0; smatrix[ 2] = 0; smatrix[ 3] = 0;
|
|
263
|
+ smatrix[ 4] = 0; smatrix[ 5] = sy; smatrix[ 6] = 0; smatrix[ 7] = 0;
|
|
264
|
+ smatrix[ 8] = 0; smatrix[ 9] = 0; smatrix[10] = sz; smatrix[11] = 0;
|
|
265
|
+ smatrix[12] = 0; smatrix[13] = 0; smatrix[14] = 0; smatrix[15] = 1;
|
|
266
|
+
|
|
267
|
+ copy(mMatrix, tmp);
|
|
268
|
+ multiply(tmp, smatrix, mMatrix);
|
|
269
|
+}
|
|
270
|
+
|
|
271
|
+void Matrix::rotate(const vec_t *xyz) {
|
|
272
|
+ rotate(xyz[0], xyz[1], xyz[2]);
|
|
273
|
+}
|
|
274
|
+
|
|
275
|
+void Matrix::rotate(vec_t ax, vec_t ay, vec_t az) {
|
|
276
|
+ matrix_t xmat, ymat, zmat, tmp, tmp2;
|
|
277
|
+
|
|
278
|
+ xmat[ 0]=1; xmat[ 1]=0; xmat[ 2]=0; xmat[ 3]=0;
|
|
279
|
+ xmat[ 4]=0; xmat[ 5]=cosf(ax); xmat[ 6]=sinf(ax); xmat[ 7]=0;
|
|
280
|
+ xmat[ 8]=0; xmat[ 9]=-sinf(ax); xmat[10]=cosf(ax); xmat[11]=0;
|
|
281
|
+ xmat[12]=0; xmat[13]=0; xmat[14]=0; xmat[15]=1;
|
|
282
|
+
|
|
283
|
+ ymat[ 0]=cosf(ay); ymat[ 1]=0; ymat[ 2]=-sinf(ay); ymat[ 3]=0;
|
|
284
|
+ ymat[ 4]=0; ymat[ 5]=1; ymat[ 6]=0; ymat[ 7]=0;
|
|
285
|
+ ymat[ 8]=sinf(ay); ymat[ 9]=0; ymat[10]=cosf(ay); ymat[11]=0;
|
|
286
|
+ ymat[12]=0; ymat[13]=0; ymat[14]=0; ymat[15]=1;
|
|
287
|
+
|
|
288
|
+ zmat[ 0]=cosf(az); zmat[ 1]=sinf(az); zmat[ 2]=0; zmat[ 3]=0;
|
|
289
|
+ zmat[ 4]=-sinf(az); zmat[ 5]=cosf(az); zmat[ 6]=0; zmat[ 7]=0;
|
|
290
|
+ zmat[ 8]=0; zmat[ 9]=0; zmat[10]=1; zmat[11]=0;
|
|
291
|
+ zmat[12]=0; zmat[13]=0; zmat[14]=0; zmat[15]=1;
|
|
292
|
+
|
|
293
|
+ multiply(mMatrix, ymat, tmp);
|
|
294
|
+ multiply(tmp, xmat, tmp2);
|
|
295
|
+ multiply(tmp2, zmat, mMatrix);
|
|
296
|
+}
|
|
297
|
+
|
|
298
|
+void Matrix::translate(const vec_t *xyz) {
|
|
299
|
+ translate(xyz[0], xyz[1], xyz[2]);
|
|
300
|
+}
|
|
301
|
+
|
|
302
|
+void Matrix::translate(vec_t tx, vec_t ty, vec_t tz) {
|
|
303
|
+ matrix_t tmat, tmp;
|
|
304
|
+
|
|
305
|
+ tmat[ 0]=1; tmat[ 1]=0; tmat[ 2]=0; tmat[ 3]=0;
|
|
306
|
+ tmat[ 4]=0; tmat[ 5]=1; tmat[ 6]=0; tmat[ 7]=0;
|
|
307
|
+ tmat[ 8]=0; tmat[ 9]=0; tmat[10]=1; tmat[11]=0;
|
|
308
|
+ tmat[12]=tx; tmat[13]=ty; tmat[14]=tz; tmat[15]=1;
|
|
309
|
+
|
|
310
|
+ copy(mMatrix, tmp);
|
|
311
|
+ multiply(tmp, tmat, mMatrix);
|
|
312
|
+}
|
|
313
|
+
|
|
314
|
+void Matrix::copy(matrix_t source, matrix_t dest) {
|
|
315
|
+ for (int i = 0; i < 16; i++)
|
|
316
|
+ dest[i] = source[i];
|
|
317
|
+}
|
|
318
|
+
|
|
319
|
+void Matrix::multiply(const matrix_t a, const matrix_t b, matrix_t result) {
|
|
320
|
+ /* Generated code for matrix mult
|
|
321
|
+ * Code used:
|
|
322
|
+
|
|
323
|
+ // char order is argument
|
|
324
|
+ int i, j, k;
|
|
325
|
+ if (order == 'r') {
|
|
326
|
+ printf("// Row order\n");
|
|
327
|
+ } else {
|
|
328
|
+ printf("// Column order\n");
|
|
329
|
+ }
|
|
330
|
+ for (i = 0; i < 4; ++i) {
|
|
331
|
+ for (j = 0; j < 4; ++j) {
|
|
332
|
+ if (order == 'r') {
|
|
333
|
+ printf("result[%2i] = ", j+i*4);
|
|
334
|
+ } else {
|
|
335
|
+ printf("result[%2i] = ", j+i*4);
|
|
336
|
+ }
|
|
337
|
+ for (k = 0; k < 4; ++k) {
|
|
338
|
+ if (order == 'r') {
|
|
339
|
+ printf("a[%2i] * b[%2i]%s",
|
|
340
|
+ k+i*4, j+k*4, (k == 3) ? ";\n" : " + ");
|
|
341
|
+ } else {
|
|
342
|
+ printf("a[%2i] * b[%2i]%s",
|
|
343
|
+ i+k*4, k+j*4, (k == 3) ? ";\n" : " + ");
|
|
344
|
+ }
|
|
345
|
+ //sum+=(elements[i+k*4]*m.elements[k+j*4]);
|
|
346
|
+ }
|
|
347
|
+ //result.elements[i+j*4]=sum;
|
|
348
|
+ }
|
|
349
|
+ printf("\n");
|
|
350
|
+ }
|
|
351
|
+ printf("\n");
|
|
352
|
+ printf("// Transpose\n");
|
|
353
|
+ for(i = 0; i < 4; ++i) {
|
|
354
|
+ for (j = 0; j < 4; ++j) {
|
|
355
|
+ printf("a[%2i] = b[%2i]%s",
|
|
356
|
+ j+i*4, i+j*4, (j == 3) ? ";\n" : "; ");
|
|
357
|
+ }
|
|
358
|
+ }
|
|
359
|
+
|
|
360
|
+ * was in test/Matrix.cpp
|
|
361
|
+ */
|
|
362
|
+#ifdef COLUMN_ORDER
|
|
363
|
+ /* Column order */
|
|
364
|
+ result[ 0] = a[ 0] * b[ 0] + a[ 4] * b[ 1] + a[ 8] * b[ 2] + a[12] * b[ 3];
|
|
365
|
+ result[ 1] = a[ 0] * b[ 4] + a[ 4] * b[ 5] + a[ 8] * b[ 6] + a[12] * b[ 7];
|
|
366
|
+ result[ 2] = a[ 0] * b[ 8] + a[ 4] * b[ 9] + a[ 8] * b[10] + a[12] * b[11];
|
|
367
|
+ result[ 3] = a[ 0] * b[12] + a[ 4] * b[13] + a[ 8] * b[14] + a[12] * b[15];
|
|
368
|
+
|
|
369
|
+ result[ 4] = a[ 1] * b[ 0] + a[ 5] * b[ 1] + a[ 9] * b[ 2] + a[13] * b[ 3];
|
|
370
|
+ result[ 5] = a[ 1] * b[ 4] + a[ 5] * b[ 5] + a[ 9] * b[ 6] + a[13] * b[ 7];
|
|
371
|
+ result[ 6] = a[ 1] * b[ 8] + a[ 5] * b[ 9] + a[ 9] * b[10] + a[13] * b[11];
|
|
372
|
+ result[ 7] = a[ 1] * b[12] + a[ 5] * b[13] + a[ 9] * b[14] + a[13] * b[15];
|
|
373
|
+
|
|
374
|
+ result[ 8] = a[ 2] * b[ 0] + a[ 6] * b[ 1] + a[10] * b[ 2] + a[14] * b[ 3];
|
|
375
|
+ result[ 9] = a[ 2] * b[ 4] + a[ 6] * b[ 5] + a[10] * b[ 6] + a[14] * b[ 7];
|
|
376
|
+ result[10] = a[ 2] * b[ 8] + a[ 6] * b[ 9] + a[10] * b[10] + a[14] * b[11];
|
|
377
|
+ result[11] = a[ 2] * b[12] + a[ 6] * b[13] + a[10] * b[14] + a[14] * b[15];
|
|
378
|
+
|
|
379
|
+ result[12] = a[ 3] * b[ 0] + a[ 7] * b[ 1] + a[11] * b[ 2] + a[15] * b[ 3];
|
|
380
|
+ result[13] = a[ 3] * b[ 4] + a[ 7] * b[ 5] + a[11] * b[ 6] + a[15] * b[ 7];
|
|
381
|
+ result[14] = a[ 3] * b[ 8] + a[ 7] * b[ 9] + a[11] * b[10] + a[15] * b[11];
|
|
382
|
+ result[15] = a[ 3] * b[12] + a[ 7] * b[13] + a[11] * b[14] + a[15] * b[15];
|
|
383
|
+#else
|
|
384
|
+ /* Row order */
|
|
385
|
+ result[ 0] = a[ 0] * b[ 0] + a[ 1] * b[ 4] + a[ 2] * b[ 8] + a[ 3] * b[12];
|
|
386
|
+ result[ 1] = a[ 0] * b[ 1] + a[ 1] * b[ 5] + a[ 2] * b[ 9] + a[ 3] * b[13];
|
|
387
|
+ result[ 2] = a[ 0] * b[ 2] + a[ 1] * b[ 6] + a[ 2] * b[10] + a[ 3] * b[14];
|
|
388
|
+ result[ 3] = a[ 0] * b[ 3] + a[ 1] * b[ 7] + a[ 2] * b[11] + a[ 3] * b[15];
|
|
389
|
+
|
|
390
|
+ result[ 4] = a[ 4] * b[ 0] + a[ 5] * b[ 4] + a[ 6] * b[ 8] + a[ 7] * b[12];
|
|
391
|
+ result[ 5] = a[ 4] * b[ 1] + a[ 5] * b[ 5] + a[ 6] * b[ 9] + a[ 7] * b[13];
|
|
392
|
+ result[ 6] = a[ 4] * b[ 2] + a[ 5] * b[ 6] + a[ 6] * b[10] + a[ 7] * b[14];
|
|
393
|
+ result[ 7] = a[ 4] * b[ 3] + a[ 5] * b[ 7] + a[ 6] * b[11] + a[ 7] * b[15];
|
|
394
|
+
|
|
395
|
+ result[ 8] = a[ 8] * b[ 0] + a[ 9] * b[ 4] + a[10] * b[ 8] + a[11] * b[12];
|
|
396
|
+ result[ 9] = a[ 8] * b[ 1] + a[ 9] * b[ 5] + a[10] * b[ 9] + a[11] * b[13];
|
|
397
|
+ result[10] = a[ 8] * b[ 2] + a[ 9] * b[ 6] + a[10] * b[10] + a[11] * b[14];
|
|
398
|
+ result[11] = a[ 8] * b[ 3] + a[ 9] * b[ 7] + a[10] * b[11] + a[11] * b[15];
|
|
399
|
+
|
|
400
|
+ result[12] = a[12] * b[ 0] + a[13] * b[ 4] + a[14] * b[ 8] + a[15] * b[12];
|
|
401
|
+ result[13] = a[12] * b[ 1] + a[13] * b[ 5] + a[14] * b[ 9] + a[15] * b[13];
|
|
402
|
+ result[14] = a[12] * b[ 2] + a[13] * b[ 6] + a[14] * b[10] + a[15] * b[14];
|
|
403
|
+ result[15] = a[12] * b[ 3] + a[13] * b[ 7] + a[14] * b[11] + a[15] * b[15];
|
|
404
|
+#endif
|
|
405
|
+}
|
|
406
|
+
|