diff --git a/puremd_rc_1003/sPuReMD/GMRES.c b/puremd_rc_1003/sPuReMD/GMRES.c
index ab790fe81952d106039286f318e2d7ef5b523515..8488ac69ff5a0180314ab28eb3513347824a9f53 100644
--- a/puremd_rc_1003/sPuReMD/GMRES.c
+++ b/puremd_rc_1003/sPuReMD/GMRES.c
@@ -25,8 +25,13 @@
#include "vector.h"
-/* sparse matrix-vector product */
-void Sparse_MatVec( sparse_matrix *A, real *x, real *b )
+/* sparse matrix-vector product Ax=b
+ * where:
+ * A: lower triangular matrix
+ * x: vector
+ * b: vector (result) */
+static void Sparse_MatVec( const sparse_matrix * const A,
+ const real * const x, real * const b )
{
int i, j, k, n, si, ei;
real H;
@@ -51,7 +56,47 @@ void Sparse_MatVec( sparse_matrix *A, real *x, real *b )
}
// the diagonal entry is the last one in
- b[i] += A->entries[k].val * x[i];
+ b[i] += A->entries[A->start[i+1]-1].val * x[i];
+ }
+}
+
+
+/* sparse matrix-vector product Ax=b
+ * ommitting diagonal of A (for Jacobi iteration),
+ * A: triangular matrix
+ * x: vector
+ * b: vector (result)
+ * tri: triangularity of A (lower/upper) */
+static void Sparse_MatVec2( const sparse_matrix * const A,
+ const TRIANGULARITY tri, const real * const x, real * const b)
+{
+ int i, k, n, si = 0, ei = 0;
+
+ n = A->n;
+ for( i = 0; i < n; ++i )
+ {
+ b[i] = 0;
+ }
+
+ for( i = 0; i < n; ++i )
+ {
+ if(tri == LOWER)
+ {
+ si = A->start[i];
+ ei = A->start[i+1]-1;
+ }
+ else if(tri == UPPER)
+ {
+
+ si = A->start[i]+1;
+ ei = A->start[i+1];
+ }
+
+ for( k = si; k < ei; ++k )
+ {
+ b[A->entries[k].j] += A->entries[k].j * x[i];
+ b[i] += A->entries[k].j * x[A->entries[k].j];
+ }
}
}
@@ -113,62 +158,29 @@ void Backward_Subs( sparse_matrix *U, real *y, real *x )
*
* Note: Newmann series arises from series expansion of the inverse of
* the coefficient matrix in the triangular system. */
-void Jacobi_Iter( const sparse_matrix *G, const sparse_matrix *Dinv,
- const real *b, const real *x_0, real *x, const int maxiter )
+static void Jacobi_Iter( const sparse_matrix * const G, const TRIANGULARITY tri,
+ const sparse_matrix * const Dinv, const real * const b,
+ real * const x, const unsigned int maxiter )
{
- unsigned int i, iter = 0;
+ unsigned int iter = 0;
real *Dinv_b;
- if( (Dinv_b = (real*) malloc(sizeof(real)*(Dinv->m))) == NULL )
+ if( (Dinv_b = (real*) malloc(sizeof(real)*(Dinv->n))) == NULL )
{
fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
exit(INSUFFICIENT_SPACE);
}
/* precompute and cache, as invariant in loop below */
- Sparse_MatVec( (sparse_matrix*) Dinv, (real*) b, Dinv_b );
+ Sparse_MatVec( (sparse_matrix*)Dinv, (real*)b, Dinv_b );
do
{
// x_{k+1} = G*x_{k} + Dinv*b;
- Sparse_MatVec( (sparse_matrix*) &G, x, x );
- //#pragma omp parallel for
+ Sparse_MatVec2( (sparse_matrix*)G, tri, x, x );
Vector_Add2( x, Dinv_b, Dinv->n );
++iter;
} while( iter < maxiter );
-
-// Dinv_c = Dinv*c;
-// if ((nargin < 4) || isempty(x0))
-// x0 = zeros(size(c));
-// end
-// if ((nargin < 5) || isempty(tol))
-// tol = Inf;
-// end
-// if ((nargin == 7) && ~isempty(fp))
-// write = true;
-// else
-// write = false;
-// end
-// x = x0;
-// x_old = x0;
-// relres = Inf;
-// relres_old = Inf;
-// resvec = [];
-// iters = 0;
-// first = true;
-//
-// % while first || (relres >= tol && relres < relres_old && iters < maxiters)
-// % while first || (relres >= tol && iters < maxiters)
-// while first || (iters < maxiters)
-// x_old = x;
-// relres_old = relres;
-// x = G*x_old + Dinv_c;
-// iters = iters + 1;
-//
-// first = false;
-// relres = norm(x-x_old);
-// resvec(iters) = relres;
-// end
free(Dinv_b);
}
@@ -603,92 +615,39 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
{
int i, j, k, itr, N, pj, si, ei;
real cc, tmp1, tmp2, temp, bnorm;
- sparse_matrix *G_L, *G_U, *Dinv_L, *Dinv_U;
+ sparse_matrix *Dinv_L, *Dinv_U;
N = H->n;
bnorm = Norm( b, N );
- // TODO: compute Jacobi iteration matrices
- if( Allocate_Matrix( &G_L, L->n, (L->m - L->n) ) == 0 ||
- Allocate_Matrix( &G_U, U->n, (U->m - U->n) ) == 0 ||
- Allocate_Matrix( &Dinv_L, L->n, L->n ) == 0 ||
- Allocate_Matrix( &Dinv_U, U->n, U->n ) == 0 ){
- fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
- exit(INSUFFICIENT_SPACE);
- }
- /* truncated Newmann series: x_{k+1} = Gx_k + D^{-1}b
+ /* Compute Jacobi iteration matrices from
+ * truncated Newmann series: x_{k+1} = Gx_k + D^{-1}b
* where:
* G = I - D^{-1}R
* R = triangular matrix
* D = diagonal matrix, diagonals from R */
-
- /* construct D^{-1}_L and D^{-1}_U */
- for( i = 0; i < N; ++i )
+ if( Allocate_Matrix( &Dinv_L, L->n, L->n ) == 0 ||
+ Allocate_Matrix( &Dinv_U, U->n, U->n ) == 0 )
{
- si = L->start[i];
- ei = L->start[i+1];
- for( pj = si; pj < ei; ++pj )
- {
- if( L->entries[pj].j == i)
- {
- Dinv_L->start[i] = i;
- Dinv_L->entries[i].j = i;
- Dinv_L->entries[i].val = 1.0 / L->entries[pj].val;
- break;
- }
- }
-
- si = U->start[i];
- ei = U->start[i+1];
- for( pj = si; pj < ei; ++pj )
- {
- if( U->entries[pj].j == i)
- {
- Dinv_U->start[i] = i;
- Dinv_U->entries[i].j = i;
- Dinv_U->entries[i].val = 1.0 / U->entries[pj].val;
- break;
- }
- }
+ fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
+ exit(INSUFFICIENT_SPACE);
}
- /* construct G_L and G_U */
- G_L->n = L->n;
- G_L->m = (L->m-L->n);
- G_U->n = U->n;
- G_U->m = (U->m-U->n);
- j = 0;
- k = 0;
+ /* construct D^{-1}_L and D^{-1}_U */
for( i = 0; i < N; ++i )
{
- si = L->start[i];
- ei = L->start[i+1];
-
- for( pj = si; pj < ei; ++pj )
- {
- /* zeros along diagonal, so skip */
- if( L->entries[pj].j != i )
- {
- G_L->entries[j].j = L->entries[pj].j;
- G_L->entries[j].val = -(Dinv_L->entries[i].val * L->entries[j].val);
- ++j;
- }
- }
+ si = L->start[i+1]-1;
+ Dinv_L->start[i] = i;
+ Dinv_L->entries[i].j = i;
+ Dinv_L->entries[i].val = 1. / L->entries[si].val;
si = U->start[i];
- ei = U->start[i+1];
-
- for( pj = si; pj < ei; ++pj )
- {
- /* zeros along diagonal, so skip */
- if( U->entries[pj].j != i )
- {
- G_U->entries[k].j = U->entries[pj].j;
- G_U->entries[k].val = -(Dinv_U->entries[i].val * U->entries[k].val);
- ++k;
- }
- }
+ Dinv_U->start[i] = i;
+ Dinv_U->entries[i].j = i;
+ Dinv_U->entries[i].val = 1. / U->entries[si].val;
}
+ Dinv_L->start[N] = N;
+ Dinv_U->start[N] = N;
/* GMRES outer-loop */
for( itr = 0; itr < MAX_ITR; ++itr ) {
@@ -697,12 +656,10 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
Vector_Sum( workspace->v[0], 1., b, -1., workspace->b_prm, N );
// TODO: add parameters to config file
// TODO: cache results and use for inital guess?
- Jacobi_Iter( (const sparse_matrix*)G_L, (const sparse_matrix*)Dinv_L,
- (const real*)workspace->v[0], (const real*)workspace->v[0],
- workspace->v[0], 100 );
- Jacobi_Iter( (const sparse_matrix*)G_U, (const sparse_matrix*)Dinv_U,
- (const real*)workspace->v[0], (const real*)workspace->v[0],
- workspace->v[0], 100 );
+ Jacobi_Iter( (const sparse_matrix*)L, LOWER, (const sparse_matrix*)Dinv_L,
+ (const real*)workspace->v[0], workspace->v[0], 100 );
+ Jacobi_Iter( (const sparse_matrix*)U, UPPER, (const sparse_matrix*)Dinv_U,
+ (const real*)workspace->v[0], workspace->v[0], 100 );
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] );
@@ -713,12 +670,10 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
Sparse_MatVec( H, workspace->v[j], workspace->v[j+1] );
// TODO: add parameters to config file
// TODO: cache results and use for inital guess?
- Jacobi_Iter( (const sparse_matrix*)G_L, (const sparse_matrix*)Dinv_L,
- (const real*)workspace->v[j+1], (const real*)workspace->v[j+1],
- workspace->v[j+1], 100 );
- Jacobi_Iter( (const sparse_matrix*)G_U, (const sparse_matrix*)Dinv_U,
- (const real*)workspace->v[j+1], (const real*)workspace->v[j+1],
- workspace->v[j+1], 100 );
+ Jacobi_Iter( (const sparse_matrix*)L, LOWER, (const sparse_matrix*)Dinv_L,
+ (const real*)workspace->v[j+1], workspace->v[j+1], 100 );
+ Jacobi_Iter( (const sparse_matrix*)U, UPPER, (const sparse_matrix*)Dinv_U,
+ (const real*)workspace->v[j+1], workspace->v[j+1], 100 );
/* apply modified Gram-Schmidt to orthogonalize the new residual */
for( i = 0; i < j-1; i++ ) workspace->h[i][j] = 0;
@@ -799,7 +754,11 @@ int PGMRES_Jacobi( static_storage *workspace, sparse_matrix *H, real *b, real to
// itr, j, fabs( workspace->g[j] ) / bnorm );
// data->timing.matvec += itr * RESTART + j;
- if( itr >= MAX_ITR ) {
+ Deallocate_Matrix( Dinv_U );
+ Deallocate_Matrix( Dinv_L );
+
+ if( itr >= MAX_ITR )
+ {
fprintf( stderr, "GMRES convergence failed\n" );
// return -1;
return itr * (RESTART+1) + j + 1;
diff --git a/puremd_rc_1003/sPuReMD/GMRES.h b/puremd_rc_1003/sPuReMD/GMRES.h
index b624c4f6aac827ccdb9a50e1536cab3c252d0757..c5db381f95387c99140f6b9a402390dd52c3131b 100644
--- a/puremd_rc_1003/sPuReMD/GMRES.h
+++ b/puremd_rc_1003/sPuReMD/GMRES.h
@@ -26,6 +26,12 @@
#include "mytypes.h"
+typedef enum
+{
+ LOWER = 0,
+ UPPER = 1,
+} TRIANGULARITY;
+
int GMRES( static_storage*, sparse_matrix*,
real*, real, real*, FILE* );
diff --git a/puremd_rc_1003/sPuReMD/QEq.c b/puremd_rc_1003/sPuReMD/QEq.c
index d3c2e6c4f6ce2db3878c22933a4842669c586a51..ad1c0cbc91e001bb8a46ef0f0eea49017806dc39 100644
--- a/puremd_rc_1003/sPuReMD/QEq.c
+++ b/puremd_rc_1003/sPuReMD/QEq.c
@@ -104,6 +104,7 @@ int Estimate_LU_Fill( sparse_matrix *A, real *droptol )
}
+/* Incomplete Cholesky factorization with thresholding */
void ICHOLT( sparse_matrix *A, real *droptol,
sparse_matrix *L, sparse_matrix *U )
{
@@ -194,18 +195,25 @@ void ICHOLT( sparse_matrix *A, real *droptol,
L->start[i] = Ltop;
//fprintf( stderr, "nnz(L): %d, max: %d\n", Ltop, L->n * 50 );
+ /* U = L^T (Cholesky factorization) */
for( i = 1; i <= U->n; ++i )
+ {
Utop[i] = U->start[i] = U->start[i] + U->start[i-1] + 1;
-
+ }
for( i = 0; i < L->n; ++i )
- for( pj = L->start[i]; pj < L->start[i+1]; ++pj ){
+ {
+ for( pj = L->start[i]; pj < L->start[i+1]; ++pj )
+ {
j = L->entries[pj].j;
U->entries[Utop[j]].j = i;
U->entries[Utop[j]].val = L->entries[pj].val;
Utop[j]++;
}
+ }
//fprintf( stderr, "nnz(U): %d, max: %d\n", Utop[U->n], U->n * 50 );
+
+ free(Utop);
}
@@ -324,10 +332,10 @@ void QEq( reax_system *system, control_params *control, simulation_data *data,
// Print_Linear_System( system, control, workspace, far_nbrs, data->step );
//TODO: add parameters in control file for solver choice and options
- //matvecs = GMRES( workspace, workspace->H,
- // workspace->b_s, control->q_err, workspace->s[0], out_control->log );
- //matvecs += GMRES( workspace, workspace->H,
- // workspace->b_t, control->q_err, workspace->t[0], out_control->log );
+// matvecs = GMRES( workspace, workspace->H,
+// workspace->b_s, control->q_err, workspace->s[0], out_control->log );
+// matvecs += GMRES( workspace, workspace->H,
+// workspace->b_t, control->q_err, workspace->t[0], out_control->log );
//matvecs = GMRES_HouseHolder( workspace, workspace->H,
// workspace->b_s, control->q_err, workspace->s[0], out_control->log );
diff --git a/puremd_rc_1003/sPuReMD/allocate.h b/puremd_rc_1003/sPuReMD/allocate.h
index 60f67b0383824ae1de1c0e9c75d3c3eca4cef859..c7f9bf15dd432fe931025b8fd39815bea577f9f3 100644
--- a/puremd_rc_1003/sPuReMD/allocate.h
+++ b/puremd_rc_1003/sPuReMD/allocate.h
@@ -28,6 +28,8 @@ void Reallocate( reax_system*, static_storage*, list**, int );
int Allocate_Matrix( sparse_matrix**, int, int );
+void Deallocate_Matrix( sparse_matrix* );
+
int Allocate_HBond_List( int, int, int*, int*, list* );
int Allocate_Bond_List( int, int*, list* );
diff --git a/puremd_rc_1003/sPuReMD/analyze.c b/puremd_rc_1003/sPuReMD/analyze.c
index 95d14ee0f0baf090998184f49eda681b5fde8126..ac5cfe3054aa500a9a0dc80354626cdd84cad8a9 100644
--- a/puremd_rc_1003/sPuReMD/analyze.c
+++ b/puremd_rc_1003/sPuReMD/analyze.c
@@ -728,8 +728,9 @@ void Calculate_Drift( reax_system *system, control_params *control,
/* the atom has moved almost half the box size.
exclude it from further drift computations as it might have an
improper contribution due to periodic boudnaries. */
- workspace->x_old[i][0] = workspace->x_old[i][2] =
- workspace->x_old[i][2] = -999999999999.0;
+ workspace->x_old[i][0] = -999999999999.0;
+ workspace->x_old[i][2] = -999999999999.0;
+ workspace->x_old[i][2] = -999999999999.0;
continue;
}
diff --git a/puremd_rc_1003/sPuReMD/param.c b/puremd_rc_1003/sPuReMD/param.c
index 38b80efca82891f75d51966b66ad966de9e4ac64..045e97bc7e5f020d86eb311ef8a0184ea9f14122 100644
--- a/puremd_rc_1003/sPuReMD/param.c
+++ b/puremd_rc_1003/sPuReMD/param.c
@@ -880,8 +880,9 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
tmp[i] = (char*) malloc(sizeof(char)*MAX_LINE);
/* read control parameters file */
- while (!feof(fp)) {
- fgets(s,MAX_LINE,fp);
+// while (!feof(fp)) {
+ while(fgets(s,MAX_LINE,fp))
+ {
c=Tokenize(s,&tmp);
if( strcmp(tmp[0], "simulation_name") == 0 ) {
@@ -1163,6 +1164,11 @@ char Read_Control_File( FILE* fp, reax_system *system, control_params* control,
}
}
+ if(ferror(fp))
+ {
+ fprintf(stderr, "Error reading control file. Terminating.\n");
+ exit( INVALID_INPUT );
+ }
/* determine target T */
if( control->T_mode == 0 )
diff --git a/puremd_rc_1003/sPuReMD/pdb_tools.c b/puremd_rc_1003/sPuReMD/pdb_tools.c
index 2986b7583f1992d99d399792bdceedb108eb02d7..0cd0331bed28c7d4f31c9b7af2c3f026befcf3fb 100644
--- a/puremd_rc_1003/sPuReMD/pdb_tools.c
+++ b/puremd_rc_1003/sPuReMD/pdb_tools.c
@@ -100,9 +100,9 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
/* count number of atoms in the pdb file */
system->N = 0;
- while (!feof(pdb)) {
- s[0] = 0;
- fgets( s, MAX_LINE, pdb );
+ s[0] = 0;
+ while (fgets( s, MAX_LINE, pdb ))
+ {
tmp[0][0] = 0;
c = Tokenize( s, &tmp );
@@ -110,7 +110,14 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
if( strncmp( tmp[0], "ATOM", 4 ) == 0 ||
strncmp( tmp[0], "HETATM", 6 ) == 0 )
(system->N)++;
+ s[0] = 0;
}
+ if(ferror(pdb))
+ {
+ fprintf(stderr, "Error reading PDB file. Terminating.\n");
+ exit( INVALID_INPUT );
+ }
+
fclose(pdb);
#if defined(DEBUG_FOCUS)
fprintf( stderr, "system->N: %d\n", system->N );
@@ -131,7 +138,11 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
/* start reading and processing pdb file */
- pdb = fopen(pdb_file,"r");
+ if ( (pdb = fopen(pdb_file,"r")) == NULL )
+ {
+ fprintf( stderr, "Error opening the pdb file!\n" );
+ exit( FILE_NOT_FOUND_ERR );
+ }
c = 0;
c1 = 0;
@@ -282,7 +293,11 @@ char Write_PDB( reax_system* system, control_params *control,
/* open output pdb file */
sprintf( temp, "%s%d.pdb", control->sim_name, data->step );
- out_control->pdb = fopen( temp, "w" );
+ if ( (out_control->pdb = fopen( temp, "w" )) == NULL )
+ {
+ fprintf( stderr, "Error opening the pdb output file!\n" );
+ exit( FILE_NOT_FOUND_ERR );
+ }
/* Writing Box information */
@@ -426,8 +441,12 @@ char Read_BGF( char* bgf_file, reax_system* system, control_params *control,
workspace->restricted_list[i] = (int*) calloc( MAX_RESTRICT, sizeof(int) );
- /* start reading and processing pdb file */
- bgf = fopen( bgf_file, "r" );
+ /* start reading and processing bgf file */
+ if ( (bgf = fopen( bgf_file, "r" )) == NULL )
+ {
+ fprintf( stderr, "Error opening the bgf file!\n" );
+ exit( FILE_NOT_FOUND_ERR );
+ }
atom_cnt = 0;
token_cnt = 0;
diff --git a/puremd_rc_1003/sPuReMD/testmd.c b/puremd_rc_1003/sPuReMD/testmd.c
index dbe87b36258ba6c5df022f50df754813c8164447..527adca07f3fde49e0a9d5129e75372a4da3e572 100644
--- a/puremd_rc_1003/sPuReMD/testmd.c
+++ b/puremd_rc_1003/sPuReMD/testmd.c
@@ -82,8 +82,16 @@ void Read_System( char *geof, char *ff, char *ctrlf,
{
FILE *ffield, *ctrl;
- ffield = fopen( ff, "r" );
- ctrl = fopen( ctrlf, "r" );
+ if ( (ffield = fopen( ff, "r" )) == NULL )
+ {
+ fprintf( stderr, "Error opening the ffield file!\n" );
+ exit( FILE_NOT_FOUND_ERR );
+ }
+ if ( (ctrl = fopen( ctrlf, "r" )) == NULL )
+ {
+ fprintf( stderr, "Error opening the ffield file!\n" );
+ exit( FILE_NOT_FOUND_ERR );
+ }
/* ffield file */
Read_Force_Field( ffield, &(system->reaxprm) );
@@ -120,6 +128,12 @@ void Read_System( char *geof, char *ff, char *ctrlf,
}
+static void usage(char* argv[])
+{
+ fprintf(stderr, "usage: ./%s geometry ffield control\n", argv[0]);
+}
+
+
int main(int argc, char* argv[])
{
reax_system system;
@@ -130,6 +144,12 @@ int main(int argc, char* argv[])
output_controls out_control;
evolve_function Evolve;
int steps;
+
+ if( argc != 4 )
+ {
+ usage(argv);
+ exit( INVALID_INPUT );
+ }
lists = (list*) malloc( sizeof(list) * LIST_N );
diff --git a/puremd_rc_1003/sPuReMD/vector.c b/puremd_rc_1003/sPuReMD/vector.c
index be34e137423797d090b31e7eb4027a7cc7db8f2a..d6bed7ae87081656283b860184b00e2c9e8e37c4 100644
--- a/puremd_rc_1003/sPuReMD/vector.c
+++ b/puremd_rc_1003/sPuReMD/vector.c
@@ -22,115 +22,155 @@
#include "vector.h"
-inline int Vector_isZero( real* v, int k )
+inline int Vector_isZero( const real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- if( fabs( v[k] ) > ALMOST_ZERO )
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ if( fabs( v[i] ) > ALMOST_ZERO )
+ {
return 0;
+ }
+ }
return 1;
}
-inline void Vector_MakeZero( real *v, int k )
+inline void Vector_MakeZero( real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- v[k] = 0;
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ v[i] = ZERO;
+ }
}
-inline void Vector_Copy( real* dest, real* v, int k )
+inline void Vector_Copy( real * const dest, const real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- dest[k] = v[k];
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ dest[i] = v[i];
+ }
}
-inline void Vector_Scale( real* dest, real c, real* v, int k )
+inline void Vector_Scale( real * const dest, const real c, const real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- dest[k] = c * v[k];
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ dest[i] = c * v[i];
+ }
}
-inline void Vector_Sum( real* dest, real c, real* v, real d, real* y, int k )
+inline void Vector_Sum( real * const dest, const real c, const real * const v, const real d,
+ const real * const y, const unsigned int k )
{
- for( --k; k>=0; --k )
- dest[k] = c * v[k] + d * y[k];
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ dest[i] = c * v[i] + d * y[i];
+ }
}
-inline void Vector_Add( real* dest, real c, real* v, int k )
+inline void Vector_Add( real * const dest, const real c, const real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- dest[k] += c * v[k];
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ dest[i] += c * v[i];
+ }
}
-inline void Vector_Add2( real* dest, real* v, int k )
+inline void Vector_Add2( real * const dest, const real * const v, const unsigned int k )
{
- for( --k; k>=0; --k )
- dest[k] += v[k];
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ dest[i] += v[i];
+ }
}
-void Vector_Print( FILE *fout, char *vname, real *v, int k )
+void Vector_Print( FILE * const fout, const char * const vname, const real * const v,
+ const unsigned int k )
{
- int i;
+ unsigned int i;
fprintf( fout, "%s:\n", vname );
for( i = 0; i < k; ++i )
+ {
fprintf( fout, "%24.15e\n", v[i] );
+ }
fprintf( fout, "\n" );
}
-inline real Dot( real* v1, real* v2, int k )
+inline real Dot( const real * const v1, const real * const v2, const unsigned int k )
{
- real ret = 0;
-
- for( --k; k>=0; --k )
- ret += v1[k] * v2[k];
+ real ret = ZERO;
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ ret += v1[i] * v2[i];
+ }
return ret;
}
-inline real Norm( real* v1, int k )
+inline real Norm( const real * const v1, const unsigned int k )
{
- real ret = 0;
-
- for( --k; k>=0; --k )
- ret += SQR( v1[k] );
+ real ret = ZERO;
+ unsigned int i;
+
+ for( i=0; i<k; ++i )
+ {
+ ret += SQR( v1[i] );
+ }
return SQRT( ret );
}
-inline void rvec_Copy( rvec dest, rvec src )
+inline void rvec_Copy( rvec dest, const rvec src )
{
dest[0] = src[0], dest[1] = src[1], dest[2] = src[2];
}
-inline void rvec_Scale( rvec ret, real c, rvec v )
+inline void rvec_Scale( rvec ret, const real c, const rvec v )
{
ret[0] = c * v[0], ret[1] = c * v[1], ret[2] = c * v[2];
}
-inline void rvec_Add( rvec ret, rvec v )
+inline void rvec_Add( rvec ret, const rvec v )
{
ret[0] += v[0], ret[1] += v[1], ret[2] += v[2];
}
-inline void rvec_ScaledAdd( rvec ret, real c, rvec v )
+inline void rvec_ScaledAdd( rvec ret, const real c, const rvec v )
{
ret[0] += c * v[0], ret[1] += c * v[1], ret[2] += c * v[2];
}
-inline void rvec_Sum( rvec ret, rvec v1 ,rvec v2 )
+inline void rvec_Sum( rvec ret, const rvec v1 , const rvec v2 )
{
ret[0] = v1[0] + v2[0];
ret[1] = v1[1] + v2[1];
@@ -138,7 +178,8 @@ inline void rvec_Sum( rvec ret, rvec v1 ,rvec v2 )
}
-inline void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 )
+inline void rvec_ScaledSum( rvec ret, const real c1, const rvec v1,
+ const real c2, const rvec v2 )
{
ret[0] = c1 * v1[0] + c2 * v2[0];
ret[1] = c1 * v1[1] + c2 * v2[1];
@@ -146,19 +187,20 @@ inline void rvec_ScaledSum( rvec ret, real c1, rvec v1 ,real c2, rvec v2 )
}
-inline real rvec_Dot( rvec v1, rvec v2 )
+inline real rvec_Dot( const rvec v1, const rvec v2 )
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
-inline real rvec_ScaledDot( real c1, rvec v1, real c2, rvec v2 )
+inline real rvec_ScaledDot( const real c1, const rvec v1,
+ const real c2, const rvec v2 )
{
return (c1*c2) * (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
}
-inline void rvec_Multiply( rvec r, rvec v1, rvec v2 )
+inline void rvec_Multiply( rvec r, const rvec v1, const rvec v2 )
{
r[0] = v1[0] * v2[0];
r[1] = v1[1] * v2[1];
@@ -166,7 +208,7 @@ inline void rvec_Multiply( rvec r, rvec v1, rvec v2 )
}
-inline void rvec_iMultiply( rvec r, ivec v1, rvec v2 )
+inline void rvec_iMultiply( rvec r, const ivec v1, const rvec v2 )
{
r[0] = v1[0] * v2[0];
r[1] = v1[1] * v2[1];
@@ -174,7 +216,7 @@ inline void rvec_iMultiply( rvec r, ivec v1, rvec v2 )
}
-inline void rvec_Divide( rvec r, rvec v1, rvec v2 )
+inline void rvec_Divide( rvec r, const rvec v1, const rvec v2 )
{
r[0] = v1[0] / v2[0];
r[1] = v1[1] / v2[1];
@@ -182,7 +224,7 @@ inline void rvec_Divide( rvec r, rvec v1, rvec v2 )
}
-inline void rvec_iDivide( rvec r, rvec v1, ivec v2 )
+inline void rvec_iDivide( rvec r, const rvec v1, const ivec v2 )
{
r[0] = v1[0] / v2[0];
r[1] = v1[1] / v2[1];
@@ -190,7 +232,7 @@ inline void rvec_iDivide( rvec r, rvec v1, ivec v2 )
}
-inline void rvec_Invert( rvec r, rvec v )
+inline void rvec_Invert( rvec r, const rvec v )
{
r[0] = 1. / v[0];
r[1] = 1. / v[1];
@@ -198,7 +240,7 @@ inline void rvec_Invert( rvec r, rvec v )
}
-inline void rvec_Cross( rvec ret, rvec v1, rvec v2 )
+inline void rvec_Cross( rvec ret, const rvec v1, const rvec v2 )
{
ret[0] = v1[1] * v2[2] - v1[2] * v2[1];
ret[1] = v1[2] * v2[0] - v1[0] * v2[2];
@@ -206,34 +248,40 @@ inline void rvec_Cross( rvec ret, rvec v1, rvec v2 )
}
-inline void rvec_OuterProduct( rtensor r, rvec v1, rvec v2 )
+inline void rvec_OuterProduct( rtensor r, const rvec v1, const rvec v2 )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
r[i][j] = v1[i] * v2[j];
+ }
+ }
}
-inline real rvec_Norm_Sqr( rvec v )
+inline real rvec_Norm_Sqr( const rvec v )
{
return SQR(v[0]) + SQR(v[1]) + SQR(v[2]);
}
-inline real rvec_Norm( rvec v )
+inline real rvec_Norm( const rvec v )
{
return SQRT( SQR(v[0]) + SQR(v[1]) + SQR(v[2]) );
}
-inline int rvec_isZero( rvec v )
+inline int rvec_isZero( const rvec v )
{
if( fabs(v[0]) > ALMOST_ZERO ||
fabs(v[1]) > ALMOST_ZERO ||
fabs(v[2]) > ALMOST_ZERO )
+ {
return 0;
+ }
return 1;
}
@@ -254,7 +302,7 @@ inline void rvec_Random( rvec v )
inline void rtensor_Multiply( rtensor ret, rtensor m1, rtensor m2 )
{
- int i, j, k;
+ unsigned int i, j, k;
rtensor temp;
// check if the result matrix is the same as one of m1, m2.
@@ -287,36 +335,46 @@ inline void rtensor_Multiply( rtensor ret, rtensor m1, rtensor m2 )
}
-inline void rtensor_MatVec( rvec ret, rtensor m, rvec v )
+inline void rtensor_MatVec( rvec ret, rtensor m, const rvec v )
{
- int i;
+ unsigned int i;
rvec temp;
// if ret is the same vector as v, we cannot modify the
// contents of v until all computation is finished.
if( ret == v )
+ {
+ for( i = 0; i < 3; ++i )
{
- for( i = 0; i < 3; ++i )
- temp[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
+ temp[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
+ }
- for( i = 0; i < 3; ++i )
- ret[i] = temp[i];
+ for( i = 0; i < 3; ++i )
+ {
+ ret[i] = temp[i];
}
+ }
else
+ {
+ for( i = 0; i < 3; ++i )
{
- for( i = 0; i < 3; ++i )
- ret[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
+ ret[i] = m[i][0] * v[0] + m[i][1] * v[1] + m[i][2] * v[2];
}
+ }
}
-inline void rtensor_Scale( rtensor ret, real c, rtensor m )
+inline void rtensor_Scale( rtensor ret, const real c, rtensor m )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] = c * m[i][j];
+ }
+ }
}
@@ -325,49 +383,69 @@ inline void rtensor_Add( rtensor ret, rtensor t )
int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] += t[i][j];
+ }
+ }
}
-inline void rtensor_ScaledAdd( rtensor ret, real c, rtensor t )
+inline void rtensor_ScaledAdd( rtensor ret, const real c, rtensor t )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] += c * t[i][j];
+ }
+ }
}
inline void rtensor_Sum( rtensor ret, rtensor t1, rtensor t2 )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] = t1[i][j] + t2[i][j];
+ }
+ }
}
-inline void rtensor_ScaledSum( rtensor ret, real c1, rtensor t1,
- real c2, rtensor t2 )
+inline void rtensor_ScaledSum( rtensor ret, const real c1, rtensor t1,
+ const real c2, rtensor t2 )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] = c1 * t1[i][j] + c2 * t2[i][j];
+ }
+ }
}
inline void rtensor_Copy( rtensor ret, rtensor t )
{
- int i, j;
+ unsigned int i, j;
for( i = 0; i < 3; ++i )
+ {
for( j = 0; j < 3; ++j )
+ {
ret[i][j] = t[i][j];
+ }
+ }
}
@@ -409,9 +487,9 @@ inline real rtensor_Trace( rtensor t )
}
-void Print_rTensor(FILE* fp, rtensor t)
+void Print_rTensor(FILE * const fp, rtensor t)
{
- int i, j;
+ unsigned int i, j;
for (i=0; i < 3; i++)
{
@@ -429,13 +507,13 @@ inline void ivec_MakeZero( ivec v )
}
-inline void ivec_Copy( ivec dest, ivec src )
+inline void ivec_Copy( ivec dest, const ivec src )
{
dest[0] = src[0], dest[1] = src[1], dest[2] = src[2];
}
-inline void ivec_Scale( ivec dest, real C, ivec src )
+inline void ivec_Scale( ivec dest, const real C, const ivec src )
{
dest[0] = C * src[0];
dest[1] = C * src[1];
@@ -443,7 +521,7 @@ inline void ivec_Scale( ivec dest, real C, ivec src )
}
-inline void ivec_rScale( ivec dest, real C, rvec src )
+inline void ivec_rScale( ivec dest, const real C, const rvec src )
{
dest[0] = (int)(C * src[0]);
dest[1] = (int)(C * src[1]);
@@ -451,24 +529,28 @@ inline void ivec_rScale( ivec dest, real C, rvec src )
}
-inline int ivec_isZero( ivec v )
+inline int ivec_isZero( const ivec v )
{
if( v[0]==0 && v[1]==0 && v[2]==0 )
+ {
return 1;
+ }
return 0;
}
-inline int ivec_isEqual( ivec v1, ivec v2 )
+inline int ivec_isEqual( const ivec v1, const ivec v2 )
{
if( v1[0]==v2[0] && v1[1]==v2[1] && v1[2]==v2[2] )
+ {
return 1;
+ }
return 0;
}
-inline void ivec_Sum( ivec dest, ivec v1, ivec v2 )
+inline void ivec_Sum( ivec dest, const ivec v1, const ivec v2 )
{
dest[0] = v1[0] + v2[0];
dest[1] = v1[1] + v2[1];
diff --git a/puremd_rc_1003/sPuReMD/vector.h b/puremd_rc_1003/sPuReMD/vector.h
index 09be45dc9347f22fbcc33682f834c3fc5179cd1f..58b94180ec55fc380e24aa95af1c25de6f5cbbd0 100644
--- a/puremd_rc_1003/sPuReMD/vector.h
+++ b/puremd_rc_1003/sPuReMD/vector.h
@@ -24,61 +24,61 @@
#include "mytypes.h"
-int Vector_isZero( real*, int );
-void Vector_MakeZero( real*, int );
-void Vector_Copy( real*, real*, int );
-void Vector_Scale( real*, real, real*, int );
-void Vector_Sum( real*, real, real*, real, real*, int );
-void Vector_Add( real*, real, real*, int );
-void Vector_Add2( real*, real*, int );
-void Vector_Print( FILE*, char*, real*, int );
-real Dot( real*, real*, int );
-real Norm( real*, int );
+int Vector_isZero( const real * const, const unsigned int );
+void Vector_MakeZero( real * const, const unsigned int );
+void Vector_Copy( real * const, const real * const, const unsigned int );
+void Vector_Scale( real * const, const real, const real * const, const unsigned int );
+void Vector_Sum( real * const, const real, const real * const, const real, const real * const, const unsigned int );
+void Vector_Add( real * const, const real, const real * const, const unsigned int );
+void Vector_Add2( real * const, const real * const, const unsigned int );
+void Vector_Print( FILE * const, const char * const, const real * const, const unsigned int );
+real Dot( const real * const, const real * const, const unsigned int );
+real Norm( const real * const, const unsigned int );
-void rvec_Copy( rvec, rvec );
-void rvec_Scale( rvec, real, rvec );
-void rvec_Add( rvec, rvec );
-void rvec_ScaledAdd( rvec, real, rvec );
-void rvec_Sum( rvec, rvec, rvec );
-void rvec_ScaledSum( rvec, real, rvec, real, rvec );
-real rvec_Dot( rvec, rvec );
-real rvec_ScaledDot( real, rvec, real, rvec );
-void rvec_Multiply( rvec, rvec, rvec );
-void rvec_iMultiply( rvec, ivec, rvec );
-void rvec_Divide( rvec, rvec, rvec );
-void rvec_iDivide( rvec, rvec, ivec );
-void rvec_Invert( rvec, rvec );
-void rvec_Cross( rvec, rvec, rvec );
-void rvec_OuterProduct( rtensor, rvec, rvec );
-real rvec_Norm_Sqr( rvec );
-real rvec_Norm( rvec );
-int rvec_isZero( rvec );
+void rvec_Copy( rvec, const rvec );
+void rvec_Scale( rvec, const real, const rvec );
+void rvec_Add( rvec, const rvec );
+void rvec_ScaledAdd( rvec, const real, const rvec );
+void rvec_Sum( rvec, const rvec, const rvec );
+void rvec_ScaledSum( rvec, const real, const rvec, const real, const rvec );
+real rvec_Dot( const rvec, const rvec );
+real rvec_ScaledDot( const real, const rvec, const real, const rvec );
+void rvec_Multiply( rvec, const rvec, const rvec );
+void rvec_iMultiply( rvec, const ivec, const rvec );
+void rvec_Divide( rvec, const rvec, const rvec );
+void rvec_iDivide( rvec, const rvec, const ivec );
+void rvec_Invert( rvec, const rvec );
+void rvec_Cross( rvec, const rvec, const rvec );
+void rvec_OuterProduct( rtensor, const rvec, const rvec );
+real rvec_Norm_Sqr( const rvec );
+real rvec_Norm( const rvec );
+int rvec_isZero( const rvec );
void rvec_MakeZero( rvec );
void rvec_Random( rvec );
void rtensor_MakeZero( rtensor );
void rtensor_Multiply( rtensor, rtensor, rtensor );
-void rtensor_MatVec( rvec, rtensor, rvec );
-void rtensor_Scale( rtensor, real, rtensor );
+void rtensor_MatVec( rvec, rtensor, const rvec );
+void rtensor_Scale( rtensor, const real, rtensor );
void rtensor_Add( rtensor, rtensor );
-void rtensor_ScaledAdd( rtensor, real, rtensor );
+void rtensor_ScaledAdd( rtensor, const real, rtensor );
void rtensor_Sum( rtensor, rtensor, rtensor );
-void rtensor_ScaledSum( rtensor, real, rtensor, real, rtensor );
-void rtensor_Scale( rtensor, real, rtensor );
+void rtensor_ScaledSum( rtensor, const real, rtensor, const real, rtensor );
+void rtensor_Scale( rtensor, const real, rtensor );
void rtensor_Copy( rtensor, rtensor );
void rtensor_Identity( rtensor );
void rtensor_Transpose( rtensor, rtensor );
real rtensor_Det( rtensor );
real rtensor_Trace( rtensor );
-void Print_rTensor(FILE*,rtensor);
+void Print_rTensor(FILE * const, rtensor);
-int ivec_isZero( ivec );
-int ivec_isEqual( ivec, ivec );
+int ivec_isZero( const ivec );
+int ivec_isEqual( const ivec, const ivec );
void ivec_MakeZero( ivec );
-void ivec_Copy( ivec, ivec );
-void ivec_Scale( ivec, real, ivec );
-void ivec_rScale( ivec, real, rvec );
-void ivec_Sum( ivec, ivec, ivec );
+void ivec_Copy( ivec, const ivec );
+void ivec_Scale( ivec, const real, const ivec );
+void ivec_rScale( ivec, const real, const rvec );
+void ivec_Sum( ivec, const ivec, const ivec );
#endif