Multiscale Universal Interface  2.0
A Concurrent Framework for Coupling Heterogeneous Solvers
matrix_manipulation.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_MANIPULATION_H_
48 #define MUI_MATRIX_MANIPULATION_H_
49 
50 #include <cassert>
51 #include <numeric>
52 
53 namespace mui {
54 namespace linalg {
55 
56 // **************************************************
57 // ************ Public member functions *************
58 // **************************************************
59 
60 // Member function to resize an all-zero or null matrix
61 template<typename ITYPE, typename VTYPE>
62 void sparse_matrix<ITYPE,VTYPE>::resize(ITYPE r, ITYPE c) {
63  assert(((this->non_zero_elements_count()) == 0) &&
64  "MUI Error [matrix_manipulation.h]: resize function only works for all-zero matrix");
65  rows_ = r;
66  cols_ = c;
67 
68  if (matrix_format_ == format::CSR) {
69  matrix_csr.row_ptrs_.clear();
70  matrix_csr.row_ptrs_.resize((r+1), 0);
71  } else if (matrix_format_ == format::CSC) {
72  matrix_csc.col_ptrs_.clear();
73  matrix_csc.col_ptrs_.resize((c+1), 0);
74  }
75 }
76 
77 // Member function to copy a sparse_matrix
78 template<typename ITYPE, typename VTYPE>
79 void sparse_matrix<ITYPE,VTYPE>::copy(const sparse_matrix<ITYPE,VTYPE> &exist_mat) {
80 
81  // Copy the data from the existing matrix
82  assert(this->empty() &&
83  "MUI Error [matrix_manipulation.h]: copy function only works for empty (all zero elements) matrix");
84  assert((((rows_ == exist_mat.rows_) && (cols_ == exist_mat.cols_)) || ((rows_ == 0) && (cols_ == 0))) &&
85  "MUI Error [matrix_manipulation.h]: matrix size mismatch in copy function ");
86 
87  exist_mat.assert_valid_vector_size("matrix_manipulation.h", "copy()");
88 
89  std::string format_store = this->get_format();
90 
91  if (exist_mat.nnz_>0) {
92 
93  if ((rows_ == 0) && (cols_ == 0))
94  this->resize(exist_mat.rows_, exist_mat.cols_);
95 
96  nnz_ = exist_mat.nnz_;
97 
98  if (exist_mat.matrix_format_ == format::COO) {
99 
100  matrix_coo.values_.reserve(exist_mat.matrix_coo.values_.size());
101  matrix_coo.row_indices_.reserve(exist_mat.matrix_coo.row_indices_.size());
102  matrix_coo.col_indices_.reserve(exist_mat.matrix_coo.col_indices_.size());
103 
104  matrix_coo.values_ = std::vector<VTYPE>(exist_mat.matrix_coo.values_.begin(), exist_mat.matrix_coo.values_.end());
105  matrix_coo.row_indices_ = std::vector<ITYPE>(exist_mat.matrix_coo.row_indices_.begin(), exist_mat.matrix_coo.row_indices_.end());
106  matrix_coo.col_indices_ = std::vector<ITYPE>(exist_mat.matrix_coo.col_indices_.begin(), exist_mat.matrix_coo.col_indices_.end());
107 
108  matrix_format_ = exist_mat.matrix_format_;
109 
110  } else if (exist_mat.matrix_format_ == format::CSR) {
111 
112  matrix_csr.values_.reserve(exist_mat.matrix_csr.values_.size());
113  matrix_csr.row_ptrs_.reserve(exist_mat.matrix_csr.row_ptrs_.size());
114  matrix_csr.col_indices_.reserve(exist_mat.matrix_csr.col_indices_.size());
115 
116  matrix_csr.values_ = std::vector<VTYPE>(exist_mat.matrix_csr.values_.begin(), exist_mat.matrix_csr.values_.end());
117  matrix_csr.row_ptrs_ = std::vector<ITYPE>(exist_mat.matrix_csr.row_ptrs_.begin(), exist_mat.matrix_csr.row_ptrs_.end());
118  matrix_csr.col_indices_ = std::vector<ITYPE>(exist_mat.matrix_csr.col_indices_.begin(), exist_mat.matrix_csr.col_indices_.end());
119 
120  matrix_format_ = exist_mat.matrix_format_;
121 
122  } else if (exist_mat.matrix_format_ == format::CSC) {
123 
124  matrix_csc.values_.reserve(exist_mat.matrix_csc.values_.size());
125  matrix_csc.row_indices_.reserve(exist_mat.matrix_csc.row_indices_.size());
126  matrix_csc.col_ptrs_.reserve(exist_mat.matrix_csc.col_ptrs_.size());
127 
128  matrix_csc.values_ = std::vector<VTYPE>(exist_mat.matrix_csc.values_.begin(), exist_mat.matrix_csc.values_.end());
129  matrix_csc.row_indices_ = std::vector<ITYPE>(exist_mat.matrix_csc.row_indices_.begin(), exist_mat.matrix_csc.row_indices_.end());
130  matrix_csc.col_ptrs_ = std::vector<ITYPE>(exist_mat.matrix_csc.col_ptrs_.begin(), exist_mat.matrix_csc.col_ptrs_.end());
131 
132  matrix_format_ = exist_mat.matrix_format_;
133 
134  } else {
135  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised exist_mat format" << std::endl;
136  std::cerr << " Please set the matrix_format_ as:" << std::endl;
137  std::cerr << " format::COO: COOrdinate format" << std::endl;
138  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
139  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
140  std::abort();
141  }
142  }
143 
144  if (this->get_format() != format_store)
145  this->format_conversion(format_store, true, true, "overwrite");
146 
147  this->assert_valid_vector_size("matrix_manipulation.h", "copy()");
148 
149 }
150 
151 // Member function to get a segment of a sparse_matrix
152 template<typename ITYPE, typename VTYPE>
153 sparse_matrix<ITYPE,VTYPE> sparse_matrix<ITYPE,VTYPE>::segment(ITYPE r_start, ITYPE r_end, ITYPE c_start, ITYPE c_end, bool performSortAndUniqueCheck) {
154  // get segment data from the existing matrix
155  assert((r_end >= r_start) &&
156  "MUI Error [matrix_manipulation.h]: segment function r_end has to be larger or equals to r_start");
157  assert((c_end >= c_start) &&
158  "MUI Error [matrix_manipulation.h]: segment function c_end has to be larger or equals to c_start");
159  assert(((r_end < rows_) && (r_start >= 0) && (c_end < cols_) && (c_start >= 0)) &&
160  "MUI Error [matrix_manipulation.h]: Matrix index out of range in segment function");
161 
162  if (performSortAndUniqueCheck) {
163  if (!(this->is_sorted_unique("matrix_manipulation.h", "segment()"))){
164  if(matrix_format_ == format::COO) {
165  this->sort_coo(true, true, "overwrite");
166  } else if (matrix_format_ == format::CSR) {
167  this->sort_csr(true, "overwrite");
168  } else if (matrix_format_ == format::CSC) {
169  this->sort_csc(true, "overwrite");
170  } else {
171  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
172  std::cerr << " Please set the matrix_format_ as:" << std::endl;
173  std::cerr << " format::COO: COOrdinate format" << std::endl;
174  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
175  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
176  std::abort();
177  }
178  }
179  }
180 
181  sparse_matrix<ITYPE,VTYPE> res((r_end-r_start+1), (c_end-c_start+1), this->get_format());
182 
183  if(matrix_format_ == format::COO) {
184 
185  // Iterate over the existing non-zero elements
186  for (ITYPE i = 0; i < nnz_; ++i) {
187  ITYPE row = matrix_coo.row_indices_[i];
188  ITYPE col = matrix_coo.col_indices_[i];
189 
190  // Check if the current element is within the specified ranges
191  if (row >= r_start && row <= r_end && col >= c_start && col <= c_end) {
192  // Calculate the indices for the segment
193  ITYPE subRow = row - r_start;
194  ITYPE subCol = col - c_start;
195 
196  res.matrix_coo.values_.reserve(res.matrix_coo.values_.size()+1);
197  res.matrix_coo.row_indices_.reserve(res.matrix_coo.row_indices_.size()+1);
198  res.matrix_coo.col_indices_.reserve(res.matrix_coo.col_indices_.size()+1);
199 
200  // Add the element to the segment matrix_coo struct
201  res.matrix_coo.values_.emplace_back(matrix_coo.values_[i]);
202  res.matrix_coo.row_indices_.emplace_back(subRow);
203  res.matrix_coo.col_indices_.emplace_back(subCol);
204  res.nnz_++;
205 
206  }
207  }
208 
209  } else if (matrix_format_ == format::CSR) {
210 
211  res.matrix_csr.row_ptrs_.reserve(r_end-r_start+2);
212 
213  // Iterate over the row pointers and column indices of the existing non-zero elements
214  for (ITYPE row = r_start; row <= r_end; ++row) {
215  // Get the starting and ending indices for the current row
216  ITYPE start = matrix_csr.row_ptrs_[row];
217  ITYPE end = matrix_csr.row_ptrs_[row + 1];
218 
219  // Iterate over the non-zero elements in the current row
220  for (ITYPE j = start; j < end; ++j) {
221  // Get the column index of the current element
222  ITYPE col = matrix_csr.col_indices_[j];
223 
224  // Check if the current element is within the specified column range
225  if (col >= c_start && col <= c_end) {
226  // Calculate the indices for the sub-block
227  ITYPE subCol = col - c_start;
228 
229  res.matrix_csr.values_.reserve(res.matrix_csr.values_.size()+1);
230  res.matrix_csr.col_indices_.reserve(res.matrix_csr.col_indices_.size()+1);
231 
232  // Add the element to the segment matrix_csr struct
233  res.matrix_csr.values_.emplace_back(matrix_csr.values_[j]);
234  res.matrix_csr.col_indices_.emplace_back(subCol);
235  res.nnz_++;
236  }
237  }
238  // Update the row pointer for the segment
239  res.matrix_csr.row_ptrs_[row - r_start + 1] = res.matrix_csr.col_indices_.size();
240  }
241 
242  } else if (matrix_format_ == format::CSC) {
243 
244  res.matrix_csc.col_ptrs_.reserve(c_end-c_start+2);
245 
246  // Iterate over the column pointers and row indices of the existing non-zero elements
247  for (ITYPE col = c_start; col <= c_end; ++col) {
248  // Get the starting and ending indices for the current column
249  ITYPE start = matrix_csc.col_ptrs_[col];
250  ITYPE end = matrix_csc.col_ptrs_[col + 1];
251 
252  // Iterate over the non-zero elements in the current column
253  for (ITYPE j = start; j < end; ++j) {
254  // Get the row index of the current element
255  ITYPE row = matrix_csc.row_indices_[j];
256 
257  // Check if the current element is within the specified row range
258  if (row >= r_start && row <= r_end) {
259  // Calculate the indices for the sub-block
260  ITYPE subRow = row - r_start;
261 
262  res.matrix_csc.values_.reserve(res.matrix_csc.values_.size()+1);
263  res.matrix_csc.row_indices_.reserve(res.matrix_csc.row_indices_.size()+1);
264 
265  // Add the element to the segment matrix_csc struct
266  res.matrix_csc.values_.emplace_back(matrix_csc.values_[j]);
267  res.matrix_csc.row_indices_.emplace_back(subRow);
268  res.nnz_++;
269  }
270  }
271 
272  // Update the column pointer for the segment
273  res.matrix_csc.col_ptrs_[col - c_start + 1] = res.matrix_csc.row_indices_.size();
274  }
275 
276  } else {
277  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
278  std::cerr << " Please set the matrix_format_ as:" << std::endl;
279  std::cerr << " format::COO: COOrdinate format" << std::endl;
280  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
281  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
282  std::abort();
283  }
284 
285  res.assert_valid_vector_size("matrix_manipulation.h", "segment()");
286 
287  return res;
288 }
289 
290 // Member function to insert an element
291 template<typename ITYPE, typename VTYPE>
292 void sparse_matrix<ITYPE,VTYPE>::set_value(ITYPE r, ITYPE c, VTYPE val, bool performSortAndUniqueCheck) {
293  assert(((r < rows_) && (r >= 0) && (c < cols_) && (c >= 0)) &&
294  "MUI Error [matrix_manipulation.h]: Matrix index out of range in set_value function");
295 
296  if (performSortAndUniqueCheck) {
297  if (!(this->is_sorted_unique("matrix_manipulation.h", "set_value()"))){
298  if(matrix_format_ == format::COO) {
299  this->sort_coo(true, true, "overwrite");
300  } else if (matrix_format_ == format::CSR) {
301  this->sort_csr(true, "overwrite");
302  } else if (matrix_format_ == format::CSC) {
303  this->sort_csc(true, "overwrite");
304  } else {
305  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
306  std::cerr << " Please set the matrix_format_ as:" << std::endl;
307  std::cerr << " format::COO: COOrdinate format" << std::endl;
308  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
309  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
310  std::abort();
311  }
312  }
313  }
314 
315  if(matrix_format_ == format::COO) {
316 
317  if (performSortAndUniqueCheck) {
318  this->coo_element_operation(r, c, val, "overwrite", "matrix_manipulation.h", "set_value()");
319  } else {
320  this->coo_element_operation(r, c, val, "nonsort", "matrix_manipulation.h", "set_value()");
321  }
322 
323  } else if (matrix_format_ == format::CSR) {
324 
325  this->csr_element_operation(r, c, val, "overwrite", "matrix_manipulation.h", "set_value()");
326 
327  } else if (matrix_format_ == format::CSC) {
328 
329  this->csc_element_operation(r, c, val, "overwrite", "matrix_manipulation.h", "set_value()");
330 
331  } else {
332  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
333  std::cerr << " Please set the matrix_format_ as:" << std::endl;
334  std::cerr << " format::COO: COOrdinate format" << std::endl;
335  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
336  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
337  std::abort();
338  }
339 }
340 
341 // Member function to insert the same value to all elements
342 template<typename ITYPE, typename VTYPE>
343 void sparse_matrix<ITYPE,VTYPE>::set_value(VTYPE val) {
344 
345  this->set_zero();
346 
347  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
348 
349  if(matrix_format_ == format::COO) {
350  // Resize the vectors to hold all possible elements
351  matrix_coo.values_.reserve(rows_*cols_);
352  matrix_coo.row_indices_.reserve(rows_*cols_);
353  matrix_coo.col_indices_.reserve(rows_*cols_);
354  matrix_coo.values_.resize(rows_*cols_);
355  matrix_coo.row_indices_.resize(rows_*cols_);
356  matrix_coo.col_indices_.resize(rows_*cols_);
357  // Fill all elements with the given value
358  for (ITYPE i = 0; i < rows_; ++i) {
359  std::fill(matrix_coo.row_indices_.begin()+(i*cols_), matrix_coo.row_indices_.end()+((i+1)*cols_), i);
360  std::iota(matrix_coo.col_indices_.begin()+(i*cols_), matrix_coo.col_indices_.end()+((i+1)*cols_), 0);
361  std::fill(matrix_coo.values_.begin()+(i*cols_), matrix_coo.values_.end()+((i+1)*cols_), val);
362  }
363  } else if (matrix_format_ == format::CSR) {
364  // Resize the vectors to hold all possible elements
365  matrix_csr.values_.reserve(rows_*cols_);
366  matrix_csr.row_ptrs_.reserve(rows_+1);
367  matrix_csr.col_indices_.reserve(rows_*cols_);
368  matrix_csr.values_.resize(rows_*cols_);
369  matrix_csr.row_ptrs_.resize(rows_+1);
370  matrix_csr.col_indices_.resize(rows_*cols_);
371  // Fill all elements with the given value
372  matrix_csr.row_ptrs_[0] = 0;
373  for (ITYPE i = 0; i < rows_; ++i) {
374  std::iota(matrix_csr.col_indices_.begin()+(i*cols_), matrix_csr.col_indices_.end()+((i+1)*cols_), 0);
375  std::fill(matrix_csr.values_.begin()+(i*cols_), matrix_csr.values_.end()+((i+1)*cols_), val);
376  matrix_csr.row_ptrs_[i + 1] = (i+1)*cols_;
377  }
378  } else if (matrix_format_ == format::CSC) {
379  // Resize the vectors to hold all possible elements
380  matrix_csc.values_.reserve(rows_*cols_);
381  matrix_csc.row_indices_.reserve(rows_*cols_);
382  matrix_csc.col_ptrs_.reserve(cols_+1);
383  matrix_csc.values_.resize(rows_*cols_);
384  matrix_csc.row_indices_.resize(rows_*cols_);
385  matrix_csc.col_ptrs_.resize(cols_+1);
386  // Fill all elements with the given value
387  matrix_csc.col_ptrs_[0] = 0;
388  for (ITYPE i = 0; i < cols_; ++i) {
389  std::iota(matrix_csc.row_indices_.begin()+(i*rows_), matrix_csc.row_indices_.end()+((i+1)*rows_), 0);
390  std::fill(matrix_csc.values_.begin()+(i*rows_), matrix_csc.values_.end()+((i+1)*rows_), val);
391  matrix_csc.col_ptrs_[i + 1] = (i+1)*rows_;
392  }
393  } else {
394  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
395  std::cerr << " Please set the matrix_format_ as:" << std::endl;
396  std::cerr << " format::COO: COOrdinate format" << std::endl;
397  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
398  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
399  std::abort();
400  }
401  nnz_ = rows_*cols_;
402  }
403 }
404 
405 // Member function to swap two elements in a sparse matrix
406 template<typename ITYPE, typename VTYPE>
407 void sparse_matrix<ITYPE,VTYPE>::swap_elements(ITYPE r1, ITYPE c1, ITYPE r2, ITYPE c2) {
408  assert(((r1 < rows_) && (r1 >= 0) && (c1 < cols_) && (c1 >= 0) &&
409  (r2 < rows_) && (r2 >= 0) && (c2 < cols_) && (c2 >= 0)) &&
410  "MUI Error [matrix_manipulation.h]: Matrix index out of range in swap_elements function");
411  VTYPE temp = this->get_value(r1, c1);
412  this->set_value(r1, c1, this->get_value(r2, c2));
413  this->set_value(r2, c2, temp);
414 }
415 
416 // Member function to set all elements to zero and empty the sparse matrix
417 template<typename ITYPE, typename VTYPE>
418 void sparse_matrix<ITYPE,VTYPE>::set_zero() {
419  // Clear the existing elements
420  if(matrix_format_ == format::COO) {
421  matrix_coo.values_.clear();
422  matrix_coo.row_indices_.clear();
423  matrix_coo.col_indices_.clear();
424  } else if (matrix_format_ == format::CSR) {
425  matrix_csr.values_.clear();
426  matrix_csr.row_ptrs_.clear();
427  matrix_csr.col_indices_.clear();
428  matrix_csr.row_ptrs_.resize((rows_+1), 0);
429  } else if (matrix_format_ == format::CSC) {
430  matrix_csc.values_.clear();
431  matrix_csc.row_indices_.clear();
432  matrix_csc.col_ptrs_.clear();
433  matrix_csc.col_ptrs_.resize((cols_+1), 0);
434  } else {
435  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
436  std::cerr << " Please set the matrix_format_ as:" << std::endl;
437  std::cerr << " format::COO: COOrdinate format" << std::endl;
438  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
439  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
440  std::abort();
441  }
442  nnz_ = 0;
443 }
444 
445 // Member function to add scalar to a specific elements
446 template<typename ITYPE, typename VTYPE>
447 void sparse_matrix<ITYPE,VTYPE>::add_scalar(ITYPE r, ITYPE c, VTYPE val, bool performSortAndUniqueCheck) {
448  assert(((r < rows_) && (r >= 0) && (c < cols_) && (c >= 0)) &&
449  "MUI Error [matrix_manipulation.h]: Matrix index out of range in add_scalar function");
450 
451  if (performSortAndUniqueCheck) {
452  if (!(this->is_sorted_unique("matrix_manipulation.h", "add_scalar()"))){
453  if(matrix_format_ == format::COO) {
454  this->sort_coo(true, true, "overwrite");
455  } else if (matrix_format_ == format::CSR) {
456  this->sort_csr(true, "overwrite");
457  } else if (matrix_format_ == format::CSC) {
458  this->sort_csc(true, "overwrite");
459  } else {
460  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
461  std::cerr << " Please set the matrix_format_ as:" << std::endl;
462  std::cerr << " format::COO: COOrdinate format" << std::endl;
463  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
464  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
465  std::abort();
466  }
467  }
468  }
469 
470  if(matrix_format_ == format::COO) {
471 
472  this->coo_element_operation(r, c, val, "plus", "matrix_manipulation.h", "add_scalar()");
473 
474  } else if (matrix_format_ == format::CSR) {
475 
476  this->csr_element_operation(r, c, val, "plus", "matrix_manipulation.h", "add_scalar()");
477 
478  } else if (matrix_format_ == format::CSC) {
479 
480  this->csc_element_operation(r, c, val, "plus", "matrix_manipulation.h", "add_scalar()");
481 
482  } else {
483  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
484  std::cerr << " Please set the matrix_format_ as:" << std::endl;
485  std::cerr << " format::COO: COOrdinate format" << std::endl;
486  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
487  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
488  std::abort();
489  }
490 }
491 
492 // Member function to subtract a scalar from a specific elements
493 template<typename ITYPE, typename VTYPE>
494 void sparse_matrix<ITYPE,VTYPE>::subtract_scalar(ITYPE r, ITYPE c, VTYPE val, bool performSortAndUniqueCheck) {
495  assert(((r < rows_) && (r >= 0) && (c < cols_) && (c >= 0)) &&
496  "MUI Error [matrix_manipulation.h]: Matrix index out of range in subtract_scalar function");
497  if (performSortAndUniqueCheck) {
498  if (!(this->is_sorted_unique("matrix_manipulation.h", "subtract_scalar()"))){
499  if(matrix_format_ == format::COO) {
500  this->sort_coo(true, true, "overwrite");
501  } else if (matrix_format_ == format::CSR) {
502  this->sort_csr(true, "overwrite");
503  } else if (matrix_format_ == format::CSC) {
504  this->sort_csc(true, "overwrite");
505  } else {
506  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
507  std::cerr << " Please set the matrix_format_ as:" << std::endl;
508  std::cerr << " format::COO: COOrdinate format" << std::endl;
509  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
510  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
511  std::abort();
512  }
513  }
514  }
515 
516  if(matrix_format_ == format::COO) {
517 
518  this->coo_element_operation(r, c, val, "minus", "matrix_manipulation.h", "subtract_scalar()");
519 
520  } else if (matrix_format_ == format::CSR) {
521 
522  this->csr_element_operation(r, c, val, "minus", "matrix_manipulation.h", "subtract_scalar()");
523 
524  } else if (matrix_format_ == format::CSC) {
525 
526  this->csc_element_operation(r, c, val, "minus", "matrix_manipulation.h", "subtract_scalar()");
527 
528  } else {
529  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
530  std::cerr << " Please set the matrix_format_ as:" << std::endl;
531  std::cerr << " format::COO: COOrdinate format" << std::endl;
532  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
533  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
534  std::abort();
535  }
536 }
537 
538 // Member function to multiply a scalar from a specific elements
539 template<typename ITYPE, typename VTYPE>
540 void sparse_matrix<ITYPE,VTYPE>::multiply_scalar(ITYPE r, ITYPE c, VTYPE val, bool performSortAndUniqueCheck) {
541  assert(((r < rows_) && (r >= 0) && (c < cols_) && (c >= 0)) &&
542  "MUI Error [matrix_manipulation.h]: Matrix index out of range in multiply_scalar function");
543 
544  if (performSortAndUniqueCheck) {
545  if (!(this->is_sorted_unique("matrix_manipulation.h", "multiply_scalar()"))){
546  if(matrix_format_ == format::COO) {
547  this->sort_coo(true, true, "overwrite");
548  } else if (matrix_format_ == format::CSR) {
549  this->sort_csr(true, "overwrite");
550  } else if (matrix_format_ == format::CSC) {
551  this->sort_csc(true, "overwrite");
552  } else {
553  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
554  std::cerr << " Please set the matrix_format_ as:" << std::endl;
555  std::cerr << " format::COO: COOrdinate format" << std::endl;
556  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
557  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
558  std::abort();
559  }
560  }
561  }
562 
563  if(matrix_format_ == format::COO) {
564 
565  this->coo_element_operation(r, c, val, "multiply", "matrix_manipulation.h", "multiply_scalar()");
566 
567  } else if (matrix_format_ == format::CSR) {
568 
569  this->csr_element_operation(r, c, val, "multiply", "matrix_manipulation.h", "multiply_scalar()");
570 
571  } else if (matrix_format_ == format::CSC) {
572 
573  this->csc_element_operation(r, c, val, "multiply", "matrix_manipulation.h", "multiply_scalar()");
574 
575  } else {
576  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
577  std::cerr << " Please set the matrix_format_ as:" << std::endl;
578  std::cerr << " format::COO: COOrdinate format" << std::endl;
579  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
580  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
581  std::abort();
582  }
583 }
584 
585 // Overloaded assignment operator
586 template <typename ITYPE, typename VTYPE>
587 sparse_matrix<ITYPE,VTYPE>& sparse_matrix<ITYPE,VTYPE>::operator=(const sparse_matrix<ITYPE,VTYPE> &exist_mat) {
588  if (this != &exist_mat) { // check for self-assignment
589  // copy the values from the other matrix to this matrix
590  assert(this->empty() &&
591  "MUI Error [matrix_manipulation.h]: assignment operator '=' only works for empty (all zero elements) matrix");
592 
593  if((rows_ == 0) && (cols_ == 0)){
594  rows_ = exist_mat.rows_;
595  cols_ = exist_mat.cols_;
596  }
597 
598  assert(((rows_ == exist_mat.rows_) && (cols_ == exist_mat.cols_)) &&
599  "MUI Error [matrix_manipulation.h]: matrix size mismatch in assignment operator '='");
600 
601  std::string format_store = this->get_format();
602  if (this->get_format() != exist_mat.get_format()) {
603  this->format_conversion(exist_mat.get_format(), true, true, "overwrite");
604  }
605 
606  // Copy the data from the existing matrix
607  if (matrix_format_ == format::COO) {
608  matrix_coo.values_ = std::vector<VTYPE>(exist_mat.matrix_coo.values_.begin(), exist_mat.matrix_coo.values_.end());
609  matrix_coo.row_indices_ = std::vector<ITYPE>(exist_mat.matrix_coo.row_indices_.begin(), exist_mat.matrix_coo.row_indices_.end());
610  matrix_coo.col_indices_ = std::vector<ITYPE>(exist_mat.matrix_coo.col_indices_.begin(), exist_mat.matrix_coo.col_indices_.end());
611  } else if (matrix_format_ == format::CSR) {
612  matrix_csr.values_ = std::vector<VTYPE>(exist_mat.matrix_csr.values_.begin(), exist_mat.matrix_csr.values_.end());
613  matrix_csr.row_ptrs_ = std::vector<ITYPE>(exist_mat.matrix_csr.row_ptrs_.begin(), exist_mat.matrix_csr.row_ptrs_.end());
614  matrix_csr.col_indices_ = std::vector<ITYPE>(exist_mat.matrix_csr.col_indices_.begin(), exist_mat.matrix_csr.col_indices_.end());
615  } else if (matrix_format_ == format::CSC) {
616  matrix_csc.values_ = std::vector<VTYPE>(exist_mat.matrix_csc.values_.begin(), exist_mat.matrix_csc.values_.end());
617  matrix_csc.row_indices_ = std::vector<ITYPE>(exist_mat.matrix_csc.row_indices_.begin(), exist_mat.matrix_csc.row_indices_.end());
618  matrix_csc.col_ptrs_ = std::vector<ITYPE>(exist_mat.matrix_csc.col_ptrs_.begin(), exist_mat.matrix_csc.col_ptrs_.end());
619  } else {
620  std::cerr << "MUI Error [matrix_ctor_dtor.h]: Unrecognised matrix format for matrix constructor" << std::endl;
621  std::cerr << " Please set the matrix_format_ as:" << std::endl;
622  std::cerr << " format::COO: COOrdinate format" << std::endl;
623  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
624  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
625  std::abort();
626  }
627 
628  nnz_ = exist_mat.nnz_;
629 
630  if (this->get_format() != format_store) {
631  this->format_conversion(format_store, true, true, "overwrite");
632  }
633  }
634  return *this;
635 }
636 
637 // Member function to convert the format of the sparse matrix
638 template<typename ITYPE, typename VTYPE>
639 void sparse_matrix<ITYPE,VTYPE>::format_conversion(const std::string &format, bool performSortAndUniqueCheck, bool deduplication, const std::string &deduplication_mode) {
640 
641  std::string matrix_format = string_to_upper(trim(format));
642 
643  if (matrix_format_ == format::COO) {
644 
645  if (matrix_format == "COO") {
646 
647  std::cout << "MUI [matrix_manipulation.h]: Convert matrix format from COO to COO (do nothing)." << std::endl;
648 
649  } else if (matrix_format == "CSR") {
650 
651  if (performSortAndUniqueCheck) {
652  if (!(this->is_sorted_unique("matrix_manipulation.h", "format_conversion()"))){
653  this->sparse_matrix<ITYPE,VTYPE>::sort_coo(true, deduplication, deduplication_mode);
654  }
655  }
656 
657  this->sparse_matrix<ITYPE,VTYPE>::coo_to_csr();
658 
659  } else if (matrix_format == "CSC") {
660 
661  this->sparse_matrix<ITYPE,VTYPE>::sort_coo(false, deduplication, deduplication_mode);
662 
663  this->sparse_matrix<ITYPE,VTYPE>::coo_to_csc();
664 
665  } else {
666 
667  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised format type: " << format << " for matrix format conversion" << std::endl;
668  std::cerr << " Please set the format string as:" << std::endl;
669  std::cerr << " 'COO': COOrdinate format" << std::endl;
670  std::cerr << " 'CSR' (default): Compressed Sparse Row format" << std::endl;
671  std::cerr << " 'CSC': Compressed Sparse Column format" << std::endl;
672  std::abort();
673 
674  }
675 
676  } else if (matrix_format_ == format::CSR) {
677 
678  if (matrix_format == "COO") {
679 
680  if (performSortAndUniqueCheck) {
681  if (!(this->is_sorted_unique("matrix_manipulation.h", "format_conversion()"))){
682  this->sparse_matrix<ITYPE,VTYPE>::sort_csr(deduplication, deduplication_mode);
683  }
684  }
685 
686  this->sparse_matrix<ITYPE,VTYPE>::csr_to_coo();
687 
688  } else if (matrix_format == "CSR") {
689 
690  std::cout << "MUI [matrix_manipulation.h]: Convert matrix format from CSR to CSR (do nothing)." << std::endl;
691 
692  } else if (matrix_format == "CSC") {
693 
694  if (performSortAndUniqueCheck) {
695  if (!(this->is_sorted_unique("matrix_manipulation.h", "format_conversion()"))){
696  this->sparse_matrix<ITYPE,VTYPE>::sort_csr(deduplication, deduplication_mode);
697  }
698  }
699 
700  this->sparse_matrix<ITYPE,VTYPE>::csr_to_csc();
701 
702  } else {
703 
704  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised format type: " << format << " for matrix format conversion" << std::endl;
705  std::cerr << " Please set the format string as:" << std::endl;
706  std::cerr << " 'COO': COOrdinate format" << std::endl;
707  std::cerr << " 'CSR' (default): Compressed Sparse Row format" << std::endl;
708  std::cerr << " 'CSC': Compressed Sparse Column format" << std::endl;
709  std::abort();
710 
711  }
712 
713  } else if (matrix_format_ == format::CSC) {
714 
715  if (matrix_format == "COO") {
716 
717  if (performSortAndUniqueCheck) {
718  if (!(this->is_sorted_unique("matrix_manipulation.h", "format_conversion()"))){
719  this->sparse_matrix<ITYPE,VTYPE>::sort_csc(deduplication, deduplication_mode);
720  }
721  }
722 
723  this->sparse_matrix<ITYPE,VTYPE>::csc_to_coo();
724 
725  } else if (matrix_format == "CSR") {
726 
727  if (performSortAndUniqueCheck) {
728  if (!(this->is_sorted_unique("matrix_manipulation.h", "format_conversion()"))){
729  this->sparse_matrix<ITYPE,VTYPE>::sort_csc(deduplication, deduplication_mode);
730  }
731  }
732 
733  this->sparse_matrix<ITYPE,VTYPE>::csc_to_csr();
734 
735  } else if (matrix_format == "CSC") {
736 
737  std::cout << "MUI [matrix_manipulation.h]: Convert matrix format from CSC to CSC (do nothing)." << std::endl;
738 
739  } else {
740 
741  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised format type: " << format << " for matrix format conversion" << std::endl;
742  std::cerr << " Please set the format string as:" << std::endl;
743  std::cerr << " 'COO': COOrdinate format" << std::endl;
744  std::cerr << " 'CSR' (default): Compressed Sparse Row format" << std::endl;
745  std::cerr << " 'CSC': Compressed Sparse Column format" << std::endl;
746  std::abort();
747 
748  }
749 
750  } else {
751 
752  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format for matrix format conversion" << std::endl;
753  std::cerr << " Please set the matrix_format_ as:" << std::endl;
754  std::cerr << " format::COO: COOrdinate format" << std::endl;
755  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
756  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
757  std::abort();
758 
759  }
760 }
761 
762 // Member function to sort the entries for the sparse matrix
763 template<typename ITYPE, typename VTYPE>
764 void sparse_matrix<ITYPE,VTYPE>::sort_deduplication(bool check_sorted_unique, bool deduplication, const std::string &deduplication_mode, bool is_row_major) {
765 
766  if (check_sorted_unique) {
767  if (!(this->is_sorted_unique("matrix_manipulation.h", "sort_deduplication()"))){
768  if(matrix_format_ == format::COO) {
769  this->sort_coo(is_row_major, deduplication, deduplication_mode);
770  } else if (matrix_format_ == format::CSR) {
771  this->sort_csr(deduplication, deduplication_mode);
772  } else if (matrix_format_ == format::CSC) {
773  this->sort_csc(deduplication, deduplication_mode);
774  } else {
775  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
776  std::cerr << " Please set the matrix_format_ as:" << std::endl;
777  std::cerr << " format::COO: COOrdinate format" << std::endl;
778  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
779  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
780  std::abort();
781  }
782  }
783  } else {
784  if(matrix_format_ == format::COO) {
785  this->sort_coo(is_row_major, deduplication, deduplication_mode);
786  } else if (matrix_format_ == format::CSR) {
787  this->sort_csr(deduplication, deduplication_mode);
788  } else if (matrix_format_ == format::CSC) {
789  this->sort_csc(deduplication, deduplication_mode);
790  } else {
791  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised matrix format" << std::endl;
792  std::cerr << " Please set the matrix_format_ as:" << std::endl;
793  std::cerr << " format::COO: COOrdinate format" << std::endl;
794  std::cerr << " format::CSR (default): Compressed Sparse Row format" << std::endl;
795  std::cerr << " format::CSC: Compressed Sparse Column format" << std::endl;
796  std::abort();
797  }
798  }
799 }
800 
801 
802 // **************************************************
803 // ********** Protected member functions ************
804 // **************************************************
805 
806 // Protected member function to sort the entries by row and column for sparse matrix with COO format
807 template<typename ITYPE, typename VTYPE>
808 void sparse_matrix<ITYPE,VTYPE>::sort_coo(bool is_row_major, bool deduplication, const std::string &deduplication_mode) {
809 
810  assert((matrix_format_ == format::COO) &&
811  "MUI Error [matrix_manipulation.h]: sort_coo() can only applied to sparse matrix with COO format. Please consider to convert the format into COO before use this function.");
812 
813  if (matrix_coo.values_.size() <= 1) return;
814 
815  std::string deduplication_mode_trim = string_to_lower(trim(deduplication_mode));
816 
817  if ((deduplication_mode_trim != "plus") && (deduplication_mode_trim != "minus") && (deduplication_mode_trim != "multiply") && (deduplication_mode_trim != "overwrite")) {
818  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised deduplication mode: " << deduplication_mode << " for sort_coo() function" << std::endl;
819  std::cerr << " Please set the deduplication mode as:" << std::endl;
820  std::cerr << " 'plus': Sum up values for all duplicated elements" << std::endl;
821  std::cerr << " 'minus': Take the difference among duplicated elements according to their position (former element minus later element)" << std::endl;
822  std::cerr << " 'multiply': Take the product among all duplicated elements" << std::endl;
823  std::cerr << " 'overwrite' (default): Keeps only the value of the last duplicated element" << std::endl;
824  std::abort();
825  }
826 
827  // Create a vector of indices to hold the original positions
828  std::vector<ITYPE> sorted_indices(matrix_coo.row_indices_.size());
829  std::iota(sorted_indices.begin(), sorted_indices.end(), 0);
830 
831  // Sort the indices based on row and column indices
832  if (is_row_major) {
833  std::sort(sorted_indices.begin(), sorted_indices.end(),
834  [&](ITYPE i, ITYPE j) {
835  return (matrix_coo.row_indices_[i] == matrix_coo.row_indices_[j])
836  ? matrix_coo.col_indices_[i] < matrix_coo.col_indices_[j]
837  : matrix_coo.row_indices_[i] < matrix_coo.row_indices_[j];
838  });
839  } else {
840  std::sort(sorted_indices.begin(), sorted_indices.end(),
841  [&](ITYPE i, ITYPE j) {
842  return (matrix_coo.col_indices_[i] == matrix_coo.col_indices_[j])
843  ? matrix_coo.row_indices_[i] < matrix_coo.row_indices_[j]
844  : matrix_coo.col_indices_[i] < matrix_coo.col_indices_[j];
845  });
846  }
847 
848  // Reorder the COO entries based on the sorted indices
849  std::vector<VTYPE> sorted_values;
850  std::vector<ITYPE> sorted_row_indices;
851  std::vector<ITYPE> sorted_column_indices;
852 
853  sorted_values.reserve(sorted_indices.size());
854  sorted_row_indices.reserve(sorted_indices.size());
855  sorted_column_indices.reserve(sorted_indices.size());
856 
857  for (ITYPE i = 0; i < static_cast<ITYPE>(sorted_indices.size()); ++i) {
858  if (std::abs(matrix_coo.values_[sorted_indices[i]]) >= std::numeric_limits<VTYPE>::min()) {
859  sorted_values.emplace_back(matrix_coo.values_[sorted_indices[i]]);
860  sorted_row_indices.emplace_back(matrix_coo.row_indices_[sorted_indices[i]]);
861  sorted_column_indices.emplace_back(matrix_coo.col_indices_[sorted_indices[i]]);
862  }
863  }
864 
865  if (deduplication) {
866 
867  // Deduplicate the COO entries based on the sorted matrix
868  std::vector<VTYPE> deduplicated_values;
869  std::vector<ITYPE> deduplicated_row_indices;
870  std::vector<ITYPE> deduplicated_column_indices;
871 
872  deduplicated_values.reserve(sorted_indices.size());
873  deduplicated_row_indices.reserve(sorted_indices.size());
874  deduplicated_column_indices.reserve(sorted_indices.size());
875 
876  ITYPE prev_row = sorted_row_indices[0];
877  ITYPE prev_col = sorted_column_indices[0];
878  VTYPE sum_value = sorted_values[0];
879  VTYPE difference_value = sorted_values[0];
880  VTYPE product_value = sorted_values[0];
881  VTYPE last_value = sorted_values[0];
882  bool is_multiply_by_zero = true;
883 
884  for (ITYPE i = 1; i < static_cast<ITYPE>(sorted_indices.size()); ++i) {
885  ITYPE curr_row = sorted_row_indices[i];
886  ITYPE curr_col = sorted_column_indices[i];
887 
888  if (i != static_cast<ITYPE>((sorted_indices.size()-1))) {
889  if (i == 1) {
890  if ((curr_row == prev_row) && (curr_col == prev_col)) {
891 
892  sum_value += sorted_values[i];
893  difference_value -= sorted_values[i];
894  product_value *= sorted_values[i];
895  last_value = sorted_values[i];
896  is_multiply_by_zero = false;
897 
898  } else {
899  if ((std::abs(sorted_values[0]) >= std::numeric_limits<VTYPE>::min()) && (deduplication_mode_trim != "multiply")) {
900  deduplicated_values.emplace_back(sorted_values[0]);
901  deduplicated_row_indices.emplace_back(prev_row);
902  deduplicated_column_indices.emplace_back(prev_col);
903  }
904 
905  prev_row = curr_row;
906  prev_col = curr_col;
907  sum_value = sorted_values[i];
908  difference_value = sorted_values[i];
909  product_value = sorted_values[i];
910  last_value = sorted_values[i];
911  is_multiply_by_zero = true;
912 
913  }
914  } else {
915  if ((curr_row == prev_row) && (curr_col == prev_col)) {
916 
917  sum_value += sorted_values[i];
918  difference_value -= sorted_values[i];
919  product_value *= sorted_values[i];
920  last_value = sorted_values[i];
921  is_multiply_by_zero = false;
922 
923  } else {
924 
925  if ((deduplication_mode_trim == "plus") && (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
926  deduplicated_values.emplace_back(sum_value);
927  deduplicated_row_indices.emplace_back(prev_row);
928  deduplicated_column_indices.emplace_back(prev_col);
929  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
930  deduplicated_values.emplace_back(difference_value);
931  deduplicated_row_indices.emplace_back(prev_row);
932  deduplicated_column_indices.emplace_back(prev_col);
933  } else if ((deduplication_mode_trim == "multiply") && (!is_multiply_by_zero) && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
934  deduplicated_values.emplace_back(product_value);
935  deduplicated_row_indices.emplace_back(prev_row);
936  deduplicated_column_indices.emplace_back(prev_col);
937  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
938  deduplicated_values.emplace_back(last_value);
939  deduplicated_row_indices.emplace_back(prev_row);
940  deduplicated_column_indices.emplace_back(prev_col);
941  }
942 
943  prev_row = curr_row;
944  prev_col = curr_col;
945  sum_value = sorted_values[i];
946  difference_value = sorted_values[i];
947  product_value = sorted_values[i];
948  last_value = sorted_values[i];
949  is_multiply_by_zero = true;
950 
951  }
952  }
953  } else {
954  if ((curr_row == prev_row) && (curr_col == prev_col)) {
955 
956  sum_value += sorted_values[i];
957  difference_value -= sorted_values[i];
958  product_value *= sorted_values[i];
959  last_value = sorted_values[i];
960  is_multiply_by_zero = false;
961 
962  if ((deduplication_mode_trim == "plus") && (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
963  deduplicated_values.emplace_back(sum_value);
964  deduplicated_row_indices.emplace_back(prev_row);
965  deduplicated_column_indices.emplace_back(prev_col);
966  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
967  deduplicated_values.emplace_back(difference_value);
968  deduplicated_row_indices.emplace_back(prev_row);
969  deduplicated_column_indices.emplace_back(prev_col);
970  } else if ((deduplication_mode_trim == "multiply") && (!is_multiply_by_zero) && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
971  deduplicated_values.emplace_back(product_value);
972  deduplicated_row_indices.emplace_back(prev_row);
973  deduplicated_column_indices.emplace_back(prev_col);
974  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
975  deduplicated_values.emplace_back(last_value);
976  deduplicated_row_indices.emplace_back(prev_row);
977  deduplicated_column_indices.emplace_back(prev_col);
978  }
979 
980  } else {
981 
982  if ((deduplication_mode_trim == "plus") && (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
983  deduplicated_values.emplace_back(sum_value);
984  deduplicated_row_indices.emplace_back(prev_row);
985  deduplicated_column_indices.emplace_back(prev_col);
986  deduplicated_values.emplace_back(sorted_values[i]);
987  deduplicated_row_indices.emplace_back(curr_row);
988  deduplicated_column_indices.emplace_back(curr_col);
989  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
990  deduplicated_values.emplace_back(difference_value);
991  deduplicated_row_indices.emplace_back(prev_row);
992  deduplicated_column_indices.emplace_back(prev_col);
993  deduplicated_values.emplace_back(sorted_values[i]);
994  deduplicated_row_indices.emplace_back(curr_row);
995  deduplicated_column_indices.emplace_back(curr_col);
996  } else if ((deduplication_mode_trim == "multiply") && (!is_multiply_by_zero) && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
997  deduplicated_values.emplace_back(product_value);
998  deduplicated_row_indices.emplace_back(prev_row);
999  deduplicated_column_indices.emplace_back(prev_col);
1000  deduplicated_values.emplace_back(sorted_values[i]);
1001  deduplicated_row_indices.emplace_back(curr_row);
1002  deduplicated_column_indices.emplace_back(curr_col);
1003  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
1004  deduplicated_values.emplace_back(last_value);
1005  deduplicated_row_indices.emplace_back(prev_row);
1006  deduplicated_column_indices.emplace_back(prev_col);
1007  deduplicated_values.emplace_back(sorted_values[i]);
1008  deduplicated_row_indices.emplace_back(curr_row);
1009  deduplicated_column_indices.emplace_back(curr_col);
1010  }
1011 
1012  }
1013  }
1014  }
1015 
1016  // Update the COOData struct with the deduplicated entries
1017  std::swap(matrix_coo.values_, deduplicated_values);
1018  std::swap(matrix_coo.row_indices_, deduplicated_row_indices);
1019  std::swap(matrix_coo.col_indices_, deduplicated_column_indices);
1020 
1021  nnz_ = matrix_coo.values_.size();
1022 
1023  } else {
1024 
1025  // Update the COOData struct with the sorted entries
1026  std::swap(matrix_coo.values_, sorted_values);
1027  std::swap(matrix_coo.row_indices_, sorted_row_indices);
1028  std::swap(matrix_coo.col_indices_, sorted_column_indices);
1029 
1030  nnz_ = matrix_coo.values_.size();
1031 
1032  }
1033 }
1034 
1035 // Protected member function to sort the entries for sparse matrix with CSR format
1036 template<typename ITYPE, typename VTYPE>
1037 void sparse_matrix<ITYPE,VTYPE>::sort_csr(bool deduplication, const std::string &deduplication_mode) {
1038 
1039  assert((matrix_format_ == format::CSR) &&
1040  "MUI Error [matrix_manipulation.h]: sort_csr() can only applied to sparse matrix with CSR format. Please consider to convert the format into CSR before use this function.");
1041 
1042  if (matrix_csr.values_.size() <= 1) return;
1043 
1044  std::string deduplication_mode_trim = string_to_lower(trim(deduplication_mode));
1045 
1046  if ((deduplication_mode_trim != "plus") && (deduplication_mode_trim != "minus") && (deduplication_mode_trim != "multiply") && (deduplication_mode_trim != "overwrite")) {
1047  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised deduplication mode: " << deduplication_mode << " for sort_csr() function" << std::endl;
1048  std::cerr << " Please set the deduplication mode as:" << std::endl;
1049  std::cerr << " 'plus': Sum up values for all duplicated elements" << std::endl;
1050  std::cerr << " 'minus': Take the difference among duplicated elements according to their position (former element minus later element)" << std::endl;
1051  std::cerr << " 'multiply': Take the product among all duplicated elements" << std::endl;
1052  std::cerr << " 'overwrite' (default): Keeps only the value of the last duplicated element" << std::endl;
1053  std::abort();
1054  }
1055 
1056  // Temporary vectors for sorted and deduplicated data
1057  std::vector<ITYPE> sorted_col_indices;
1058  std::vector<VTYPE> sorted_values;
1059 
1060  // Iterate over each row
1061  for (ITYPE row = 0; row < rows_; ++row) {
1062 
1063  assert((matrix_csr.row_ptrs_[row] <= matrix_csr.row_ptrs_[row+1]) &&
1064  "MUI Error [matrix_manipulation.h]: sort_csr() unable to sort and deduplication the sparse matrix with CSR format as the row_ptrs_ vector is not in correct order.");
1065 
1066  // Sort CSR column indices in this row
1067  std::vector<std::pair<ITYPE,VTYPE>> sorted_indices_pair;
1068  sorted_indices_pair.reserve(matrix_csr.row_ptrs_[row+1] - matrix_csr.row_ptrs_[row]);
1069 
1070  for(ITYPE j = matrix_csr.row_ptrs_[row]; j < matrix_csr.row_ptrs_[row+1]; ++j){
1071  sorted_indices_pair.emplace_back(std::make_pair(matrix_csr.col_indices_[j],matrix_csr.values_[j]));
1072  }
1073 
1074  std::sort(sorted_indices_pair.begin(),sorted_indices_pair.end(),[](const std::pair<ITYPE,VTYPE>& x, const std::pair<ITYPE,VTYPE>& y) {
1075  return x.first <= y.first;
1076  });
1077 
1078  for(ITYPE j = matrix_csr.row_ptrs_[row], n = 0; j < matrix_csr.row_ptrs_[row+1]; ++j, ++n){
1079  matrix_csr.col_indices_[j] = sorted_indices_pair[n].first;
1080  matrix_csr.values_[j] = sorted_indices_pair[n].second;
1081  }
1082 
1083  // Clear the temporary vector
1084  sorted_indices_pair.clear();
1085 
1086  if (deduplication) {
1087  // Vector to track unique column indices
1088  std::vector<ITYPE> uniqueColumns;
1089  uniqueColumns.reserve(matrix_csr.row_ptrs_[row + 1]-matrix_csr.row_ptrs_[row]);
1090 
1091  // Vector to track values corresponding to unique column indices
1092  std::vector<VTYPE> uniqueValues;
1093  uniqueValues.reserve(matrix_csr.row_ptrs_[row + 1]-matrix_csr.row_ptrs_[row]);
1094 
1095  // Iterate over the elements in the row
1096  for (ITYPE i = matrix_csr.row_ptrs_[row]; i < matrix_csr.row_ptrs_[row + 1]; ++i) {
1097  ITYPE col = matrix_csr.col_indices_[i];
1098  VTYPE val = matrix_csr.values_[i];
1099 
1100  // Check if the column index already exists
1101  auto it = std::lower_bound(uniqueColumns.begin(), uniqueColumns.end(), col);
1102  ITYPE index = std::distance(uniqueColumns.begin(), it);
1103  if (it != uniqueColumns.end() && *it == col) {
1104  // Column index already exists, overwrite the value
1105  VTYPE sum_value = uniqueValues[index] + val;
1106  VTYPE difference_value = uniqueValues[index] - val;
1107  VTYPE product_value = uniqueValues[index] * val;
1108  VTYPE last_value = val;
1109  if ((deduplication_mode_trim == "plus")&& (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
1110  uniqueValues[index] = sum_value;
1111  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
1112  uniqueValues[index] = difference_value;
1113  } else if ((deduplication_mode_trim == "multiply") && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
1114  uniqueValues[index] = product_value;
1115  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
1116  uniqueValues[index] = last_value;
1117  }
1118  } else {
1119  // Column index does not exist, insert it
1120  VTYPE sum_value = val;
1121  VTYPE difference_value = val;
1122  VTYPE product_value = val;
1123  VTYPE last_value = val;
1124  if ((deduplication_mode_trim == "plus")&& (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
1125  uniqueColumns.emplace_back(col);
1126  uniqueValues.emplace_back(sum_value);
1127  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
1128  uniqueColumns.emplace_back(col);
1129  uniqueValues.emplace_back(difference_value);
1130  } else if ((deduplication_mode_trim == "multiply") && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
1131  uniqueColumns.emplace_back(col);
1132  uniqueValues.emplace_back(product_value);
1133  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
1134  uniqueColumns.emplace_back(col);
1135  uniqueValues.emplace_back(last_value);
1136  }
1137  }
1138  }
1139 
1140  // Insert the unique column indices and values into the temporary vector
1141  sorted_col_indices.insert(sorted_col_indices.end(), uniqueColumns.begin(), uniqueColumns.end());
1142  sorted_values.insert(sorted_values.end(), uniqueValues.begin(), uniqueValues.end());
1143 
1144  // Update the row_ptr vector
1145  matrix_csr.row_ptrs_[row + 1] = matrix_csr.row_ptrs_[row] + uniqueColumns.size();
1146  }
1147  }
1148 
1149  if (deduplication) {
1150  // Assign the sorted and deduplicated data back to the col_index and value vectors
1151  matrix_csr.col_indices_.clear();
1152  matrix_csr.col_indices_.shrink_to_fit();
1153  matrix_csr.col_indices_.swap(sorted_col_indices);
1154  matrix_csr.values_.clear();
1155  matrix_csr.values_.shrink_to_fit();
1156  matrix_csr.values_.swap(sorted_values);
1157  nnz_ = matrix_csr.values_.size();
1158  }
1159 
1160  // Clear the temporary vector
1161  sorted_col_indices.clear();
1162  sorted_values.clear();
1163 
1164 }
1165 
1166 // Protected member function to sort the entries for sparse matrix with CSC format
1167 template<typename ITYPE, typename VTYPE>
1168 void sparse_matrix<ITYPE,VTYPE>::sort_csc(bool deduplication, const std::string &deduplication_mode) {
1169 
1170  assert((matrix_format_ == format::CSC) &&
1171  "MUI Error [matrix_manipulation.h]: sort_csc() can only applied to sparse matrix with CSC format. Please consider to convert the format into CSC before use this function.");
1172 
1173  if (matrix_csc.values_.size() <= 1) return;
1174 
1175  std::string deduplication_mode_trim = string_to_lower(trim(deduplication_mode));
1176 
1177  if ((deduplication_mode_trim != "plus") && (deduplication_mode_trim != "minus") && (deduplication_mode_trim != "multiply") && (deduplication_mode_trim != "overwrite")) {
1178  std::cerr << "MUI Error [matrix_manipulation.h]: Unrecognised deduplication mode: " << deduplication_mode << " for sort_csc() function" << std::endl;
1179  std::cerr << " Please set the deduplication mode as:" << std::endl;
1180  std::cerr << " 'plus': Sum up values for all duplicated elements" << std::endl;
1181  std::cerr << " 'minus': Take the difference among duplicated elements according to their position (former element minus later element)" << std::endl;
1182  std::cerr << " 'multiply': Take the product among all duplicated elements" << std::endl;
1183  std::cerr << " 'overwrite' (default): Keeps only the value of the last duplicated element" << std::endl;
1184  std::abort();
1185  }
1186 
1187  // Temporary vectors for sorted and deduplicated data
1188  std::vector<ITYPE> sorted_row_indices;
1189  std::vector<VTYPE> sorted_values;
1190 
1191  // Iterate over each column
1192  for (ITYPE col = 0; col < cols_; ++col) {
1193 
1194  assert((matrix_csc.col_ptrs_[col] <= matrix_csc.col_ptrs_[col+1]) &&
1195  "MUI Error [matrix_manipulation.h]: sort_csc() unable to sort and deduplication the sparse matrix with CSC format as the col_ptrs_ vector is not in correct order.");
1196 
1197  // Sort CSC row indices in this column
1198  std::vector<std::pair<ITYPE,VTYPE>> sorted_indices_pair;
1199  sorted_indices_pair.reserve(matrix_csc.col_ptrs_[col+1] - matrix_csc.col_ptrs_[col]);
1200 
1201  for(ITYPE j = matrix_csc.col_ptrs_[col]; j < matrix_csc.col_ptrs_[col+1]; ++j){
1202  sorted_indices_pair.emplace_back(std::make_pair(matrix_csc.row_indices_[j],matrix_csc.values_[j]));
1203  }
1204 
1205  std::sort(sorted_indices_pair.begin(),sorted_indices_pair.end(),[](const std::pair<ITYPE,VTYPE>& x, const std::pair<ITYPE,VTYPE>& y) {
1206  return x.first <= y.first;
1207  });
1208 
1209  for(ITYPE j = matrix_csc.col_ptrs_[col], n = 0; j < matrix_csc.col_ptrs_[col+1]; ++j, ++n){
1210  matrix_csc.row_indices_[j] = sorted_indices_pair[n].first;
1211  matrix_csc.values_[j] = sorted_indices_pair[n].second;
1212  }
1213 
1214  // Clear the temporary vector
1215  sorted_indices_pair.clear();
1216 
1217  if (deduplication) {
1218  // Vector to track unique row indices
1219  std::vector<ITYPE> uniqueRows;
1220  uniqueRows.reserve(matrix_csc.col_ptrs_[col + 1]-matrix_csc.col_ptrs_[col]);
1221 
1222  // Vector to track values corresponding to unique row indices
1223  std::vector<VTYPE> uniqueValues;
1224  uniqueValues.reserve(matrix_csc.col_ptrs_[col + 1]-matrix_csc.col_ptrs_[col]);
1225 
1226  // Iterate over the elements in the column
1227  for (ITYPE i = matrix_csc.col_ptrs_[col]; i < matrix_csc.col_ptrs_[col + 1]; ++i) {
1228  ITYPE row = matrix_csc.row_indices_[i];
1229  VTYPE val = matrix_csc.values_[i];
1230 
1231  // Check if the row index already exists
1232  auto it = std::lower_bound(uniqueRows.begin(), uniqueRows.end(), row);
1233  ITYPE index = std::distance(uniqueRows.begin(), it);
1234  if (it != uniqueRows.end() && *it == row) {
1235  // Row index already exists, overwrite the value
1236  VTYPE sum_value = uniqueValues[index] + val;
1237  VTYPE difference_value = uniqueValues[index] - val;
1238  VTYPE product_value = uniqueValues[index] * val;
1239  VTYPE last_value = val;
1240  if ((deduplication_mode_trim == "plus")&& (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
1241  uniqueValues[index] = sum_value;
1242  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
1243  uniqueValues[index] = difference_value;
1244  } else if ((deduplication_mode_trim == "multiply") && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
1245  uniqueValues[index] = product_value;
1246  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
1247  uniqueValues[index] = last_value;
1248  }
1249  } else {
1250  // Column index does not exist, insert it
1251  VTYPE sum_value = val;
1252  VTYPE difference_value = val;
1253  VTYPE product_value = val;
1254  VTYPE last_value = val;
1255  if ((deduplication_mode_trim == "plus")&& (std::abs(sum_value) >= std::numeric_limits<VTYPE>::min())) {
1256  uniqueRows.emplace_back(row);
1257  uniqueValues.emplace_back(sum_value);
1258  } else if ((deduplication_mode_trim == "minus") && (std::abs(difference_value) >= std::numeric_limits<VTYPE>::min())) {
1259  uniqueRows.emplace_back(row);
1260  uniqueValues.emplace_back(difference_value);
1261  } else if ((deduplication_mode_trim == "multiply") && (std::abs(product_value) >= std::numeric_limits<VTYPE>::min())) {
1262  uniqueRows.emplace_back(row);
1263  uniqueValues.emplace_back(product_value);
1264  } else if ((deduplication_mode_trim == "overwrite") && (std::abs(last_value) >= std::numeric_limits<VTYPE>::min())) {
1265  uniqueRows.emplace_back(row);
1266  uniqueValues.emplace_back(last_value);
1267  }
1268  }
1269  }
1270 
1271  // Insert the unique row indices and values into the temporary vector
1272  sorted_row_indices.insert(sorted_row_indices.end(), uniqueRows.begin(), uniqueRows.end());
1273  sorted_values.insert(sorted_values.end(), uniqueValues.begin(), uniqueValues.end());
1274 
1275  // Update the col_ptrs_ vector
1276  matrix_csc.col_ptrs_[col + 1] = matrix_csc.col_ptrs_[col] + uniqueRows.size();
1277  }
1278  }
1279 
1280  if (deduplication) {
1281  // Assign the sorted and deduplicated data back to the row_indices and value vectors
1282  matrix_csc.row_indices_.clear();
1283  matrix_csc.row_indices_.shrink_to_fit();
1284  matrix_csc.row_indices_.swap(sorted_row_indices);
1285  matrix_csc.values_.clear();
1286  matrix_csc.values_.shrink_to_fit();
1287  matrix_csc.values_.swap(sorted_values);
1288  nnz_ = matrix_csc.values_.size();
1289  }
1290 
1291  // Clear the temporary vector
1292  sorted_row_indices.clear();
1293  sorted_values.clear();
1294 
1295 }
1296 
1297 // Protected member function for element operation of COO matrix
1298 template<typename ITYPE, typename VTYPE>
1299 void sparse_matrix<ITYPE,VTYPE>::coo_element_operation(ITYPE r, ITYPE c, VTYPE val, const std::string &operation_mode, const std::string &file_name_input, const std::string &function_name_input) {
1300 
1301  std::string operation_mode_trim = string_to_lower(trim(operation_mode));
1302 
1303  std::string file_name;
1304  std::string function_name;
1305 
1306  if (file_name_input.empty()) {
1307  file_name = "matrix_manipulation.h";
1308  } else {
1309  file_name = file_name_input;
1310  }
1311 
1312  if (function_name_input.empty()) {
1313  function_name = "coo_element_operation()";
1314  } else {
1315  function_name = function_name_input;
1316  }
1317 
1318  if ((operation_mode_trim != "plus") and (operation_mode_trim != "minus") and (operation_mode_trim != "multiply") and (operation_mode_trim != "overwrite") and (operation_mode_trim != "nonsort")) {
1319  std::cerr << "MUI Error [" << file_name << "]: Unrecognised COO element operation mode: " << operation_mode_trim << " for " << function_name << " function" << std::endl;
1320  std::cerr << " Please set the COO element operation mode as:" << std::endl;
1321  std::cerr << " 'plus': Sum up elemental values" << std::endl;
1322  std::cerr << " 'minus': Take the difference among elemental values" << std::endl;
1323  std::cerr << " 'multiply': Take the product among elemental values" << std::endl;
1324  std::cerr << " 'overwrite' (default): Keeps only the last elemental value" << std::endl;
1325  std::cerr << " 'nonsort': Append the element value without sort or deduplication" << std::endl;
1326  std::abort();
1327  }
1328 
1329  if (operation_mode_trim == "nonsort") {
1330  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1331  matrix_coo.values_.reserve(matrix_coo.values_.size()+1);
1332  matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size()+1);
1333  matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size()+1);
1334  matrix_coo.values_.emplace_back(val);
1335  matrix_coo.row_indices_.emplace_back(r);
1336  matrix_coo.col_indices_.emplace_back(c);
1337  nnz_++;
1338  }
1339  return;
1340  }
1341 
1342  bool isElementAdded = false;
1343 
1344  // Find the range of column indices for the given row using lower_bound and upper_bound
1345  auto lower = std::lower_bound(matrix_coo.row_indices_.begin(), matrix_coo.row_indices_.end(), r);
1346  auto upper = std::upper_bound(matrix_coo.row_indices_.begin(), matrix_coo.row_indices_.end(), r);
1347 
1348  ITYPE insert_position = std::distance(matrix_coo.row_indices_.begin(), lower);
1349 
1350  // Iterate over the range of column indices for the given row
1351  for (auto it = lower; it != upper; ++it) {
1352  insert_position = std::distance(matrix_coo.row_indices_.begin(), it);
1353  if (matrix_coo.col_indices_[insert_position] < c) {
1354  // Do nothing
1355  } else if (matrix_coo.col_indices_[insert_position] == c) {
1356  // Found an existing entry with the same row and column, update the value
1357  if (operation_mode_trim == "plus") {
1358  // Check if the sum is zero, erase the element if so, overwrite the element if not
1359  if (std::abs(matrix_coo.values_[insert_position] + val) >= std::numeric_limits<VTYPE>::min()) {
1360  matrix_coo.values_[insert_position] += val;
1361  } else {
1362  matrix_coo.row_indices_.erase(matrix_coo.row_indices_.begin() + insert_position);
1363  matrix_coo.col_indices_.erase(matrix_coo.col_indices_.begin() + insert_position);
1364  matrix_coo.values_.erase(matrix_coo.values_.begin() + insert_position);
1365  nnz_--;
1366  }
1367  isElementAdded = true;
1368  break;
1369  } else if (operation_mode_trim == "minus") {
1370  // Check if the difference is zero, erase the element if so, overwrite the element if not
1371  if (std::abs(matrix_coo.values_[insert_position] - val) >= std::numeric_limits<VTYPE>::min()) {
1372  matrix_coo.values_[insert_position] -= val;
1373  } else {
1374  matrix_coo.row_indices_.erase(matrix_coo.row_indices_.begin() + insert_position);
1375  matrix_coo.col_indices_.erase(matrix_coo.col_indices_.begin() + insert_position);
1376  matrix_coo.values_.erase(matrix_coo.values_.begin() + insert_position);
1377  nnz_--;
1378  }
1379  isElementAdded = true;
1380  break;
1381  } else if (operation_mode_trim == "multiply") {
1382  // Check if the product is zero, erase the element if so, overwrite the element if not
1383  if (std::abs(matrix_coo.values_[insert_position] * val) >= std::numeric_limits<VTYPE>::min()) {
1384  matrix_coo.values_[insert_position] *= val;
1385  } else {
1386  matrix_coo.row_indices_.erase(matrix_coo.row_indices_.begin() + insert_position);
1387  matrix_coo.col_indices_.erase(matrix_coo.col_indices_.begin() + insert_position);
1388  matrix_coo.values_.erase(matrix_coo.values_.begin() + insert_position);
1389  nnz_--;
1390  }
1391  isElementAdded = true;
1392  break;
1393  } else { // overwrite
1394  // Check if the value is zero, erase the element if so, overwrite the element if not
1395  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1396  matrix_coo.values_[insert_position] = val;
1397  } else {
1398  matrix_coo.row_indices_.erase(matrix_coo.row_indices_.begin() + insert_position);
1399  matrix_coo.col_indices_.erase(matrix_coo.col_indices_.begin() + insert_position);
1400  matrix_coo.values_.erase(matrix_coo.values_.begin() + insert_position);
1401  nnz_--;
1402  }
1403  isElementAdded = true;
1404  break;
1405  }
1406  } else {
1407  if (operation_mode_trim != "multiply") {
1408  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1409  // Insert a new entry at the correct sorted position
1410  matrix_coo.values_.reserve(matrix_coo.values_.size()+1);
1411  matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size()+1);
1412  matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size()+1);
1413 
1414  matrix_coo.row_indices_.insert(matrix_coo.row_indices_.begin() + insert_position, r);
1415  matrix_coo.col_indices_.insert(matrix_coo.col_indices_.begin() + insert_position, c);
1416  if (operation_mode_trim == "minus") {
1417  matrix_coo.values_.insert(matrix_coo.values_.begin() + insert_position, -val);
1418  } else {
1419  matrix_coo.values_.insert(matrix_coo.values_.begin() + insert_position, val);
1420  }
1421  nnz_++;
1422  }
1423  }
1424  isElementAdded = true;
1425  break;
1426  }
1427  }
1428 
1429  if (operation_mode_trim != "multiply") {
1430  if ((!isElementAdded) && (std::abs(val) >= std::numeric_limits<VTYPE>::min())) {
1431  insert_position = std::distance(matrix_coo.row_indices_.begin(), upper);
1432  // Insert a new entry at the correct sorted position
1433  matrix_coo.values_.reserve(matrix_coo.values_.size()+1);
1434  matrix_coo.row_indices_.reserve(matrix_coo.row_indices_.size()+1);
1435  matrix_coo.col_indices_.reserve(matrix_coo.col_indices_.size()+1);
1436 
1437  matrix_coo.row_indices_.insert(matrix_coo.row_indices_.begin() + insert_position, r);
1438  matrix_coo.col_indices_.insert(matrix_coo.col_indices_.begin() + insert_position, c);
1439  if (operation_mode_trim == "minus") {
1440  matrix_coo.values_.insert(matrix_coo.values_.begin() + insert_position, -val);
1441  } else {
1442  matrix_coo.values_.insert(matrix_coo.values_.begin() + insert_position, val);
1443  }
1444  nnz_++;
1445  }
1446  }
1447 }
1448 
1449 // Protected member function for element operation of CSR matrix
1450 template<typename ITYPE, typename VTYPE>
1451 void sparse_matrix<ITYPE,VTYPE>::csr_element_operation(ITYPE r, ITYPE c, VTYPE val, const std::string &operation_mode, const std::string &file_name_input, const std::string &function_name_input) {
1452 
1453  std::string operation_mode_trim = string_to_lower(trim(operation_mode));
1454 
1455  std::string file_name;
1456  std::string function_name;
1457 
1458  if (file_name_input.empty()) {
1459  file_name = "matrix_manipulation.h";
1460  } else {
1461  file_name = file_name_input;
1462  }
1463 
1464  if (function_name_input.empty()) {
1465  function_name = "csr_element_operation()";
1466  } else {
1467  function_name = function_name_input;
1468  }
1469 
1470  if ((operation_mode_trim != "plus") and (operation_mode_trim != "minus") and (operation_mode_trim != "multiply") and (operation_mode_trim != "overwrite")) {
1471  std::cerr << "MUI Error [" << file_name << "]: Unrecognised CSR element operation mode: " << operation_mode_trim << " for " << function_name << " function" << std::endl;
1472  std::cerr << " Please set the CSR element operation mode as:" << std::endl;
1473  std::cerr << " 'plus': Sum up elemental values" << std::endl;
1474  std::cerr << " 'minus': Take the difference among elemental values" << std::endl;
1475  std::cerr << " 'multiply': Take the product among elemental values" << std::endl;
1476  std::cerr << " 'overwrite' (default): Keeps only the last elemental value" << std::endl;
1477  std::abort();
1478  }
1479 
1480  // Find the range of indices for the given row
1481  ITYPE start = matrix_csr.row_ptrs_[r];
1482  ITYPE end = matrix_csr.row_ptrs_[r + 1];
1483 
1484  // Find the column index position within the range
1485  auto it = std::lower_bound(matrix_csr.col_indices_.begin()+start, matrix_csr.col_indices_.begin()+end, c);
1486 
1487  // Check if the column index already exists in the row
1488  if (it != matrix_csr.col_indices_.begin()+end && *it == c) {
1489  // Get the index of the found column index
1490  ITYPE insert_position = std::distance(matrix_csr.col_indices_.begin(), it);
1491  if (operation_mode_trim == "plus") {
1492  // Check if the sum value is zero
1493  if (std::abs(matrix_csr.values_[insert_position] + val) >= std::numeric_limits<VTYPE>::min()) {
1494  // Update the existing value with the new value
1495  matrix_csr.values_[insert_position] += val;
1496  } else {
1497  // Erase the existing entry from the vectors
1498  matrix_csr.col_indices_.erase(matrix_csr.col_indices_.begin()+insert_position);
1499  matrix_csr.values_.erase(matrix_csr.values_.begin()+insert_position);
1500 
1501  // Adjust the row pointers after the erased element
1502  for (ITYPE i = r + 1; i < static_cast<ITYPE>(matrix_csr.row_ptrs_.size()); ++i) {
1503  matrix_csr.row_ptrs_[i]--;
1504  }
1505  nnz_--;
1506  }
1507  } else if (operation_mode_trim == "minus") {
1508  // Check if the difference value is zero
1509  if (std::abs(matrix_csr.values_[insert_position] - val) >= std::numeric_limits<VTYPE>::min()) {
1510  // Update the existing value with the new value
1511  matrix_csr.values_[insert_position] -= val;
1512  } else {
1513  // Erase the existing entry from the vectors
1514  matrix_csr.col_indices_.erase(matrix_csr.col_indices_.begin()+insert_position);
1515  matrix_csr.values_.erase(matrix_csr.values_.begin()+insert_position);
1516 
1517  // Adjust the row pointers after the erased element
1518  for (ITYPE i = r + 1; i < static_cast<ITYPE>(matrix_csr.row_ptrs_.size()); ++i) {
1519  matrix_csr.row_ptrs_[i]--;
1520  }
1521  nnz_--;
1522  }
1523  } else if (operation_mode_trim == "multiply") {
1524  // Check if the product value is zero
1525  if (std::abs(matrix_csr.values_[insert_position] * val) >= std::numeric_limits<VTYPE>::min()) {
1526  // Update the existing value with the new value
1527  matrix_csr.values_[insert_position] *= val;
1528  } else {
1529  // Erase the existing entry from the vectors
1530  matrix_csr.col_indices_.erase(matrix_csr.col_indices_.begin()+insert_position);
1531  matrix_csr.values_.erase(matrix_csr.values_.begin()+insert_position);
1532 
1533  // Adjust the row pointers after the erased element
1534  for (ITYPE i = r + 1; i < static_cast<ITYPE>(matrix_csr.row_ptrs_.size()); ++i) {
1535  matrix_csr.row_ptrs_[i]--;
1536  }
1537  nnz_--;
1538  }
1539  } else { // overwrite
1540  // Check if the new value is zero
1541  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1542  // Update the existing value with the new value
1543  matrix_csr.values_[insert_position] = val;
1544  } else {
1545  // Erase the existing entry from the vectors
1546  matrix_csr.col_indices_.erase(matrix_csr.col_indices_.begin()+insert_position);
1547  matrix_csr.values_.erase(matrix_csr.values_.begin()+insert_position);
1548 
1549  // Adjust the row pointers after the erased element
1550  for (ITYPE i = r + 1; i < static_cast<ITYPE>(matrix_csr.row_ptrs_.size()); ++i) {
1551  matrix_csr.row_ptrs_[i]--;
1552  }
1553  nnz_--;
1554  }
1555  }
1556  } else {
1557  if (operation_mode_trim != "multiply") {
1558  // Check if the new value is zero
1559  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1560  // Insert the new value and column index at the determined position
1561  ITYPE insert_position = std::distance(matrix_csr.col_indices_.begin(), it);
1562  matrix_csr.col_indices_.reserve(matrix_csr.col_indices_.size()+1);
1563  matrix_csr.values_.reserve(matrix_csr.values_.size()+1);
1564  matrix_csr.col_indices_.insert(matrix_csr.col_indices_.begin()+insert_position, c);
1565  if (operation_mode_trim == "minus") {
1566  matrix_csr.values_.insert(matrix_csr.values_.begin()+insert_position, -val);
1567  } else {
1568  matrix_csr.values_.insert(matrix_csr.values_.begin()+insert_position, val);
1569  }
1570 
1571  // Adjust the row pointers after the inserted element
1572  for (ITYPE i = r + 1; i < static_cast<ITYPE>(matrix_csr.row_ptrs_.size()); ++i) {
1573  matrix_csr.row_ptrs_[i]++;
1574  }
1575  nnz_++;
1576  } else {
1577  // No existing entry to update, do nothing
1578  }
1579  }
1580  }
1581 }
1582 
1583 // Protected member function for element operation of CSC matrix
1584 template<typename ITYPE, typename VTYPE>
1585 void sparse_matrix<ITYPE,VTYPE>::csc_element_operation(ITYPE r, ITYPE c, VTYPE val, const std::string &operation_mode, const std::string &file_name_input, const std::string &function_name_input) {
1586 
1587  std::string operation_mode_trim = string_to_lower(trim(operation_mode));
1588 
1589  std::string file_name;
1590  std::string function_name;
1591 
1592  if (file_name_input.empty()) {
1593  file_name = "matrix_manipulation.h";
1594  } else {
1595  file_name = file_name_input;
1596  }
1597 
1598  if (function_name_input.empty()) {
1599  function_name = "csc_element_operation()";
1600  } else {
1601  function_name = function_name_input;
1602  }
1603 
1604  if ((operation_mode_trim != "plus") and (operation_mode_trim != "minus") and (operation_mode_trim != "multiply") and (operation_mode_trim != "overwrite")) {
1605  std::cerr << "MUI Error [" << file_name << "]: Unrecognised CSC element operation mode: " << operation_mode_trim << " for " << function_name << " function" << std::endl;
1606  std::cerr << " Please set the CSC element operation mode as:" << std::endl;
1607  std::cerr << " 'plus': Sum up elemental values" << std::endl;
1608  std::cerr << " 'minus': Take the difference among elemental values" << std::endl;
1609  std::cerr << " 'multiply': Take the product among elemental values" << std::endl;
1610  std::cerr << " 'overwrite' (default): Keeps only the last elemental value" << std::endl;
1611  std::abort();
1612  }
1613 
1614 
1615  // Find the range of column indices for the given column
1616  ITYPE start = matrix_csc.col_ptrs_[c];
1617  ITYPE end = matrix_csc.col_ptrs_[c + 1];
1618 
1619  // Find the position of the row index within the range
1620  auto it = std::lower_bound(matrix_csc.row_indices_.begin()+start, matrix_csc.row_indices_.begin()+end, r);
1621 
1622  // Check if the row index already exists
1623  if (it != matrix_csc.row_indices_.begin()+end && *it == r) {
1624  // Get the index of the found column index
1625  ITYPE insert_position = std::distance(matrix_csc.row_indices_.begin(), it);
1626  if (operation_mode_trim == "plus") {
1627  // Check if the sum value is zero
1628  if (std::abs(matrix_csc.values_[insert_position] + val) >= std::numeric_limits<VTYPE>::min()) {
1629  // Update the existing value with the new value
1630  matrix_csc.values_[insert_position] += val;
1631  } else {
1632  // Erase the existing entry from the vectors
1633  matrix_csc.row_indices_.erase(matrix_csc.row_indices_.begin()+insert_position);
1634  matrix_csc.values_.erase(matrix_csc.values_.begin()+insert_position);
1635 
1636  // Adjust the column pointers after the erased element
1637  for (ITYPE i = c + 1; i < static_cast<ITYPE>(matrix_csc.col_ptrs_.size()); ++i) {
1638  matrix_csc.col_ptrs_[i]--;
1639  }
1640  nnz_--;
1641  }
1642  } else if (operation_mode_trim == "minus") {
1643  // Check if the difference value is zero
1644  if (std::abs(matrix_csc.values_[insert_position] - val) >= std::numeric_limits<VTYPE>::min()) {
1645  // Update the existing value with the new value
1646  matrix_csc.values_[insert_position] -= val;
1647  } else {
1648  // Erase the existing entry from the vectors
1649  matrix_csc.row_indices_.erase(matrix_csc.row_indices_.begin()+insert_position);
1650  matrix_csc.values_.erase(matrix_csc.values_.begin()+insert_position);
1651 
1652  // Adjust the column pointers after the erased element
1653  for (ITYPE i = c + 1; i < static_cast<ITYPE>(matrix_csc.col_ptrs_.size()); ++i) {
1654  matrix_csc.col_ptrs_[i]--;
1655  }
1656  nnz_--;
1657  }
1658  } else if (operation_mode_trim == "multiply") {
1659  // Check if the product value is zero
1660  if (std::abs(matrix_csc.values_[insert_position] * val) >= std::numeric_limits<VTYPE>::min()) {
1661  // Update the existing value with the new value
1662  matrix_csc.values_[insert_position] *= val;
1663  } else {
1664  // Erase the existing entry from the vectors
1665  matrix_csc.row_indices_.erase(matrix_csc.row_indices_.begin()+insert_position);
1666  matrix_csc.values_.erase(matrix_csc.values_.begin()+insert_position);
1667 
1668  // Adjust the column pointers after the erased element
1669  for (ITYPE i = c + 1; i < static_cast<ITYPE>(matrix_csc.col_ptrs_.size()); ++i) {
1670  matrix_csc.col_ptrs_[i]--;
1671  }
1672  nnz_--;
1673  }
1674  } else { // overwrite
1675  // Check if the new value is zero
1676  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1677  // Update the existing value with the new value
1678  matrix_csc.values_[insert_position] = val;
1679  } else {
1680  // Erase the existing entry from the vectors
1681  matrix_csc.row_indices_.erase(matrix_csc.row_indices_.begin()+insert_position);
1682  matrix_csc.values_.erase(matrix_csc.values_.begin()+insert_position);
1683 
1684  // Adjust the column pointers after the erased element
1685  for (ITYPE i = c + 1; i < static_cast<ITYPE>(matrix_csc.col_ptrs_.size()); ++i) {
1686  matrix_csc.col_ptrs_[i]--;
1687  }
1688  nnz_--;
1689  }
1690  }
1691  } else {
1692  if (operation_mode_trim != "multiply") {
1693  // Check if the new value is zero
1694  if (std::abs(val) >= std::numeric_limits<VTYPE>::min()) {
1695  // Insert the new value and row index at the determined position
1696  ITYPE insert_position = std::distance(matrix_csc.row_indices_.begin(), it);
1697  matrix_csc.row_indices_.reserve(matrix_csc.row_indices_.size()+1);
1698  matrix_csc.values_.reserve(matrix_csc.values_.size()+1);
1699  matrix_csc.row_indices_.insert(matrix_csc.row_indices_.begin()+insert_position, r);
1700  if (operation_mode_trim == "minus") {
1701  matrix_csc.values_.insert(matrix_csc.values_.begin()+insert_position, -val);
1702  } else {
1703  matrix_csc.values_.insert(matrix_csc.values_.begin()+insert_position, val);
1704  }
1705  // Adjust the row pointers after the inserted element
1706  for (ITYPE i = c + 1; i < static_cast<ITYPE>(matrix_csc.col_ptrs_.size()); ++i) {
1707  matrix_csc.col_ptrs_[i]++;
1708  }
1709  nnz_++;
1710  } else {
1711  // No existing entry to update, do nothing
1712  }
1713  }
1714  }
1715 
1716 }
1717 
1718 // Protected member function to convert COO matrix into CSR matrix
1719 template<typename ITYPE, typename VTYPE>
1720 void sparse_matrix<ITYPE,VTYPE>::coo_to_csr() {
1721 
1722  assert((matrix_format_ == format::COO) &&
1723  "MUI Error [matrix_manipulation.h]: coo_to_csr() can only applied to sparse matrix with COO format.");
1724 
1725  if (matrix_coo.values_.size() == 0) {
1726  matrix_format_ = format::CSR;
1727  matrix_csr.row_ptrs_.resize((rows_+1), 0);
1728  return;
1729  }
1730 
1731  // Determine the number of non-zero entries in each row
1732  std::vector<ITYPE> rowPtr(rows_+1, 0);
1733 
1734  for (ITYPE i = 0; i < static_cast<ITYPE>(matrix_coo.row_indices_.size()); ++i) {
1735  if (matrix_coo.row_indices_[i] >= (rows_+1)) {
1736  std::cerr << "MUI Error [matrix_manipulation.h]: row index: " << matrix_coo.row_indices_[i] << " at the " << i << "th matrix_coo.row_indices_ is out of range (" << (rows_+1) << ") in coo_to_csr()" << std::endl;
1737  std::abort();
1738  } else {
1739  rowPtr[matrix_coo.row_indices_[i]]++;
1740  }
1741 
1742  }
1743 
1744  for(ITYPE i = 0, accumulator = 0; i < rows_; ++i){
1745  ITYPE temp = rowPtr[i];
1746  rowPtr[i] = accumulator;
1747  accumulator += temp;
1748  }
1749  rowPtr[rows_] = nnz_;
1750 
1751  // Fill the CSRData struct
1752  matrix_csr.row_ptrs_.reserve(rows_+1);
1753  matrix_csr.values_.reserve(nnz_);
1754  matrix_csr.col_indices_.reserve(nnz_);
1755 
1756  std::swap(matrix_csr.values_, matrix_coo.values_);
1757  std::swap(matrix_csr.row_ptrs_, rowPtr);
1758  std::swap(matrix_csr.col_indices_, matrix_coo.col_indices_);
1759 
1760  // Deallocate the memory
1761  matrix_coo.values_.clear();
1762  matrix_coo.row_indices_.clear();
1763  matrix_coo.col_indices_.clear();
1764  rowPtr.clear();
1765 
1766  // Reset the matrix format
1767  matrix_format_ = format::CSR;
1768 }
1769 
1770 // Protected member function to convert COO matrix into CSC matrix
1771 template<typename ITYPE, typename VTYPE>
1772 void sparse_matrix<ITYPE,VTYPE>::coo_to_csc() {
1773 
1774  assert((matrix_format_ == format::COO) &&
1775  "MUI Error [matrix_manipulation.h]: coo_to_csc() can only applied to sparse matrix with COO format.");
1776 
1777  if (matrix_coo.values_.size() == 0) {
1778  matrix_format_ = format::CSC;
1779  matrix_csc.col_ptrs_.resize((cols_+1), 0);
1780  return;
1781  }
1782 
1783  // Determine the number of non-zero entries in each column
1784  std::vector<ITYPE> colPtr(cols_+1, 0);
1785 
1786  for (ITYPE i = 0; i < static_cast<ITYPE>(matrix_coo.col_indices_.size()); ++i) {
1787  colPtr[matrix_coo.col_indices_[i]]++;
1788  }
1789 
1790  for(ITYPE i = 0, accumulator = 0; i < cols_; ++i){
1791  ITYPE temp = colPtr[i];
1792  colPtr[i] = accumulator;
1793  accumulator += temp;
1794  }
1795  colPtr[cols_] = nnz_;
1796 
1797  // Fill the CSCData struct
1798  matrix_csc.values_.reserve(nnz_);
1799  matrix_csc.row_indices_.reserve(nnz_);
1800  matrix_csc.col_ptrs_.reserve(cols_+1);
1801 
1802  std::swap(matrix_csc.values_, matrix_coo.values_);
1803  std::swap(matrix_csc.row_indices_, matrix_coo.row_indices_);
1804  std::swap(matrix_csc.col_ptrs_, colPtr);
1805 
1806  // Deallocate the memory
1807  matrix_coo.values_.clear();
1808  matrix_coo.row_indices_.clear();
1809  matrix_coo.col_indices_.clear();
1810  colPtr.clear();
1811 
1812  // Reset the matrix format
1813  matrix_format_ = format::CSC;
1814 }
1815 
1816 // Protected member function to convert CSR matrix into COO matrix
1817 template<typename ITYPE, typename VTYPE>
1818 void sparse_matrix<ITYPE,VTYPE>::csr_to_coo() {
1819 
1820  assert((matrix_format_ == format::CSR) &&
1821  "MUI Error [matrix_manipulation.h]: csr_to_coo() can only applied to sparse matrix with CSR format.");
1822 
1823  if (matrix_csr.values_.size() == 0) {
1824  matrix_format_ = format::COO;
1825  return;
1826  }
1827 
1828  matrix_coo.values_.reserve(nnz_);
1829  matrix_coo.row_indices_.reserve(nnz_);
1830  matrix_coo.col_indices_.reserve(nnz_);
1831 
1832  // Populate the matrix_coo vectors from the matrix_csr vectors
1833  for (ITYPE row = 0; row < rows_; ++row) {
1834  for (ITYPE row_nnz = matrix_csr.row_ptrs_[row]; row_nnz < matrix_csr.row_ptrs_[row + 1]; ++row_nnz) {
1835  matrix_coo.values_.emplace_back(matrix_csr.values_[row_nnz]);
1836  matrix_coo.row_indices_.emplace_back(row);
1837  matrix_coo.col_indices_.emplace_back(matrix_csr.col_indices_[row_nnz]);
1838  }
1839  }
1840 
1841  // Deallocate the memory
1842  matrix_csr.values_.clear();
1843  matrix_csr.row_ptrs_.clear();
1844  matrix_csr.col_indices_.clear();
1845 
1846  // Reset the matrix format
1847  matrix_format_ = format::COO;
1848 }
1849 
1850 // Protected member function to convert CSR matrix into CSC matrix
1851 template<typename ITYPE, typename VTYPE>
1852 void sparse_matrix<ITYPE,VTYPE>::csr_to_csc() {
1853 
1854  assert((matrix_format_ == format::CSR) &&
1855  "MUI Error [matrix_manipulation.h]: csr_to_csc() can only applied to sparse matrix with CSR format.");
1856 
1857  if (matrix_csr.values_.size() == 0) {
1858  matrix_format_ = format::CSC;
1859  matrix_csc.col_ptrs_.resize((cols_+1), 0);
1860  return;
1861  }
1862 
1863  matrix_csc.col_ptrs_.resize((cols_+1), 0);
1864  matrix_csc.values_.resize(nnz_, 0);
1865  matrix_csc.row_indices_.resize(nnz_, 0);
1866 
1867  // Determine the number of non-zero entries in each column
1868  for (ITYPE i = 0; i < cols_; ++i) {
1869  for (ITYPE j = 0; j < static_cast<ITYPE>(matrix_csr.col_indices_.size()); ++j) {
1870  if (matrix_csr.col_indices_[j] == i) {
1871  ++matrix_csc.col_ptrs_[i];
1872  }
1873  }
1874  }
1875 
1876  for(ITYPE i = 0, accumulator = 0; i < cols_; ++i){
1877  ITYPE temp = matrix_csc.col_ptrs_[i];
1878  matrix_csc.col_ptrs_[i] = accumulator;
1879  accumulator += temp;
1880  }
1881 
1882  matrix_csc.col_ptrs_[cols_] = nnz_;
1883 
1884  for (ITYPE row = 0; row < rows_; ++row) {
1885  for (ITYPE row_nnz = matrix_csr.row_ptrs_[row]; row_nnz < matrix_csr.row_ptrs_[row + 1]; ++row_nnz) {
1886  ITYPE col = matrix_csr.col_indices_[row_nnz];
1887  ITYPE dest = matrix_csc.col_ptrs_[col];
1888  matrix_csc.row_indices_[dest] = row;
1889  matrix_csc.values_[dest] = matrix_csr.values_[row_nnz];
1890  matrix_csc.col_ptrs_[col]++;
1891  }
1892  }
1893 
1894  for (ITYPE col = 0, last = 0; col <= cols_; ++col) {
1895  ITYPE temp = matrix_csc.col_ptrs_[col];
1896  matrix_csc.col_ptrs_[col] = last;
1897  last = temp;
1898  }
1899 
1900  // Deallocate the memory
1901  matrix_csr.values_.clear();
1902  matrix_csr.row_ptrs_.clear();
1903  matrix_csr.col_indices_.clear();
1904 
1905  // Reset the matrix format
1906  matrix_format_ = format::CSC;
1907 
1908 }
1909 
1910 // Protected member function to convert CSC matrix into COO matrix
1911 template<typename ITYPE, typename VTYPE>
1912 void sparse_matrix<ITYPE,VTYPE>::csc_to_coo() {
1913 
1914  assert((matrix_format_ == format::CSC) &&
1915  "MUI Error [matrix_manipulation.h]: csc_to_coo() can only applied to sparse matrix with CSC format.");
1916 
1917  if (matrix_csc.values_.size() == 0) {
1918  matrix_format_ = format::COO;
1919  return;
1920  }
1921 
1922  matrix_coo.values_.reserve(nnz_);
1923  matrix_coo.row_indices_.reserve(nnz_);
1924  matrix_coo.col_indices_.reserve(nnz_);
1925 
1926  // Populate the matrix_coo vectors from the matrix_csc vectors
1927  for (ITYPE col = 0; col < cols_; ++col) {
1928  for (ITYPE col_nnz = matrix_csc.col_ptrs_[col]; col_nnz < matrix_csc.col_ptrs_[col + 1]; ++col_nnz) {
1929  matrix_coo.values_.emplace_back(matrix_csc.values_[col_nnz]);
1930  matrix_coo.row_indices_.emplace_back(matrix_csc.row_indices_[col_nnz]);
1931  matrix_coo.col_indices_.emplace_back(col);
1932  }
1933  }
1934 
1935  // Deallocate the memory
1936  matrix_csc.values_.clear();
1937  matrix_csc.row_indices_.clear();
1938  matrix_csc.col_ptrs_.clear();
1939 
1940  // Reset the matrix format
1941  matrix_format_ = format::COO;
1942 }
1943 
1944 // Protected member function to convert CSC matrix into CSR matrix
1945 template<typename ITYPE, typename VTYPE>
1946 void sparse_matrix<ITYPE,VTYPE>::csc_to_csr() {
1947 
1948  assert((matrix_format_ == format::CSC) &&
1949  "MUI Error [matrix_manipulation.h]: csc_to_csr() can only applied to sparse matrix with CSC format.");
1950 
1951  if (matrix_csc.values_.size() == 0) {
1952  matrix_format_ = format::CSR;
1953  matrix_csr.row_ptrs_.resize((rows_+1), 0);
1954  return;
1955  }
1956 
1957  matrix_csr.values_.resize(nnz_, 0);
1958  matrix_csr.row_ptrs_.resize((rows_+1), 0);
1959  matrix_csr.col_indices_.resize(nnz_, 0);
1960 
1961  // Determine the number of non-zero entries in each row
1962  for (ITYPE i = 0; i < static_cast<ITYPE>(matrix_csc.row_indices_.size()); ++i) {
1963  matrix_csr.row_ptrs_[matrix_csc.row_indices_[i]]++;
1964  }
1965 
1966  for(ITYPE i = 0, accumulator = 0; i < rows_; ++i){
1967  ITYPE temp = matrix_csr.row_ptrs_[i];
1968  matrix_csr.row_ptrs_[i] = accumulator;
1969  accumulator += temp;
1970  }
1971 
1972  matrix_csr.row_ptrs_[rows_] = nnz_;
1973 
1974  for (ITYPE col = 0; col < cols_; ++col) {
1975  for (ITYPE col_nnz = matrix_csc.col_ptrs_[col]; col_nnz < matrix_csc.col_ptrs_[col + 1]; ++col_nnz) {
1976  ITYPE row = matrix_csc.row_indices_[col_nnz];
1977  ITYPE dest = matrix_csr.row_ptrs_[row];
1978  matrix_csr.col_indices_[dest] = col;
1979  matrix_csr.values_[dest] = matrix_csc.values_[col_nnz];
1980  matrix_csr.row_ptrs_[row]++;
1981  }
1982  }
1983 
1984  for (ITYPE row = 0, last = 0; row <= rows_; ++row) {
1985  ITYPE temp = matrix_csr.row_ptrs_[row];
1986  matrix_csr.row_ptrs_[row] = last;
1987  last = temp;
1988  }
1989 
1990  // Deallocate the memory
1991  matrix_csc.values_.clear();
1992  matrix_csc.row_indices_.clear();
1993  matrix_csc.col_ptrs_.clear();
1994 
1995  // Reset the matrix format
1996  matrix_format_ = format::CSR;
1997 
1998 }
1999 
2000 // Protected member function to clear all vectors of the sparse matrix
2001 template<typename ITYPE, typename VTYPE>
2002 void sparse_matrix<ITYPE,VTYPE>::clear_vectors() {
2003  // Clear the existing elements
2004  matrix_coo.values_.clear();
2005  matrix_coo.row_indices_.clear();
2006  matrix_coo.col_indices_.clear();
2007  matrix_csr.values_.clear();
2008  matrix_csr.row_ptrs_.clear();
2009  matrix_csr.col_indices_.clear();
2010  matrix_csc.values_.clear();
2011  matrix_csc.row_indices_.clear();
2012  matrix_csc.col_ptrs_.clear();
2013  nnz_ = 0;
2014 }
2015 
2016 } // linalg
2017 } // mui
2018 
2019 #endif /* MUI_MATRIX_MANIPULATION_H_ */
mui::linalg::sparse_matrix
Definition: matrix.h:61
mui::linalg::sparse_matrix::clear_vectors
void clear_vectors()
Definition: matrix_manipulation.h:2002
mui::linalg::sparse_matrix::sort_deduplication
void sort_deduplication(bool=true, bool=true, const std::string &="overwrite", bool=true)
Definition: matrix_manipulation.h:764
mui::linalg::sparse_matrix::csr_to_csc
void csr_to_csc()
Definition: matrix_manipulation.h:1852
mui::linalg::sparse_matrix::csc_element_operation
void csc_element_operation(ITYPE, ITYPE, VTYPE, const std::string &, const std::string &={}, const std::string &={})
Definition: matrix_manipulation.h:1585
mui::linalg::sparse_matrix::coo_element_operation
void coo_element_operation(ITYPE, ITYPE, VTYPE, const std::string &, const std::string &={}, const std::string &={})
Definition: matrix_manipulation.h:1299
mui::linalg::sparse_matrix::sort_csr
void sort_csr(bool=false, const std::string &="overwrite")
Definition: matrix_manipulation.h:1037
mui::linalg::sparse_matrix::set_zero
void set_zero()
Definition: matrix_manipulation.h:418
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=(const sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_manipulation.h:587
mui::linalg::sparse_matrix::segment
sparse_matrix< ITYPE, VTYPE > segment(ITYPE, ITYPE, ITYPE, ITYPE, bool=true)
Definition: matrix_manipulation.h:153
mui::linalg::sparse_matrix::csc_to_coo
void csc_to_coo()
Definition: matrix_manipulation.h:1912
mui::linalg::sparse_matrix::get_format
std::string get_format() const
Definition: matrix_io_info.h:681
mui::linalg::sparse_matrix::coo_to_csc
void coo_to_csc()
Definition: matrix_manipulation.h:1772
mui::linalg::sparse_matrix::assert_valid_vector_size
void assert_valid_vector_size(const std::string &={}, const std::string &={}) const
Definition: matrix_asserts.h:61
mui::linalg::sparse_matrix::resize
void resize(ITYPE, ITYPE)
Definition: matrix_manipulation.h:62
mui::linalg::sparse_matrix::add_scalar
void add_scalar(ITYPE, ITYPE, VTYPE, bool=true)
Definition: matrix_manipulation.h:447
mui::linalg::sparse_matrix::swap_elements
void swap_elements(ITYPE, ITYPE, ITYPE, ITYPE)
Definition: matrix_manipulation.h:407
mui::linalg::sparse_matrix::multiply_scalar
void multiply_scalar(ITYPE, ITYPE, VTYPE, bool=true)
Definition: matrix_manipulation.h:540
mui::linalg::sparse_matrix::csr_element_operation
void csr_element_operation(ITYPE, ITYPE, VTYPE, const std::string &, const std::string &={}, const std::string &={})
Definition: matrix_manipulation.h:1451
mui::linalg::sparse_matrix::copy
void copy(const sparse_matrix< ITYPE, VTYPE > &)
Definition: matrix_manipulation.h:79
mui::linalg::sparse_matrix::coo_to_csr
void coo_to_csr()
Definition: matrix_manipulation.h:1720
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::csc_to_csr
void csc_to_csr()
Definition: matrix_manipulation.h:1946
mui::linalg::sparse_matrix::csr_to_coo
void csr_to_coo()
Definition: matrix_manipulation.h:1818
mui::linalg::sparse_matrix::subtract_scalar
void subtract_scalar(ITYPE, ITYPE, VTYPE, bool=true)
Definition: matrix_manipulation.h:494
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::string_to_lower
std::string string_to_lower(const std::string &s)
Definition: linalg_util.h:78
mui::linalg::trim
std::string trim(const std::string &s)
Definition: linalg_util.h:73
mui::linalg::string_to_upper
std::string string_to_upper(const std::string &s)
Definition: linalg_util.h:85
mui
Definition: comm.h:54
mui::swap
void swap(storage< Args... > &lhs, storage< Args... > &rhs)
Definition: dynstorage.h:234