From 716ebd3fb2e5144416bfd25dadd929a9957e94e5 Mon Sep 17 00:00:00 2001
From: Abdullah Alperen <alperena@msu.edu>
Date: Thu, 13 Sep 2018 15:07:45 -0400
Subject: [PATCH] PuReMD: add sai_thres variable
---
PuReMD/src/control.c | 7 +++
PuReMD/src/geo_tools.c | 22 ++++-----
PuReMD/src/geo_tools.h | 8 ++--
PuReMD/src/linear_solvers.c | 89 ++++++++++++++++---------------------
PuReMD/src/linear_solvers.h | 2 +-
PuReMD/src/parallelreax.c | 10 ++++-
PuReMD/src/qEq.c | 8 ++--
PuReMD/src/reax_types.h | 5 ++-
8 files changed, 80 insertions(+), 71 deletions(-)
diff --git a/PuReMD/src/control.c b/PuReMD/src/control.c
index cd95c246..cbdd4147 100644
--- a/PuReMD/src/control.c
+++ b/PuReMD/src/control.c
@@ -75,6 +75,7 @@ char Read_Control_File( char *control_file, control_params* control,
control->qeq_freq = 1;
control->q_err = 1e-6;
+ control->sai_thres = 0.1;
control->refactor = 100;
control->droptol = 1e-2;;
@@ -203,6 +204,7 @@ char Read_Control_File( char *control_file, control_params* control,
{
ival = atoi( tmp[1] );
control->reneighbor = ival;
+ control->refactor = ival;
}
else if ( strcmp(tmp[0], "vlist_buffer") == 0 )
{
@@ -249,6 +251,11 @@ char Read_Control_File( char *control_file, control_params* control,
val = atof( tmp[1] );
control->q_err = val;
}
+ else if ( strcmp(tmp[0], "sai_thres") == 0 )
+ {
+ val = atof( tmp[1] );
+ control->sai_thres = val;
+ }
else if ( strcmp(tmp[0], "ilu_refactor") == 0 )
{
ival = atoi( tmp[1] );
diff --git a/PuReMD/src/geo_tools.c b/PuReMD/src/geo_tools.c
index 480b9ae7..c01f65cb 100644
--- a/PuReMD/src/geo_tools.c
+++ b/PuReMD/src/geo_tools.c
@@ -57,7 +57,7 @@ void Count_Geo_Atoms( FILE *geo, reax_system *system )
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d@count atoms:\n", system->my_rank );
- fprintf( stderr, "p%d: bigNNN = %d\n", system->my_rank, system->bigN );
+ fprintf( stderr, "p%d: bigN = %d\n", system->my_rank, system->bigN );
fprintf( stderr, "p%d: n = %d\n", system->my_rank, system->n );
fprintf( stderr, "p%d: N = %d\n\n", system->my_rank, system->N );
#endif
@@ -239,12 +239,12 @@ void Count_PDB_Atoms( FILE *geo, reax_system *system )
system->N = system->n;
- //#if defined(DEBUG)
+#if defined(DEBUG)
fprintf( stderr, "p%d@count atoms:\n", system->my_rank );
- fprintf( stderr, "p%d: bigNNN = %d\n", system->my_rank, system->bigN );
+ fprintf( stderr, "p%d: bigN = %d\n", system->my_rank, system->bigN );
fprintf( stderr, "p%d: n = %d\n", system->my_rank, system->n );
fprintf( stderr, "p%d: N = %d\n\n", system->my_rank, system->N );
- //#endif
+#endif
}
@@ -497,7 +497,7 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
Also, we do not write connect lines yet.
*/
-char Write_PDB(reax_system* system, reax_list* bonds, simulation_data *data,
+char Write_PDB(reax_system* system, reax_list** bonds, simulation_data *data,
control_params *control, mpi_datatypes *mpi_data,
output_controls *out_control)
{
@@ -551,13 +551,13 @@ char Write_PDB(reax_system* system, reax_list* bonds, simulation_data *data,
sprintf(fname, "%s-%d.pdb", control->sim_name, data->step);
pdb = fopen(fname, "w");
- fprintf( pdb, PDB_CRYST1_FORMAT_O,
+ /*fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
system->big_box.box_norms[0], system->big_box.box_norms[1],
system->big_box.box_norms[2],
RAD2DEG(alpha), RAD2DEG(beta), RAD2DEG(gamma), " ", 0 );
fprintf( out_control->log, "Box written\n" );
- fflush( out_control->log );
+ fflush( out_control->log );*/
}
/*write atom lines to buffer*/
@@ -566,11 +566,13 @@ char Write_PDB(reax_system* system, reax_list* bonds, simulation_data *data,
p_atom = &(system->my_atoms[i]);
strncpy(name, p_atom->name, 8);
Trim_Spaces(name);
- sprintf( line, PDB_ATOM_FORMAT_O,
+ /*sprintf( line, PDB_ATOM_FORMAT_O,
"ATOM ", p_atom->orig_id, p_atom->name, ' ', "REX", ' ', 1, ' ',
p_atom->x[0], p_atom->x[1], p_atom->x[2],
- 1.0, 0.0, "0", name, " " );
- fprintf(stderr, "PDB NAME <%s>\n", p_atom->name);
+ 1.0, 0.0, "0", name, " " );*/
+ sprintf( line, PDB_ATOM_FORMAT_O,
+ p_atom->orig_id, p_atom->x[0], p_atom->x[1], p_atom->x[2] );
+ //fprintf( stderr, "PDB NAME <%s>\n", p_atom->name);
strncpy( buffer + i * PDB_ATOM_FORMAT_O_LENGTH, line,
PDB_ATOM_FORMAT_O_LENGTH );
}
diff --git a/PuReMD/src/geo_tools.h b/PuReMD/src/geo_tools.h
index 80786856..d34caaff 100644
--- a/PuReMD/src/geo_tools.h
+++ b/PuReMD/src/geo_tools.h
@@ -110,14 +110,16 @@ COLUMNS DATA TYPE FIELD DEFINITION
#define PDB_CONECT_FORMAT "%6s%5d%5d%5d%5d%5d\n"
#define PDB_CRYST1_FORMAT "%6s%9s%9s%9s%7s%7s%7s%11s%4s\n"
-#define PDB_ATOM_FORMAT_O "%6s%5d %4s%c%3s %c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %-4s%2s%2s\n"
-#define PDB_ATOM_FORMAT_O_LENGTH 81
+//#define PDB_ATOM_FORMAT_O "%6s%5d %4s%c%3s %c%4d%c %8.3f%8.3f%8.3f%6.2f%6.2f %-4s%2s%2s\n"
+#define PDB_ATOM_FORMAT_O "%5d%8.3f%8.3f%8.3f\n"
+//#define PDB_ATOM_FORMAT_O_LENGTH 81
+#define PDB_ATOM_FORMAT_O_LENGTH 30
#define PDB_CRYST1_FORMAT_O "%6s%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f%11s%4d\n"
char Read_PDB( char*, reax_system*, control_params*,
simulation_data*, storage*, mpi_datatypes* );
-char Write_PDB( reax_system*, reax_list*, simulation_data*,
+char Write_PDB( reax_system*, reax_list**, simulation_data*,
control_params*, mpi_datatypes*, output_controls* );
#endif
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index 9ec9f524..0bd9faeb 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -104,10 +104,8 @@ int find_bucket( double *list, int len, double a )
void setup_sparse_approx_inverse( reax_system *system, storage *workspace,
mpi_datatypes *mpi_data, sparse_matrix *A, sparse_matrix **A_spar_patt,
- int nprocs, double filter )
+ int nprocs, real filter )
{
- /* Print_Sparse_Matrix2( system, A, "Charge_Matrix_MPI_Step0.txt" ); */
-
int i, bin, total, pos, push;
int n, n_gather, m, s_local, s, n_local;
int target_proc;
@@ -146,8 +144,6 @@ void setup_sparse_approx_inverse( reax_system *system, storage *workspace,
comm = mpi_data->world;
- // fprintf("%d\n",);
-
if ( *A_spar_patt == NULL )
{
Allocate_Matrix2(A_spar_patt, A->n, system->local_cap, A->m, comm );
@@ -492,8 +488,6 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
/* use a Dist-like approach to send the row information */
for ( d = 0; d < 3; ++d)
{
- //fprintf( stderr, "p%d, d = %d\n", system->my_rank, d);
-
flag1 = 0;
flag2 = 0;
cnt = 0;
@@ -681,14 +675,11 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
}
}
- fprintf( stderr,"p%d, After Manual Dist Call\n", system->my_rank );
-
X = (int *) malloc( sizeof(int) * (system->bigN + 1) );
pos_x = (int *) malloc( sizeof(int) * (system->bigN + 1) );
for ( i = 0; i < A_spar_patt->n; ++i )
{
- //fprintf( stderr, "p%d, column %d\n", system->my_rank, i);
N = 0;
M = 0;
for ( k = 0; k <= system->bigN; ++k )
@@ -764,8 +755,9 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
local_pos = A_spar_patt->entries[ A_spar_patt->start[i] + d_j ].j;
if( local_pos < 0 || local_pos >= system->N )
{
- fprintf( stdout, "THE LOCAL POSITION OF THE ATOM IS NOT VALID, STOP THE EXECUTION\n");
- fflush( stdout );
+ fprintf( stderr, "THE LOCAL POSITION OF THE ATOM IS NOT VALID, STOP THE EXECUTION\n");
+ fflush( stderr );
+
}
if( local_pos < A->n )
{
@@ -774,8 +766,8 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
atom = &system->my_atoms[ A->entries[d_i].j ];
if (pos_x[ atom->orig_id ] >= M || d_j >= N )
{
- fprintf( stdout, "CANNOT MAP IT TO THE DENSE MATRIX, STOP THE EXECUTION, orig_id = %d, i = %d, j = %d, M = %d N = %d\n", atom->orig_id, pos_x[ atom->orig_id ], d_j, M, N );
- fflush( stdout );
+ fprintf( stderr, "CANNOT MAP IT TO THE DENSE MATRIX, STOP THE EXECUTION, orig_id = %d, i = %d, j = %d, M = %d N = %d\n", atom->orig_id, pos_x[ atom->orig_id ], d_j, M, N );
+ fflush( stderr );
}
if ( X[ atom->orig_id ] == 1 )
{
@@ -789,8 +781,8 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
{
if (pos_x[ j_list[local_pos][d_i] ] >= M || d_j >= N )
{
- fprintf( stdout, "CANNOT MAP IT TO THE DENSE MATRIX, STOP THE EXECUTION, %d %d\n", pos_x[ j_list[local_pos][d_i] ], d_j);
- fflush( stdout );
+ fprintf( stderr, "CANNOT MAP IT TO THE DENSE MATRIX, STOP THE EXECUTION, %d %d\n", pos_x[ j_list[local_pos][d_i] ], d_j);
+ fflush( stderr );
}
if ( X[ j_list[local_pos][d_i] ] == 1 )
{
@@ -811,7 +803,7 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
e_j[identity_pos] = 1.0;
/* Solve the overdetermined system AX = B through the least-squares problem:
- * * * min ||B - AX||_2 */
+ * min ||B - AX||_2 */
m = M;
n = N;
nrhs = 1;
@@ -842,9 +834,6 @@ void sparse_approx_inverse(reax_system *system, storage *workspace, mpi_datatype
free( e_j );
}
-
- fprintf( stderr,"p%d, After Building SAI Matrix\n", system->my_rank );
-
free( pos_x);
free( X );
@@ -1047,7 +1036,7 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H,
workspace->t[j] = workspace->x[j][1];
}
matvecs = CG( system, workspace, H, workspace->b_t, tol,
- workspace->t,mpi_data, fout );
+ workspace->t,mpi_data, fout, -1 );
for ( j = 0; j < n; ++j )
{
workspace->x[j][1] = workspace->t[j];
@@ -1060,7 +1049,7 @@ int dual_CG( reax_system *system, storage *workspace, sparse_matrix *H,
workspace->s[j] = workspace->x[j][0];
}
matvecs = CG( system, workspace, H, workspace->b_s, tol, workspace->s,
- mpi_data, fout );
+ mpi_data, fout, -1 );
for ( j = 0; j < system->n; ++j )
{
workspace->x[j][0] = workspace->s[j];
@@ -1126,34 +1115,34 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
* where:
* A: matrix, stored in CSR format
* x: vector
- * b: vector (result) */
+ * b: vector (result)
static void Sparse_MatVec_full( const sparse_matrix * const A,
const real * const x, real * const b )
{
- //TODO: implement full SpMV in MPI
- // int i, pj;
- //
- // Vector_MakeZero( b, A->n );
- //
- //#ifdef _OPENMP
- // #pragma omp for schedule(static)
- //#endif
- // for ( i = 0; i < A->n; ++i )
- // {
- // for ( pj = A->start[i]; pj < A->start[i + 1]; ++pj )
- // {
- // b[i] += A->val[pj] * x[A->j[pj]];
- // }
- // }
-}
+ TODO: implement full SpMV in MPI
+ int i, pj;
+
+ Vector_MakeZero( b, A->n );
+
+ #ifdef _OPENMP
+ #pragma omp for schedule(static)
+ #endif
+ for ( i = 0; i < A->n; ++i )
+ {
+ for ( pj = A->start[i]; pj < A->start[i + 1]; ++pj )
+ {
+ b[i] += A->val[pj] * x[A->j[pj]];
+ }
+ }
+}*/
int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
- real tol, real *x, mpi_datatypes* mpi_data, FILE *fout )
+ real tol, real *x, mpi_datatypes* mpi_data, FILE *fout, int step )
{
- int i, j, scale;
+ int i, scale;
real tmp, alpha, beta, b_norm;
- real sig_old, sig_new, sig0;
+ real sig_old, sig_new;
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -1167,7 +1156,6 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
Dist( system, mpi_data, x, MPI_DOUBLE, scale, real_packer );
Sparse_MatVec( H, x, workspace->q, system->N );
#if defined(HALF_LIST)
- // tryQEq
Coll( system, mpi_data, workspace->q, MPI_DOUBLE, scale, real_unpacker );
#endif
@@ -1193,7 +1181,6 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
b_norm = Parallel_Norm( b, system->n, mpi_data->world );
sig_new = Parallel_Dot(workspace->r, workspace->d, system->n, mpi_data->world);
- sig0 = sig_new;
#if defined(CG_PERFORMANCE)
if ( system->my_rank == MASTER_NODE )
@@ -1248,7 +1235,7 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
if ( i >= 300 )
{
- fprintf( stderr, "CG convergence failed!\n" );
+ fprintf( stderr, "CG convergence failed at step %d!\n", step );
return i;
}
@@ -1269,7 +1256,7 @@ int CG_test( reax_system *system, storage *workspace, sparse_matrix *H,
{
int i, j, scale;
real tmp, alpha, beta, b_norm;
- real sig_old, sig_new, sig0;
+ real sig_old, sig_new;
scale = sizeof(real) / sizeof(void);
b_norm = Parallel_Norm( b, system->n, mpi_data->world );
@@ -1293,7 +1280,6 @@ int CG_test( reax_system *system, storage *workspace, sparse_matrix *H,
sig_new = Parallel_Dot( workspace->r, workspace->d, system->n,
mpi_data->world );
- sig0 = sig_new;
#if defined(DEBUG)
//if( system->my_rank == MASTER_NODE ) {
fprintf( stderr, "p%d CG:sig_new=%24.15e,d_norm=%24.15e,q_norm=%24.15e\n",
@@ -1437,11 +1423,15 @@ int PCG( reax_system *system, storage *workspace,
sparse_matrix *L, sparse_matrix *U, real *x,
mpi_datatypes* mpi_data, FILE *fout )
{
- int i, me, n, N, scale;
+ int i, n, N, scale;
real tmp, alpha, beta, b_norm, r_norm, sig_old, sig_new;
MPI_Comm world;
+#if defined(DEBUG_FOCUS)
+ int me;
me = system->my_rank;
+#endif
+
n = system->n;
N = system->N;
world = mpi_data->world;
@@ -1525,7 +1515,7 @@ int sCG( reax_system *system, storage *workspace, sparse_matrix *H,
{
int i, j;
real tmp, alpha, beta, b_norm;
- real sig_old, sig_new, sig0;
+ real sig_old, sig_new;
b_norm = Norm( b, system->n );
#if defined(DEBUG)
@@ -1547,7 +1537,6 @@ int sCG( reax_system *system, storage *workspace, sparse_matrix *H,
workspace->d[j] = workspace->r[j] * workspace->Hdia_inv[j]; //pre-condition
sig_new = Dot( workspace->r, workspace->d, system->n );
- sig0 = sig_new;
#if defined(DEBUG)
if ( system->my_rank == MASTER_NODE )
{
diff --git a/PuReMD/src/linear_solvers.h b/PuReMD/src/linear_solvers.h
index 292a02ca..3b6fc1a0 100644
--- a/PuReMD/src/linear_solvers.h
+++ b/PuReMD/src/linear_solvers.h
@@ -37,7 +37,7 @@ int GMRES_HouseHolder( reax_system*, storage*, sparse_matrix*,
int dual_CG( reax_system*, storage*, sparse_matrix*,
rvec2*, real, rvec2*, mpi_datatypes*, FILE* );
int CG( reax_system*, storage*, sparse_matrix*,
- real*, real, real*, mpi_datatypes*, FILE* );
+ real*, real, real*, mpi_datatypes*, FILE*, int );
int PCG( reax_system*, storage*, sparse_matrix*, real*, real,
sparse_matrix*, sparse_matrix*, real*, mpi_datatypes*, FILE* );
int sCG( reax_system*, storage*, sparse_matrix*,
diff --git a/PuReMD/src/parallelreax.c b/PuReMD/src/parallelreax.c
index 4b401ad6..ea30ceed 100644
--- a/PuReMD/src/parallelreax.c
+++ b/PuReMD/src/parallelreax.c
@@ -211,6 +211,7 @@ int main( int argc, char* argv[] )
#endif
/* start the simulation */
+ int total_itr = data->timing.s_matvecs + data->timing.t_matvecs;
for ( ++data->step; data->step <= control->nsteps; data->step++ )
{
if ( control->T_mode )
@@ -218,6 +219,10 @@ int main( int argc, char* argv[] )
Evolve( system, control, data, workspace, lists, out_control, mpi_data );
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;
+ }
Output_Results( system, control, data, lists, out_control, mpi_data );
//Analysis(system, control, data, workspace, lists, out_control, mpi_data);
@@ -230,6 +235,8 @@ int main( int argc, char* argv[] )
else if ( out_control->restart_format == WRITE_BINARY )
Write_Binary_Restart( system, control, data, out_control, mpi_data );
}
+ /*if(data->step == 1 || data->step == control->nsteps)
+ Write_PDB( system, lists, data, control, mpi_data, out_control );*/
#if defined(DEBUG)
fprintf( stderr, "p%d: step%d completed\n", system->my_rank, data->step );
@@ -242,9 +249,10 @@ int main( int argc, char* argv[] )
{
t_elapsed = Get_Timing_Info( t_start );
fprintf( out_control->out, "Total Simulation Time: %.2f secs\n", t_elapsed );
+ fprintf( out_control->log, "Avg. # of Solver Itrs: %.2f\n", total_itr/((double)control->nsteps) );
}
- // Write_PDB( &system, &(lists[BOND]), &out_control );
+ //Write_PDB( system, lists, data, control, mpi_data, out_control );
Close_Output_Files( system, control, out_control, mpi_data );
MPI_Finalize();
diff --git a/PuReMD/src/qEq.c b/PuReMD/src/qEq.c
index bd22b54f..87b45bdd 100644
--- a/PuReMD/src/qEq.c
+++ b/PuReMD/src/qEq.c
@@ -373,7 +373,7 @@ static void Setup_Preconditioner_QEq( reax_system *system, control_params *contr
//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,
- control->nprocs, 0.1 );
+ control->nprocs, control->sai_thres );
}
static void Compute_Preconditioner_QEq( reax_system *system, control_params *control,
@@ -406,7 +406,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
//t_matvecs = 0;
#if defined(SAI_PRECONDITIONER)
- //if( control->refactor > 0 && ((data->step - data->prev_steps) % control->refactor == 0))
+ if( control->refactor > 0 && ((data->step - data->prev_steps) % control->refactor == 0))
{
Setup_Preconditioner_QEq( system, control, data, workspace, mpi_data );
@@ -418,7 +418,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
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 );
+ control->q_err, workspace->s, mpi_data, out_control->log , data->step );
for ( j = 0; j < system->n; ++j )
workspace->x[j][0] = workspace->s[j];
@@ -432,7 +432,7 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
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 );
+ control->q_err, workspace->t, mpi_data, out_control->log, data->step );
for ( j = 0; j < system->n; ++j )
workspace->x[j][1] = workspace->t[j];
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 30a3daca..58ac8498 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -41,8 +41,8 @@
#define PURE_REAX
//#define LAMMPS_REAX
-#define SAI_PRECONDITIONER
-//#define HALF_LIST
+//#define SAI_PRECONDITIONER
+#define HALF_LIST
//#define DEBUG
//#define DEBUG_FOCUS
//#define TEST_ENERGY
@@ -535,6 +535,7 @@ typedef struct
int qeq_freq;
real q_err;
+ real sai_thres;
int refactor;
real droptol;
--
GitLab