PuReMD: cleand code and add timing performance

PuReMD/src/control.c

@@ -74,6 +74,7 @@ char Read_Control_File( char *control_file, control_params* control,
     control->tabulate = 0;
 
     control->qeq_freq = 1;
+    control->cm_solver_max_iters = 100;
     control->q_err = 1e-6;
     control->sai_thres = 0.1;
     control->refactor = 100;
@@ -246,11 +247,21 @@ char Read_Control_File( char *control_file, control_params* control,
             ival = atoi( tmp[1] );
             control->qeq_freq = ival;
         }
+        else if ( strcmp(tmp[0], "cm_solver_max_iters") == 0 )
+        {
+            ival = atoi( tmp[1] );
+            control->cm_solver_max_iters = ival;
+        }
         else if ( strcmp(tmp[0], "q_err") == 0 )
         {
             val = atof( tmp[1] );
             control->q_err = val;
         }
+        else if ( strcmp(tmp[0], "cm_solver_pre_comp_type") == 0 )
+        {
+            val = atoi( tmp[1] );
+            control->cm_solver_pre_comp_type = val;
+        }
         else if ( strcmp(tmp[0], "sai_thres") == 0 )
         {
             val = atof( tmp[1] );

PuReMD/src/forces.c

@@ -968,7 +968,7 @@ void Compute_Forces( reax_system *system, control_params *control,
     MPI_Comm comm;
     int qeq_flag;
 #if defined(LOG_PERFORMANCE)
-    real t_start = 0;
+    real t_start = 0, t_elapsed;
 
     //MPI_Barrier( mpi_data->world );
     if ( system->my_rank == MASTER_NODE )
@@ -1005,7 +1005,9 @@ void Compute_Forces( reax_system *system, control_params *control,
 #if defined(LOG_PERFORMANCE)
     //MPI_Barrier( mpi_data->world );
     if ( system->my_rank == MASTER_NODE )
-        Update_Timing_Info( &t_start, &(data->timing.bonded) );
+    {
+        t_start = Get_Time( );
+    }
 #endif
 #if defined(DEBUG_FOCUS)
     fprintf( stderr, "p%d @ step%d: completed bonded\n",
@@ -1022,7 +1024,10 @@ void Compute_Forces( reax_system *system, control_params *control,
 #if defined(LOG_PERFORMANCE)
     //MPI_Barrier( mpi_data->world );
     if ( system->my_rank == MASTER_NODE )
-        Update_Timing_Info( &t_start, &(data->timing.qEq) );
+    {
+        t_elapsed = Get_Timing_Info( t_start );
+        data->timing.cm += t_elapsed;
+    }
 #endif
 #if defined(DEBUG_FOCUS)
     fprintf(stderr, "p%d @ step%d: qeq completed\n", system->my_rank, data->step);

PuReMD/src/linear_solvers.h

@@ -24,11 +24,11 @@
 
 #include "reax_types.h"
 
-void setup_sparse_approx_inverse( reax_system*, storage*, mpi_datatypes*, 
+real setup_sparse_approx_inverse( reax_system*, simulation_data*, storage*, mpi_datatypes*, 
         sparse_matrix *, sparse_matrix **, int, double );
 
-void sparse_approx_inverse( reax_system*, storage*, mpi_datatypes*, 
-        sparse_matrix*, sparse_matrix*, sparse_matrix** );
+real sparse_approx_inverse( reax_system*, simulation_data*, storage*, mpi_datatypes*, 
+        sparse_matrix*, sparse_matrix*, sparse_matrix**, int );
 
 int GMRES( reax_system*, storage*, sparse_matrix*,
            real*, real, real*, mpi_datatypes*, FILE* );
@@ -36,7 +36,7 @@ int GMRES_HouseHolder( reax_system*, storage*, sparse_matrix*,
                        real*, real, real*, mpi_datatypes*, FILE* );
 int dual_CG( reax_system*, storage*, sparse_matrix*,
              rvec2*, real, rvec2*, mpi_datatypes*, FILE* );
-int CG( reax_system*, storage*, sparse_matrix*,
+int CG( reax_system*, control_params*, simulation_data*, storage*, sparse_matrix*,
         real*, real, real*, mpi_datatypes*, FILE*, int );
 int PCG( reax_system*, storage*, sparse_matrix*, real*, real,
          sparse_matrix*, sparse_matrix*, real*, mpi_datatypes*, FILE* );

PuReMD/src/parallelreax.c

@@ -211,7 +211,7 @@ int main( int argc, char* argv[] )
 #endif
 
     /* start the simulation */
-    int total_itr = data->timing.s_matvecs + data->timing.t_matvecs;
+    int total_itr = data->timing.cm_solver_iters;
     for ( ++data->step; data->step <= control->nsteps; data->step++ )
     {
         if ( control->T_mode )
@@ -221,7 +221,7 @@ int main( int argc, char* argv[] )
         Post_Evolve(system, control, data, workspace, lists, out_control, mpi_data);
         if( system->my_rank == MASTER_NODE )
         {
-            total_itr += data->timing.s_matvecs + data->timing.t_matvecs;
+            total_itr += data->timing.cm_solver_iters;
         }
         Output_Results( system, control, data, lists, out_control, mpi_data );
         //Analysis(system, control, data, workspace, lists, out_control, mpi_data);

PuReMD/src/qEq.c

@@ -368,20 +368,41 @@ void Calculate_Charges( reax_system *system, storage *workspace,
 static void Setup_Preconditioner_QEq( reax_system *system, control_params *control, 
         simulation_data *data, storage *workspace, mpi_datatypes *mpi_data )
 {
+    real time, t_sort, t_pc, total_sort, total_pc;
+
     /* sort H needed for SpMV's in linear solver, H or H_sp needed for preconditioning */
+    time = Get_Time( );
     Sort_Matrix_Rows( workspace->H );
+    t_sort = Get_Timing_Info( time );
 
-    //TODO: add sai filtering value, which will be passed as the last parameter
-    setup_sparse_approx_inverse( system, workspace, mpi_data, workspace->H, &workspace->H_spar_patt, 
+    t_pc = setup_sparse_approx_inverse( system, data, workspace, mpi_data, workspace->H, &workspace->H_spar_patt, 
             control->nprocs, control->sai_thres );
+
+
+    MPI_Reduce(&t_sort, &total_sort, 1, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world);
+    MPI_Reduce(&t_pc, &total_pc, 1, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world);
+
+    if( system->my_rank == MASTER_NODE )
+    {
+        data->timing.cm_sort_mat_rows += total_sort / control->nprocs;
+        data->timing.cm_solver_pre_comp += total_pc / control->nprocs;
+    }
 }
 
 static void Compute_Preconditioner_QEq( reax_system *system, control_params *control, 
         simulation_data *data, storage *workspace, mpi_datatypes *mpi_data )
 {
+    real t_pc, total_pc;
 #if defined(HAVE_LAPACKE) || defined(HAVE_LAPACKE_MKL)
-    sparse_approx_inverse( system, workspace, mpi_data,
-            workspace->H, workspace->H_spar_patt, &workspace->H_app_inv );
+    t_pc = sparse_approx_inverse( system, data, workspace, mpi_data,
+            workspace->H, workspace->H_spar_patt, &workspace->H_app_inv, control->nprocs );
+
+    MPI_Reduce(&t_pc, &total_pc, 1, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world);
+
+    if( system->my_rank == MASTER_NODE )
+    {
+        data->timing.cm_solver_pre_comp += total_pc / control->nprocs;
+    }
 #else
     fprintf( stderr, "[ERROR] LAPACKE support disabled. Re-compile before enabling. Terminating...\n" );
     exit( INVALID_INPUT );
@@ -392,7 +413,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
         storage *workspace, output_controls *out_control,
         mpi_datatypes *mpi_data )
 {
-    int j, s_matvecs, t_matvecs;
+    int j, iters;
 
     Init_MatVec( system, data, control, workspace, mpi_data );
 
@@ -401,44 +422,35 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
     //Print_Linear_System( system, control, workspace, data->step );
 #endif
 
-    //s_matvecs = dual_CG(system, workspace, workspace->H, workspace->b,
-    //          control->q_err, workspace->x, mpi_data, out_control->log);
-    //t_matvecs = 0;
-
-#if defined(SAI_PRECONDITIONER) && !defined(HALF_LIST)
-    if( control->refactor > 0 && ((data->step - data->prev_steps) % control->refactor == 0))
+    if( control->cm_solver_pre_comp_type == SAI_PC )
     {
-        Setup_Preconditioner_QEq( system, control, data, workspace, mpi_data );
+        if( control->refactor > 0 && ((data->step - data->prev_steps) % control->refactor == 0))
+        {
+            Setup_Preconditioner_QEq( system, control, data, workspace, mpi_data );
 
-        Compute_Preconditioner_QEq( system, control, data, workspace, mpi_data );
+            Compute_Preconditioner_QEq( system, control, data, workspace, mpi_data );
+        }
     }
-#endif
 
 
     for ( j = 0; j < system->n; ++j )
         workspace->s[j] = workspace->x[j][0];
-    s_matvecs = CG(system, workspace, workspace->H, workspace->b_s,//newQEq sCG
-            control->q_err, workspace->s, mpi_data, out_control->log , data->step );
+    iters = CG(system, control, data, workspace, workspace->H, workspace->b_s,
+            control->q_err, workspace->s, mpi_data, out_control->log , control->nprocs );
     for ( j = 0; j < system->n; ++j )
         workspace->x[j][0] = workspace->s[j];
 
-    //s_matvecs = PCG( system, workspace, workspace->H, workspace->b_s,
-    //   control->q_err, workspace->L, workspace->U, workspace->s,
-    //   mpi_data, out_control->log );
 #if defined(DEBUG)
     fprintf( stderr, "p%d: first CG completed\n", system->my_rank );
 #endif
 
     for ( j = 0; j < system->n; ++j )
         workspace->t[j] = workspace->x[j][1];
-    t_matvecs = CG(system, workspace, workspace->H, workspace->b_t,//newQEq sCG
-            control->q_err, workspace->t, mpi_data, out_control->log, data->step );
+    iters += CG(system, control, data, workspace, workspace->H, workspace->b_t,//newQEq sCG
+            control->q_err, workspace->t, mpi_data, out_control->log, control->nprocs );
     for ( j = 0; j < system->n; ++j )
         workspace->x[j][1] = workspace->t[j];
 
-    //t_matvecs = PCG( system, workspace, workspace->H, workspace->b_t,
-    //   control->q_err, workspace->L, workspace->U, workspace->t,
-    //   mpi_data, out_control->log );
 #if defined(DEBUG)
     fprintf( stderr, "p%d: second CG completed\n", system->my_rank );
 #endif
@@ -452,8 +464,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
 #if defined(LOG_PERFORMANCE)
     if ( system->my_rank == MASTER_NODE )
     {
-        data->timing.s_matvecs += s_matvecs;
-        data->timing.t_matvecs += t_matvecs;
+        data->timing.cm_solver_iters += iters;
     }
 #endif
 }

PuReMD/src/reax_types.h

@@ -41,18 +41,18 @@
 
 #define PURE_REAX
 //#define LAMMPS_REAX
-#define SAI_PRECONDITIONER
 //#define HALF_LIST
 //#define DEBUG
 //#define DEBUG_FOCUS
 //#define TEST_ENERGY
 //#define TEST_FORCES
-#define CG_PERFORMANCE
+//#define CG_PERFORMANCE
 #define LOG_PERFORMANCE
 #define STANDARD_BOUNDARIES
 //#define OLD_BOUNDARIES
 //#define MIDPOINT_BOUNDARIES
 
+#define ZERO                    (0.000000000000000e+00)
 #define REAX_MAX_STR            1024
 #define REAX_MAX_NBRS           6
 #define REAX_MAX_3BODY_PARAM    5
@@ -69,6 +69,20 @@
 #define NAN_REAL(a) (0)
 #endif
 
+
+/* preconditioner computation type for charge method linear solver */
+enum pre_comp
+{
+    NONE_PC = 0,
+    DIAG_PC = 1,
+    ICHOLT_PC = 2,
+    ILU_PAR_PC = 3,
+    ILUT_PAR_PC = 4,
+    ILU_SUPERLU_MT_PC = 5,
+    SAI_PC = 6,
+};
+
+
 /********************** TYPE DEFINITIONS ********************/
 typedef int  ivec[3];
 typedef double real;
@@ -534,7 +548,9 @@ typedef struct
     int tabulate;
 
     int qeq_freq;
+    int cm_solver_max_iters;
     real q_err;
+    int cm_solver_pre_comp_type;
     real sai_thres;
     int refactor;
     real droptol;
@@ -603,19 +619,44 @@ typedef struct
 
 typedef struct
 {
+    /* start time of event */
     real start;
+    /* end time of event */
     real end;
+    /* total elapsed time of event */
     real elapsed;
-
+    /* total simulation time */
     real total;
+    /* communication time */
     real comm;
+    /* neighbor list generation time */
     real nbrs;
+    /* force initialization time */
     real init_forces;
+    /* bonded force calculation time */
     real bonded;
+    /* non-bonded force calculation time */
     real nonb;
-    real qEq;
-    int  s_matvecs;
-    int  t_matvecs;
+    /* atomic charge distribution calculation time */
+    real cm;
+    /**/
+    real cm_sort_mat_rows;
+    /**/
+    real cm_solver_comm;
+    /**/
+    real cm_solver_pre_comp;
+    /**/
+    real cm_solver_pre_app; // update CG()
+    /* num. of steps in iterative linear solver for charge distribution */
+    int cm_solver_iters;
+    /**/
+    real cm_solver_spmv; // update CG()
+    /**/
+    real cm_solver_vector_ops; // update CG()
+    /**/
+    real cm_solver_orthog;
+    /**/
+    real cm_solver_tri_solve;
 } reax_timing;
 
 

PuReMD/src/reset_tools.c

@@ -106,15 +106,14 @@ void Reset_Simulation_Data( simulation_data* data, int virial )
 
 void Reset_Timing( reax_timing *rt )
 {
-    rt->total = Get_Time();
+    /*rt->total = Get_Time();
     rt->comm = 0;
     rt->nbrs = 0;
     rt->init_forces = 0;
     rt->bonded = 0;
     rt->nonb = 0;
     rt->qEq = 0;
-    rt->s_matvecs = 0;
-    rt->t_matvecs = 0;
+    rt->cm_solver_iters = 0;*/
 }
 
 #ifdef TEST_FORCES