Skip to content

Commit

Permalink
fixed small bug in classical interpolation
Browse files Browse the repository at this point in the history
  • Loading branch information
Amanda Bienz committed Oct 29, 2018
1 parent 4094260 commit 0a99e70
Show file tree
Hide file tree
Showing 6 changed files with 272 additions and 176 deletions.
9 changes: 7 additions & 2 deletions examples/benchmark_tap_amg.cpp
View file
Expand Up @@ -147,7 +147,7 @@ int main(int argc, char* argv[])
A = readParMatrix(file);
}

if (system != 2)
/* if (system != 2)
{
A->tap_comm = new TAPComm(A->partition, A->off_proc_column_map,
A->on_proc_column_map);
Expand Down Expand Up @@ -326,11 +326,16 @@ int main(int argc, char* argv[])
}
delete ml;
delete A;
HYPRE_IJMatrixDestroy(A_h_ij);
HYPRE_IJVectorDestroy(x_h_ij);
HYPRE_IJVectorDestroy(b_h_ij);
*/
printf("%d\n", A->global_num_rows);
// write_par_mm(A, "grad_div.mtx");
delete A;


MPI_Finalize();
return 0;
}
Expand Down
18 changes: 13 additions & 5 deletions raptor/gallery/external/hypre_wrapper.cpp
View file
Expand Up @@ -57,6 +57,14 @@ HYPRE_IJMatrix convert(raptor::ParCSRMatrix* A_rap, MPI_Comm comm_mat)
local_col_start += col_sizes[i];
}

// Send new rows for off_proc_cols
aligned_vector<int> rows(A_rap->local_num_rows);
for (int i = 0; i < A_rap->local_num_rows; i++)
{
rows[i] = i + local_row_start;
}
aligned_vector<int>& off_rows = A_rap->comm->communicate(rows);

n_rows = A_rap->local_num_rows;
if (n_rows)
{
Expand All @@ -78,10 +86,10 @@ HYPRE_IJMatrix convert(raptor::ParCSRMatrix* A_rap, MPI_Comm comm_mat)
{
row_start = A_rap->on_proc->idx1[i];
row_end = A_rap->on_proc->idx1[i+1];
global_row = A_rap->local_row_map[i];
global_row = i + local_row_start;
for (int j = row_start; j < row_end; j++)
{
global_col = A_rap->on_proc_column_map[A_rap->on_proc->idx2[j]];
global_col = A_rap->on_proc->idx2[j] + local_col_start;
value = A_rap->on_proc->vals[j];
HYPRE_IJMatrixSetValues(A, 1, &one, &global_row, &global_col, &value);
}
Expand All @@ -90,10 +98,10 @@ HYPRE_IJMatrix convert(raptor::ParCSRMatrix* A_rap, MPI_Comm comm_mat)
{
row_start = A_rap->off_proc->idx1[i];
row_end = A_rap->off_proc->idx1[i+1];
global_row = A_rap->local_row_map[i];
global_row = i + local_row_start;
for (int j = row_start; j < row_end; j++)
{
global_col = A_rap->off_proc_column_map[A_rap->off_proc->idx2[j]];
global_col = off_rows[A_rap->off_proc->idx2[j]];
value = A_rap->off_proc->vals[j];
HYPRE_IJMatrixSetValues(A, 1, &one, &global_row, &global_col, &value);
}
Expand Down Expand Up @@ -202,7 +210,7 @@ HYPRE_Solver hypre_create_hierarchy(hypre_ParCSRMatrix* A,
HYPRE_BoomerAMGSetRelaxOrder(amg_data, 0);
// HYPRE_BoomerAMGSetRelaxWt(amg_data, 3.0/4); // set omega for SOR
HYPRE_BoomerAMGSetNumFunctions(amg_data, num_functions);
HYPRE_BoomerAMGSetRAP2(amg_data, 1);
//HYPRE_BoomerAMGSetRAP2(amg_data, 1);

HYPRE_BoomerAMGSetPrintLevel(amg_data, 0);
HYPRE_BoomerAMGSetMaxIter(amg_data, 1000);
Expand Down
2 changes: 2 additions & 0 deletions raptor/raptor.hpp
View file
Expand Up @@ -38,8 +38,10 @@

// Matrix IO
#include "gallery/matrix_IO.hpp"
#include "gallery/matrix_market.hpp"
#ifndef NO_MPI
#include "gallery/par_matrix_IO.hpp"
#include "gallery/par_matrix_market.hpp"
#endif

// Matrix external gallery
Expand Down
17 changes: 13 additions & 4 deletions raptor/ruge_stuben/par_interpolation.cpp
View file
Expand Up @@ -904,6 +904,9 @@ ParCSRMatrix* mod_classical_interpolation(ParCSRMatrix* A,
bool tap_interp, int num_variables, int* variables,
data_t* comm_t, data_t* comm_mat_t)
{
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);

int start, end;
int start_k, end_k;
int end_S;
Expand Down Expand Up @@ -1215,9 +1218,12 @@ ParCSRMatrix* mod_classical_interpolation(ParCSRMatrix* A,
}
}
}
else if (num_variables == 1 || variables[i] == variables[col])
else if (states[col] != NoNeighbors)
{
weak_sum += val;
if (num_variables == 1 || variables[i] == variables[col])
{
weak_sum += val;
}
}
}

Expand Down Expand Up @@ -1294,9 +1300,12 @@ ParCSRMatrix* mod_classical_interpolation(ParCSRMatrix* A,
}
}
}
else if (num_variables == 1 || variables[i] == off_variables[col])
else if (off_proc_states[col] != NoNeighbors)
{
weak_sum += val;
if (num_variables == 1 || variables[i] == off_variables[col])
{
weak_sum += val;
}
}
}

