PuReMD: beginning setting up SAI preconditioning.

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,

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);

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] );

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;