57729e5b
David Mayerich
changed directory...
|
1
2
3
4
5
6
|
#ifndef RTS_MATRIX_H
#define RTS_MATRIX_H
//#include "rts/vector.h"
#include <string.h>
#include <iostream>
|
cce7daf9
Pavel Govyadinov
added glObj.h to ...
|
7
|
#include <stim/math/vector.h>
|
308a743c
David Mayerich
fixed class compa...
|
8
|
#include <stim/math/vec3.h>
|
1c92dffe
David Mayerich
simplified CUDA i...
|
9
|
//#include <stim/cuda/cudatools/callable.h>
|
57729e5b
David Mayerich
changed directory...
|
10
|
|
8a86bd56
David Mayerich
changed rts names...
|
11
|
namespace stim{
|
57729e5b
David Mayerich
changed directory...
|
12
|
|
6c4afcac
David Mayerich
introduced a gene...
|
13
14
|
template <class T>
class matrix {
|
57729e5b
David Mayerich
changed directory...
|
15
|
//the matrix will be stored in column-major order (compatible with OpenGL)
|
6c4afcac
David Mayerich
introduced a gene...
|
16
17
18
19
|
T* M; //pointer to the matrix data
size_t R; //number of rows
size_t C; //number of colums
|
65c1aac0
David Mayerich
fixed memory leak...
|
20
|
/*void init(size_t rows, size_t cols){
|
6c4afcac
David Mayerich
introduced a gene...
|
21
22
|
R = rows;
C = cols;
|
65c1aac0
David Mayerich
fixed memory leak...
|
23
24
25
26
27
28
29
30
31
32
33
|
if (R == 0 || C == 0) M = NULL;
else
M = (T*)malloc(R * C * sizeof(T)); //allocate space for the matrix
}*/
T get(const size_t row, const size_t col) const {
if (row >= R || col >= C) {
std::cout << "ERROR: row or column out of range." << std::endl;
exit(1);
}
return M[col * R + row];
|
6c4afcac
David Mayerich
introduced a gene...
|
34
|
}
|
57729e5b
David Mayerich
changed directory...
|
35
|
|
6c4afcac
David Mayerich
introduced a gene...
|
36
|
T& at(size_t row, size_t col){
|
65c1aac0
David Mayerich
fixed memory leak...
|
37
38
39
40
|
if (row >= R || col >= C) {
std::cout << "ERROR: row or column out of range." << std::endl;
exit(1);
}
|
6c4afcac
David Mayerich
introduced a gene...
|
41
|
return M[col * R + row];
|
57729e5b
David Mayerich
changed directory...
|
42
43
|
}
|
6c4afcac
David Mayerich
introduced a gene...
|
44
|
public:
|
65c1aac0
David Mayerich
fixed memory leak...
|
45
46
47
48
49
50
|
matrix() {
R = 0;
C = 0;
M = NULL;
}
|
6c4afcac
David Mayerich
introduced a gene...
|
51
|
matrix(size_t rows, size_t cols) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
52
53
|
R = rows;
C = cols;
|
ef1f6218
David Mayerich
fixed some memory...
|
54
55
56
|
M = NULL;
if (R * C > 0)
M = (T*) malloc(R * C * sizeof(T));
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
57
58
|
}
|
ef1f6218
David Mayerich
fixed some memory...
|
59
|
matrix(size_t rows, size_t cols, const T* data) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
60
61
|
R = rows;
C = cols;
|
ef1f6218
David Mayerich
fixed some memory...
|
62
63
64
|
M = NULL;
if (R * C > 0)
M = (T*)malloc(R * C * sizeof(T));
|
dd5aab2f
David Mayerich
fixed errors in s...
|
65
|
memcpy(M, data, R * C * sizeof(T));
|
8e56a0a7
Pavel Govyadinov
Added the propose...
|
66
|
}
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
67
|
|
6c4afcac
David Mayerich
introduced a gene...
|
68
|
matrix(const matrix<T>& cpy){
|
ef1f6218
David Mayerich
fixed some memory...
|
69
70
71
|
M = NULL;
if (cpy.R * cpy.C > 0)
M = (T*)malloc(cpy.R * cpy.C * sizeof(T));
|
65c1aac0
David Mayerich
fixed memory leak...
|
72
73
74
75
|
memcpy(M, cpy.M, cpy.R * cpy.C * sizeof(T));
R = cpy.R;
C = cpy.C;
|
6c4afcac
David Mayerich
introduced a gene...
|
76
|
}
|
19a3973b
David Mayerich
updates to matrix...
|
77
|
|
6c4afcac
David Mayerich
introduced a gene...
|
78
|
~matrix() {
|
6c4afcac
David Mayerich
introduced a gene...
|
79
|
if(M) free(M);
|
65c1aac0
David Mayerich
fixed memory leak...
|
80
81
|
M = NULL;
R = C = 0;
|
19a3973b
David Mayerich
updates to matrix...
|
82
83
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
84
|
size_t rows() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
85
|
return R;
|
57729e5b
David Mayerich
changed directory...
|
86
87
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
88
|
size_t cols() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
89
90
|
return C;
}
|
57729e5b
David Mayerich
changed directory...
|
91
|
|
1c92dffe
David Mayerich
simplified CUDA i...
|
92
|
T& operator()(size_t row, size_t col) {
|
6c4afcac
David Mayerich
introduced a gene...
|
93
94
95
|
return at(row, col);
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
96
97
|
matrix<T>& operator=(const T rhs) {
//init(R, C);
|
6c4afcac
David Mayerich
introduced a gene...
|
98
99
100
101
102
103
104
|
size_t N = R * C;
for(size_t n=0; n<N; n++)
M[n] = rhs;
return *this;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
105
106
107
108
109
110
111
112
113
114
|
matrix<T>& operator=(const matrix<T>& rhs){
if (this != &rhs) { //if the matrix isn't self-assigned
T* new_matrix = new T[rhs.R * rhs.C]; //allocate new resources
memcpy(new_matrix, rhs.M, rhs.R * rhs.C * sizeof(T)); //copy the matrix
delete[] M; //delete the previous array
M = new_matrix;
R = rhs.R;
C = rhs.C;
}
|
57729e5b
David Mayerich
changed directory...
|
115
116
|
return *this;
}
|
42bb075c
Jiaming Guo
last changes for ...
|
117
|
|
6c4afcac
David Mayerich
introduced a gene...
|
118
|
//element-wise operations
|
65c1aac0
David Mayerich
fixed memory leak...
|
119
|
matrix<T> operator+(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
120
121
122
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
|
42bb075c
Jiaming Guo
last changes for ...
|
123
|
for(int i=0; i<N; i++)
|
6c4afcac
David Mayerich
introduced a gene...
|
124
125
126
|
result.M[i] = M[i] + rhs; //calculate the operation and assign to result
return result;
|
42bb075c
Jiaming Guo
last changes for ...
|
127
|
}
|
57729e5b
David Mayerich
changed directory...
|
128
|
|
65c1aac0
David Mayerich
fixed memory leak...
|
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
matrix<T> operator+(const matrix<T> rhs) const {
if (R != rhs.R || C != rhs.C) {
std::cout << "ERROR: addition is only defined for matrices that are the same size." << std::endl;
exit(1);
}
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for (int i = 0; i < N; i++)
result.M[i] = M[i] + rhs.M[i]; //calculate the operation and assign to result
return result;
}
matrix<T> operator-(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
144
145
146
|
return operator+(-rhs); //add the negative of rhs
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
matrix<T> operator-(const matrix<T> rhs) const {
return operator+(-rhs);
}
matrix<T> operator-() const {
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for (int i = 0; i < N; i++)
result.M[i] = -M[i]; //calculate the operation and assign to result
return result;
}
matrix<T> operator*(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
162
163
164
165
166
167
168
169
170
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for(int i=0; i<N; i++)
result.M[i] = M[i] * rhs; //calculate the operation and assign to result
return result;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
171
|
matrix<T> operator/(const T rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
172
173
174
175
176
|
matrix<T> result(R, C); //create a result matrix
size_t N = R * C;
for(int i=0; i<N; i++)
result.M[i] = M[i] / rhs; //calculate the operation and assign to result
|
487a9b49
David Mayerich
added the ability...
|
177
|
|
6c4afcac
David Mayerich
introduced a gene...
|
178
179
|
return result;
}
|
57729e5b
David Mayerich
changed directory...
|
180
|
|
6c4afcac
David Mayerich
introduced a gene...
|
181
|
//matrix multiplication
|
65c1aac0
David Mayerich
fixed memory leak...
|
182
|
matrix<T> operator*(const matrix<T> rhs) const {
|
6c4afcac
David Mayerich
introduced a gene...
|
183
184
185
186
187
|
if(C != rhs.R){
std::cout<<"ERROR: matrix multiplication is undefined for matrices of size ";
std::cout<<"[ "<<R<<" x "<<C<<" ] and [ "<<rhs.R<<" x "<<rhs.C<<"]"<<std::endl;
exit(1);
}
|
57729e5b
David Mayerich
changed directory...
|
188
|
|
6c4afcac
David Mayerich
introduced a gene...
|
189
190
|
matrix<T> result(R, rhs.C); //create the output matrix
T inner; //stores the running inner product
|
1c92dffe
David Mayerich
simplified CUDA i...
|
191
192
193
|
size_t c, r, i;
for(c = 0; c < rhs.C; c++){
for(r = 0; r < R; r++){
|
6c4afcac
David Mayerich
introduced a gene...
|
194
|
inner = (T)0;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
195
|
for(i = 0; i < C; i++){
|
65c1aac0
David Mayerich
fixed memory leak...
|
196
|
inner += get(r, i) * rhs.get(i, c);
|
6c4afcac
David Mayerich
introduced a gene...
|
197
198
199
200
|
}
result.M[c * R + r] = inner;
}
}
|
57729e5b
David Mayerich
changed directory...
|
201
202
203
|
return result;
}
|
6c4afcac
David Mayerich
introduced a gene...
|
204
205
206
207
|
//returns a pointer to the raw matrix data (in column major format)
T* data(){
return M;
}
|
8e4f8364
David Mayerich
started a new opt...
|
208
|
|
6c4afcac
David Mayerich
introduced a gene...
|
209
|
//return a transposed matrix
|
65c1aac0
David Mayerich
fixed memory leak...
|
210
|
matrix<T> transpose() const {
|
6c4afcac
David Mayerich
introduced a gene...
|
211
212
213
214
215
216
217
|
matrix<T> result(C, R);
size_t c, r;
for(c = 0; c < C; c++){
for(r = 0; r < R; r++){
result.M[r * C + c] = M[c * R + r];
}
}
|
8e4f8364
David Mayerich
started a new opt...
|
218
219
220
|
return result;
}
|
ef1f6218
David Mayerich
fixed some memory...
|
221
222
223
224
225
226
|
// Reshapes the matrix in place
void reshape(size_t rows, size_t cols) {
R = rows;
C = cols;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
227
228
229
230
231
232
233
234
235
236
237
|
///Calculate and return the determinant of the matrix
T det() const {
if (R != C) {
std::cout << "ERROR: a determinant can only be calculated for a square matrix." << std::endl;
exit(1);
}
if (R == 1) return M[0]; //if the matrix only contains one value, return it
int r, c, ri, cia, cib;
T a = 0;
T b = 0;
|
f7cf19e4
David Mayerich
fixed warnings, t...
|
238
|
for (c = 0; c < (int)C; c++) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
239
240
|
for (r = 0; r < R; r++) {
ri = r;
|
f7cf19e4
David Mayerich
fixed warnings, t...
|
241
242
|
cia = (r + c) % (int)C;
cib = ((int)C - 1 - r) % (int)C;
|
65c1aac0
David Mayerich
fixed memory leak...
|
243
244
245
246
247
248
249
|
a += get(ri, cia);
b += get(ri, cib);
}
}
return a - b;
}
|
d299c720
David Mayerich
added row and col...
|
250
251
252
253
254
255
256
257
|
/// Sum all elements in the matrix
T sum() const {
size_t N = R * C; //calculate the number of elements in the matrix
T s = (T)0; //allocate a register to store the sum
for (size_t n = 0; n < N; n++) s += M[n]; //perform the summation
return s;
}
|
1c92dffe
David Mayerich
simplified CUDA i...
|
258
|
/// Sort rows of the matrix by the specified indices
|
65c1aac0
David Mayerich
fixed memory leak...
|
259
|
matrix<T> sort_rows(size_t* idx) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
260
261
262
263
264
265
266
267
268
269
270
|
matrix<T> result(C, R); //create the output matrix
size_t r, c;
for (c = 0; c < C; c++) { //for each column
for (r = 0; r < R; r++) { //for each row element
result.M[c * R + r] = M[c * R + idx[r]]; //copy each element of the row into its new position
}
}
return result;
}
/// Sort columns of the matrix by the specified indices
|
65c1aac0
David Mayerich
fixed memory leak...
|
271
|
matrix<T> sort_cols(size_t* idx) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
272
273
274
275
276
277
278
279
|
matrix<T> result(C, R);
size_t c;
for (c = 0; c < C; c++) { //for each column
memcpy(&result.M[c * R], &M[idx[c] * R], sizeof(T) * R); //copy the entire column from this matrix to the appropriate location
}
return result;
}
|
d299c720
David Mayerich
added row and col...
|
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
|
/// Return the column specified by index i
matrix<T> col(size_t i) {
matrix<T> c(R, 1); //create a single column matrix
memcpy(c.data(), &data()[R*i], C * sizeof(T)); //copy the column
return c;
}
/// Return the row specified by index i
matrix<T> row(size_t i) {
matrix<T> r(1, C); //create a single row matrix
for (size_t c = 0; c < C; c++)
r(0, c) = at(i, c);
return r;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
295
|
std::string toStr() const {
|
57729e5b
David Mayerich
changed directory...
|
296
297
|
std::stringstream ss;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
298
|
for(int r = 0; r < R; r++) {
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
299
|
ss << "| ";
|
1c92dffe
David Mayerich
simplified CUDA i...
|
300
|
for(int c=0; c<C; c++) {
|
6c4afcac
David Mayerich
introduced a gene...
|
301
|
ss << M[c * R + r] << " ";
|
57729e5b
David Mayerich
changed directory...
|
302
|
}
|
b2d10ff1
Pavel Govyadinov
final edits, expa...
|
303
|
ss << "|" << std::endl;
|
57729e5b
David Mayerich
changed directory...
|
304
|
}
|
57729e5b
David Mayerich
changed directory...
|
305
306
|
return ss.str();
}
|
19a3973b
David Mayerich
updates to matrix...
|
307
|
|
65c1aac0
David Mayerich
fixed memory leak...
|
308
309
|
void csv(std::ostream& out) const {
//std::stringstream csvss;
|
c495d065
David Mayerich
Fixed stim::matri...
|
310
|
for (size_t i = 0; i < R; i++) {
|
65c1aac0
David Mayerich
fixed memory leak...
|
311
312
313
314
|
out << std::fixed << M[i];
for (size_t j = 1; j < C; j++)
out << ", " << std::fixed << M[j * R + i];
out << std::endl;
|
1c92dffe
David Mayerich
simplified CUDA i...
|
315
|
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
316
317
318
319
320
321
322
323
|
//return csvss.str();
}
std::string csv() const {
std::stringstream csvss;
int digits = std::numeric_limits<double>::max_digits10;
csvss.precision(digits);
csv(csvss);
|
1c92dffe
David Mayerich
simplified CUDA i...
|
324
325
326
|
return csvss.str();
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
327
328
|
|
1c92dffe
David Mayerich
simplified CUDA i...
|
329
|
//save the data as a CSV file
|
65c1aac0
David Mayerich
fixed memory leak...
|
330
|
void csv(std::string filename) const {
|
1c92dffe
David Mayerich
simplified CUDA i...
|
331
332
333
334
335
|
ofstream basisfile(filename.c_str());
basisfile << csv();
basisfile.close();
}
|
dd5aab2f
David Mayerich
fixed errors in s...
|
336
|
static matrix<T> I(size_t N) {
|
dd5aab2f
David Mayerich
fixed errors in s...
|
337
338
339
340
341
342
343
344
|
matrix<T> result(N, N); //create the identity matrix
memset(result.M, 0, N * N * sizeof(T)); //set the entire matrix to zero
for (size_t n = 0; n < N; n++) {
result(n, n) = (T)1; //set the diagonal component to 1
}
return result;
}
|
65c1aac0
David Mayerich
fixed memory leak...
|
345
346
347
348
349
350
351
352
353
354
355
356
357
|
//loads a matrix from a stream in CSV format
void csv(std::istream& in) {
size_t c, r;
T v;
for (r = 0; r < R; r++) {
for (c = 0; c < C; c++) {
in >> v;
if (in.peek() == ',') in.seekg(1, std::ios::cur);
at(r, c) = v;;
}
}
}
|
57729e5b
David Mayerich
changed directory...
|
358
359
360
361
|
};
} //end namespace rts
|
57729e5b
David Mayerich
changed directory...
|
362
363
|
#endif
|