Expand Down
179 changes: 95 additions & 84 deletions raptor/ruge_stuben/tests/test_hypre.cpp
View file
Expand Up @@ -4,13 +4,7 @@

#include "gtest/gtest.h"
#include "mpi.h"
#include "core/types.hpp"
#include "core/par_matrix.hpp"
#include "gallery/par_stencil.hpp"
#include "gallery/laplacian27pt.hpp"
#include "gallery/external/hypre_wrapper.hpp"
#include "ruge_stuben/par_cf_splitting.hpp"
#include "ruge_stuben/par_interpolation.hpp"
#include "raptor.hpp"
#include "tests/hypre_compare.hpp"
#include <iostream>
#include <fstream>
Expand Down Expand Up @@ -44,108 +38,125 @@ int main(int argc, char** argv)
return temp;
} // end of main() //

TEST(TestHypreAgg, TestsInRuge_Stuben)
TEST(TestHypre, TestsInRuge_Stuben)
{
int rank, num_procs;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

FILE* f;
double* weights;
int cf;

int n = 25;
aligned_vector<int> grid(3, n);
double* stencil = laplace_stencil_27pt();

std::vector<ParCSRMatrix*> A_array;
std::vector<ParCSRMatrix*> P_array;

ParCSRMatrix* A = par_stencil_grid(stencil, grid.data(), 3);
A_array.push_back(A);
delete[] stencil;
form_hypre_weights(&weights, A->local_num_rows);

double strong_threshold = 0.25;
int hyp_coarsen_type = 0; // CLJP
int hyp_interp_type = 0; // ModClassical
int p_max_elmts = 0;
int agg_num_levels = 0;
ParCSRMatrix* A;
ParVector x, b;
ParMultilevel* ml;
HYPRE_IJMatrix Aij;
hypre_ParCSRMatrix* A_hyp;
int n;
aligned_vector<int> grid;
double* stencil;


/************************************
**** Test Laplacian
***********************************/
n = 25;
grid.resize(3);
std::fill(grid.begin(), grid.end(), n);
stencil = laplace_stencil_27pt();
A = par_stencil_grid(stencil, grid.data(), 3);
x = ParVector(A->global_num_rows, A->local_num_rows, A->partition->first_local_row);
b = ParVector(A->global_num_rows, A->local_num_rows, A->partition->first_local_row);
delete[] stencil;

ml = new ParRugeStubenSolver(strong_threshold, CLJP, ModClassical, Classical, SOR);
form_hypre_weights(&ml->weights, A->local_num_rows);
ml->setup(A);

Aij = convert(A);
HYPRE_IJMatrixGetObject(Aij, (void**) &A_hyp);
compare(A, A_hyp);

hypre_ParCSRMatrix* S_hyp;
hypre_ParCSRMatrix* P_hyp;
hypre_ParCSRMatrix* Ac_hyp;
int* states_hypre;
int* coarse_dof_func;
int* coarse_pnts_gbl;
// Convert vectors... needed for Hypre Setup
hypre_ParVector* x_hyp;
hypre_ParVector* b_hyp;
HYPRE_IJVector x_h_ij = convert(x);
HYPRE_IJVector b_h_ij = convert(b);
HYPRE_IJVectorGetObject(x_h_ij, (void **) &x_hyp);
HYPRE_IJVectorGetObject(b_h_ij, (void **) &b_hyp);

aligned_vector<int> states;
aligned_vector<int> off_proc_states;
// Setup Hypre Hierarchy
HYPRE_Solver solver_data = hypre_create_hierarchy(A_hyp, x_hyp, b_hyp,
hyp_coarsen_type, hyp_interp_type, p_max_elmts, agg_num_levels,
strong_threshold);

int nrows = A_array[0]->global_num_rows;
int level = 0;
while (nrows > 50)
hypre_ParCSRMatrix** A_array = hypre_ParAMGDataAArray((hypre_ParAMGData*) solver_data);
hypre_ParCSRMatrix** P_array = hypre_ParAMGDataPArray((hypre_ParAMGData*) solver_data);

for (int level = 0; level < ml->num_levels - 1; level++)
{
ParCSRMatrix* Al = A_array[level];

// Create Strength of Connection Matrix
ParCSRMatrix* Sl = Al->strength(Classical, 0.25);
hypre_BoomerAMGCreateS(A_hyp, 0.25, 1.0, 1, NULL, &S_hyp);
compareS(Sl, S_hyp);

// C/F Splitting (PMIS)
split_cljp(Sl, states, off_proc_states, false, weights);
hypre_BoomerAMGCoarsen(S_hyp, A_hyp, 0, 0, &states_hypre);
compare_states(Al->local_num_rows, states, states_hypre);

// Extended Interpolation
ParCSRMatrix* Pl = mod_classical_interpolation(Al, Sl, states, off_proc_states, false);
P_array.push_back(Pl);
hypre_BoomerAMGCoarseParms(MPI_COMM_WORLD, Al->local_num_rows, 1, NULL, states_hypre,
&coarse_dof_func, &coarse_pnts_gbl);
hypre_BoomerAMGBuildInterp(A_hyp, states_hypre, S_hyp, coarse_pnts_gbl, 1, NULL,
0, 0.0, 0.0, NULL, &P_hyp);
compare(Pl, P_hyp);

// SpGEMM (Form Ac)
ParCSRMatrix* APl = Al->mult(Pl);
ParCSCMatrix* Pcsc = Pl->to_ParCSC();
APl->comm = new ParComm(APl->partition, APl->off_proc_column_map, APl->on_proc_column_map);
ParCSRMatrix* Ac = APl->mult_T(Pcsc);
Ac->comm = new ParComm(Ac->partition, Ac->off_proc_column_map, Ac->on_proc_column_map);
A_array.push_back(Ac);
hypre_BoomerAMGBuildCoarseOperator(P_hyp, A_hyp, P_hyp, &Ac_hyp);
compare(Ac, Ac_hyp);

if (level > 0)
hypre_ParCSRMatrixDestroy(A_hyp);
A_hyp = Ac_hyp;

nrows = Ac->global_num_rows;
level++;

hypre_TFree(states_hypre, HYPRE_MEMORY_HOST);
hypre_ParCSRMatrixDestroy(P_hyp);
hypre_ParCSRMatrixDestroy(S_hyp);

delete APl;
delete Pcsc;
delete Sl;
compare(ml->levels[level]->P, P_array[level]);
compare(ml->levels[level+1]->A, A_array[level+1]);
}
hypre_ParCSRMatrixDestroy(Ac_hyp);

for (std::vector<ParCSRMatrix*>::iterator it = A_array.begin(); it != A_array.end(); ++it)
delete *it;
hypre_BoomerAMGDestroy(solver_data);
HYPRE_IJMatrixDestroy(Aij);
delete ml;
delete A;


/************************************
**** Test Anisotropic Diffusion
***********************************/
n = 100;
grid.resize(2);
std::fill(grid.begin(), grid.end(), n);
stencil = diffusion_stencil_2d(0.001, M_PI/4.0);
A = par_stencil_grid(stencil, grid.data(), 2);
x = ParVector(A->global_num_rows, A->local_num_rows, A->partition->first_local_row);
b = ParVector(A->global_num_rows, A->local_num_rows, A->partition->first_local_row);
delete[] stencil;

ml = new ParRugeStubenSolver(strong_threshold, CLJP, ModClassical, Classical, SOR);
form_hypre_weights(&ml->weights, A->local_num_rows);
ml->setup(A);

Aij = convert(A);
HYPRE_IJMatrixGetObject(Aij, (void**) &A_hyp);
compare(A, A_hyp);

// Convert vectors... needed for Hypre Setup
x_h_ij = convert(x);
b_h_ij = convert(b);
HYPRE_IJVectorGetObject(x_h_ij, (void **) &x_hyp);
HYPRE_IJVectorGetObject(b_h_ij, (void **) &b_hyp);

// Setup Hypre Hierarchy
solver_data = hypre_create_hierarchy(A_hyp, x_hyp, b_hyp,
hyp_coarsen_type, hyp_interp_type, p_max_elmts, agg_num_levels,
strong_threshold);

for (std::vector<ParCSRMatrix*>::iterator it = P_array.begin(); it != P_array.end(); ++it)
delete *it;
A_array = hypre_ParAMGDataAArray((hypre_ParAMGData*) solver_data);
P_array = hypre_ParAMGDataPArray((hypre_ParAMGData*) solver_data);

delete[] weights;
for (int level = 0; level < ml->num_levels - 1; level++)
{
compare(ml->levels[level]->P, P_array[level]);
compare(ml->levels[level+1]->A, A_array[level+1]);
}

hypre_BoomerAMGDestroy(solver_data);
HYPRE_IJMatrixDestroy(Aij);
delete ml;
delete A;


} // end of TEST(TestParSplitting, TestsInRuge_Stuben) //
} // end of TEST(TestHypre, TestsInRuge_Stuben) //



Expand Down

0 comments on commit 0a99e70

Please sign in to comment.