From 09df8f0820a7591eefd7380cd3ec525803ec50d1 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Wed, 24 Feb 2016 12:59:39 -0500
Subject: [PATCH] sPuReMD: fix diagonal scaling in factorization.
---
.gitignore | 5 +++++
puremd_rc_1003/sPuReMD/QEq.c | 27 +++++++++++++++------------
puremd_rc_1003/sPuReMD/forces.c | 1 +
3 files changed, 21 insertions(+), 12 deletions(-)
diff --git a/.gitignore b/.gitignore
index be8bc24c..5ad3a172 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,8 @@
+*.log
+*.pdb
+*.pot
+*.trj
+
# TeX
*.aux
*.log
diff --git a/puremd_rc_1003/sPuReMD/QEq.c b/puremd_rc_1003/sPuReMD/QEq.c
index 8ef7ddae..6c268de2 100644
--- a/puremd_rc_1003/sPuReMD/QEq.c
+++ b/puremd_rc_1003/sPuReMD/QEq.c
@@ -333,7 +333,7 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
sum = DAD->entries[j].val - sum;
/* diagonal entries */
- if( (U_t->start[k] - 1) == j )
+ if( (k - 1) == U_t->entries[j].j )
{
/* sanity check */
if( sum < ZERO )
@@ -352,8 +352,9 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
}
}
- /* apply inverse transformation D^{-1}DADD^{-1} = A \approx D^{-1}U^{T}UD^{-1},
- * i.e., D^{-1}U^{T} */
+ /* apply inverse transformation D^{-1}U^{T},
+ * since DAD \approx U^{T}U, so
+ * D^{-1}DADD^{-1} = A \approx D^{-1}U^{T}UD^{-1} */
for ( i = 0; i < A->n; ++i )
{
for ( pj = A->start[i]; pj < A->start[i + 1]; ++pj )
@@ -362,10 +363,10 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
}
}
- /* transpose U^T and copy into U */
+ /* transpose U^{T} and copy into U */
for ( i = 0; i < U_t->n; ++i )
{
- /* count nonzeros in each row of U (column-wise),
+ /* count nonzeros in each row of U (columns of U^{T}),
* store in U->start */
for ( pj = U_t->start[i]; pj < U_t->start[i + 1]; ++pj )
{
@@ -375,7 +376,7 @@ static void ICHOL_PAR( const sparse_matrix * const A, const unsigned int sweeps,
/* set correct row pointer in U */
for ( i = 1; i <= U->n; ++i )
{
- Utop[i] = U->start[i] = U->start[i] + U->start[i - 1] + 1;
+ Utop[i] = U->start[i] = U->start[i] + U->start[i - 1];
}
/* for each row */
for ( i = 0; i < U_t->n; ++i )
@@ -424,8 +425,9 @@ void Init_MatVec( reax_system *system, control_params *control,
// fprintf( stderr, "not enough memory for LU matrices. terminating.\n" );
// exit(INSUFFICIENT_SPACE);
// }
- if ( Allocate_Matrix( &(workspace->L), far_nbrs->n, workspace->H->m ) == 0 ||
- Allocate_Matrix( &(workspace->U), far_nbrs->n, workspace->H->m ) == 0 )
+ /* 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);
@@ -439,13 +441,14 @@ void Init_MatVec( reax_system *system, control_params *control,
// ICHOLT( workspace->H, workspace->droptol, workspace->L, workspace->U );
// TODO: add parameters for sweeps to control file
- ICHOL_PAR( workspace->H, 1, workspace->L, workspace->U );
+ ICHOL_PAR( workspace->H, 1, workspace->L, workspace->U );
+#if defined(DEBUG_FOCUS)
sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
Print_Sparse_Matrix2( workspace->L, fname );
sprintf( fname, "%s.U%d.out", control->sim_name, data->step );
Print_Sparse_Matrix2( workspace->U, fname );
-#if defined(DEBUG_FOCUS)
+
fprintf( stderr, "icholt-" );
//sprintf( fname, "%s.L%d.out", control->sim_name, data->step );
//Print_Sparse_Matrix2( workspace->L, fname );
@@ -522,8 +525,8 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
Init_MatVec( system, control, data, workspace, far_nbrs );
- //if( data->step % 10 == 0 )
- // Print_Linear_System( system, control, workspace, far_nbrs, data->step );
+// if( data->step % 10 == 0 )
+// Print_Linear_System( system, control, workspace, data->step );
//TODO: add parameters in control file for solver choice and options
// matvecs = GMRES( workspace, workspace->H,
diff --git a/puremd_rc_1003/sPuReMD/forces.c b/puremd_rc_1003/sPuReMD/forces.c
index ae6128c6..d9348b33 100644
--- a/puremd_rc_1003/sPuReMD/forces.c
+++ b/puremd_rc_1003/sPuReMD/forces.c
@@ -318,6 +318,7 @@ void Init_Forces( reax_system *system, control_params *control,
atom_j = &(system->atoms[j]);
flag = 0;
+ /*TODO: build H with smaller cutoff, use for preconditioning*/
if ((data->step - data->prev_steps) % control->reneighbor == 0)
{
if ( nbr_pj->d <= control->r_cut)
--
GitLab