suanPan
sort_rcm.h
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
40#ifndef RCM_H
41#define RCM_H
42
45#include <Toolbox/container.h>
46
47uvec sort_rcm(const std::vector<uvec>&, const uvec&);
48
50
51template<sp_d eT> uvec sort_rcm(const SpMat<eT>& MEAT) {
52 suanpan_assert([&] { if(!MEAT.is_square()) throw logic_error("can only be applied to square matrix"); });
53
54 wall_clock TM;
55 TM.tic();
56
58 auto S = MEAT.n_cols;
59
61 uvec E(S, fill::none);
62 suanpan::for_each(S, [&](const uword I) { E(I) = MEAT.col(I).n_nonzero; });
63
64 std::vector<uvec> A(S);
65 suanpan::for_each(S, [&](const uword K) {
66 unsigned J = 0;
67 uvec IDX(E(K));
68 for(auto L = MEAT.begin_col(K); L != MEAT.end_col(K); ++L) IDX(J++) = L.row();
69 A[K] = IDX(sort_index(E(IDX)));
70 });
71
73 uvec G = sort_index(E);
74
76 std::vector<bool> M(S, false);
78 uvec R(S, fill::zeros);
79
84
87 uword IDXA = 0, IDXB = S - 1, IDXC = S - 1;
88
91
97 while(IDXA < S) {
98 if(IDXB == IDXC) {
100 while(IDXA < S && M[G(IDXA)]) ++IDXA;
103 if(IDXA == S) break;
105 R(IDXC--) = G(IDXA);
107 M[G(IDXA++)] = true;
108 }
114 for(const auto& IDX : A[R(IDXB--)]) if(!M[IDX]) M[R(IDXC--) = IDX] = true;
115 }
116
117 suanpan_debug("RCM algorithm takes {:.5E} seconds.\n", TM.toc());
118
119 return R;
120}
121
122template<sp_d eT> uvec sort_rcm(const Mat<eT>& MEAT) { return sort_rcm(SpMat<eT>(MEAT)); }
123
124template<sp_d dt> uvec sort_rcm(const csc_form<dt, uword>& csc_mat) {
126 const auto S = csc_mat.n_cols;
127
129 uvec E(S, fill::none);
130 suanpan::for_each(S, [&](const uword I) { E(I) = csc_mat.col(I + 1) - csc_mat.col(I); });
131 std::vector<uvec> A(S);
132 suanpan::for_each(S, [&](const uword I) {
133 const uvec IDX(csc_mat.row_mem() + csc_mat.col(I), E(I));
134 A[I] = IDX(sort_index(E(IDX)));
135 });
136
137 return sort_rcm(A, E);
138}
139
140template<sp_d dt, sp_i it> uvec sort_rcm(triplet_form<dt, it>& triplet_mat) { return sort_rcm(csc_form<dt, uword>(triplet_mat)); }
141
142#endif
143
Definition: csc_form.hpp:25
const index_t n_cols
Definition: csc_form.hpp:51
const index_t * row_mem() const
Definition: csc_form.hpp:61
index_t col(const index_t I) const
Definition: csc_form.hpp:75
Definition: triplet_form.hpp:62
uvec sort_rcm(const std::vector< uvec > &, const uvec &)
Definition: sort_rcm.cpp:78
std::vector< T > vector
Definition: container.h:53
std::unordered_set< T > unordered_set
Definition: container.h:55
void for_each(const IT start, const IT end, F &&FN)
Definition: utility.h:28
#define suanpan_debug(...)
Definition: suanPan.h:307
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296