Multiscale Universal Interface  2.0
A Concurrent Framework for Coupling Heterogeneous Solvers
matrix_arithmetic.h
Go to the documentation of this file.
1 /*****************************************************************************
2 * Multiscale Universal Interface Code Coupling Library *
3 * *
4 * Copyright (C) 2023 W. Liu *
5 * *
6 * This software is jointly licensed under the Apache License, Version 2.0 *
7 * and the GNU General Public License version 3, you may use it according *
8 * to either. *
9 * *
10 * ** Apache License, version 2.0 ** *
11 * *
12 * Licensed under the Apache License, Version 2.0 (the "License"); *
13 * you may not use this file except in compliance with the License. *
14 * You may obtain a copy of the License at *
15 * *
16 * http://www.apache.org/licenses/LICENSE-2.0 *
17 * *
18 * Unless required by applicable law or agreed to in writing, software *
19 * distributed under the License is distributed on an "AS IS" BASIS, *
20 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
21 * See the License for the specific language governing permissions and *
22 * limitations under the License. *
23 * *
24 * ** GNU General Public License, version 3 ** *
25 * *
26 * This program is free software: you can redistribute it and/or modify *
27 * it under the terms of the GNU General Public License as published by *
28 * the Free Software Foundation, either version 3 of the License, or *
29 * (at your option) any later version. *
30 * *
31 * This program is distributed in the hope that it will be useful, *
32 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
33 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
34 * GNU General Public License for more details. *
35 * *
36 * You should have received a copy of the GNU General Public License *
37 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
38 *****************************************************************************/
39 
47 #ifndef MUI_MATRIX_ARITHMETIC_H_
48 #define MUI_MATRIX_ARITHMETIC_H_
49 
50 #include <cassert>
51 #include <math.h>
52 
53 namespace mui {
54 namespace linalg {
55 
56 // **************************************************
57 // ************ Public member functions *************
58 // **************************************************
59 
60 // Overload addition operator to perform sparse matrix addition
61 template<typename ITYPE, typename VTYPE>
62 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::operator+(sparse_matrix<ITYPE,VTYPE> &addend) {
63 
64  if (rows_ != addend.rows_ || cols_ != addend.cols_) {
65  std::cerr << "MUI Error [matrix_arithmetic.h]: matrix size mismatch during matrix addition" << std::endl;
66  std::abort();
67  }
68 
69  if (addend.matrix_format_ != matrix_format_) {
70  addend.format_conversion(this->get_format(), true, true, "overwrite");
71  } else {
72  if (!addend.is_sorted_unique("matrix_arithmetic.h", "operator+()")){
73  if (addend.matrix_format_ == format::COO) {
74  addend.sort_coo(true, true, "overwrite");
75  } else if (addend.matrix_format_ == format::CSR) {
76  addend.sort_csr(true, "overwrite");
77  } else if (addend.matrix_format_ == format::CSC) {
78  addend.sort_csc(true, "overwrite");
79  } else {
80  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised addend matrix format for matrix operator+()" << std::endl;
81  std::cerr << " Please set the addend matrix_format_ as:" << std::endl;
82  std::cerr << " format::COO: COOrdinate format" << std::endl;
83  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
84  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
85  std::abort();
86  }
87  }
88  }
89 
90  if (!this->is_sorted_unique("matrix_arithmetic.h", "operator+()")){
91  if (matrix_format_ == format::COO) {
92  this->sort_coo(true, true, "overwrite");
93  } else if (matrix_format_ == format::CSR) {
94  this->sort_csr(true, "overwrite");
95  } else if (matrix_format_ == format::CSC) {
96  this->sort_csc(true, "overwrite");
97  } else {
98  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator+()" << std::endl;
99  std::cerr << " Please set the matrix_format_ as:" << std::endl;
100  std::cerr << " format::COO: COOrdinate format" << std::endl;
101  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
102  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
103  std::abort();
104  }
105  }
106 
107  // Create a new sparse matrix object for the result
108  sparse_matrix<ITYPE,VTYPE> res(rows_, cols_, this->get_format());
109 
110  if (matrix_format_ == format::COO) {
111 
112  // Perform element-wise addition of the COO vectors
113  res.matrix_coo.values_.reserve(matrix_coo.values_.size() + addend.matrix_coo.values_.size());
114  res.matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size() + addend.matrix_coo.row_indices_.size());
115  res.matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size() + addend.matrix_coo.col_indices_.size());
116 
117  // Insert the COO vectors of the initial sparse matrix to the result sparse matrix
118  res.matrix_coo.values_ = std::vector<VTYPE>(matrix_coo.values_.begin(), matrix_coo.values_.end());
119  res.matrix_coo.row_indices_ = std::vector<ITYPE>(matrix_coo.row_indices_.begin(), matrix_coo.row_indices_.end());
120  res.matrix_coo.col_indices_ = std::vector<ITYPE>(matrix_coo.col_indices_.begin(), matrix_coo.col_indices_.end());
121 
122  // Append the addend COO vectors to the result sparse matrix
123  res.matrix_coo.values_.insert(res.matrix_coo.values_.end(), addend.matrix_coo.values_.begin(), addend.matrix_coo.values_.end());
124  res.matrix_coo.row_indices_.insert(res.matrix_coo.row_indices_.end(), addend.matrix_coo.row_indices_.begin(), addend.matrix_coo.row_indices_.end());
125  res.matrix_coo.col_indices_.insert(res.matrix_coo.col_indices_.end(), addend.matrix_coo.col_indices_.begin(), addend.matrix_coo.col_indices_.end());
126 
127  // Sort and deduplicate the result
128  res.sort_coo(true, true, "plus");
129  res.nnz_ = res.matrix_coo.values_.size();
130 
131  } else if (matrix_format_ == format::CSR) {
132 
133  // Perform element-wise addition of the CSR vectors
134  res.matrix_csr.values_.reserve(matrix_csr.values_.size() + addend.matrix_csr.values_.size());
135  res.matrix_csr.row_ptrs_.resize(rows_ + 1);
136  res.matrix_csr.col_indices_.reserve(matrix_csr.col_indices_.size() + addend.matrix_csr.col_indices_.size());
137 
138  ITYPE row = 0;
139  while (row < rows_) {
140  ITYPE start = matrix_csr.row_ptrs_[row];
141  ITYPE end = matrix_csr.row_ptrs_[row + 1];
142 
143  ITYPE addend_start = addend.matrix_csr.row_ptrs_[row];
144  ITYPE addend_end = addend.matrix_csr.row_ptrs_[row + 1];
145 
146  res.matrix_csr.row_ptrs_[0] = 0;
147 
148  // Merge the values and column indices of the two rows
149  ITYPE i = start;
150  ITYPE j = addend_start;
151  while (i < end && j < addend_end) {
152  ITYPE col = matrix_csr.col_indices_[i];
153  ITYPE addend_col = addend.matrix_csr.col_indices_[j];
154 
155  if (col == addend_col) {
156  // Add the corresponding values if the columns match
157  if (std::abs(matrix_csr.values_[i] + addend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()){
158  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i] + addend.matrix_csr.values_[j]);
159  res.matrix_csr.col_indices_.emplace_back(col);
160  }
161  i++;
162  j++;
163  } else if (col < addend_col) {
164  // Add the current value from the initial matrix
165  if (std::abs(matrix_csr.values_[i]) >= std::numeric_limits<VTYPE>::min()){
166  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i]);
167  res.matrix_csr.col_indices_.emplace_back(col);
168  }
169  i++;
170  } else {
171  // Add the current value from the addend matrix
172  if (std::abs(addend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()){
173  res.matrix_csr.values_.emplace_back(addend.matrix_csr.values_[j]);
174  res.matrix_csr.col_indices_.emplace_back(addend_col);
175  }
176  j++;
177  }
178  }
179 
180  // Add any remaining elements from the initial matrix
181  for (; i < end; i++) {
182  if (std::abs(matrix_csr.values_[i]) >= std::numeric_limits<VTYPE>::min()){
183  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i]);
184  res.matrix_csr.col_indices_.emplace_back(matrix_csr.col_indices_[i]);
185  }
186  }
187 
188  // Add any remaining elements from the addend matrix
189  for (; j < addend_end; j++) {
190  if (std::abs(addend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()){
191  res.matrix_csr.values_.emplace_back(addend.matrix_csr.values_[j]);
192  res.matrix_csr.col_indices_.emplace_back(addend.matrix_csr.col_indices_[j]);
193  }
194  }
195 
196  // Update the row pointer
197  res.nnz_ = res.matrix_csr.col_indices_.size();
198  res.matrix_csr.row_ptrs_[row + 1] = res.nnz_;
199 
200  row++;
201  }
202 
203  } else if (matrix_format_ == format::CSC) {
204 
205  // Perform element-wise addition of the CSC vectors
206  res.matrix_csc.values_.reserve(matrix_csc.values_.size() + addend.matrix_csc.values_.size());
207  res.matrix_csc.row_indices_.reserve(matrix_csc.row_indices_.size() + addend.matrix_csc.row_indices_.size());
208  res.matrix_csc.col_ptrs_.resize(cols_ + 1);
209 
210  ITYPE column = 0;
211  while (column < cols_) {
212  ITYPE start = matrix_csc.col_ptrs_[column];
213  ITYPE end = matrix_csc.col_ptrs_[column + 1];
214 
215  ITYPE addend_start = addend.matrix_csc.col_ptrs_[column];
216  ITYPE addend_end = addend.matrix_csc.col_ptrs_[column + 1];
217 
218  res.matrix_csc.col_ptrs_[0] = 0;
219 
220  // Merge the values and row indices of the two columns
221  ITYPE i = start;
222  ITYPE j = addend_start;
223  while (i < end && j < addend_end) {
224  ITYPE row = matrix_csc.row_indices_[i];
225  ITYPE addend_row = addend.matrix_csc.row_indices_[j];
226 
227  if (row == addend_row) {
228  // Add the corresponding values if the columns match
229  if (std::abs(matrix_csc.values_[i] + addend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
230  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i] + addend.matrix_csc.values_[j]);
231  res.matrix_csc.row_indices_.emplace_back(row);
232  }
233  i++;
234  j++;
235  } else if (row < addend_row) {
236  // Add the current value from the initial matrix
237  if (std::abs(matrix_csc.values_[i]) >= std::numeric_limits<VTYPE>::min()) {
238  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i]);
239  res.matrix_csc.row_indices_.emplace_back(row);
240  }
241  i++;
242  } else {
243  // Add the current value from the addend matrix
244  if (std::abs(addend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
245  res.matrix_csc.values_.emplace_back(addend.matrix_csc.values_[j]);
246  res.matrix_csc.row_indices_.emplace_back(addend_row);
247  }
248  j++;
249  }
250  }
251 
252  // Add any remaining elements from the initial matrix
253  for (; i < end; i++) {
254  if (std::abs(matrix_csc.values_[i]) >= std::numeric_limits<VTYPE>::min()){
255  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i]);
256  res.matrix_csc.row_indices_.emplace_back(matrix_csc.row_indices_[i]);
257  }
258  }
259 
260  // Add any remaining elements from the addend matrix
261  for (; j < addend_end; j++) {
262  if (std::abs(addend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()){
263  res.matrix_csc.values_.emplace_back(addend.matrix_csc.values_[j]);
264  res.matrix_csc.row_indices_.emplace_back(addend.matrix_csc.row_indices_[j]);
265  }
266  }
267 
268  // Update the column pointer
269  res.nnz_ = res.matrix_csc.row_indices_.size();
270  res.matrix_csc.col_ptrs_[column + 1] = res.nnz_;
271 
272  column++;
273  }
274 
275  } else {
276  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator+()" << std::endl;
277  std::cerr << " Please set the matrix_format_ as:" << std::endl;
278  std::cerr << " format::COO: COOrdinate format" << std::endl;
279  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
280  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
281  std::abort();
282  }
283 
284  return res;
285 }
286 
287 // Overload subtraction operator to perform sparse matrix subtraction
288 template<typename ITYPE, typename VTYPE>
289 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::operator-(sparse_matrix<ITYPE,VTYPE> &subtrahend) {
290  if (rows_ != subtrahend.rows_ || cols_ != subtrahend.cols_) {
291  std::cerr << "MUI Error [matrix_arithmetic.h]: matrix size mismatch during matrix subtraction" << std::endl;
292  std::abort();
293  }
294 
295  if (subtrahend.matrix_format_ != matrix_format_) {
296  subtrahend.format_conversion(this->get_format(), true, true, "overwrite");
297  } else {
298  if (!subtrahend.is_sorted_unique("matrix_arithmetic.h", "operator-()")){
299  if (subtrahend.matrix_format_ == format::COO) {
300  subtrahend.sort_coo(true, true, "overwrite");
301  } else if (subtrahend.matrix_format_ == format::CSR) {
302  subtrahend.sort_csr(true, "overwrite");
303  } else if (subtrahend.matrix_format_ == format::CSC) {
304  subtrahend.sort_csc(true, "overwrite");
305  } else {
306  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised subtrahend matrix format for matrix operator-()" << std::endl;
307  std::cerr << " Please set the subtrahend matrix_format_ as:" << std::endl;
308  std::cerr << " format::COO: COOrdinate format" << std::endl;
309  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
310  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
311  std::abort();
312  }
313  }
314  }
315 
316  if (!this->is_sorted_unique("matrix_arithmetic.h", "operator-()")){
317  if (matrix_format_ == format::COO) {
318  this->sort_coo(true, true, "overwrite");
319  } else if (matrix_format_ == format::CSR) {
320  this->sort_csr(true, "overwrite");
321  } else if (matrix_format_ == format::CSC) {
322  this->sort_csc(true, "overwrite");
323  } else {
324  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator-()" << std::endl;
325  std::cerr << " Please set the matrix_format_ as:" << std::endl;
326  std::cerr << " format::COO: COOrdinate format" << std::endl;
327  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
328  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
329  std::abort();
330  }
331  }
332 
333  // Create a new sparse matrix object for the result
334  sparse_matrix<ITYPE,VTYPE> res(rows_, cols_, this->get_format());
335 
336  if (matrix_format_ == format::COO) {
337 
338  std::vector<VTYPE> subtrahend_value;
339  subtrahend_value.reserve(subtrahend.matrix_coo.values_.size());
340 
341  for (VTYPE &element : subtrahend.matrix_coo.values_) {
342  subtrahend_value.emplace_back(element*(-1));
343  }
344 
345  // Perform element-wise subtrahend of the COO vectors
346  res.matrix_coo.values_.reserve(matrix_coo.values_.size() + subtrahend.matrix_coo.values_.size());
347  res.matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size() + subtrahend.matrix_coo.row_indices_.size());
348  res.matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size() + subtrahend.matrix_coo.col_indices_.size());
349 
350  // Insert the COO vectors of the initial sparse matrix to the result sparse matrix
351  res.matrix_coo.values_ = std::vector<VTYPE>(matrix_coo.values_.begin(), matrix_coo.values_.end());
352  res.matrix_coo.row_indices_ = std::vector<ITYPE>(matrix_coo.row_indices_.begin(), matrix_coo.row_indices_.end());
353  res.matrix_coo.col_indices_ = std::vector<ITYPE>(matrix_coo.col_indices_.begin(), matrix_coo.col_indices_.end());
354 
355  // Append the subtrahend COO vectors to the result sparse matrix
356  res.matrix_coo.values_.insert(res.matrix_coo.values_.end(), subtrahend_value.begin(), subtrahend_value.end());
357  res.matrix_coo.row_indices_.insert(res.matrix_coo.row_indices_.end(), subtrahend.matrix_coo.row_indices_.begin(), subtrahend.matrix_coo.row_indices_.end());
358  res.matrix_coo.col_indices_.insert(res.matrix_coo.col_indices_.end(), subtrahend.matrix_coo.col_indices_.begin(), subtrahend.matrix_coo.col_indices_.end());
359 
360  // Sort and deduplicate the result
361  res.sort_coo(true, true, "plus");
362 
363  } else if (matrix_format_ == format::CSR) {
364 
365  // Perform element-wise subtraction of the CSR vectors
366  res.matrix_csr.values_.reserve(matrix_csr.values_.size() + subtrahend.matrix_csr.values_.size());
367  res.matrix_csr.row_ptrs_.resize(rows_ + 1);
368  res.matrix_csr.col_indices_.reserve(matrix_csr.col_indices_.size() + subtrahend.matrix_csr.col_indices_.size());
369 
370  ITYPE row = 0;
371  while (row < rows_) {
372  ITYPE start = matrix_csr.row_ptrs_[row];
373  ITYPE end = matrix_csr.row_ptrs_[row + 1];
374 
375  ITYPE subtrahend_start = subtrahend.matrix_csr.row_ptrs_[row];
376  ITYPE subtrahend_end = subtrahend.matrix_csr.row_ptrs_[row + 1];
377 
378  res.matrix_csr.row_ptrs_[0] = 0;
379 
380  // Merge the values and column indices of the two rows
381  ITYPE i = start;
382  ITYPE j = subtrahend_start;
383  while (i < end && j < subtrahend_end) {
384  ITYPE col = matrix_csr.col_indices_[i];
385  ITYPE subtrahend_col = subtrahend.matrix_csr.col_indices_[j];
386 
387  if (col == subtrahend_col) {
388  // Add the corresponding values if the columns match
389  if (std::abs(matrix_csr.values_[i] - subtrahend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
390  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i] - subtrahend.matrix_csr.values_[j]);
391  res.matrix_csr.col_indices_.emplace_back(col);
392  }
393  i++;
394  j++;
395  } else if (col < subtrahend_col) {
396  // Add the current value from the initial matrix
397  if (std::abs(matrix_csr.values_[i]) >= std::numeric_limits<VTYPE>::min()) {
398  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i]);
399  res.matrix_csr.col_indices_.emplace_back(col);
400  }
401  i++;
402  } else {
403  // Add the current value from the subtrahend matrix
404  if (std::abs(-subtrahend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
405  res.matrix_csr.values_.emplace_back(-subtrahend.matrix_csr.values_[j]);
406  res.matrix_csr.col_indices_.emplace_back(subtrahend_col);
407  }
408  j++;
409  }
410  }
411 
412  // Add any remaining elements from the initial matrix
413  for (; i < end; i++) {
414  if (std::abs(matrix_csr.values_[i]) >= std::numeric_limits<VTYPE>::min()) {
415  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i]);
416  res.matrix_csr.col_indices_.emplace_back(matrix_csr.col_indices_[i]);
417  }
418  }
419 
420  // Add any remaining elements from the subtrahend matrix
421  for (; j < subtrahend_end; j++) {
422  if (std::abs(-subtrahend.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
423  res.matrix_csr.values_.emplace_back(-subtrahend.matrix_csr.values_[j]);
424  res.matrix_csr.col_indices_.emplace_back(subtrahend.matrix_csr.col_indices_[j]);
425  }
426  }
427 
428  // Update the row pointer
429  res.nnz_ = res.matrix_csr.col_indices_.size();
430  res.matrix_csr.row_ptrs_[row + 1] = res.nnz_;
431 
432  row++;
433  }
434 
435  } else if (matrix_format_ == format::CSC) {
436 
437  // Perform element-wise subtraction of the CSC vectors
438  res.matrix_csc.values_.reserve(matrix_csc.values_.size() + subtrahend.matrix_csc.values_.size());
439  res.matrix_csc.row_indices_.reserve(matrix_csc.row_indices_.size() + subtrahend.matrix_csc.row_indices_.size());
440  res.matrix_csc.col_ptrs_.resize(cols_ + 1);
441 
442  ITYPE column = 0;
443  while (column < cols_) {
444  ITYPE start = matrix_csc.col_ptrs_[column];
445  ITYPE end = matrix_csc.col_ptrs_[column + 1];
446 
447  ITYPE subtrahend_start = subtrahend.matrix_csc.col_ptrs_[column];
448  ITYPE subtrahend_end = subtrahend.matrix_csc.col_ptrs_[column + 1];
449 
450  res.matrix_csc.col_ptrs_[0] = 0;
451 
452  // Merge the values and row indices of the two columns
453  ITYPE i = start;
454  ITYPE j = subtrahend_start;
455  while (i < end && j < subtrahend_end) {
456  ITYPE row = matrix_csc.row_indices_[i];
457  ITYPE subtrahend_row = subtrahend.matrix_csc.row_indices_[j];
458 
459  if (row == subtrahend_row) {
460  // Add the corresponding values if the columns match
461  if (std::abs(matrix_csc.values_[i] - subtrahend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
462  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i] - subtrahend.matrix_csc.values_[j]);
463  res.matrix_csc.row_indices_.emplace_back(row);
464  }
465  i++;
466  j++;
467  } else if (row < subtrahend_row) {
468  // Add the current value from the initial matrix
469  if (std::abs(matrix_csc.values_[i]) >= std::numeric_limits<VTYPE>::min()) {
470  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i]);
471  res.matrix_csc.row_indices_.emplace_back(row);
472  }
473  i++;
474  } else {
475  // Add the current value from the subtrahend matrix
476  if (std::abs(-subtrahend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
477  res.matrix_csc.values_.emplace_back(-subtrahend.matrix_csc.values_[j]);
478  res.matrix_csc.row_indices_.emplace_back(subtrahend_row);
479  }
480  j++;
481  }
482  }
483 
484  // Add any remaining elements from the initial matrix
485  for (; i < end; i++) {
486  if (std::abs(matrix_csc.values_[i]) >= std::numeric_limits<VTYPE>::min()){
487  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i]);
488  res.matrix_csc.row_indices_.emplace_back(matrix_csc.row_indices_[i]);
489  }
490  }
491 
492  // Add any remaining elements from the subtrahend matrix
493  for (; j < subtrahend_end; j++) {
494  if (std::abs(-subtrahend.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()){
495  res.matrix_csc.values_.emplace_back(-subtrahend.matrix_csc.values_[j]);
496  res.matrix_csc.row_indices_.emplace_back(subtrahend.matrix_csc.row_indices_[j]);
497  }
498  }
499 
500  // Update the column pointer
501  res.nnz_ = res.matrix_csc.row_indices_.size();
502  res.matrix_csc.col_ptrs_[column + 1] = res.nnz_;
503 
504  column++;
505  }
506 
507  } else {
508  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator-()" << std::endl;
509  std::cerr << " Please set the matrix_format_ as:" << std::endl;
510  std::cerr << " format::COO: COOrdinate format" << std::endl;
511  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
512  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
513  std::abort();
514  }
515 
516  return res;
517 }
518 
519 // Overload multiplication operator to perform sparse matrix multiplication
520 template<typename ITYPE, typename VTYPE>
521 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::operator*(sparse_matrix<ITYPE,VTYPE> &multiplicand) {
522 
523  if (cols_ != multiplicand.rows_) {
524  std::cerr << "MUI Error [matrix_arithmetic.h]: matrix size mismatch during matrix multiplication" << std::endl;
525  std::abort();
526  }
527 
528  if (multiplicand.matrix_format_ != matrix_format_) {
529  multiplicand.format_conversion(this->get_format(), true, true, "overwrite");
530  } else {
531  if (!multiplicand.is_sorted_unique("matrix_arithmetic.h", "operator*()")){
532  if (multiplicand.matrix_format_ == format::COO) {
533  multiplicand.sort_coo(true, true, "overwrite");
534  } else if (multiplicand.matrix_format_ == format::CSR) {
535  multiplicand.sort_csr(true, "overwrite");
536  } else if (multiplicand.matrix_format_ == format::CSC) {
537  multiplicand.sort_csc(true, "overwrite");
538  } else {
539  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised multiplicand matrix format for matrix operator*()" << std::endl;
540  std::cerr << " Please set the multiplicand matrix_format_ as:" << std::endl;
541  std::cerr << " format::COO: COOrdinate format" << std::endl;
542  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
543  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
544  std::abort();
545  }
546  }
547  }
548 
549  if (!this->is_sorted_unique("matrix_arithmetic.h", "operator*()")){
550  if (matrix_format_ == format::COO) {
551  this->sort_coo(true, true, "overwrite");
552  } else if (matrix_format_ == format::CSR) {
553  this->sort_csr(true, "overwrite");
554  } else if (matrix_format_ == format::CSC) {
555  this->sort_csc(true, "overwrite");
556  } else {
557  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator*()" << std::endl;
558  std::cerr << " Please set the matrix_format_ as:" << std::endl;
559  std::cerr << " format::COO: COOrdinate format" << std::endl;
560  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
561  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
562  std::abort();
563  }
564  }
565 
566  // Create a new sparse matrix object for the result
567  sparse_matrix<ITYPE,VTYPE> res(rows_, multiplicand.cols_, this->get_format());
568 
569  if (matrix_format_ == format::COO) {
570 
571  // Perform element-wise multiplication of the COO vectors
572  res.matrix_coo.values_.reserve((matrix_coo.values_.size() <= multiplicand.matrix_coo.values_.size()) ? multiplicand.matrix_coo.values_.size() : matrix_coo.values_.size());
573  res.matrix_coo.row_indices_.reserve((matrix_coo.row_indices_.size() <= multiplicand.matrix_coo.row_indices_.size()) ? multiplicand.matrix_coo.row_indices_.size() : matrix_coo.row_indices_.size());
574  res.matrix_coo.col_indices_.reserve((matrix_coo.col_indices_.size() <= multiplicand.matrix_coo.col_indices_.size()) ? multiplicand.matrix_coo.col_indices_.size() : matrix_coo.col_indices_.size());
575 
576  for (ITYPE i = 0; i < static_cast<ITYPE>(matrix_coo.row_indices_.size()); ++i) {
577  for (ITYPE j = 0; j < static_cast<ITYPE>(multiplicand.matrix_coo.col_indices_.size()); ++j) {
578  if (matrix_coo.col_indices_[i] == multiplicand.matrix_coo.row_indices_[j]) {
579  // Multiply the corresponding values if the columns match
580  VTYPE value = matrix_coo.values_[i] * multiplicand.matrix_coo.values_[j];
581  if (std::abs(value) >= std::numeric_limits<VTYPE>::min()) {
582  res.matrix_coo.values_.emplace_back(value);
583  res.matrix_coo.row_indices_.emplace_back(matrix_coo.row_indices_[i]);
584  res.matrix_coo.col_indices_.emplace_back(multiplicand.matrix_coo.col_indices_[j]);
585  }
586  }
587  }
588  }
589 
590  // Sort and deduplicate the result
591  res.sort_coo(true, true, "plus");
592 
593  res.nnz_ = res.matrix_coo.values_.size();
594 
595  } else if (matrix_format_ == format::CSR) {
596 
597  // Perform element-wise multiplication of the CSR vectors
598  res.matrix_csr.values_.reserve((matrix_csr.values_.size() <= multiplicand.matrix_csr.values_.size()) ? multiplicand.matrix_csr.values_.size() : matrix_csr.values_.size());
599  res.matrix_csr.row_ptrs_.resize(rows_+1);
600  res.matrix_csr.col_indices_.reserve((matrix_csr.col_indices_.size() <= multiplicand.matrix_csr.col_indices_.size()) ? multiplicand.matrix_csr.col_indices_.size() : matrix_csr.col_indices_.size());
601 
602  // Initialize a vector to store the intermediate results
603  std::vector<VTYPE> intermediate(multiplicand.cols_, 0.0);
604 
605  res.matrix_csr.row_ptrs_[0] = 0;
606 
607  // Iterate over each row of the initial matrix
608  for (ITYPE i = 0; i < rows_; ++i) {
609  // Clear the intermediate results vector for each row
610  std::fill(intermediate.begin(), intermediate.end(), 0.0);
611 
612  ITYPE start = matrix_csr.row_ptrs_[i];
613  ITYPE end = matrix_csr.row_ptrs_[i + 1];
614 
615  // Iterate over the non-zero elements of the row
616  for (ITYPE j = start; j < end; ++j) {
617  // Get the column index and value of the element
618  ITYPE col = matrix_csr.col_indices_[j];
619  VTYPE value = matrix_csr.values_[j];
620 
621  ITYPE multiplicand_start = multiplicand.matrix_csr.row_ptrs_[col];
622  ITYPE multiplicand_end = multiplicand.matrix_csr.row_ptrs_[col + 1];
623 
624  // Multiply the element with the corresponding column of the other matrix
625  for (ITYPE k = multiplicand_start; k < multiplicand_end; ++k) {
626  ITYPE multiplicand_col = multiplicand.matrix_csr.col_indices_[k];
627  VTYPE multiplicand_value = multiplicand.matrix_csr.values_[k];
628  intermediate[multiplicand_col] += value * multiplicand_value;
629  }
630  }
631 
632  // Add the intermediate results to the result vectors
633  for (ITYPE j = 0; j < multiplicand.cols_; ++j) {
634  VTYPE result_value = intermediate[j];
635  if (std::abs(result_value) >= std::numeric_limits<VTYPE>::min()) {
636  res.matrix_csr.values_.emplace_back(result_value);
637  res.matrix_csr.col_indices_.emplace_back(j);
638  }
639  }
640  res.matrix_csr.row_ptrs_[i+1]=res.matrix_csr.values_.size();
641  }
642 
643  res.nnz_ = res.matrix_csr.values_.size();
644 
645  } else if (matrix_format_ == format::CSC) {
646 
647  // Perform element-wise multiplication of the CSC vectors
648  res.matrix_csc.values_.reserve((matrix_csc.values_.size() <= multiplicand.matrix_csc.values_.size()) ? multiplicand.matrix_csc.values_.size() : matrix_csc.values_.size());
649  res.matrix_csc.row_indices_.reserve((matrix_csc.row_indices_.size() <= multiplicand.matrix_csc.row_indices_.size()) ? multiplicand.matrix_csc.row_indices_.size() : matrix_csc.row_indices_.size());
650  res.matrix_csc.col_ptrs_.resize(cols_+1);
651 
652  // Initialize a vector to store the intermediate results
653  std::vector<VTYPE> intermediate(rows_, 0.0);
654 
655  res.matrix_csc.col_ptrs_[0] = 0;
656 
657  // Iterate over each column of the initial matrix
658  for (ITYPE j = 0; j < multiplicand.cols_; ++j) {
659  // Clear the intermediate results vector for each row
660  std::fill(intermediate.begin(), intermediate.end(), 0.0);
661 
662  ITYPE multiplicand_start = multiplicand.matrix_csc.col_ptrs_[j];
663  ITYPE multiplicand_end = multiplicand.matrix_csc.col_ptrs_[j + 1];
664 
665  // Iterate over the non-zero elements of the cloumn
666  for (ITYPE k = multiplicand_start; k < multiplicand_end; ++k) {
667  // Get the row index and value of the element
668  ITYPE multiplicand_row = multiplicand.matrix_csc.row_indices_[k];
669  VTYPE multiplicand_value = multiplicand.matrix_csc.values_[k];
670 
671  ITYPE start = matrix_csc.col_ptrs_[multiplicand_row];
672  ITYPE end = matrix_csc.col_ptrs_[multiplicand_row + 1];
673 
674  // Multiply the element with the corresponding column of the other matrix
675  for (ITYPE i = start; i < end; ++i) {
676  ITYPE row = matrix_csc.row_indices_[i];
677  VTYPE value = matrix_csc.values_[i];
678  intermediate[row] += value * multiplicand_value;
679  }
680  }
681 
682  // Add the intermediate results to the result vectors
683  for (ITYPE i = 0; i < multiplicand.rows_; ++i) {
684  VTYPE result_value = intermediate[i];
685  if (std::abs(result_value) >= std::numeric_limits<VTYPE>::min()) {
686  res.matrix_csc.values_.emplace_back(result_value);
687  res.matrix_csc.row_indices_.emplace_back(i);
688  }
689  }
690  res.matrix_csc.col_ptrs_[j+1]=res.matrix_csc.values_.size();
691  }
692 
693  res.nnz_ = res.matrix_csc.values_.size();
694 
695  } else {
696  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix operator*()" << std::endl;
697  std::cerr << " Please set the matrix_format_ as:" << std::endl;
698  std::cerr << " format::COO: COOrdinate format" << std::endl;
699  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
700  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
701  std::abort();
702  }
703 
704  return res;
705 }
706 
707 // Overload multiplication operator to perform scalar multiplication A*x
708 template <typename ITYPE, typename VTYPE>
709 template <typename STYPE>
710 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::operator*(const STYPE &scalar) const {
711  static_assert(std::is_convertible<STYPE, VTYPE>::value,
712  "MUI Error [matrix_arithmetic.h]: scalar type cannot be converted to matrix element type in scalar multiplication");
713 
714  // Create a new sparse matrix object for the result
715  sparse_matrix<ITYPE,VTYPE> res(*this);
716 
717  if (matrix_format_ == format::COO) {
718 
719  for (VTYPE &element : res.matrix_coo.values_) {
720  if (static_cast<VTYPE>(scalar) >= std::numeric_limits<VTYPE>::min())
721  element *= scalar;
722  }
723 
724  } else if (matrix_format_ == format::CSR) {
725 
726  for (VTYPE &element : res.matrix_csr.values_) {
727  if (static_cast<VTYPE>(scalar) >= std::numeric_limits<VTYPE>::min())
728  element *= scalar;
729  }
730 
731  } else if (matrix_format_ == format::CSC) {
732 
733  for (VTYPE &element : res.matrix_csc.values_) {
734  if (static_cast<VTYPE>(scalar) >= std::numeric_limits<VTYPE>::min())
735  element *= scalar;
736  }
737 
738  } else {
739  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix scalar operator*()" << std::endl;
740  std::cerr << " Please set the matrix_format_ as:" << std::endl;
741  std::cerr << " format::COO: COOrdinate format" << std::endl;
742  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
743  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
744  std::abort();
745  }
746 
747  return res;
748 
749 }
750 
751 // Overload multiplication operator to perform scalar multiplication x*A
752 template<typename ITYPE, typename VTYPE, typename STYPE>
753 sparse_matrix<ITYPE,VTYPE> operator*(const STYPE &scalar, const sparse_matrix<ITYPE,VTYPE> &exist_mat) {
754  return exist_mat * scalar;
755 }
756 
757 // Member function of dot product
758 template <typename ITYPE, typename VTYPE>
759 VTYPE sparse_matrix<ITYPE,VTYPE>::dot_product(sparse_matrix<ITYPE,VTYPE> &exist_mat) const {
760  assert(((cols_ == 1)&&(exist_mat.cols_ == 1)) &&
761  "MUI Error [matrix_arithmetic.h]: dot_product function only works for column vectors");
762  sparse_matrix<ITYPE,VTYPE> tempThis(*this);
763  sparse_matrix<ITYPE,VTYPE> thisT(tempThis.transpose());
764  sparse_matrix<ITYPE,VTYPE> tempMat(thisT * exist_mat);
765  assert(((tempMat.get_rows() == 1)&&(tempMat.get_cols() == 1)) &&
766  "MUI Error [matrix_arithmetic.h]: result of dot_product function should be a scalar");
767  return (tempMat.get_value(0,0));
768 }
769 
770 // Member function of Hadamard product
771 template <typename ITYPE, typename VTYPE>
772 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::hadamard_product(sparse_matrix<ITYPE,VTYPE> &exist_mat) {
773  if (rows_ != exist_mat.rows_ || cols_ != exist_mat.cols_) {
774  std::cerr << "MUI Error [matrix_arithmetic.h]: matrix size mismatch during matrix Hadamard product" << std::endl;
775  std::abort();
776  }
777 
778  if (exist_mat.matrix_format_ != matrix_format_) {
779  exist_mat.format_conversion(this->get_format(), true, true, "overwrite");
780  } else {
781  if (!exist_mat.is_sorted_unique("matrix_arithmetic.h", "hadamard_product()")){
782  if (exist_mat.matrix_format_ == format::COO) {
783  exist_mat.sort_coo(true, true, "overwrite");
784  } else if (exist_mat.matrix_format_ == format::CSR) {
785  exist_mat.sort_csr(true, "overwrite");
786  } else if (exist_mat.matrix_format_ == format::CSC) {
787  exist_mat.sort_csc(true, "overwrite");
788  } else {
789  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised exist_mat matrix format for matrix hadamard_product()" << std::endl;
790  std::cerr << " Please set the exist_mat matrix_format_ as:" << std::endl;
791  std::cerr << " format::COO: COOrdinate format" << std::endl;
792  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
793  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
794  std::abort();
795  }
796  }
797  }
798 
799  if (!this->is_sorted_unique("matrix_arithmetic.h", "hadamard_product()")){
800  if (matrix_format_ == format::COO) {
801  this->sort_coo(true, true, "overwrite");
802  } else if (matrix_format_ == format::CSR) {
803  this->sort_csr(true, "overwrite");
804  } else if (matrix_format_ == format::CSC) {
805  this->sort_csc(true, "overwrite");
806  } else {
807  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix hadamard_product()" << std::endl;
808  std::cerr << " Please set the matrix_format_ as:" << std::endl;
809  std::cerr << " format::COO: COOrdinate format" << std::endl;
810  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
811  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
812  std::abort();
813  }
814  }
815 
816  // Create a new sparse matrix object for the result
817  sparse_matrix<ITYPE,VTYPE> res(rows_, cols_, this->get_format());
818 
819  if (matrix_format_ == format::COO) {
820 
821  // Perform hadamard product of the COO vectors
822  res.matrix_coo.values_.reserve(matrix_coo.values_.size() + exist_mat.matrix_coo.values_.size());
823  res.matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size() + exist_mat.matrix_coo.row_indices_.size());
824  res.matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size() + exist_mat.matrix_coo.col_indices_.size());
825 
826  // Insert the COO vectors of the initial sparse matrix to the result sparse matrix
827  res.matrix_coo.values_ = std::vector<VTYPE>(matrix_coo.values_.begin(), matrix_coo.values_.end());
828  res.matrix_coo.row_indices_ = std::vector<ITYPE>(matrix_coo.row_indices_.begin(), matrix_coo.row_indices_.end());
829  res.matrix_coo.col_indices_ = std::vector<ITYPE>(matrix_coo.col_indices_.begin(), matrix_coo.col_indices_.end());
830 
831  // Append the exist_mat COO vectors to the result sparse matrix
832  res.matrix_coo.values_.insert(res.matrix_coo.values_.end(), exist_mat.matrix_coo.values_.begin(), exist_mat.matrix_coo.values_.end());
833  res.matrix_coo.row_indices_.insert(res.matrix_coo.row_indices_.end(), exist_mat.matrix_coo.row_indices_.begin(), exist_mat.matrix_coo.row_indices_.end());
834  res.matrix_coo.col_indices_.insert(res.matrix_coo.col_indices_.end(), exist_mat.matrix_coo.col_indices_.begin(), exist_mat.matrix_coo.col_indices_.end());
835 
836  // Sort and deduplicate the result
837  res.sort_coo(true, true, "multiply");
838 
839  } else if (matrix_format_ == format::CSR) {
840 
841  // Perform element-wise hadamard product of the CSR vectors
842  res.matrix_csr.values_.reserve(matrix_csr.values_.size() + exist_mat.matrix_csr.values_.size());
843  res.matrix_csr.row_ptrs_.resize(rows_ + 1);
844  res.matrix_csr.col_indices_.reserve(matrix_csr.col_indices_.size() + exist_mat.matrix_csr.col_indices_.size());
845 
846  res.matrix_csr.row_ptrs_[0] = 0;
847 
848  ITYPE row = 0;
849  while (row < rows_) {
850  ITYPE start = matrix_csr.row_ptrs_[row];
851  ITYPE end = matrix_csr.row_ptrs_[row + 1];
852 
853  ITYPE exist_mat_start = exist_mat.matrix_csr.row_ptrs_[row];
854  ITYPE exist_mat_end = exist_mat.matrix_csr.row_ptrs_[row + 1];
855 
856  // Merge the values and column indices of the two rows
857  ITYPE i = start;
858  ITYPE j = exist_mat_start;
859  while (i < end && j < exist_mat_end) {
860  ITYPE col = matrix_csr.col_indices_[i];
861  ITYPE exist_mat_col = exist_mat.matrix_csr.col_indices_[j];
862 
863  if (col == exist_mat_col) {
864  // Add the corresponding values if the columns match
865  if (std::abs(matrix_csr.values_[i] * exist_mat.matrix_csr.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
866  res.matrix_csr.values_.emplace_back(matrix_csr.values_[i] * exist_mat.matrix_csr.values_[j]);
867  res.matrix_csr.col_indices_.emplace_back(col);
868  }
869  i++;
870  j++;
871  } else if (col < exist_mat_col) {
872  i++;
873  } else {
874  j++;
875  }
876  }
877 
878  // Update the row pointer
879  res.nnz_ = res.matrix_csr.col_indices_.size();
880  res.matrix_csr.row_ptrs_[row + 1] = res.nnz_;
881 
882  row++;
883  }
884 
885  } else if (matrix_format_ == format::CSC) {
886 
887  // Perform element-wise hadamard product of the CSC vectors
888  res.matrix_csc.values_.reserve(matrix_csc.values_.size() + exist_mat.matrix_csc.values_.size());
889  res.matrix_csc.row_indices_.reserve(matrix_csc.row_indices_.size() + exist_mat.matrix_csc.row_indices_.size());
890  res.matrix_csc.col_ptrs_.resize(cols_ + 1);
891 
892  res.matrix_csc.col_ptrs_[0] = 0;
893 
894  ITYPE column = 0;
895  while (column < cols_) {
896  ITYPE start = matrix_csc.col_ptrs_[column];
897  ITYPE end = matrix_csc.col_ptrs_[column + 1];
898 
899  ITYPE exist_mat_start = exist_mat.matrix_csc.col_ptrs_[column];
900  ITYPE exist_mat_end = exist_mat.matrix_csc.col_ptrs_[column + 1];
901 
902  // Merge the values and row indices of the two columns
903  ITYPE i = start;
904  ITYPE j = exist_mat_start;
905  while (i < end && j < exist_mat_end) {
906  ITYPE row = matrix_csc.row_indices_[i];
907  ITYPE exist_mat_row = exist_mat.matrix_csc.row_indices_[j];
908 
909  if ((row == exist_mat_row) && std::abs(matrix_csc.values_[i] * exist_mat.matrix_csc.values_[j]) >= std::numeric_limits<VTYPE>::min()) {
910  // Add the corresponding values if the columns match
911  res.matrix_csc.values_.emplace_back(matrix_csc.values_[i] * exist_mat.matrix_csc.values_[j]);
912  res.matrix_csc.row_indices_.emplace_back(row);
913  i++;
914  j++;
915  } else if (row < exist_mat_row) {
916  i++;
917  } else {
918  j++;
919  }
920  }
921 
922  // Update the column pointer
923  res.nnz_ = res.matrix_csc.row_indices_.size();
924  res.matrix_csc.col_ptrs_[column + 1] = res.nnz_;
925 
926  column++;
927  }
928 
929  } else {
930  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix hadamard_product()" << std::endl;
931  std::cerr << " Please set the matrix_format_ as:" << std::endl;
932  std::cerr << " format::COO: COOrdinate format" << std::endl;
933  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
934  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
935  std::abort();
936  }
937 
938  return res;
939 }
940 
941 // Member function to get transpose of matrix
942 template <typename ITYPE, typename VTYPE>
943 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::transpose(bool performSortAndUniqueCheck) const {
944 
945  sparse_matrix<ITYPE,VTYPE> res(*this);
946 
947  if (matrix_format_ == format::COO) {
948 
949  if (performSortAndUniqueCheck){
950  if (!res.is_sorted_unique("matrix_arithmetic.h", "transpose()")){
951  res.sort_coo(true, true, "overwrite");
952  }
953  }
954 
955  res.index_reinterpretation();
956 
957  } else if (matrix_format_ == format::CSR) {
958 
959  res.format_conversion("CSC", performSortAndUniqueCheck, performSortAndUniqueCheck, "overwrite");
960 
961  res.format_reinterpretation();
962 
963  } else if (matrix_format_ == format::CSC) {
964 
965  res.format_conversion("CSR", performSortAndUniqueCheck, performSortAndUniqueCheck, "overwrite");
966 
967  res.format_reinterpretation();
968 
969  } else {
970  std::cerr << "MUI Error [matrix_arithmetic.h]: Unrecognised matrix format for matrix transpose()" << std::endl;
971  std::cerr << " Please set the matrix_format_ as:" << std::endl;
972  std::cerr << " format::COO: COOrdinate format" << std::endl;
973  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
974  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
975  std::abort();
976  }
977 
978  return res;
979 
980 }
981 
982 // Member function to perform LU decomposition
983 template <typename ITYPE, typename VTYPE>
984 void sparse_matrix<ITYPE,VTYPE>::lu_decomposition(sparse_matrix<ITYPE,VTYPE> &L, sparse_matrix<ITYPE,VTYPE> &U) const {
985  if (((L.get_rows() != 0) && (L.get_rows() != rows_)) ||
986  ((U.get_rows() != 0) && (U.get_rows() != rows_)) ||
987  ((L.get_cols() != 0) && (L.get_cols() != cols_)) ||
988  ((U.get_cols() != 0) && (U.get_cols() != cols_))) {
989  std::cerr << "MUI Error [matrix_arithmetic.h]: L & U Matrices must be null or same size of initial matrix in LU decomposition" << std::endl;
990  std::abort();
991  }
992 
993  if ((!L.empty()) || (!U.empty())) {
994  std::cerr << "MUI Error [matrix_arithmetic.h]: L & U Matrices must be empty in LU decomposition" << std::endl;
995  std::abort();
996  }
997 
998  if (rows_ != cols_) {
999  std::cerr << "MUI Error [matrix_arithmetic.h]: Only square matrix can perform LU decomposition" << std::endl;
1000  std::abort();
1001  }
1002 
1003  if ((L.get_rows() != rows_) || (L.get_cols() != cols_)) {
1004  // Resize the lower triangular matrix
1005  L.resize(rows_, cols_);
1006  }
1007 
1008  if ((U.get_rows() != rows_) || (U.get_cols() != cols_)) {
1009  // Resize the upper triangular matrix
1010  U.resize(rows_, cols_);
1011  }
1012 
1013  ITYPE n = rows_;
1014  for (ITYPE i = 0; i < rows_; ++i) {
1015  // Calculate the upper triangular matrix
1016  for (ITYPE k = i; k < cols_; ++k) {
1017  VTYPE sum = 0.0;
1018  for (ITYPE j = 0; j < i; ++j) {
1019  sum += L.get_value(i, j) * U.get_value(j, k);
1020  }
1021  U.set_value(i, k, (this->get_value(i, k) - sum));
1022  }
1023 
1024  // Calculate the lower triangular matrix
1025  for (ITYPE k = i; k < rows_; k++) {
1026  if (i == k) {
1027  L.set_value(i, i, static_cast<VTYPE>(1.0));
1028  } else {
1029  VTYPE sum = 0.0;
1030  for (ITYPE j = 0; j < i; ++j) {
1031  sum += L.get_value(k, j) * U.get_value(j, i);
1032  }
1033  assert(std::abs(U.get_value(i, i)) >= std::numeric_limits<VTYPE>::min() &&
1034  "MUI Error [matrix_arithmetic.h]: Divide by zero assert for U.get_value(i, i)");
1035  L.set_value(k, i, (this->get_value(k, i) - sum) / U.get_value(i, i));
1036  }
1037  }
1038  }
1039 }
1040 
1041 // Member function to perform QR decomposition
1042 template <typename ITYPE, typename VTYPE>
1043 void sparse_matrix<ITYPE,VTYPE>::qr_decomposition(sparse_matrix<ITYPE,VTYPE> &Q, sparse_matrix<ITYPE,VTYPE> &R) const {
1044  if (((Q.get_rows() != 0) && (Q.get_rows() != rows_)) ||
1045  ((R.get_rows() != 0) && (R.get_rows() != rows_)) ||
1046  ((Q.get_cols() != 0) && (Q.get_cols() != cols_)) ||
1047  ((R.get_cols() != 0) && (R.get_cols() != cols_))) {
1048  std::cerr << "MUI Error [matrix_arithmetic.h]: Q & R Matrices must be null in QR decomposition" << std::endl;
1049  std::abort();
1050  }
1051  if ((!Q.empty()) || (!R.empty())) {
1052  std::cerr << "MUI Error [matrix_arithmetic.h]: Q & R Matrices must be empty in QR decomposition" << std::endl;
1053  std::abort();
1054  }
1055  assert((rows_ >= cols_) &&
1056  "MUI Error [matrix_arithmetic.h]: number of rows of matrix should larger or equals to number of columns in QR decomposition");
1057 
1058  if ((Q.get_rows() != rows_) || (Q.get_cols() != cols_)) {
1059  // Resize the orthogonal matrix
1060  Q.resize(rows_, cols_);
1061  }
1062  if ((R.get_rows() != rows_) || (R.get_cols() != cols_)) {
1063  // Resize the upper triangular matrix
1064  R.resize(rows_, cols_);
1065  }
1066 
1067  // Get a copy of the matrix
1068  sparse_matrix<ITYPE,VTYPE> mat_copy (*this);
1069  // Diagonal elements
1070  std::vector<VTYPE> r_diag (cols_);
1071 
1072  // Calculate the diagonal element values
1073  for (ITYPE c = 0; c <cols_; ++c) {
1074  VTYPE nrm (0.0);
1075 
1076  // Compute 2-norm of k-th column without under/overflow.
1077  for (ITYPE r = c; r < rows_; ++r)
1078  nrm = std::sqrt((nrm * nrm) + (mat_copy.get_value(r, c) * mat_copy.get_value(r, c)));
1079 
1080  if (nrm != static_cast<VTYPE>(0.0)) {
1081 
1082  // Form k-th Householder vector.
1083  if (mat_copy.get_value(c, c) < static_cast<VTYPE>(0.0))
1084  nrm = -nrm;
1085 
1086  for (ITYPE r = c; r < rows_; ++r)
1087  mat_copy.set_value(r, c, (mat_copy.get_value(r, c)/nrm));
1088 
1089  mat_copy.set_value(c, c, (mat_copy.get_value(c, c) + static_cast<VTYPE>(1.0)));
1090 
1091  // Apply transformation to remaining columns.
1092  for (ITYPE j = c + 1; j < cols_; ++j) {
1093  VTYPE s = 0.0;
1094 
1095  for (ITYPE r = c; r < rows_; ++r)
1096  s += mat_copy.get_value(r, c) * mat_copy.get_value(r, j);
1097 
1098  s /= -mat_copy.get_value(c, c);
1099  for (ITYPE r = c; r < rows_; ++r)
1100  mat_copy.set_value(r, j, (mat_copy.get_value(r, j) + s * mat_copy.get_value(r, c)));
1101  }
1102  }
1103  r_diag[c] = -nrm;
1104  }
1105 
1106  // Calculate the orthogonal matrix
1107  for (ITYPE c = cols_ - 1; c >= 0; --c) {
1108  Q.set_value(c, c, static_cast<VTYPE>(1.0));
1109 
1110  for (ITYPE cc = c; cc < cols_; ++cc)
1111  if (mat_copy.get_value(c, c) != static_cast<VTYPE>(0.0)) {
1112  VTYPE s=0.0;
1113 
1114  for (ITYPE r = c; r < rows_; ++r)
1115  s += mat_copy.get_value(r, c) * Q.get_value(r, cc);
1116 
1117  s /= -mat_copy.get_value(c, c);
1118  for (ITYPE r = c; r < rows_; ++r)
1119  Q.set_value(r, cc, (Q.get_value(r, cc) + s * mat_copy.get_value(r, c)));
1120  }
1121  }
1122 
1123  // Calculate the upper triangular matrix
1124  for (ITYPE c = 0; c < cols_; ++c)
1125  for (ITYPE r = 0; r < rows_; ++r)
1126  if (c < r)
1127  R.set_value(c, r, mat_copy.get_value(c, r));
1128  else if (c == r)
1129  R.set_value(c, r, r_diag[c]);
1130 }
1131 
1132 // Member function to get the inverse of matrix by using Gaussian elimination
1133 template <typename ITYPE, typename VTYPE>
1134 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::inverse() const {
1135  if (rows_ != cols_) {
1136  std::cerr << "MUI Error [matrix_arithmetic.h]: Matrix must be square to find its inverse" << std::endl;
1137  std::abort();
1138  }
1139 
1140  sparse_matrix<ITYPE,VTYPE> mat_copy (*this);
1141  sparse_matrix<ITYPE,VTYPE> inverse_mat (rows_,"identity", this->get_format());
1142 
1143  for (ITYPE r = 0; r < rows_; ++r) {
1144 
1145  ITYPE max_row = r;
1146  VTYPE max_value= static_cast<VTYPE>(-1.0);
1147 
1148  // Partial pivoting for Gaussian elimination
1149  ITYPE ppivot;
1150  for (ITYPE rb = r; rb < rows_; ++rb) {
1151  const VTYPE tmp = std::abs(mat_copy.get_value(rb, r));
1152 
1153  if ((tmp > max_value) && (std::abs(tmp) >= std::numeric_limits<VTYPE>::min())) {
1154  max_value = tmp;
1155  max_row = rb;
1156  }
1157  }
1158 
1159  assert(std::abs(mat_copy.get_value(max_row, r)) >= std::numeric_limits<VTYPE>::min() &&
1160  "MUI Error [matrix_arithmetic.h]: Divide by zero assert for mat_copy.get_value(max_row, r). Cannot perform matrix invert due to singular matrix.");
1161 
1162  if (max_row != r) {
1163  for (ITYPE c = 0; c < cols_; ++c)
1164  mat_copy.swap_elements(r, c, max_row, c);
1165  ppivot = max_row;
1166  } else {
1167  ppivot = 0;
1168  }
1169 
1170  const ITYPE indx = ppivot;
1171 
1172  if (indx != 0)
1173  for (ITYPE c = 0; c < cols_; ++c)
1174  inverse_mat.swap_elements(r, c, indx, c);
1175 
1176  const VTYPE diag = mat_copy.get_value(r, r);
1177 
1178  for (ITYPE c = 0; c < cols_; ++c) {
1179  mat_copy.set_value(r, c, (mat_copy.get_value(r, c) / diag));
1180  inverse_mat.set_value(r, c, (inverse_mat.get_value(r, c) / diag));
1181  }
1182 
1183  for (ITYPE rr = 0; rr < rows_; ++rr)
1184  if (rr != r) {
1185  const VTYPE off_diag = mat_copy.get_value(rr, r);
1186 
1187  for (ITYPE c = 0; c < cols_; ++c) {
1188  mat_copy.set_value(rr, c, (mat_copy.get_value(rr, c) - off_diag * mat_copy.get_value(r, c)));
1189  inverse_mat.set_value(rr, c, (inverse_mat.get_value(rr, c) - off_diag * inverse_mat.get_value(r, c)));
1190  }
1191  }
1192  }
1193  return inverse_mat;
1194 }
1195 
1196 // **************************************************
1197 // ********** Protected member functions ************
1198 // **************************************************
1199 
1200 // Protected member function to reinterpret the row and column indexes for sparse matrix with COO format - helper function on matrix transpose
1201 template<typename ITYPE, typename VTYPE>
1202 void sparse_matrix<ITYPE,VTYPE>::index_reinterpretation() {
1203  assert((matrix_format_ == format::COO) &&
1204  "MUI Error [matrix_arithmetic.h]: index_reinterpretation() is for COO format only.");
1205 
1206  std::swap(matrix_coo.row_indices_, matrix_coo.col_indices_);
1207  ITYPE temp_index = rows_;
1208  rows_ = cols_;
1209  cols_ = temp_index;
1210 
1211 }
1212 
1213 // Protected member function to reinterpret the format of sparse matrix between CSR format and CSC format - helper function on matrix transpose
1214 template<typename ITYPE, typename VTYPE>
1215 void sparse_matrix<ITYPE,VTYPE>::format_reinterpretation() {
1216  assert(((matrix_format_ == format::CSR) || (matrix_format_ == format::CSC)) &&
1217  "MUI Error [matrix_arithmetic.h]: format_reinterpretation() is for CSR or CSC format.");
1218 
1219  if (matrix_format_ == format::CSR) {
1220 
1221  matrix_csc.col_ptrs_.swap(matrix_csr.row_ptrs_);
1222  matrix_csc.row_indices_.swap(matrix_csr.col_indices_);
1223  matrix_csc.values_.swap(matrix_csr.values_);
1224 
1225  ITYPE temp_index = rows_;
1226  rows_ = cols_;
1227  cols_ = temp_index;
1228 
1229  matrix_format_ = format::CSC;
1230 
1231  matrix_csr.row_ptrs_.clear();
1232  matrix_csr.col_indices_.clear();
1233  matrix_csr.values_.clear();
1234 
1235  } else if (matrix_format_ == format::CSC) {
1236 
1237  matrix_csr.row_ptrs_.swap(matrix_csc.col_ptrs_);
1238  matrix_csr.col_indices_.swap(matrix_csc.row_indices_);
1239  matrix_csr.values_.swap(matrix_csc.values_);
1240 
1241  ITYPE temp_index = rows_;
1242  rows_ = cols_;
1243  cols_ = temp_index;
1244 
1245  matrix_format_ = format::CSR;
1246 
1247  matrix_csc.col_ptrs_.clear();
1248  matrix_csc.row_indices_.clear();
1249  matrix_csc.values_.clear();
1250 
1251  }
1252 
1253 }
1254 
1255 } // linalg
1256 } // mui
1257 
1258 #endif /* MUI_MATRIX_ARITHMETIC_H_ */
mui::linalg::sparse_matrix
Definition: matrix.h:61
mui::linalg::sparse_matrix::index_reinterpretation
void index_reinterpretation()
Definition: matrix_arithmetic.h:1202
mui::linalg::sparse_matrix::get_rows
ITYPE get_rows() const
Definition: matrix_io_info.h:579
mui::linalg::sparse_matrix::is_sorted_unique
bool is_sorted_unique(const std::string &={}, const std::string &={}) const
Definition: matrix_io_info.h:701
mui::linalg::sparse_matrix::operator*
sparse_matrix< ITYPE, VTYPE > operator*(sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_arithmetic.h:521
mui::linalg::sparse_matrix::dot_product
VTYPE dot_product(sparse_matrix< ITYPE, VTYPE > &) const
Definition: matrix_arithmetic.h:759
mui::linalg::sparse_matrix::sort_csr
void sort_csr(bool=false, const std::string &="overwrite")
Definition: matrix_manipulation.h:1037
mui::linalg::sparse_matrix::transpose
sparse_matrix< ITYPE, VTYPE > transpose(bool=true) const
Definition: matrix_arithmetic.h:943
mui::linalg::sparse_matrix::qr_decomposition
void qr_decomposition(sparse_matrix< ITYPE, VTYPE > &, sparse_matrix< ITYPE, VTYPE > &) const
Definition: matrix_arithmetic.h:1043
mui::linalg::sparse_matrix::set_value
void set_value(ITYPE, ITYPE, VTYPE, bool=true)
Definition: matrix_manipulation.h:292
mui::linalg::sparse_matrix::operator+
sparse_matrix< ITYPE, VTYPE > operator+(sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_arithmetic.h:62
mui::linalg::sparse_matrix::get_value
VTYPE get_value(ITYPE, ITYPE) const
Definition: matrix_io_info.h:523
mui::linalg::sparse_matrix::hadamard_product
sparse_matrix< ITYPE, VTYPE > hadamard_product(sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_arithmetic.h:772
mui::linalg::sparse_matrix::resize
void resize(ITYPE, ITYPE)
Definition: matrix_manipulation.h:62
mui::linalg::sparse_matrix::get_cols
ITYPE get_cols() const
Definition: matrix_io_info.h:585
mui::linalg::sparse_matrix::inverse
sparse_matrix< ITYPE, VTYPE > inverse() const
Definition: matrix_arithmetic.h:1134
mui::linalg::sparse_matrix::swap_elements
void swap_elements(ITYPE, ITYPE, ITYPE, ITYPE)
Definition: matrix_manipulation.h:407
mui::linalg::sparse_matrix::sort_coo
void sort_coo(bool=true, bool=false, const std::string &="overwrite")
Definition: matrix_manipulation.h:808
mui::linalg::sparse_matrix::format_conversion
void format_conversion(const std::string &="COO", bool=true, bool=false, const std::string &="overwrite")
Definition: matrix_manipulation.h:639
mui::linalg::sparse_matrix::format_reinterpretation
void format_reinterpretation()
Definition: matrix_arithmetic.h:1215
mui::linalg::sparse_matrix::operator-
sparse_matrix< ITYPE, VTYPE > operator-(sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_arithmetic.h:289
mui::linalg::sparse_matrix::empty
bool empty() const
Definition: matrix_io_info.h:665
mui::linalg::sparse_matrix::lu_decomposition
void lu_decomposition(sparse_matrix< ITYPE, VTYPE > &, sparse_matrix< ITYPE, VTYPE > &) const
Definition: matrix_arithmetic.h:984
mui::linalg::sparse_matrix::sort_csc
void sort_csc(bool=false, const std::string &="overwrite")
Definition: matrix_manipulation.h:1168
mui::dim::min
u u u u u u min
Definition: dim.h:289
mui::linalg::operator*
sparse_matrix< ITYPE, VTYPE > operator*(const STYPE &scalar, const sparse_matrix< ITYPE, VTYPE > &exist_mat)
Definition: matrix_arithmetic.h:753
mui
Definition: comm.h:54
mui::sum
SCALAR sum(vexpr< E, SCALAR, D > const &u)
Definition: point.h:362
mui::swap
void swap(storage< Args... > &lhs, storage< Args... > &rhs)
Definition: dynstorage.h:234