matrix.h
3.43 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
#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, rhs, 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()(int row, int col) {
return at(row, col);
}
matrix<T> operator=(T rhs) {
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
for(size_t c = 0; c < rhs.C; c++){
for(size_t r = 0; r < R; r++){
inner = (T)0;
for(size_t 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;
}
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();
}
};
} //end namespace rts
#endif