openGPMP
Open Source Mathematics Package
mtx_naive.cpp
Go to the documentation of this file.
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 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  int rows = A.size();
58  int cols = rows > 0 ? A.size() / rows : 0;
59 
60  for (int i = 0; i < rows; ++i) {
61  for (int j = 0; j < cols; ++j) {
62  // perform matrix addition
63  C[i * cols + j] = A[i * cols + j] + B[i * cols + j];
64  }
65  }
66 }
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 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  const int size = A.size();
93 
94  for (int i = 0; i < size; ++i) {
95  for (int j = 0; j < size; ++j) {
96  // perform matrix addition
97  C[i][j] = A[i][j] + B[i][j];
98  }
99  }
100 }
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 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  const int size = A.size();
124 
125  for (int i = 0; i < size; ++i) {
126  for (int j = 0; j < size; ++j) {
127  // Perform matrix subtraction
128  C[i][j] = A[i][j] - B[i][j];
129  }
130  }
131 }
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 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  int64_t nrows = A.size();
158  int64_t ncols = A[0].size();
159 
160  for (int64_t i = 0; i < nrows; ++i) {
161  for (int64_t j = 0; j < ncols; ++j) {
162  C[i][j] = 0.0;
163  for (int64_t k = 0; k < ncols; ++k) {
164  C[i][j] += A[i][k] * B[k][j];
165  }
166  }
167  }
168 }
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 }*/
gpmp::linalg::Mtx::std_mtx_mult
void std_mtx_mult(const T *A, const T *B, U *C, int rows_a, int cols_a, int cols_b)
Definition: mtx.hpp:540
gpmp::linalg::Mtx::std_mtx_add
void std_mtx_add(const T *A, const T *B, T *C, int rows, int cols)
Perform matrix addition on two matrices as flat arrays.
Definition: mtx.hpp:510
gpmp::linalg::Mtx::std_mtx_sub
void std_mtx_sub(const T *A, const T *B, T *C, int rows, int cols)
Perform matrix subtraction on two matrices as flat arrays.
Definition: mtx.hpp:529
python.linalg.C
list C
Definition: linalg.py:24
python.linalg.A
list A
Definition: linalg.py:22
python.linalg.B
list B
Definition: linalg.py:23
test_linalg.rows
int rows
Definition: test_linalg.py:13
test_linalg.cols
int cols
Definition: test_linalg.py:14
Generated by doxygen 1.9.1