sort_rcm.h

Go to the documentation of this file.

/*******************************************************************************
 * Copyright (C) 2017-2023 Theodore Chang
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ******************************************************************************/
#ifndef RCM_H
#define RCM_H
 
#include <Domain/MetaMat/triplet_form.hpp>
#include <Domain/MetaMat/csc_form.hpp>
 
uvec sort_rcm(const std::vector<uvec>&, const uvec&);
 
template<typename eT> uvec sort_rcm(const SpMat<eT>& MEAT) {
    suanpan_assert([&] { if(!MEAT.is_square()) throw logic_error("can only be applied to square matrix"); });
 
    wall_clock TM;
    TM.tic();
 
    auto S = MEAT.n_cols;
 
    uvec E(S, fill::none);
    suanpan_for(0llu, S, [&](const uword I) { E(I) = MEAT.col(I).n_nonzero; });
 
    std::vector<uvec> A(S);
    suanpan_for(0llu, S, [&](const uword K) {
        unsigned J = 0;
        uvec IDX(E(K));
        for(auto L = MEAT.begin_col(K); L != MEAT.end_col(K); ++L) IDX(J++) = L.row();
        A[K] = IDX(sort_index(E(IDX)));
    });
 
    uvec G = sort_index(E);
 
    std::vector<bool> M(S, false);
    uvec R(S, fill::zeros);
 
 
    uword IDXA = 0, IDXB = S - 1, IDXC = S - 1;
 
 
    while(IDXA < S) {
        if(IDXB == IDXC) {
            while(IDXA < S && M[G(IDXA)]) ++IDXA;
            if(IDXA == S) break;
            R(IDXC--) = G(IDXA);
            M[G(IDXA++)] = true;
        }
        for(const auto& IDX : A[R(IDXB--)]) if(!M[IDX]) M[R(IDXC--) = IDX] = true;
    }
 
    suanpan_debug("RCM algorithm takes {:.5E} seconds.\n", TM.toc());
 
    return R;
}
 
template<typename eT> uvec sort_rcm(const Mat<eT>& MEAT) { return sort_rcm(SpMat<eT>(MEAT)); }
 
template<typename dt> uvec sort_rcm(const csc_form<dt, uword>& csc_mat) {
    auto S = csc_mat.n_cols;
 
    uvec E(S, fill::none);
    suanpan_for(0llu, S, [&](const uword I) { E(I) = csc_mat.col(I + 1) - csc_mat.col(I); });
    std::vector<uvec> A(S);
    suanpan_for(0llu, S, [&](const uword I) {
        const uvec IDX(csc_mat.row_mem() + csc_mat.col(I), E(I));
        A[I] = IDX(sort_index(E(IDX)));
    });
 
    return sort_rcm(A, E);
}
 
template<typename dt, typename it> uvec sort_rcm(triplet_form<dt, it>& triplet_mat) {
    csc_form<dt, uword> csc_mat(triplet_mat);
 
    return sort_rcm(csc_mat);
}
 
#endif