diff --git a/puremd_rc_1003/environ/param.gpu.water b/puremd_rc_1003/environ/param.gpu.water
index 061aad32e64e6d53543fa3ff073067352ce10b93..a1366285d5e1ba655f66d9aea163c5ff3e5dc18e 100644
--- a/puremd_rc_1003/environ/param.gpu.water
+++ b/puremd_rc_1003/environ/param.gpu.water
@@ -16,6 +16,7 @@ hbond_cutoff 7.50 ! cutoff distance for hydrogen bond interactions
q_err 1e-6 ! relative residual norm threshold used in GMRES
solver 2 ! method used to solve charge equilibration phase (QEq)
pre_comp 1 ! method used to compute QEq preconditioner, if applicable
+pre_j_iters 50 ! number of Jacobi iterations used for applying QEq precondition, if applicable
pre_refactor 100 ! nsteps to recompute preconditioner
pre_droptol 0.01 ! threshold tolerance for dropping values in preconditioner computation, if applicable
diff --git a/puremd_rc_1003/sPuReMD/GMRES.c b/puremd_rc_1003/sPuReMD/GMRES.c
index 2d4150ac86ddf5ff2f1ab3c9f2f6cb03c3592c14..a1d8c66e7fc1b47ca7a16c1c01e79eb4d6c899dc 100644
--- a/puremd_rc_1003/sPuReMD/GMRES.c
+++ b/puremd_rc_1003/sPuReMD/GMRES.c
@@ -883,7 +883,8 @@ int PGMRES( static_storage *workspace, sparse_matrix *H, real *b, real tol,
* with preconditioner using factors LU \approx H
* and Jacobi iteration for approximate factor application */
int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real tol,
- sparse_matrix *L, sparse_matrix *U, real *x, FILE *fout, real *time, real *spmv_time )
+ sparse_matrix *L, sparse_matrix *U, real *x, unsigned int iters,
+ FILE *fout, real *time, real *spmv_time )
{
int i, j, k, itr, N, si;
real cc, tmp1, tmp2, temp, bnorm;
@@ -926,10 +927,9 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
*spmv_time += (stop.tv_sec + stop.tv_usec / 1000000.0)
- (start.tv_sec + start.tv_usec / 1000000.0);
Vector_Sum( workspace->v[0], 1., b, -1., workspace->b_prm, N );
- // TODO: add parameters to config file
gettimeofday( &start, NULL );
- Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[0], workspace->v[0], 50 );
- Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[0], workspace->v[0], 50 );
+ Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[0], workspace->v[0], iters );
+ Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[0], workspace->v[0], iters );
gettimeofday( &stop, NULL );
*time += (stop.tv_sec + stop.tv_usec / 1000000.0)
- (start.tv_sec + start.tv_usec / 1000000.0);
@@ -946,10 +946,9 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
gettimeofday( &stop, NULL );
*spmv_time += (stop.tv_sec + stop.tv_usec / 1000000.0)
- (start.tv_sec + start.tv_usec / 1000000.0);
- // TODO: add parameters to config file
gettimeofday( &start, NULL );
- Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[j + 1], workspace->v[j + 1], 50 );
- Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[j + 1], workspace->v[j + 1], 50 );
+ Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[j + 1], workspace->v[j + 1], iters );
+ Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[j + 1], workspace->v[j + 1], iters );
gettimeofday( &stop, NULL );
*time += (stop.tv_sec + stop.tv_usec / 1000000.0)
- (start.tv_sec + start.tv_usec / 1000000.0);
diff --git a/puremd_rc_1003/sPuReMD/GMRES.h b/puremd_rc_1003/sPuReMD/GMRES.h
index 7961802963b0afa10862fb16bcbc8ad7bf876acb..391c4777b95131e55fbcfac8393638f06a3686de 100644
--- a/puremd_rc_1003/sPuReMD/GMRES.h
+++ b/puremd_rc_1003/sPuReMD/GMRES.h
@@ -42,7 +42,7 @@ int PGMRES( static_storage*, sparse_matrix*, real*, real,
sparse_matrix*, sparse_matrix*, real*, FILE*, real*, real* );
int PGMRES_Jacobi( static_storage*, sparse_matrix*, real*, real,
- sparse_matrix*, sparse_matrix*, real*, FILE*, real*, real* );
+ sparse_matrix*, sparse_matrix*, real*, unsigned int, FILE*, real*, real* );
int PCG( static_storage*, sparse_matrix*, real*, real,
sparse_matrix*, sparse_matrix*, real*, FILE* );
diff --git a/puremd_rc_1003/sPuReMD/QEq.c b/puremd_rc_1003/sPuReMD/QEq.c
index 04c9a0cc568c4c4dc59311a7acb87889c4564400..807715e8cbf1d5549959955e17e2ae7854aeb4d3 100644
--- a/puremd_rc_1003/sPuReMD/QEq.c
+++ b/puremd_rc_1003/sPuReMD/QEq.c
@@ -241,13 +241,16 @@ real ICHOLT( sparse_matrix *A, real *droptol,
* Edmond Chow and Aftab Patel
* Fine-Grained Parallel Incomplete LU Factorization
* SIAM J. Sci. Comp. */
-static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
+static real ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
sparse_matrix * const U_t, sparse_matrix * const U )
{
unsigned int i, j, k, pj, x = 0, y = 0, ei_x, ei_y;
real *D, *D_inv, sum;
sparse_matrix *DAD;
int *Utop;
+ struct timeval start, stop;
+
+ gettimeofday( &start, NULL );
if ( Allocate_Matrix( &DAD, A->n, A->m ) == 0 )
{
@@ -297,7 +300,7 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
{
/* for each nonzero */
#pragma omp parallel for schedule(guided) \
- default(none) shared(DAD) private(sum, ei_x, ei_y, k) firstprivate(x, y)
+ default(none) shared(DAD, stderr) private(sum, ei_x, ei_y, k) firstprivate(x, y)
for ( j = 0; j < A->start[A->n]; ++j )
{
sum = ZERO;
@@ -419,6 +422,10 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
free(D_inv);
free(D);
free(Utop);
+
+ gettimeofday( &stop, NULL );
+ return (stop.tv_sec + stop.tv_usec / 1000000.0)
+ - (start.tv_sec + start.tv_usec / 1000000.0);
}
@@ -437,24 +444,36 @@ void Init_MatVec( reax_system *system, control_params *control,
Sort_Matrix_Rows( workspace->H );
Sort_Matrix_Rows( workspace->H_sp );
//fprintf( stderr, "H matrix sorted\n" );
- Calculate_Droptol( workspace->H, workspace->droptol, control->droptol );
- //fprintf( stderr, "drop tolerances calculated\n" );
+ if ( control->pre_comp == ICHOLT_PC )
+ {
+ Calculate_Droptol( workspace->H, workspace->droptol, control->droptol );
+ //fprintf( stderr, "drop tolerances calculated\n" );
+ }
if ( workspace->L == NULL )
{
- fillin = Estimate_LU_Fill( workspace->H, workspace->droptol );
- if ( Allocate_Matrix( &(workspace->L), far_nbrs->n, fillin ) == 0 ||
- Allocate_Matrix( &(workspace->U), far_nbrs->n, fillin ) == 0 )
- {
- fprintf( stderr, "not enough memory for LU matrices. terminating.\n" );
- exit(INSUFFICIENT_SPACE);
- }
- /* factors have sparsity pattern as H */
-// if ( Allocate_Matrix( &(workspace->L), workspace->H->n, workspace->H->m ) == 0 ||
-// Allocate_Matrix( &(workspace->U), workspace->H->n, workspace->H->m ) == 0 )
-// {
-// fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
-// exit(INSUFFICIENT_SPACE);
-// }
+ switch ( control->pre_comp )
+ {
+ case ICHOLT_PC:
+ fillin = Estimate_LU_Fill( workspace->H, workspace->droptol );
+ if ( Allocate_Matrix( &(workspace->L), far_nbrs->n, fillin ) == 0 ||
+ Allocate_Matrix( &(workspace->U), far_nbrs->n, fillin ) == 0 )
+ {
+ fprintf( stderr, "not enough memory for LU matrices. terminating.\n" );
+ exit(INSUFFICIENT_SPACE);
+ }
+ break;
+ case ICHOL_PAR_PC:
+ /* ICHOL_PAR: factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), workspace->H->n, workspace->H->m ) == 0 ||
+ Allocate_Matrix( &(workspace->U), workspace->H->n, workspace->H->m ) == 0 )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit(INSUFFICIENT_SPACE);
+ }
+ break;
+ default:
+ break;
+ }
#if defined(DEBUG_FOCUS)
fprintf( stderr, "fillin = %d\n", fillin );
fprintf( stderr, "allocated memory: L = U = %ldMB\n",
@@ -462,10 +481,18 @@ void Init_MatVec( reax_system *system, control_params *control,
#endif
}
- data->timing.pre_comp += ICHOLT( workspace->H, workspace->droptol, workspace->L, workspace->U );
-// data->timing.pre_comp += ICHOLT( workspace->H_sp, workspace->droptol, workspace->L, workspace->U );
- // TODO: add parameters for sweeps to control file
-// ICHOL_PAR( workspace->H, 1, workspace->L, workspace->U );
+ switch ( control->pre_comp )
+ {
+ case ICHOLT_PC:
+ data->timing.pre_comp += ICHOLT( workspace->H, workspace->droptol, workspace->L, workspace->U );
+// data->timing.pre_comp += ICHOLT( workspace->H_sp, workspace->droptol, workspace->L, workspace->U );
+ break;
+ case ICHOL_PAR_PC:
+ data->timing.pre_comp += ICHOL_PAR( workspace->H, 1, workspace->L, workspace->U );
+ break;
+ default:
+ break;
+ }
// fprintf( stderr, "condest = %f\n", condest(workspace->L, workspace->U) );
@@ -554,41 +581,61 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
// if( data->step == 0 || data->step == 100 )
// Print_Linear_System( system, control, workspace, data->step );
- //TODO: add parameters in control file for solver choice and options
-// matvecs = GMRES( workspace, workspace->H,
-// workspace->b_s, control->q_err, workspace->s[0], out_control->log, &(data->timing.pre_app) );
-// matvecs += GMRES( workspace, workspace->H,
-// workspace->b_t, control->q_err, workspace->t[0], out_control->log, &(data->timing.pre_app) );
-
-// matvecs = GMRES_HouseHolder( workspace, workspace->H,
-// workspace->b_s, control->q_err, workspace->s[0], out_control->log );
-// matvecs += GMRES_HouseHolder( workspace, workspace->H,
-// workspace->b_t, control->q_err, workspace->t[0], out_control->log );
-
- matvecs = PGMRES( workspace, workspace->H, workspace->b_s, control->q_err,
- workspace->L, workspace->U, workspace->s[0], out_control->log, &(data->timing.pre_app) );
- matvecs += PGMRES( workspace, workspace->H, workspace->b_t, control->q_err,
- workspace->L, workspace->U, workspace->t[0], out_control->log, &(data->timing.pre_app) );
-
-// matvecs = PGMRES_Jacobi( workspace, workspace->H, workspace->b_s, control->q_err,
-// workspace->L, workspace->U, workspace->s[0], out_control->log, &(data->timing.pre_app) );
-// matvecs += PGMRES_Jacobi( workspace, workspace->H, workspace->b_t, control->q_err,
-// workspace->L, workspace->U, workspace->t[0], out_control->log, &(data->timing.pre_app) );
-
- //matvecs=PCG( workspace, workspace->H, workspace->b_s, control->q_err,
- // workspace->L, workspace->U, workspace->s[0], out_control->log ) + 1;
- ///matvecs+=PCG( workspace, workspace->H, workspace->b_t, control->q_err,
- // workspace->L, workspace->U, workspace->t[0], out_control->log ) + 1;
-
- //matvecs = CG( workspace, workspace->H,
- // workspace->b_s, control->q_err, workspace->s[0], out_control->log ) + 1;
- //matvecs += CG( workspace, workspace->H,
- // workspace->b_t, control->q_err, workspace->t[0], out_control->log ) + 1;
-
- //matvecs = SDM( workspace, workspace->H,
- // workspace->b_s, control->q_err, workspace->s[0], out_control->log ) + 1;
- //matvecs += SDM( workspace, workspace->H,
- // workspace->b_t, control->q_err, workspace->t[0], out_control->log ) + 1;
+ switch ( control->solver )
+ {
+ case GMRES_S:
+ matvecs = GMRES( workspace, workspace->H,
+ workspace->b_s, control->q_err, workspace->s[0], out_control->log,
+ &(data->timing.pre_app), &(data->timing.spmv) );
+ matvecs += GMRES( workspace, workspace->H,
+ workspace->b_t, control->q_err, workspace->t[0], out_control->log,
+ &(data->timing.pre_app), &(data->timing.spmv) );
+ break;
+ case GMRES_H_S:
+ matvecs = GMRES_HouseHolder( workspace, workspace->H,
+ workspace->b_s, control->q_err, workspace->s[0], out_control->log );
+ matvecs += GMRES_HouseHolder( workspace, workspace->H,
+ workspace->b_t, control->q_err, workspace->t[0], out_control->log );
+ break;
+ case PGMRES_S:
+ matvecs = PGMRES( workspace, workspace->H, workspace->b_s, control->q_err,
+ workspace->L, workspace->U, workspace->s[0], out_control->log,
+ &(data->timing.pre_app), &(data->timing.spmv) );
+ matvecs += PGMRES( workspace, workspace->H, workspace->b_t, control->q_err,
+ workspace->L, workspace->U, workspace->t[0], out_control->log,
+ &(data->timing.pre_app), &(data->timing.spmv) );
+ break;
+ case PGMRES_J_S:
+ matvecs = PGMRES_Jacobi( workspace, workspace->H, workspace->b_s, control->q_err,
+ workspace->L, workspace->U, workspace->s[0], control->jacobi_iters,
+ out_control->log, &(data->timing.pre_app), &(data->timing.spmv) );
+ matvecs += PGMRES_Jacobi( workspace, workspace->H, workspace->b_t, control->q_err,
+ workspace->L, workspace->U, workspace->t[0], control->jacobi_iters,
+ out_control->log, &(data->timing.pre_app), &(data->timing.spmv) );
+ break;
+ case CG_S:
+ matvecs = CG( workspace, workspace->H,
+ workspace->b_s, control->q_err, workspace->s[0], out_control->log ) + 1;
+ matvecs += CG( workspace, workspace->H,
+ workspace->b_t, control->q_err, workspace->t[0], out_control->log ) + 1;
+ break;
+ case PCG_S:
+ matvecs = PCG( workspace, workspace->H, workspace->b_s, control->q_err,
+ workspace->L, workspace->U, workspace->s[0], out_control->log ) + 1;
+ matvecs += PCG( workspace, workspace->H, workspace->b_t, control->q_err,
+ workspace->L, workspace->U, workspace->t[0], out_control->log ) + 1;
+ break;
+ case SDM_S:
+ matvecs = SDM( workspace, workspace->H,
+ workspace->b_s, control->q_err, workspace->s[0], out_control->log ) + 1;
+ matvecs += SDM( workspace, workspace->H,
+ workspace->b_t, control->q_err, workspace->t[0], out_control->log ) + 1;
+ break;
+ default:
+ fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
data->timing.matvecs += matvecs;
#if defined(DEBUG_FOCUS)
diff --git a/puremd_rc_1003/sPuReMD/mytypes.h b/puremd_rc_1003/sPuReMD/mytypes.h
index 32d95d6e27af5404e614172c81b12a36bcefdce8..3f98f74b51101c718135d0279d4fc08e3d816c3c 100644
--- a/puremd_rc_1003/sPuReMD/mytypes.h
+++ b/puremd_rc_1003/sPuReMD/mytypes.h
@@ -128,18 +128,37 @@ typedef real rvec[3];
typedef int ivec[3];
typedef real rtensor[3][3];
-enum {NVE = 0, NVT = 1, NPT = 2, sNPT = 3, iNPT = 4, ensNR = 5, bNVT = 6};
-enum {FAR_NBRS = 0, NEAR_NBRS = 1, THREE_BODIES = 2, BONDS = 3, OLD_BONDS = 4,
- HBONDS = 5, DBO = 6, DDELTA = 7, LIST_N = 8
- };
-enum {TYP_VOID = 0, TYP_THREE_BODY = 1, TYP_BOND = 2, TYP_HBOND = 3, TYP_DBO = 4,
- TYP_DDELTA = 5, TYP_FAR_NEIGHBOR = 6, TYP_NEAR_NEIGHBOR = 7, TYP_N = 8
- };
-enum {UNKNOWN = 0, WATER = 1};
-enum {NO_ANALYSIS = 0, FRAGMENTS = 1, REACTIONS = 2, NUM_ANALYSIS = 3};
-enum {WRITE_ASCII = 0, WRITE_BINARY = 1, RF_N = 2};
-enum {XYZ = 0, PDB = 1, BGF = 2, ASCII_RESTART = 3, BINARY_RESTART = 4, GF_N = 5};
-
+/* config params */
+enum ensemble {
+ NVE = 0, NVT = 1, NPT = 2, sNPT = 3, iNPT = 4, ensNR = 5, bNVT = 6
+};
+enum interaction_list_offets {
+ FAR_NBRS = 0, NEAR_NBRS = 1, THREE_BODIES = 2, BONDS = 3, OLD_BONDS = 4,
+ HBONDS = 5, DBO = 6, DDELTA = 7, LIST_N = 8
+};
+enum interaction_type {
+ TYP_VOID = 0, TYP_THREE_BODY = 1, TYP_BOND = 2, TYP_HBOND = 3, TYP_DBO = 4,
+ TYP_DDELTA = 5, TYP_FAR_NEIGHBOR = 6, TYP_NEAR_NEIGHBOR = 7, TYP_N = 8
+};
+enum molecule_type {
+ UNKNOWN = 0, WATER = 1
+};
+enum molecular_analysis_type {
+ NO_ANALYSIS = 0, FRAGMENTS = 1, REACTIONS = 2, NUM_ANALYSIS = 3
+};
+enum restart_format {
+ WRITE_ASCII = 0, WRITE_BINARY = 1, RF_N = 2
+};
+enum geometry {
+ XYZ = 0, PDB = 1, BGF = 2, ASCII_RESTART = 3, BINARY_RESTART = 4, GF_N = 5
+};
+enum solver {
+ GMRES_S = 0, GMRES_H_S = 1, PGMRES_S = 2,
+ PGMRES_J_S = 3, CG_S = 4, PCG_S = 5, SDM_S = 6,
+};
+enum pre_comp {
+ DIAG_PC = 0, ICHOLT_PC = 1, ICHOL_PAR_PC = 2,
+};
/* Global params mapping */
@@ -443,8 +462,12 @@ typedef struct
int freq_diffusion_coef;
int restrict_type;
+ unsigned int solver;
+ unsigned int pre_comp;
int refactor;
real droptol;
+ unsigned int jacobi_iters;
+
int molec_anal;
int freq_molec_anal;
@@ -751,6 +774,7 @@ typedef struct
} static_storage;
+/* interaction lists */
typedef struct
{
int n;
diff --git a/puremd_rc_1003/sPuReMD/param.c b/puremd_rc_1003/sPuReMD/param.c
index 89bcd550608f7466d54ab1d5c60135830fb248c6..406a86935441302cd318b4ccf24af04ede5dfa0d 100644
--- a/puremd_rc_1003/sPuReMD/param.c
+++ b/puremd_rc_1003/sPuReMD/param.c
@@ -839,7 +839,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
strcpy( control->sim_name, "default.sim" );
control->restart = 0;
- out_control->restart_format = 1;
+ out_control->restart_format = WRITE_BINARY;
out_control->restart_freq = 0;
strcpy( control->restart_from, "default.res" );
out_control->restart_freq = 0;
@@ -847,11 +847,11 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->reposition_atoms = 0;
- control->ensemble = 0;
+ control->ensemble = NVE;
control->nsteps = 0;
control->dt = 0.25;
- control->geo_format = 1;
+ control->geo_format = PDB;
control->restrict_bonds = 0;
control->periodic_boundaries = 1;
@@ -870,8 +870,11 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
control->q_err = 0.000001;
control->tabulate = 0;
+ control->solver = PGMRES_S;
+ control->pre_comp = ICHOLT_PC;
control->refactor = 100;
control->droptol = 0.01;
+ control->jacobi_iters = 0;
control->T_init = 0.;
control->T_final = 300.;
@@ -911,7 +914,7 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
out_control->bond_info = 0;
out_control->angle_info = 0;
- control->molec_anal = 0;
+ control->molec_anal = NO_ANALYSIS;
control->freq_molec_anal = 0;
control->bg_cut = 0.3;
control->num_ignored = 0;
@@ -1037,16 +1040,31 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
val = atof( tmp[1] );
control->q_err = val;
}
- else if ( strcmp(tmp[0], "ilu_refactor") == 0 )
+ else if ( strcmp(tmp[0], "solver") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->solver = ival;
+ }
+ else if ( strcmp(tmp[0], "pre_comp") == 0 )
+ {
+ ival = atoi( tmp[1] );
+ control->pre_comp = ival;
+ }
+ else if ( strcmp(tmp[0], "pre_refactor") == 0 )
{
ival = atoi( tmp[1] );
control->refactor = ival;
}
- else if ( strcmp(tmp[0], "ilu_droptol") == 0 )
+ else if ( strcmp(tmp[0], "pre_droptol") == 0 )
{
val = atof( tmp[1] );
control->droptol = val;
}
+ else if ( strcmp(tmp[0], "pre_j_iters") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->jacobi_iters = val;
+ }
else if ( strcmp(tmp[0], "temp_init") == 0 )
{
val = atof(tmp[1]);