matrix.h
4.91 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
#ifndef RTS_MATRIX_H
#define RTS_MATRIX_H
//#include "rts/vector.h"
#include <string.h>
#include <iostream>
#include <stim/math/vector.h>
#include <stim/math/vec3.h>
//#include <stim/cuda/cudatools/callable.h>
namespace stim{
template <class T>
class matrix {
//the matrix will be stored in column-major order (compatible with OpenGL)
T* M; //pointer to the matrix data
size_t R; //number of rows
size_t C; //number of colums
void init(size_t rows, size_t cols){
R = rows;
C = cols;
M = (T*) malloc (R * C * sizeof(T)); //allocate space for the matrix
}
T& at(size_t row, size_t col){
return M[col * R + row];
}
public:
matrix(size_t rows, size_t cols) {
init(rows, cols); //initialize memory
}
matrix(size_t rows, size_t cols, T* data) {
init(rows, cols);
memcpy(M, data, R * C * sizeof(T));
}
matrix(const matrix<T>& cpy){
init(cpy.R, cpy.C);
memcpy(M, cpy.M, R * C * sizeof(T));
}
~matrix() {
R = C = 0;
if(M) free(M);
}
size_t rows(){
return R;
}
size_t cols(){
return C;
}
T& operator()(size_t row, size_t col) {
return at(row, col);
}
matrix<T> operator=(T rhs) {
if (&rhs == this)
return *this;
init(R, C);
size_t N = R * C;
for(size_t n=0; n<N; n++)
M[n] = rhs;
return *this;
}
matrix<T>& operator=(matrix<T> rhs){
init(rhs.R, rhs.C);
memcpy(M, rhs.M, R * C * sizeof(T));
return *this;
}
//element-wise operations
matrix<T> operator+(T rhs) {
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;
}
matrix<T> operator-(T rhs) {
return operator+(-rhs); //add the negative of rhs
}
matrix<T> operator*(T rhs) {
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;
}
matrix<T> operator/(T rhs) {
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;
}
//matrix multiplication
matrix<T> operator*(matrix<T> rhs){
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);
}
matrix<T> result(R, rhs.C); //create the output matrix
T inner; //stores the running inner product
size_t c, r, i;
for(c = 0; c < rhs.C; c++){
for(r = 0; r < R; r++){
inner = (T)0;
for(i = 0; i < C; i++){
inner += at(r, i) * rhs.at(i, c);
}
result.M[c * R + r] = inner;
}
}
return result;
}
//returns a pointer to the raw matrix data (in column major format)
T* data(){
return M;
}
//return a transposed matrix
matrix<T> transpose(){
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];
}
}
return result;
}
/// Sort rows of the matrix by the specified indices
matrix<T> sort_rows(size_t* idx) {
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
matrix<T> sort_cols(size_t* idx) {
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;
}
std::string toStr() {
std::stringstream ss;
for(int r = 0; r < R; r++) {
ss << "| ";
for(int c=0; c<C; c++) {
ss << M[c * R + r] << " ";
}
ss << "|" << std::endl;
}
return ss.str();
}
std::string csv() {
std::stringstream csvss;
for (size_t i = 0; i < R; i++) {
csvss << M[i];
for (size_t j = 0; j < C; j++) csvss << ", " << M[j * R + i];
csvss << std::endl;
}
return csvss.str();
}
//save the data as a CSV file
void csv(std::string filename) {
ofstream basisfile(filename.c_str());
basisfile << csv();
basisfile.close();
}
static matrix<T> I(size_t N) {
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;
}
};
} //end namespace rts
#endif