Line data Source code
1 : /*************************************************************************
2 : *
3 : * Project
4 : * _____ _____ __ __ _____
5 : * / ____| __ \| \/ | __ \
6 : * ___ _ __ ___ _ __ | | __| |__) | \ / | |__) |
7 : * / _ \| '_ \ / _ \ '_ \| | |_ | ___/| |\/| | ___/
8 : *| (_) | |_) | __/ | | | |__| | | | | | | |
9 : * \___/| .__/ \___|_| |_|\_____|_| |_| |_|_|
10 : * | |
11 : * |_|
12 : *
13 : * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
14 : *
15 : * This software is licensed as described in the file LICENSE, which
16 : * you should have received as part of this distribution. The terms
17 : * among other details are referenced in the official documentation
18 : * seen here : https://akielaries.github.io/openGPMP/ along with
19 : * important files seen in this project.
20 : *
21 : * You may opt to use, copy, modify, merge, publish, distribute
22 : * and/or sell copies of the Software, and permit persons to whom
23 : * the Software is furnished to do so, under the terms of the
24 : * LICENSE file. As this is an Open Source effort, all implementations
25 : * must be of the same methodology.
26 : *
27 : *
28 : *
29 : * This software is distributed on an AS IS basis, WITHOUT
30 : * WARRANTY OF ANY KIND, either express or implied.
31 : *
32 : ************************************************************************/
33 : #include <cassert>
34 : #include <cstddef>
35 : #include <cstdint>
36 : #include <iostream>
37 : #include <openGPMP/linalg/mtx.hpp>
38 : #include <vector>
39 :
40 : /************************************************************************
41 : *
42 : * Standard/Naive Matrix Operations on Arrays
43 : *
44 : ************************************************************************/
45 :
46 : /************************************************************************
47 : *
48 : * Standard/Naive Matrix Operations on vector<>
49 : *
50 : ************************************************************************/
51 : // naive matrix addition algorithm on vector<>
52 : template <typename T>
53 0 : void gpmp::linalg::Mtx::std_mtx_add(const std::vector<T> &A,
54 : const std::vector<T> &B,
55 : std::vector<T> &C) {
56 : // MTX A AND B MUST BE SAME SIZE
57 0 : int rows = A.size();
58 0 : int cols = rows > 0 ? A.size() / rows : 0;
59 :
60 0 : for (int i = 0; i < rows; ++i) {
61 0 : for (int j = 0; j < cols; ++j) {
62 : // perform matrix addition
63 0 : C[i * cols + j] = A[i * cols + j] + B[i * cols + j];
64 : }
65 : }
66 0 : }
67 :
68 : // instantiations for types accepted by templated std_mtx_add function for
69 : // flat vectors
70 : template void gpmp::linalg::Mtx::std_mtx_add(const std::vector<int> &A,
71 : const std::vector<int> &B,
72 : std::vector<int> &C);
73 :
74 : template void gpmp::linalg::Mtx::std_mtx_add(const std::vector<double> &A,
75 : const std::vector<double> &B,
76 : std::vector<double> &C);
77 :
78 : template void gpmp::linalg::Mtx::std_mtx_add(const std::vector<float> &A,
79 : const std::vector<float> &B,
80 : std::vector<float> &C);
81 :
82 : /************************************************************************
83 : *
84 : * Standard/Naive Matrix Operations on vector<vector>
85 : *
86 : ************************************************************************/
87 : // naive matrix addition algorithm on vector<vector>
88 : template <typename T>
89 6 : void gpmp::linalg::Mtx::std_mtx_add(const std::vector<std::vector<T>> &A,
90 : const std::vector<std::vector<T>> &B,
91 : std::vector<std::vector<T>> &C) {
92 6 : const int size = A.size();
93 :
94 4230 : for (int i = 0; i < size; ++i) {
95 4206720 : for (int j = 0; j < size; ++j) {
96 : // perform matrix addition
97 4202496 : C[i][j] = A[i][j] + B[i][j];
98 : }
99 : }
100 6 : }
101 :
102 : // instantiations for types accepted by templated std_mtx_add function
103 : template void
104 : gpmp::linalg::Mtx::std_mtx_add(const std::vector<std::vector<int>> &A,
105 : const std::vector<std::vector<int>> &B,
106 : std::vector<std::vector<int>> &C);
107 :
108 : template void
109 : gpmp::linalg::Mtx::std_mtx_add(const std::vector<std::vector<double>> &A,
110 : const std::vector<std::vector<double>> &B,
111 : std::vector<std::vector<double>> &C);
112 :
113 : template void
114 : gpmp::linalg::Mtx::std_mtx_add(const std::vector<std::vector<float>> &A,
115 : const std::vector<std::vector<float>> &B,
116 : std::vector<std::vector<float>> &C);
117 :
118 : // naive matrix subtraction algorithm
119 : template <typename T>
120 0 : void gpmp::linalg::Mtx::std_mtx_sub(const std::vector<std::vector<T>> &A,
121 : const std::vector<std::vector<T>> &B,
122 : std::vector<std::vector<T>> &C) {
123 0 : const int size = A.size();
124 :
125 0 : for (int i = 0; i < size; ++i) {
126 0 : for (int j = 0; j < size; ++j) {
127 : // Perform matrix subtraction
128 0 : C[i][j] = A[i][j] - B[i][j];
129 : }
130 : }
131 0 : }
132 :
133 : // instantiations for types accepted by templated std_mtx_sub function
134 : template void
135 : gpmp::linalg::Mtx::std_mtx_sub(const std::vector<std::vector<int>> &A,
136 : const std::vector<std::vector<int>> &B,
137 : std::vector<std::vector<int>> &C);
138 :
139 : template void
140 : gpmp::linalg::Mtx::std_mtx_sub(const std::vector<std::vector<double>> &A,
141 : const std::vector<std::vector<double>> &B,
142 : std::vector<std::vector<double>> &C);
143 :
144 : template void
145 : gpmp::linalg::Mtx::std_mtx_sub(const std::vector<std::vector<float>> &A,
146 : const std::vector<std::vector<float>> &B,
147 : std::vector<std::vector<float>> &C);
148 :
149 : // naive matrix multiplication algorithm
150 : template <typename T>
151 1 : void gpmp::linalg::Mtx::std_mtx_mult(const std::vector<std::vector<T>> &A,
152 : const std::vector<std::vector<T>> &B,
153 : std::vector<std::vector<T>> &C) {
154 : assert(A.size() == B.size());
155 : assert(A[0].size() == B[0].size());
156 :
157 1 : int64_t nrows = A.size();
158 1 : int64_t ncols = A[0].size();
159 :
160 1025 : for (int64_t i = 0; i < nrows; ++i) {
161 1049600 : for (int64_t j = 0; j < ncols; ++j) {
162 1048576 : C[i][j] = 0.0;
163 1074790400 : for (int64_t k = 0; k < ncols; ++k) {
164 1073741824 : C[i][j] += A[i][k] * B[k][j];
165 : }
166 : }
167 : }
168 1 : }
169 :
170 : // instantiations for types accepted by templated std_mtx_mult function
171 : template void
172 : gpmp::linalg::Mtx::std_mtx_mult(const std::vector<std::vector<int>> &A,
173 : const std::vector<std::vector<int>> &B,
174 : std::vector<std::vector<int>> &C);
175 :
176 : template void
177 : gpmp::linalg::Mtx::std_mtx_mult(const std::vector<std::vector<double>> &A,
178 : const std::vector<std::vector<double>> &B,
179 : std::vector<std::vector<double>> &C);
180 :
181 : template void
182 : gpmp::linalg::Mtx::std_mtx_mult(const std::vector<std::vector<float>> &A,
183 : const std::vector<std::vector<float>> &B,
184 : std::vector<std::vector<float>> &C);
185 :
186 : /*
187 : // naive implementation of Strassen matrix multiplication algorithm
188 : void gpmp::linalg::Mtx::std_mtx_mult_strass(
189 : const std::vector<std::vector<int>> &A,
190 : const std::vector<std::vector<int>> &B, std::vector<std::vector<int>> &C) {
191 : int n = A.size();
192 :
193 : // base case: If the matrix size is 1x1, perform regular multiplication
194 : if (n == 1) {
195 : C[0][0] = A[0][0] * B[0][0];
196 : }
197 :
198 : // splitting the matrices into quadrants
199 : int half = n / 2;
200 : std::vector<std::vector<int>> A11(half, std::vector<int>(half));
201 : std::vector<std::vector<int>> A12(half, std::vector<int>(half));
202 : std::vector<std::vector<int>> A21(half, std::vector<int>(half));
203 : std::vector<std::vector<int>> A22(half, std::vector<int>(half));
204 : std::vector<std::vector<int>> B11(half, std::vector<int>(half));
205 : std::vector<std::vector<int>> B12(half, std::vector<int>(half));
206 : std::vector<std::vector<int>> B21(half, std::vector<int>(half));
207 : std::vector<std::vector<int>> B22(half, std::vector<int>(half));
208 : std::vector<std::vector<int>> C11(half, std::vector<int>(half));
209 : std::vector<std::vector<int>> C12(half, std::vector<int>(half));
210 : std::vector<std::vector<int>> C21(half, std::vector<int>(half));
211 : std::vector<std::vector<int>> C22(half, std::vector<int>(half));
212 :
213 : // dividing the input matrices into quadrants
214 : for (int i = 0; i < half; i++) {
215 : for (int j = 0; j < half; j++) {
216 : A11[i][j] = A[i][j];
217 : A12[i][j] = A[i][j + half];
218 : A21[i][j] = A[i + half][j];
219 : A22[i][j] = A[i + half][j + half];
220 :
221 : B11[i][j] = B[i][j];
222 : B12[i][j] = B[i][j + half];
223 : B21[i][j] = B[i + half][j];
224 : B22[i][j] = B[i + half][j + half];
225 : }
226 : }
227 :
228 : // recursive calls for sub-matrix multiplication
229 : std::vector<std::vector<int>> M1(half, std::vector<int>(half));
230 : std::vector<std::vector<int>> M2(half, std::vector<int>(half));
231 : std::vector<std::vector<int>> M3(half, std::vector<int>(half));
232 : std::vector<std::vector<int>> M4(half, std::vector<int>(half));
233 : std::vector<std::vector<int>> M5(half, std::vector<int>(half));
234 : std::vector<std::vector<int>> M6(half, std::vector<int>(half));
235 : std::vector<std::vector<int>> M7(half, std::vector<int>(half));
236 :
237 : // M1 = (A11 + A22) * (B11 + B22)
238 : std_mtx_mult_strass(std_mtx_add(A11, A22), std_mtx_add(B11, B22), M1);
239 : // M2 = (A21 + A22) * B11
240 : std_mtx_mult_strass(std_mtx_add(A21, A22), B11, M2);
241 : // M3 = A11 * (B12 - B22)
242 : std_mtx_mult_strass(A11, std_mtx_sub(B12, B22), M3);
243 : // M4 = A22 * (B21 - B11)
244 : std_mtx_mult_strass(A22, std_mtx_sub(B21, B11), M4);
245 : // M5 = (A11 + A12) * B22
246 : std_mtx_mult_strass(std_mtx_add(A11, A12), B22, M5);
247 : // M6 = (A21 - A11) * (B11 + B12)
248 : std_mtx_mult_strass(std_mtx_sub(A21, A11), std_mtx_add(B11, B12), M6);
249 : // M7 = (A12 - A22) * (B21 + B22)
250 : std_mtx_mult_strass(std_mtx_sub(A12, A22), std_mtx_add(B21, B22), M7);
251 :
252 : // Computing the sub-matrices of the result matrix
253 : std::vector<std::vector<int>> C11_temp =
254 :
255 : std_mtx_add(std_mtx_sub(std_mtx_add(M1, M4), M5), M7);
256 :
257 : std::vector<std::vector<int>> C12_temp = std_mtx_add(M3, M5);
258 : std::vector<std::vector<int>> C21_temp = std_mtx_add(M2, M4);
259 : std::vector<std::vector<int>> C22_temp =
260 : std_mtx_add(std_mtx_sub(std_mtx_add(M1, M3), M2), M6);
261 :
262 : // Combining the sub-matrices to form the resulting matrix
263 : for (int i = 0; i < half; i++) {
264 : for (int j = 0; j < half; j++) {
265 : C[i][j] = C11_temp[i][j];
266 : C[i][j + half] = C12_temp[i][j];
267 : C[i + half][j] = C21_temp[i][j];
268 : C[i + half][j + half] = C22_temp[i][j];
269 : }
270 : }
271 : }*/
|
