From 4be2c517a3d32bb40b808557bd52e9d9d2b0d2fe Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Thu, 19 Apr 2018 15:37:42 -0400
Subject: [PATCH] PuReMD: beginning setting up SAI preconditioning.
---
PuReMD/src/init_md.c | 6 ++++++
PuReMD/src/linear_solvers.c | 32 ++++++++++++++++++++++++++++++++
PuReMD/src/qEq.c | 6 ++++++
PuReMD/src/reax_types.h | 8 +++++++-
4 files changed, 51 insertions(+), 1 deletion(-)
diff --git a/PuReMD/src/init_md.c b/PuReMD/src/init_md.c
index caad5420..da661aa4 100644
--- a/PuReMD/src/init_md.c
+++ b/PuReMD/src/init_md.c
@@ -587,6 +587,12 @@ int Init_Lists( reax_system *system, control_params *control,
Allocate_Matrix( &(workspace->H), system->local_cap, Htop, comm );
workspace->L = NULL;
workspace->U = NULL;
+
+ //TODO: uncomment for SAI
+// Allocate_Matrix( &(workspace->H_spar_patt), workspace->H->n, workspace->H->m );
+// Allocate_Matrix( &(workspace->H_spar_patt_full), workspace->H->n, 2 * workspace->H->m - workspace->H->n );
+// Allocate_Matrix( &(workspace->H_app_inv), workspace->H->n, 2 * workspace->H->m - workspace->H->n );
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "p%d: allocated H matrix: Htop=%d, space=%dMB\n",
system->my_rank, Htop,
diff --git a/PuReMD/src/linear_solvers.c b/PuReMD/src/linear_solvers.c
index d08dfe1e..541a132b 100644
--- a/PuReMD/src/linear_solvers.c
+++ b/PuReMD/src/linear_solvers.c
@@ -278,6 +278,32 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b, int N )
}
+/* sparse matrix-vector product Ax = b
+ * where:
+ * A: matrix, stored in CSR format
+ * x: vector
+ * 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]];
+// }
+// }
+}
+
+
int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
real tol, real *x, mpi_datatypes* mpi_data, FILE *fout )
{
@@ -307,10 +333,13 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
#endif
Vector_Sum( workspace->r , 1., b, -1., workspace->q, system->n );
+
for ( j = 0; j < system->n; ++j )
{
workspace->d[j] = workspace->r[j] * workspace->Hdia_inv[j]; //pre-condition
}
+ //TODO: apply SAI preconditioner here, comment out diagonal preconditioning above
+// Sparse_MatVec_full( workspace->H_app_inv, workspace->r, workspace->d );
b_norm = Parallel_Norm( b, system->n, mpi_data->world );
sig_new = Parallel_Dot(workspace->r, workspace->d, system->n, mpi_data->world);
@@ -341,11 +370,14 @@ int CG( reax_system *system, storage *workspace, sparse_matrix *H, real *b,
alpha = sig_new / tmp;
Vector_Add( x, alpha, workspace->d, system->n );
Vector_Add( workspace->r, -alpha, workspace->q, system->n );
+
/* pre-conditioning */
for ( j = 0; j < system->n; ++j )
{
workspace->p[j] = workspace->r[j] * workspace->Hdia_inv[j];
}
+ //TODO: apply SAI preconditioner here, comment out diagonal preconditioning above
+// Sparse_MatVec_full( workspace->H_app_inv, workspace->r, workspace->d );
sig_old = sig_new;
sig_new = Parallel_Dot(workspace->r, workspace->p, system->n, mpi_data->world);
diff --git a/PuReMD/src/qEq.c b/PuReMD/src/qEq.c
index ed5c276b..c8a348c7 100644
--- a/PuReMD/src/qEq.c
+++ b/PuReMD/src/qEq.c
@@ -272,6 +272,12 @@ void Init_MatVec( reax_system *system, simulation_data *data,
#endif
}*/
+ //TODO: fill in code for setting up and computing SAI, see sPuReMD code,
+ // and remove diagonal preconditioner computation below (workspace->Hdia_inv)
+// setup_sparse_approx_inverse( Hptr, &workspace->H_full, &workspace->H_spar_patt,
+// &workspace->H_spar_patt_full, &workspace->H_app_inv,
+// control->cm_solver_pre_comp_sai_thres );
+
for ( i = 0; i < system->n; ++i )
{
atom = &( system->my_atoms[i] );
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 3ebcd73c..152c0c54 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -800,7 +800,13 @@ typedef struct
int *bond_mark, *done_after;
/* QEq storage */
- sparse_matrix *H, *L, *U;
+ sparse_matrix *H;
+ sparse_matrix *L;
+ sparse_matrix *U;
+ sparse_matrix *H_full;
+ sparse_matrix *H_spar_patt;
+ sparse_matrix *H_spar_patt_full;
+ sparse_matrix *H_app_inv;
real *Hdia_inv, *b_s, *b_t, *b_prc, *b_prm, *s, *t;
real *droptol;
rvec2 *b, *x;
--
GitLab