Skip to content
Snippets Groups Projects
Commit eff2a773 authored by Kurt A. O'Hearn's avatar Kurt A. O'Hearn
Browse files

sPuReMD: refactor and clean up preconditioner application code (level... 

sPuReMD: refactor and clean up preconditioner application code (level scheduling, graph coloring, Jacobi iteration). Conditionally change charge matrix for setting up preconditioner for ICHOLT and ILU(T).
parent a466098d
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
sPuReMD/src/charges.c
+ 26
− 72
View file @ eff2a773
......@@ -771,10 +771,6 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
}
data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
#if defined(DEBUG)
fprintf( stderr, "H matrix sorted\n" );
#endif
switch ( control->cm_solver_pre_comp_type )
{
case NONE_PC:
......@@ -791,18 +787,8 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
case ICHOLT_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
fprintf( stderr, "fillin = %d\n", fillin );
fprintf( stderr, "allocated memory: L = U = %ldMB\n",
fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
if ( workspace->L == NULL )
{
if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
......@@ -839,10 +825,6 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
case ILUT_PAR_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
if ( workspace->L == NULL )
{
/* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
......@@ -919,11 +901,12 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
}
data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
#if defined(DEBUG)
fprintf( stderr, "H matrix sorted\n" );
#endif
if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
control->cm_solver_pre_comp_type == ILU_PAR_PC ||
control->cm_solver_pre_comp_type == ILUT_PAR_PC )
{
Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
}
switch ( control->cm_solver_pre_comp_type )
{
......@@ -941,18 +924,8 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
case ICHOLT_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
fprintf( stderr, "fillin = %d\n", fillin );
fprintf( stderr, "allocated memory: L = U = %ldMB\n",
fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
if ( workspace->L == NULL )
{
if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
......@@ -989,10 +962,6 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
case ILUT_PAR_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
if ( workspace->L == NULL )
{
/* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
......@@ -1038,7 +1007,12 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
control->cm_solver_pre_comp_type == ILU_PAR_PC ||
control->cm_solver_pre_comp_type == ILUT_PAR_PC )
{
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
}
}
......@@ -1071,12 +1045,13 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
}
data->timing.cm_sort_mat_rows += Get_Timing_Info( time );
Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
Hptr->val[Hptr->start[system->N_cm] - 1] = 1.0;
#if defined(DEBUG)
fprintf( stderr, "H matrix sorted\n" );
#endif
if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
control->cm_solver_pre_comp_type == ILU_PAR_PC ||
control->cm_solver_pre_comp_type == ILUT_PAR_PC )
{
Hptr->val[Hptr->start[system->N + 1] - 1] = 1.0;
Hptr->val[Hptr->start[system->N_cm] - 1] = 1.0;
}
switch ( control->cm_solver_pre_comp_type )
{
......@@ -1094,18 +1069,8 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
case ICHOLT_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
fprintf( stderr, "fillin = %d\n", fillin );
fprintf( stderr, "allocated memory: L = U = %ldMB\n",
fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
if ( workspace->L == NULL )
{
if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
......@@ -1141,10 +1106,6 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
case ILUT_PAR_PC:
Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
fprintf( stderr, "drop tolerances calculated\n" );
#endif
if ( workspace->L == NULL )
{
/* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
......@@ -1190,8 +1151,13 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
Hptr->val[Hptr->start[system->N_cm] - 1] = 0.0;
if ( control->cm_solver_pre_comp_type == ICHOLT_PC ||
control->cm_solver_pre_comp_type == ILU_PAR_PC ||
control->cm_solver_pre_comp_type == ILUT_PAR_PC )
{
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
Hptr->val[Hptr->start[system->N_cm] - 1] = 0.0;
}
}
......@@ -1315,10 +1281,6 @@ static void QEq( reax_system * const system, control_params * const control,
data->timing.cm_solver_iters += iters;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "linsolve-" );
#endif
Calculate_Charges_QEq( system, workspace );
#if defined(DEBUG_FOCUS)
......@@ -1395,10 +1357,6 @@ static void EE( reax_system * const system, control_params * const control,
data->timing.cm_solver_iters += iters;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "linsolve-" );
#endif
Calculate_Charges_EE( system, workspace );
// if( data->step == control->nsteps )
......@@ -1469,10 +1427,6 @@ static void ACKS2( reax_system * const system, control_params * const control,
data->timing.cm_solver_iters += iters;
#if defined(DEBUG_FOCUS)
fprintf( stderr, "linsolve-" );
#endif
Calculate_Charges_EE( system, workspace );
}
......
sPuReMD/src/init_md.c
+ 100
− 12
View file @ eff2a773
......@@ -501,27 +501,66 @@ void Init_Workspace( reax_system *system, control_params *control,
break;
}
/* SpMV related */
#ifdef _OPENMP
workspace->b_local = (real*) smalloc( control->num_threads * system->N_cm * sizeof(real),
"Init_Workspace::b_local" );
#endif
/* level scheduling related */
workspace->levels_L = 1;
workspace->levels_U = 1;
if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
workspace->row_levels_L = (unsigned int*) smalloc( system->N_cm * sizeof(unsigned int),
"Init_Workspace::row_levels_L" );
workspace->level_rows_L = (unsigned int*) smalloc( system->N_cm * sizeof(unsigned int),
"Init_Workspace::level_rows_L" );
workspace->level_rows_cnt_L = (unsigned int*) smalloc( (system->N_cm + 1) * sizeof(unsigned int),
"Init_Workspace::level_rows_cnt_L" );
workspace->row_levels_U = (unsigned int*) smalloc( system->N_cm * sizeof(unsigned int),
"Init_Workspace::row_levels_U" );
workspace->level_rows_U = (unsigned int*) smalloc( system->N_cm * sizeof(unsigned int),
"Init_Workspace::level_rows_U" );
workspace->level_rows_cnt_U = (unsigned int*) smalloc( (system->N_cm + 1) * sizeof(unsigned int),
"Init_Workspace::level_rows_cnt_U" );
workspace->top = (unsigned int*) smalloc( (system->N_cm + 1) * sizeof(unsigned int),
"Init_Workspace::top" );
}
else
{
workspace->row_levels_L = NULL;
workspace->level_rows_L = NULL;
workspace->level_rows_cnt_L = NULL;
workspace->row_levels_U = NULL;
workspace->level_rows_U = NULL;
workspace->level_rows_cnt_U = NULL;
workspace->top = NULL;
}
/* graph coloring related */
workspace->recolor_cnt = 0;
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
workspace->color = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::color" );
"Init_Workspace::color" );
workspace->to_color = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::to_color" );
"Init_Workspace::to_color" );
workspace->conflict = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::conflict" );
"setup_graph_coloring::conflict" );
workspace->conflict_cnt = (unsigned int*) smalloc( sizeof(unsigned int) * (control->num_threads + 1),
"setup_graph_coloring::conflict_cnt" );
"Init_Workspace::conflict_cnt" );
workspace->recolor = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::recolor" );
"Init_Workspace::recolor" );
workspace->color_top = (unsigned int*) smalloc( sizeof(unsigned int) * (system->N_cm + 1),
"setup_graph_coloring::color_top" );
"Init_Workspace::color_top" );
workspace->permuted_row_col = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::premuted_row_col" );
"Init_Workspace::premuted_row_col" );
workspace->permuted_row_col_inv = (unsigned int*) smalloc( sizeof(unsigned int) * system->N_cm,
"setup_graph_coloring::premuted_row_col_inv" );
"Init_Workspace::premuted_row_col_inv" );
workspace->y_p = (real*) smalloc( sizeof(real) * system->N_cm,
"setup_graph_coloring::y_p" );
"Init_Workspace::y_p" );
workspace->x_p = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::x_p" );
}
......@@ -531,9 +570,7 @@ void Init_Workspace( reax_system *system, control_params *control,
workspace->to_color = NULL;
workspace->conflict = NULL;
workspace->conflict_cnt = NULL;
workspace->temp_ptr = NULL;
workspace->recolor = NULL;
workspace->recolor_cnt = 0;
workspace->color_top = NULL;
workspace->permuted_row_col = NULL;
workspace->permuted_row_col_inv = NULL;
......@@ -541,6 +578,29 @@ void Init_Workspace( reax_system *system, control_params *control,
workspace->x_p = NULL;
}
/* Jacobi iteration related */
if ( control->cm_solver_pre_app_type == JACOBI_ITER_PA )
{
workspace->Dinv_L = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::Dinv_L" );
workspace->Dinv_U = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::Dinv_U" );
workspace->Dinv_b = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::Dinv_b" );
workspace->rp = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::rp" );
workspace->rp2 = (real*) smalloc( sizeof(real) * system->N_cm,
"Init_Workspace::rp2" );
}
else
{
workspace->Dinv_L = NULL;
workspace->Dinv_U = NULL;
workspace->Dinv_b = NULL;
workspace->rp = NULL;
workspace->rp2 = NULL;
}
/* integrator storage */
workspace->a = (rvec *) smalloc( system->N * sizeof( rvec ),
"Init_Workspace::workspace->a" );
......@@ -993,7 +1053,7 @@ void Initialize( reax_system *system, control_params *control,
#endif
#ifdef _OPENMP
#pragma omp parallel default(shared)
#pragma omp parallel default(none) shared(control)
{
#pragma omp single
control->num_threads = omp_get_num_threads( );
......@@ -1211,6 +1271,24 @@ void Finalize_Workspace( reax_system *system, control_params *control,
break;
}
/* SpMV related */
#ifdef _OPENMP
sfree( workspace->b_local, "Finalize_Workspace::b_local" );
#endif
/* level scheduling related */
if ( control->cm_solver_pre_app_type == TRI_SOLVE_LEVEL_SCHED_PA ||
control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
sfree( workspace->row_levels_L, "Finalize_Workspace::row_levels_L" );
sfree( workspace->level_rows_L, "Finalize_Workspace::level_rows_L" );
sfree( workspace->level_rows_cnt_L, "Finalize_Workspace::level_rows_cnt_L" );
sfree( workspace->row_levels_U, "Finalize_Workspace::row_levels_U" );
sfree( workspace->level_rows_U, "Finalize_Workspace::level_rows_U" );
sfree( workspace->level_rows_cnt_U, "Finalize_Workspace::level_rows_cnt_U" );
sfree( workspace->top, "Finalize_Workspace::top" );
}
/* graph coloring related */
if ( control->cm_solver_pre_app_type == TRI_SOLVE_GC_PA )
{
......@@ -1226,6 +1304,16 @@ void Finalize_Workspace( reax_system *system, control_params *control,
sfree( workspace->x_p, "Finalize_Workspace::workspace->x_p" );
}
/* Jacobi iteration related */
if ( control->cm_solver_pre_app_type == JACOBI_ITER_PA )
{
sfree( workspace->Dinv_L, "Finalize_Workspace::Dinv_L" );
sfree( workspace->Dinv_U, "Finalize_Workspace::Dinv_U" );
sfree( workspace->Dinv_b, "Finalize_Workspace::Dinv_b" );
sfree( workspace->rp, "Finalize_Workspace::rp" );
sfree( workspace->rp2, "Finalize_Workspace::rp2" );
}
/* integrator storage */
sfree( workspace->a, "Finalize_Workspace::workspace->a" );
sfree( workspace->f_old, "Finalize_Workspace::workspace->f_old" );
......
sPuReMD/src/lin_alg.c
+ 137
− 215
View file @ eff2a773
This diff is collapsed.
sPuReMD/src/lin_alg.h
+ 5
− 3
View file @ eff2a773
......@@ -73,16 +73,18 @@ void Transpose_I( sparse_matrix * const );
void tri_solve( const sparse_matrix * const, const real * const,
real * const, const int, const TRIANGULARITY );
void tri_solve_level_sched( const sparse_matrix * const,
void tri_solve_level_sched( static_storage *,
const sparse_matrix * const,
const real * const, real * const, const int,
const TRIANGULARITY, int );
void jacobi_iter( const sparse_matrix * const, const real * const,
void jacobi_iter( const static_storage * const,
const sparse_matrix * const, const real * const,
const real * const, real * const, const TRIANGULARITY,
const unsigned int );
void setup_graph_coloring( const control_params * const,
static_storage * const, const sparse_matrix * const,
const static_storage * const, const sparse_matrix * const,
sparse_matrix **, sparse_matrix ** );
int GMRES( const static_storage * const, const control_params * const,
......
sPuReMD/src/mytypes.h
+ 25
− 1
View file @ eff2a773
......@@ -977,13 +977,29 @@ typedef struct
real *q;
real *p;
/* SpMV related storage */
#ifdef _OPENMP
real *b_local;
#endif
/* Level scheduling related storage for applying, e.g. ICHOLT and ILU(T),
* preconditioners */
int levels_L;
int levels_U;
unsigned int *row_levels_L;
unsigned int *level_rows_L;
unsigned int *level_rows_cnt_L;
unsigned int *row_levels_U;
unsigned int *level_rows_U;
unsigned int *level_rows_cnt_U;
unsigned int *top;
/* Graph coloring related storage for applying, e.g. ICHOLT and ILU(T),
* preconditioners */
unsigned int *color;
unsigned int *to_color;
unsigned int *conflict;
unsigned int *conflict_cnt;
unsigned int *temp_ptr;
unsigned int *recolor;
unsigned int recolor_cnt;
unsigned int *color_top;
......@@ -992,6 +1008,14 @@ typedef struct
real *y_p;
real *x_p;
/* Jacobi iteration related storage for applying, e.g. ICHOLT and ILU(T),
* preconditioners */
real *Dinv_L;
real *Dinv_U;
real *Dinv_b;
real *rp;
real *rp2;
int num_H;
int *hbond_index; // for hydrogen bonds
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment