Skip to content

Commit

Permalink
Cleanup.
Browse files Browse the repository at this point in the history
  • Loading branch information
@andrewreisner
andrewreisner committed Sep 12, 2024
1 parent 8972fa6 commit 7c76e36
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 96 deletions.
20 changes: 10 additions & 10 deletions raptor/ruge_stuben/par_air_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
25 changes: 7 additions & 18 deletions raptor/ruge_stuben/par_interpolation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,17 +20,6 @@ CSRMatrix* communicate(ParCSRMatrix* A, ParCSRMatrix* S, const std::vector<int>&
CSRMatrix* communicate(ParCSRMatrix* A, const std::vector<int>& states,
const std::vector<int>& off_proc_states, CommPkg* comm);
void filter_interp(ParCSRMatrix* P, const double filter_threshold);
ParCSRMatrix* extended_interpolation(ParCSRMatrix* A,
ParCSRMatrix* S, const std::vector<int>& states,
const std::vector<int>& 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<int>& states,
const std::vector<int>& off_proc_states,
bool tap_interp, int num_variables, int* variables);
ParCSRMatrix* direct_interpolation(ParCSRMatrix* A,
ParCSRMatrix* S, const std::vector<int>& states,
const std::vector<int>& off_proc_states, bool tap_interp);



Expand Down Expand Up @@ -2140,7 +2129,7 @@ T * create_R(const T & A,
std::vector<int> 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]);
}
Expand Down Expand Up @@ -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<decltype(colmap)>(colmap)(d2point);
}
}
Expand Down Expand Up @@ -2418,7 +2407,7 @@ void fill_data<ParBSRMatrix>(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<std::size_t>(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];
Expand Down Expand Up @@ -2485,11 +2474,11 @@ void fill_data<ParBSRMatrix>(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];
Expand Down
72 changes: 4 additions & 68 deletions raptor/ruge_stuben/tests/test_air.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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];
}
};
Expand All @@ -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<int> grid;

grid.resize(1, n);

std::vector<double> 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<BSRMatrix &>(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<int> & colmap) {
std::vector<int> 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<int> & colmap) {
const auto & rowmap = A.local_row_map;
for (int j = mat.idx1[i]; j < mat.idx1[i+1]; ++j) {
Expand Down Expand Up @@ -258,17 +201,10 @@ ParCSRMatrix * gen(std::size_t n) {
return A;
}
TEST(UpwindAdvection, TestsInRuge_Stuben) {
// constexpr std::size_t n{100};
// std::vector<int> grid;
// grid.resize(1, n);
// std::vector<double> 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);
}

0 comments on commit 7c76e36

Please sign in to comment.