From 583dead34502ff4847d0e5267be4bdb05fc366d1 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Wed, 17 Feb 2016 14:48:25 -0500
Subject: [PATCH] sPuReMD: code cleanup.
---
puremd_rc_1003/sPuReMD/GMRES.c | 20 +++++++++++++++-----
1 file changed, 15 insertions(+), 5 deletions(-)
diff --git a/puremd_rc_1003/sPuReMD/GMRES.c b/puremd_rc_1003/sPuReMD/GMRES.c
index 6bdde28f..aa40cf4a 100644
--- a/puremd_rc_1003/sPuReMD/GMRES.c
+++ b/puremd_rc_1003/sPuReMD/GMRES.c
@@ -682,8 +682,8 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
Sparse_MatVec( H, x, workspace->b_prm );
Vector_Sum( workspace->v[0], 1., b, -1., workspace->b_prm, N );
// TODO: add parameters to config file
- Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[0], workspace->v[0], 40 );
- Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[0], workspace->v[0], 40 );
+ Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[0], workspace->v[0], 30 );
+ Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[0], workspace->v[0], 30 );
workspace->g[0] = Norm( workspace->v[0], N );
Vector_Scale( workspace->v[0], 1. / workspace->g[0], workspace->v[0], N );
//fprintf( stderr, "res: %.15e\n", workspace->g[0] );
@@ -694,11 +694,14 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
/* matvec */
Sparse_MatVec( H, workspace->v[j], workspace->v[j + 1] );
// TODO: add parameters to config file
- Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[j + 1], workspace->v[j + 1], 40 );
- Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[j + 1], workspace->v[j + 1], 40 );
+ Jacobi_Iter( L, LOWER, Dinv_L, N, workspace->v[j + 1], workspace->v[j + 1], 30 );
+ Jacobi_Iter( U, UPPER, Dinv_U, N, workspace->v[j + 1], workspace->v[j + 1], 30 );
/* apply modified Gram-Schmidt to orthogonalize the new residual */
- for ( i = 0; i < j - 1; i++ ) workspace->h[i][j] = 0;
+ for ( i = 0; i < j - 1; i++ )
+ {
+ workspace->h[i][j] = 0;
+ }
//for( i = 0; i <= j; i++ ) {
for ( i = MAX(j - 1, 0); i <= j; i++ )
@@ -745,12 +748,15 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
}
+ /* TODO: solve using Jacobi iteration? */
/* solve Hy = g: H is now upper-triangular, do back-substitution */
for ( i = j - 1; i >= 0; i-- )
{
temp = workspace->g[i];
for ( k = j - 1; k > i; k-- )
+ {
temp -= workspace->h[i][k] * workspace->y[k];
+ }
workspace->y[i] = temp / workspace->h[i][i];
}
@@ -758,14 +764,18 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
/* update x = x_0 + Vy */
Vector_MakeZero( workspace->p, N );
for ( i = 0; i < j; i++ )
+ {
Vector_Add( workspace->p, workspace->y[i], workspace->v[i], N );
+ }
//Backward_Subs( U, workspace->p, workspace->p );
//Forward_Subs( L, workspace->p, workspace->p );
Vector_Add( x, 1., workspace->p, N );
/* stopping condition */
if ( fabs(workspace->g[j]) / bnorm <= tol )
+ {
break;
+ }
}
// Sparse_MatVec( H, x, workspace->b_prm );
--
GitLab