sort_rcm.h

Go to the documentation of this file.

/*******************************************************************************
 * Copyright (C) 2017-2022 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>

uvec sort_rcm(const std::vector<uvec>&, const uvec&);

template<typename eT> uvec sort_rcm(const SpMat<eT>& MEAT) {
#ifdef SUANPAN_DEBUG
    if(!MEAT.is_square()) throw logic_error("RCM() can only be applied to square matrix.\n");

    wall_clock TM;
    TM.tic();
#endif

    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;
    }

#ifdef SUANPAN_DEBUG
    suanpan_debug("RCM algorithm takes %.5E seconds.\n", TM.toc());
#endif

    return R;
}

template<typename eT> uvec sort_rcm(const Mat<eT>& MEAT) { return sort_rcm(SpMat<eT>(MEAT)); }

template<typename dt, typename it> uvec sort_rcm(const triplet_form<dt, it>& triplet_mat) {
    csc_form<dt, uword> csc_mat(triplet_mat);

    const uvec row_idx(csc_mat.row_idx, csc_mat.n_elem, false, false);
    const uvec col_ptr(csc_mat.col_ptr, csc_mat.n_cols + 1, false, false);
    const Col<dt> val_idx(csc_mat.val_idx, csc_mat.n_elem, false, false);

    return sort_rcm({row_idx, col_ptr, val_idx, csc_mat.n_rows, csc_mat.n_cols});
}

#endif