From 7c76e36d7df5b84c2ebb3c8f0e17f5367a9d1cf6 Mon Sep 17 00:00:00 2001 From: Andrew Reisner Date: Thu, 12 Sep 2024 10:38:19 -0600 Subject: [PATCH] Cleanup. --- raptor/ruge_stuben/par_air_solver.hpp | 20 +++---- raptor/ruge_stuben/par_interpolation.cpp | 25 +++----- raptor/ruge_stuben/tests/test_air.cpp | 72 ++---------------------- 3 files changed, 21 insertions(+), 96 deletions(-) diff --git a/raptor/ruge_stuben/par_air_solver.hpp b/raptor/ruge_stuben/par_air_solver.hpp index b1b72883..0c5fa164 100644 --- a/raptor/ruge_stuben/par_air_solver.hpp +++ b/raptor/ruge_stuben/par_air_solver.hpp @@ -10,16 +10,16 @@ namespace raptor { struct ParAIRSolver : ParMultilevel { - ParAIRSolver(double strong_threshold = 0.0, - coarsen_t coarsen_type = RS, - interp_t interp_type = OnePoint, - strength_t strength_type = Classical, - restrict_t restrict_type = AIR, - relax_t relax_type = FCJacobi) : - ParMultilevel(strong_threshold, strength_type, relax_type), - coarsen_type(coarsen_type), - interp_type(interp_type), - restrict_type(restrict_type), + ParAIRSolver(double strong_threshold_ = 0.0, + coarsen_t coarsen_type_ = RS, + interp_t interp_type_ = OnePoint, + strength_t strength_type_ = Classical, + restrict_t restrict_type_ = AIR, + relax_t relax_type_ = FCJacobi) : + ParMultilevel(strong_threshold_, strength_type_, relax_type_), + coarsen_type(coarsen_type_), + interp_type(interp_type_), + restrict_type(restrict_type_), variables(nullptr) { num_variables = 1; diff --git a/raptor/ruge_stuben/par_interpolation.cpp b/raptor/ruge_stuben/par_interpolation.cpp index 08d2c883..42908392 100644 --- a/raptor/ruge_stuben/par_interpolation.cpp +++ b/raptor/ruge_stuben/par_interpolation.cpp @@ -20,17 +20,6 @@ CSRMatrix* communicate(ParCSRMatrix* A, ParCSRMatrix* S, const std::vector& CSRMatrix* communicate(ParCSRMatrix* A, const std::vector& states, const std::vector& off_proc_states, CommPkg* comm); void filter_interp(ParCSRMatrix* P, const double filter_threshold); -ParCSRMatrix* extended_interpolation(ParCSRMatrix* A, - ParCSRMatrix* S, const std::vector& states, - const std::vector& off_proc_states, const double filter_threshold, - bool tap_interp, int num_variables, int* variables); -ParCSRMatrix* mod_classical_interpolation(ParCSRMatrix* A, - ParCSRMatrix* S, const std::vector& states, - const std::vector& off_proc_states, - bool tap_interp, int num_variables, int* variables); -ParCSRMatrix* direct_interpolation(ParCSRMatrix* A, - ParCSRMatrix* S, const std::vector& states, - const std::vector& off_proc_states, bool tap_interp); @@ -2140,7 +2129,7 @@ T * create_R(const T & A, std::vector rowmap; rowmap.reserve(local_rows); - for (int i = 0; i < splitting.on_proc.size(); ++i) { + for (std::size_t i = 0; i < splitting.on_proc.size(); ++i) { if (splitting.on_proc[i] == Selected) { rowmap.push_back(A.local_row_map[i]); } @@ -2209,14 +2198,14 @@ auto fill_colind(std::size_t row, fpoint_distance distance, ParCSRMatrix & R) { // Note: uses global indices for off_proc columns - auto expand = [](const int row, + auto expand = [](const int i, const Matrix & soc, const auto & split, Matrix & rmat, int & ind, auto && colmap) { - for (int off = soc.idx1[row]; off < soc.idx1[row+1]; ++off) { + for (int off = soc.idx1[i]; off < soc.idx1[i + 1]; ++off) { auto d2point = soc.idx2[off]; - if (split[d2point] == Unselected && d2point != row) { + if (split[d2point] == Unselected && d2point != i) { rmat.idx2[ind++] = std::forward(colmap)(d2point); } } @@ -2418,7 +2407,7 @@ void fill_data(std::size_t row, int cpoint, int ind, int ind_off, for (int block_row = 0; block_row < blocksize; ++block_row) { auto row_maj_ind = (this_row + block_row) * num_dofs + this_col; for(int block_col = 0; block_col < blocksize; ++block_col) { - if ((row_maj_ind + block_col) > A0.size()) { + if (static_cast(row_maj_ind + block_col) > A0.size()) { std::cerr << "Warning block local_air fill_data: Accessing out of bounds index building A0.\n"; } A0[row_maj_ind + block_col] = vals[block_row * blocksize + block_col]; @@ -2485,11 +2474,11 @@ void fill_data(std::size_t row, int cpoint, int ind, int ind_off, // Add solution for each block row to data array. See section on RHS for // mapping between bsr data array and row-major array solution stored in - auto get_vals = [&](int block_ind) { + auto get_vals = [&](int block) { double * vals = new double[blocksize*blocksize]; for (int this_row = 0; this_row < blocksize; ++this_row) { for (int this_col = 0; this_col < blocksize; ++this_col) { - int row_ind = num_dofs * this_row + block_ind * blocksize + this_col; + int row_ind = num_dofs * this_row + block * blocksize + this_col; int bsr_ind = this_row * blocksize + this_col; if (std::abs(b0[row_ind]) > 1e-15) vals[bsr_ind] = b0[row_ind]; diff --git a/raptor/ruge_stuben/tests/test_air.cpp b/raptor/ruge_stuben/tests/test_air.cpp index 09ffcd56..ce5f9a4f 100644 --- a/raptor/ruge_stuben/tests/test_air.cpp +++ b/raptor/ruge_stuben/tests/test_air.cpp @@ -15,7 +15,6 @@ int main(int argc, char** argv) return ret; } -#if 0 TEST(TestOnePointInterp, TestsInRuge_Stuben) { int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); constexpr std::size_t n{16}; @@ -133,7 +132,7 @@ TEST(TestLocalAIR, TestsInRuge_Stuben) { auto expected = get_expected(row); expect_t cols; auto add_cols = [&](const Matrix & mat, const auto & colmap) { - for (std::size_t off = mat.idx1[i]; off < mat.idx1[i + 1]; ++off) { + for (int off = mat.idx1[i]; off < mat.idx1[i + 1]; ++off) { cols[colmap[mat.idx2[off]]] = mat.vals[off]; } }; @@ -143,66 +142,10 @@ TEST(TestLocalAIR, TestsInRuge_Stuben) { } } -TEST(TestLocalAIRBSR, TestsInRuge_Stuben) { - int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); - constexpr std::size_t n{16}; - // constexpr std::size_t n{5}; - std::vector grid; - - grid.resize(1, n); - - std::vector stencil{{-1., 2, -1}}; - - auto A = par_stencil_grid(stencil.data(), grid.data(), 1); - auto Absr = new raptor::ParBSRMatrix(A->partition, 3, 3); - [](auto & ... mats) { - ((mats.idx1[0] = 0),...); - }(*Absr->on_proc, *Absr->off_proc); - for (int i = 0; i < A->local_num_rows; ++i) { - auto copy = [=](const Matrix & src, Matrix & dst, auto && colmap) { - auto & bsr = dynamic_cast(dst); - for (int j = src.idx1[i]; j < src.idx1[i+1]; ++j) { - dst.idx2.push_back(colmap(src.idx2[j])); - bsr.block_vals.push_back(new double[bsr.b_size]()); - for (int k = 0; k < bsr.b_size; ++k) - bsr.block_vals.back()[k] = src.vals[j]; - } - dst.idx1[i+1] = dst.idx2.size(); - }; - copy(*A->on_proc, *Absr->on_proc, [](auto c) { return c; }); - copy(*A->off_proc, *Absr->off_proc, [&](auto c) { return A->off_proc_column_map[c]; }); - } - [](auto & ... mats) { - ((mats.nnz = mats.idx2.size()),...); - }(*Absr->on_proc, *Absr->off_proc); - Absr->finalize(); - - auto get_split = [](const Matrix &, const std::vector & colmap) { - std::vector split(colmap.size(), 0); - - std::size_t i{0}; - for (auto col : colmap) split[i++] = (col % 2 == 0) ? Selected : Unselected; - - return split; - }; - splitting_t split{get_split(*A->on_proc, A->on_proc_column_map), - get_split(*A->off_proc, A->off_proc_column_map)}; - - { - int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); - if (rank == 1) - split.on_proc[2] = Unselected; - } - - auto S = A->copy(); - - auto R = local_air(*Absr, *S, split); -} -#endif -void dump(const ParCSRMatrix & A) { +void dump(std::string prefix, const ParCSRMatrix & A) { int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); - std::ofstream ofile("coarse-" + std::to_string(rank)); + std::ofstream ofile(prefix + "-" + std::to_string(rank)); auto write_row = [&](int i, const Matrix & mat, const std::vector & colmap) { const auto & rowmap = A.local_row_map; for (int j = mat.idx1[i]; j < mat.idx1[i+1]; ++j) { @@ -258,17 +201,10 @@ ParCSRMatrix * gen(std::size_t n) { return A; } TEST(UpwindAdvection, TestsInRuge_Stuben) { - // constexpr std::size_t n{100}; - // std::vector grid; - // grid.resize(1, n); - // std::vector stencil{{-1, 1, 0}}; - // auto A = par_stencil_grid(stencil.data(), grid.data(), 1); auto A = gen(100); - // dump(*A); ParAIRSolver ml; - // ParRugeStubenSolver ml; ml.max_levels = 2; ml.setup(A); - dump(*ml.levels[1]->A); + // dump("pre", *ml.levels[1]->A); }