From c0f8499b0411ba8770eda055730cca1437980bf5 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Mon, 21 Aug 2017 22:04:00 -0400
Subject: [PATCH] sPuReMD: fix linear solver to work with no preconditioner.
Fix issue when precondition recomputation frequency is zero. Other misc.
changes.
---
sPuReMD/configure.ac | 4 +-
sPuReMD/src/charges.c | 699 +++++++++++++++++++++-------------------
sPuReMD/src/geo_tools.c | 4 +
sPuReMD/src/lin_alg.c | 231 ++++++-------
sPuReMD/src/lin_alg.h | 2 +-
sPuReMD/src/mytypes.h | 4 +-
sPuReMD/src/tool_box.c | 8 +
sPuReMD/src/tool_box.h | 3 -
8 files changed, 495 insertions(+), 460 deletions(-)
diff --git a/sPuReMD/configure.ac b/sPuReMD/configure.ac
index 578c00d4..57b6b8ea 100644
--- a/sPuReMD/configure.ac
+++ b/sPuReMD/configure.ac
@@ -60,7 +60,7 @@ if test "x$ax_cv_c_compiler_vendor" = "xintel"; then
fi
# Check for OpenMP support.
-if test "x$BUILD_OPENMP" = "xyes"; then
+if test "x${BUILD_OPENMP}" = "xyes"; then
AC_OPENMP
if test "x${OPENMP_CFLAGS}" = "x"; then
AC_MSG_WARN([
@@ -70,7 +70,7 @@ if test "x$BUILD_OPENMP" = "xyes"; then
-----------------------------------------------])
else
# bug due to recent Intel compiler change (?)
- if test "x$ax_cv_c_compiler_vendor" = "xintel"; then
+ if test "x${ax_cv_c_compiler_vendor}" = "xintel"; then
OPENMP_CFLAGS="-qopenmp"
fi
AC_SUBST(AM_CFLAGS, "$OPENMP_CFLAGS")
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 01279242..6e663d91 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -32,6 +32,13 @@
#endif
+typedef struct
+{
+ unsigned int j;
+ real val;
+} sparse_matrix_entry;
+
+
#if defined(TEST_MAT)
static sparse_matrix * create_test_mat( void )
{
@@ -85,7 +92,7 @@ static sparse_matrix * create_test_mat( void )
static int compare_matrix_entry(const void *v1, const void *v2)
{
/* larger element has larger column index */
- return *(unsigned int *)v1 - *(unsigned int *)v2;
+ return ((sparse_matrix_entry *)v1)->j - ((sparse_matrix_entry *)v2)->j;
}
@@ -97,51 +104,41 @@ static int compare_matrix_entry(const void *v1, const void *v2)
*/
static void Sort_Matrix_Rows( sparse_matrix * const A )
{
- unsigned int i, j, k, si, ei, *temp_j;
- real *temp_val;
+ unsigned int i, j, si, ei;
+ sparse_matrix_entry *temp;
- #pragma omp parallel default(none) private(i, j, k, si, ei, temp_j, temp_val) shared(stderr)
+// #pragma omp parallel default(none) private(i, j, si, ei, temp) shared(stderr)
{
- if ( ( temp_j = (unsigned int*) malloc( A->n * sizeof(unsigned int)) ) == NULL
- || ( temp_val = (real*) malloc( A->n * sizeof(real)) ) == NULL )
+ if ( ( temp = (sparse_matrix_entry *) malloc( A->n * sizeof(sparse_matrix_entry)) ) == NULL )
{
fprintf( stderr, "Not enough space for matrix row sort. Terminating...\n" );
exit( INSUFFICIENT_MEMORY );
}
/* sort each row of A using column indices */
- #pragma omp for schedule(guided)
+// #pragma omp for schedule(guided)
for ( i = 0; i < A->n; ++i )
{
si = A->start[i];
ei = A->start[i + 1];
- memcpy( temp_j, A->j + si, sizeof(unsigned int) * (ei - si) );
- memcpy( temp_val, A->val + si, sizeof(real) * (ei - si) );
- //TODO: consider implementing single custom one-pass sort instead of using qsort + manual sort
+ for ( j = 0; j < (ei - si); ++j )
+ {
+ (temp + j)->j = A->j[si + j];
+ (temp + j)->val = A->val[si + j];
+ }
+
/* polymorphic sort in standard C library using column indices */
- qsort( temp_j, ei - si, sizeof(unsigned int), compare_matrix_entry );
+ qsort( temp, ei - si, sizeof(sparse_matrix_entry), compare_matrix_entry );
- /* manually sort vals */
for ( j = 0; j < (ei - si); ++j )
{
- for ( k = 0; k < (ei - si); ++k )
- {
- if ( A->j[si + j] == temp_j[k] )
- {
- A->val[si + k] = temp_val[j];
- break;
- }
-
- }
+ A->j[si + j] = (temp + j)->j;
+ A->val[si + j] = (temp + j)->val;
}
-
- /* copy sorted column indices */
- memcpy( A->j + si, temp_j, sizeof(unsigned int) * (ei - si) );
}
- free( temp_val );
- free( temp_j );
+ free( temp );
}
}
@@ -1500,45 +1497,49 @@ static void Compute_Preconditioner_QEq( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- data->timing.cm_solver_pre_comp +=
- diag_pre_comp( Hptr, workspace->Hdia_inv );
- break;
+ case NONE_PC:
+ break;
- case ICHOLT_PC:
- data->timing.cm_solver_pre_comp +=
- ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
- break;
+ case DIAG_PC:
+ data->timing.cm_solver_pre_comp +=
+ diag_pre_comp( Hptr, workspace->Hdia_inv );
+ break;
- case ILU_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
- break;
+ case ICHOLT_PC:
+ data->timing.cm_solver_pre_comp +=
+ ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
+ break;
- case ILUT_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
- workspace->L, workspace->U );
- break;
+ case ILU_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
+ break;
- case ILU_SUPERLU_MT_PC:
+ case ILUT_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
+ workspace->L, workspace->U );
+ break;
+
+ case ILU_SUPERLU_MT_PC:
#if defined(HAVE_SUPERLU_MT)
- data->timing.cm_solver_pre_comp +=
- SuperLU_Factorize( Hptr, workspace->L, workspace->U );
+ data->timing.cm_solver_pre_comp +=
+ SuperLU_Factorize( Hptr, workspace->L, workspace->U );
#else
- fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
- exit( INVALID_INPUT );
+ fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
+ exit( INVALID_INPUT );
#endif
- break;
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
#if defined(DEBUG)
- if ( control->cm_solver_pre_comp_type != DIAG_PC )
+ if ( control->cm_solver_pre_comp_type != NONE_PC &&
+ control->cm_solver_pre_comp_type != DIAG_PC )
{
fprintf( stderr, "condest = %f\n", condest(workspace->L, workspace->U) );
@@ -1840,47 +1841,51 @@ static void Compute_Preconditioner_EE( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- data->timing.cm_solver_pre_comp +=
- diag_pre_comp( Hptr, workspace->Hdia_inv );
- break;
+ case NONE_PC:
+ break;
- case ICHOLT_PC:
- data->timing.cm_solver_pre_comp +=
- ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
- break;
+ case DIAG_PC:
+ data->timing.cm_solver_pre_comp +=
+ diag_pre_comp( Hptr, workspace->Hdia_inv );
+ break;
- case ILU_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
- break;
+ case ICHOLT_PC:
+ data->timing.cm_solver_pre_comp +=
+ ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
+ break;
- case ILUT_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
- workspace->L, workspace->U );
- break;
+ case ILU_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
+ break;
- case ILU_SUPERLU_MT_PC:
+ case ILUT_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
+ workspace->L, workspace->U );
+ break;
+
+ case ILU_SUPERLU_MT_PC:
#if defined(HAVE_SUPERLU_MT)
- data->timing.cm_solver_pre_comp +=
- SuperLU_Factorize( Hptr, workspace->L, workspace->U );
+ data->timing.cm_solver_pre_comp +=
+ SuperLU_Factorize( Hptr, workspace->L, workspace->U );
#else
- fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
- exit( INVALID_INPUT );
+ fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
+ exit( INVALID_INPUT );
#endif
- break;
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
#if defined(DEBUG)
- if ( control->cm_solver_pre_comp_type != DIAG_PC )
+ if ( control->cm_solver_pre_comp_type != NONE_PC &&
+ control->cm_solver_pre_comp_type != DIAG_PC )
{
fprintf( stderr, "condest = %f\n", condest(workspace->L, workspace->U) );
@@ -1939,48 +1944,52 @@ static void Compute_Preconditioner_ACKS2( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- data->timing.cm_solver_pre_comp +=
- diag_pre_comp( Hptr, workspace->Hdia_inv );
- break;
+ case NONE_PC:
+ break;
- case ICHOLT_PC:
- data->timing.cm_solver_pre_comp +=
- ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
- break;
+ case DIAG_PC:
+ data->timing.cm_solver_pre_comp +=
+ diag_pre_comp( Hptr, workspace->Hdia_inv );
+ break;
- case ILU_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
- break;
+ case ICHOLT_PC:
+ data->timing.cm_solver_pre_comp +=
+ ICHOLT( Hptr, workspace->droptol, workspace->L, workspace->U );
+ break;
- case ILUT_PAR_PC:
- data->timing.cm_solver_pre_comp +=
- ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
- workspace->L, workspace->U );
- break;
+ case ILU_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILU_PAR( Hptr, control->cm_solver_pre_comp_sweeps, workspace->L, workspace->U );
+ break;
+
+ case ILUT_PAR_PC:
+ data->timing.cm_solver_pre_comp +=
+ ILUT_PAR( Hptr, workspace->droptol, control->cm_solver_pre_comp_sweeps,
+ workspace->L, workspace->U );
+ break;
- case ILU_SUPERLU_MT_PC:
+ case ILU_SUPERLU_MT_PC:
#if defined(HAVE_SUPERLU_MT)
- data->timing.cm_solver_pre_comp +=
- SuperLU_Factorize( Hptr, workspace->L, workspace->U );
+ data->timing.cm_solver_pre_comp +=
+ SuperLU_Factorize( Hptr, workspace->L, workspace->U );
#else
- fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
- exit( INVALID_INPUT );
+ fprintf( stderr, "SuperLU MT support disabled. Re-compile before enabling. Terminating...\n" );
+ exit( INVALID_INPUT );
#endif
- break;
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
Hptr->val[Hptr->start[system->N_cm] - 1] = 0.0;
#if defined(DEBUG)
- if ( control->cm_solver_pre_comp_type != DIAG_PC )
+ if ( control->cm_solver_pre_comp_type != NONE_PC ||
+ control->cm_solver_pre_comp_type != DIAG_PC )
{
fprintf( stderr, "condest = %f\n", condest(workspace->L, workspace->U) );
@@ -2030,109 +2039,112 @@ static void Setup_Preconditioner_QEq( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- if ( workspace->Hdia_inv == NULL )
- {
- if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
+ case NONE_PC:
+ break;
+
+ case DIAG_PC:
+ if ( workspace->Hdia_inv == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- break;
+ break;
- case ICHOLT_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ICHOLT_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
+ fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
- fprintf( stderr, "fillin = %d\n", fillin );
- fprintf( stderr, "allocated memory: L = U = %ldMB\n",
- fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
+ fprintf( stderr, "fillin = %d\n", fillin );
+ fprintf( stderr, "allocated memory: L = U = %ldMB\n",
+ fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
- if ( workspace->L == NULL )
- {
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ }
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILU_PAR_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_PAR_PC:
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILUT_PAR_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ILUT_PAR_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- if ( workspace->L == NULL )
- {
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILU_SUPERLU_MT_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_SUPERLU_MT_PC:
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
}
@@ -2174,109 +2186,112 @@ static void Setup_Preconditioner_EE( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- if ( workspace->Hdia_inv == NULL )
- {
- if ( ( workspace->Hdia_inv = (real *) calloc( system->N_cm, sizeof( real ) ) ) == NULL )
+ case NONE_PC:
+ break;
+
+ case DIAG_PC:
+ if ( workspace->Hdia_inv == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( ( workspace->Hdia_inv = (real *) calloc( system->N_cm, sizeof( real ) ) ) == NULL )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- break;
+ break;
- case ICHOLT_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ICHOLT_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
+ fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
- fprintf( stderr, "fillin = %d\n", fillin );
- fprintf( stderr, "allocated memory: L = U = %ldMB\n",
- fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
+ fprintf( stderr, "fillin = %d\n", fillin );
+ fprintf( stderr, "allocated memory: L = U = %ldMB\n",
+ fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
- if ( workspace->L == NULL )
- {
- if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
- }
+ if ( Allocate_Matrix( &(workspace->L), system->N_cm, fillin + system->N_cm ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), system->N_cm, fillin + system->N_cm ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ }
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILU_PAR_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_PAR_PC:
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILUT_PAR_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ILUT_PAR_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- if ( workspace->L == NULL )
- {
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ if ( workspace->L == NULL )
+ {
+ /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
+ }
+ else
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ //TODO: reallocate
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ break;
- case ILU_SUPERLU_MT_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_SUPERLU_MT_PC:
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
@@ -2321,108 +2336,111 @@ static void Setup_Preconditioner_ACKS2( const reax_system * const system,
switch ( control->cm_solver_pre_comp_type )
{
- case DIAG_PC:
- if ( workspace->Hdia_inv == NULL )
- {
- if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
+ case NONE_PC:
+ break;
+
+ case DIAG_PC:
+ if ( workspace->Hdia_inv == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( ( workspace->Hdia_inv = (real *) calloc( Hptr->n, sizeof( real ) ) ) == NULL )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- break;
+ break;
- case ICHOLT_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ICHOLT_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
+ fillin = Estimate_LU_Fill( Hptr, workspace->droptol );
#if defined(DEBUG)
- fprintf( stderr, "fillin = %d\n", fillin );
- fprintf( stderr, "allocated memory: L = U = %ldMB\n",
- fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
+ fprintf( stderr, "fillin = %d\n", fillin );
+ fprintf( stderr, "allocated memory: L = U = %ldMB\n",
+ fillin * (sizeof(real) + sizeof(unsigned int)) / (1024 * 1024) );
#endif
- if ( workspace->L == NULL )
- {
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, fillin ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, fillin ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILU_PAR_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_PAR_PC:
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILUT_PAR_PC:
- Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
+ case ILUT_PAR_PC:
+ Calculate_Droptol( Hptr, workspace->droptol, control->cm_solver_pre_comp_droptol );
#if defined(DEBUG_FOCUS)
- fprintf( stderr, "drop tolerances calculated\n" );
+ fprintf( stderr, "drop tolerances calculated\n" );
#endif
- if ( workspace->L == NULL )
- {
- /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ if ( workspace->L == NULL )
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ /* TODO: safest storage estimate is ILU(0) (same as lower triangular portion of H), could improve later */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ else
+ {
+ //TODO: reallocate
+ }
+ break;
- case ILU_SUPERLU_MT_PC:
- if ( workspace->L == NULL )
- {
- /* factors have sparsity pattern as H */
- if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
- Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ case ILU_SUPERLU_MT_PC:
+ if ( workspace->L == NULL )
+ {
+ /* factors have sparsity pattern as H */
+ if ( Allocate_Matrix( &(workspace->L), Hptr->n, Hptr->m ) == FAILURE ||
+ Allocate_Matrix( &(workspace->U), Hptr->n, Hptr->m ) == FAILURE )
+ {
+ fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
+ }
+ else
{
- fprintf( stderr, "not enough memory for preconditioning matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ //TODO: reallocate
}
- }
- else
- {
- //TODO: reallocate
- }
- break;
+ break;
- default:
- fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner computation method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
Hptr->val[Hptr->start[system->N + 1] - 1] = 0.0;
@@ -2505,7 +2523,8 @@ static void QEq( reax_system * const system, control_params * const control,
iters = GMRES( workspace, control, data, workspace->H,
workspace->b_s, control->cm_solver_q_err, workspace->s[0],
out_control->log,
- ((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) ? TRUE : FALSE );
+ (control->cm_solver_pre_comp_refactor > 0 &&
+ (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) ? TRUE : FALSE );
iters += GMRES( workspace, control, data, workspace->H,
workspace->b_t, control->cm_solver_q_err, workspace->t[0],
out_control->log, FALSE );
@@ -2514,7 +2533,9 @@ static void QEq( reax_system * const system, control_params * const control,
case GMRES_H_S:
iters = GMRES_HouseHolder( workspace, control, data, workspace->H,
workspace->b_s, control->cm_solver_q_err, workspace->s[0],
- out_control->log, (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0 );
+ out_control->log,
+ (control->cm_solver_pre_comp_refactor > 0 &&
+ (data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) ? TRUE : FALSE );
iters += GMRES_HouseHolder( workspace, control, data, workspace->H,
workspace->b_t, control->cm_solver_q_err, workspace->t[0],
out_control->log, 0 );
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 3d7d4ce3..f996af5e 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -503,7 +503,11 @@ char Write_PDB( reax_system* system, list* bonds, simulation_data *data,
"ATOM ", workspace->orig_id[i], p_atom->name, ' ', "REX", ' ', 1, ' ',
p_atom->x[0], p_atom->x[1], p_atom->x[2],
1.0, 0.0, "0", name, " " );
+
+#if defined(DEBUG)
fprintf( stderr, "PDB NAME <%s>\n", p_atom->name );
+#endif
+
strncpy( buffer + i * PDB_ATOM_FORMAT_O_LENGTH, line,
PDB_ATOM_FORMAT_O_LENGTH );
}
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 31bd16e1..118654b7 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -146,7 +146,7 @@ static void Sparse_MatVec( const sparse_matrix * const A,
* A: stored in CSR
* A_t: stored in CSR
*/
-void Transpose( const sparse_matrix const *A, sparse_matrix const *A_t )
+void Transpose( const sparse_matrix * const A, sparse_matrix const *A_t )
{
unsigned int i, j, pj, *A_t_top;
@@ -1022,156 +1022,160 @@ void jacobi_iter( const sparse_matrix * const R, const real * const Dinv,
* Matrices have non-zero diagonals
* Each row of a matrix has at least one non-zero (i.e., no rows with all zeros) */
static void apply_preconditioner( const static_storage * const workspace, const control_params * const control,
- const real * const y, real * const x, const int fresh_pre )
+ const real * const y, real * const x, const int fresh_pre )
{
int i, si;
- switch ( control->cm_solver_pre_app_type )
+ /* no preconditioning */
+ if ( control->cm_solver_pre_comp_type == NONE_PC )
{
- case NONE_PA:
- break;
- case TRI_SOLVE_PA:
- switch ( control->cm_solver_pre_comp_type )
- {
- case DIAG_PC:
- diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
- break;
- case ICHOLT_PC:
- case ILU_PAR_PC:
- case ILUT_PAR_PC:
- tri_solve( workspace->L, y, x, workspace->L->n, LOWER );
- tri_solve( workspace->U, x, x, workspace->U->n, UPPER );
- break;
- default:
- fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
- case TRI_SOLVE_LEVEL_SCHED_PA:
- switch ( control->cm_solver_pre_comp_type )
- {
- case DIAG_PC:
- diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
- break;
- case ICHOLT_PC:
- case ILU_PAR_PC:
- case ILUT_PAR_PC:
- tri_solve_level_sched( workspace->L, y, x, workspace->L->n, LOWER, fresh_pre );
- tri_solve_level_sched( workspace->U, x, x, workspace->U->n, UPPER, fresh_pre );
- break;
- default:
- fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
- }
- break;
- case TRI_SOLVE_GC_PA:
- switch ( control->cm_solver_pre_comp_type )
+ Vector_Copy( x, y, workspace->H->n );
+ }
+ else
+ {
+ switch ( control->cm_solver_pre_app_type )
{
- case DIAG_PC:
- fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
- exit( INVALID_INPUT );
+ case TRI_SOLVE_PA:
+ switch ( control->cm_solver_pre_comp_type )
+ {
+ case DIAG_PC:
+ diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
+ break;
+ case ICHOLT_PC:
+ case ILU_PAR_PC:
+ case ILUT_PAR_PC:
+ tri_solve( workspace->L, y, x, workspace->L->n, LOWER );
+ tri_solve( workspace->U, x, x, workspace->U->n, UPPER );
+ break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
break;
- case ICHOLT_PC:
- case ILU_PAR_PC:
- case ILUT_PAR_PC:
- #pragma omp master
+ case TRI_SOLVE_LEVEL_SCHED_PA:
+ switch ( control->cm_solver_pre_comp_type )
{
- memcpy( y_p, y, sizeof(real) * workspace->H->n );
+ case DIAG_PC:
+ diag_pre_app( workspace->Hdia_inv, y, x, workspace->H->n );
+ break;
+ case ICHOLT_PC:
+ case ILU_PAR_PC:
+ case ILUT_PAR_PC:
+ tri_solve_level_sched( workspace->L, y, x, workspace->L->n, LOWER, fresh_pre );
+ tri_solve_level_sched( workspace->U, x, x, workspace->U->n, UPPER, fresh_pre );
+ break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
}
+ break;
+ case TRI_SOLVE_GC_PA:
+ switch ( control->cm_solver_pre_comp_type )
+ {
+ case DIAG_PC:
+ fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ case ICHOLT_PC:
+ case ILU_PAR_PC:
+ case ILUT_PAR_PC:
+ #pragma omp master
+ {
+ memcpy( y_p, y, sizeof(real) * workspace->H->n );
+ }
- #pragma omp barrier
+ #pragma omp barrier
- permute_vector( y_p, workspace->H->n, FALSE, LOWER );
- tri_solve_level_sched( workspace->L, y_p, x, workspace->L->n, LOWER, fresh_pre );
- tri_solve_level_sched( workspace->U, x, x, workspace->U->n, UPPER, fresh_pre );
- permute_vector( x, workspace->H->n, TRUE, UPPER );
- break;
- default:
- fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
- exit( INVALID_INPUT );
+ permute_vector( y_p, workspace->H->n, FALSE, LOWER );
+ tri_solve_level_sched( workspace->L, y_p, x, workspace->L->n, LOWER, fresh_pre );
+ tri_solve_level_sched( workspace->U, x, x, workspace->U->n, UPPER, fresh_pre );
+ permute_vector( x, workspace->H->n, TRUE, UPPER );
break;
- }
- break;
- case JACOBI_ITER_PA:
- switch ( control->cm_solver_pre_comp_type )
- {
- case DIAG_PC:
- fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
- exit( INVALID_INPUT );
+ default:
+ fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
break;
- case ICHOLT_PC:
- case ILU_PAR_PC:
- case ILUT_PAR_PC:
- #pragma omp master
+ case JACOBI_ITER_PA:
+ switch ( control->cm_solver_pre_comp_type )
{
- if ( Dinv_L == NULL )
+ case DIAG_PC:
+ fprintf( stderr, "Unsupported preconditioner computation/application method combination. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ case ICHOLT_PC:
+ case ILU_PAR_PC:
+ case ILUT_PAR_PC:
+ #pragma omp master
{
- if ( (Dinv_L = (real*) malloc(sizeof(real) * workspace->L->n)) == NULL )
+ if ( Dinv_L == NULL )
{
- fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( (Dinv_L = (real*) malloc(sizeof(real) * workspace->L->n)) == NULL )
+ {
+ fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
}
- }
- #pragma omp barrier
+ #pragma omp barrier
- /* construct D^{-1}_L */
- if ( fresh_pre == TRUE )
- {
- #pragma omp for schedule(static)
- for ( i = 0; i < workspace->L->n; ++i )
+ /* construct D^{-1}_L */
+ if ( fresh_pre == TRUE )
{
- si = workspace->L->start[i + 1] - 1;
- Dinv_L[i] = 1. / workspace->L->val[si];
+ #pragma omp for schedule(static)
+ for ( i = 0; i < workspace->L->n; ++i )
+ {
+ si = workspace->L->start[i + 1] - 1;
+ Dinv_L[i] = 1. / workspace->L->val[si];
+ }
}
- }
- jacobi_iter( workspace->L, Dinv_L, y, x, LOWER, control->cm_solver_pre_app_jacobi_iters );
+ jacobi_iter( workspace->L, Dinv_L, y, x, LOWER, control->cm_solver_pre_app_jacobi_iters );
- #pragma omp master
- {
- if ( Dinv_U == NULL )
+ #pragma omp master
{
- if ( (Dinv_U = (real*) malloc(sizeof(real) * workspace->U->n)) == NULL )
+ if ( Dinv_U == NULL )
{
- fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
- exit( INSUFFICIENT_MEMORY );
+ if ( (Dinv_U = (real*) malloc(sizeof(real) * workspace->U->n)) == NULL )
+ {
+ fprintf( stderr, "not enough memory for Jacobi iteration matrices. terminating.\n" );
+ exit( INSUFFICIENT_MEMORY );
+ }
}
}
- }
- #pragma omp barrier
+ #pragma omp barrier
- /* construct D^{-1}_U */
- if ( fresh_pre == TRUE )
- {
- #pragma omp for schedule(static)
- for ( i = 0; i < workspace->U->n; ++i )
+ /* construct D^{-1}_U */
+ if ( fresh_pre == TRUE )
{
- si = workspace->U->start[i];
- Dinv_U[i] = 1. / workspace->U->val[si];
+ #pragma omp for schedule(static)
+ for ( i = 0; i < workspace->U->n; ++i )
+ {
+ si = workspace->U->start[i];
+ Dinv_U[i] = 1. / workspace->U->val[si];
+ }
}
- }
- jacobi_iter( workspace->U, Dinv_U, y, x, UPPER, control->cm_solver_pre_app_jacobi_iters );
+ jacobi_iter( workspace->U, Dinv_U, y, x, UPPER, control->cm_solver_pre_app_jacobi_iters );
+ break;
+ default:
+ fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
+ exit( INVALID_INPUT );
+ break;
+ }
break;
default:
fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
exit( INVALID_INPUT );
break;
- }
- break;
- default:
- fprintf( stderr, "Unrecognized preconditioner application method. Terminating...\n" );
- exit( INVALID_INPUT );
- break;
+ }
}
-
- return;
}
@@ -1384,7 +1388,8 @@ int GMRES( const static_storage * const workspace, const control_params * const
#pragma omp master
{
time_start = Get_Time( );
- if ( control->cm_solver_pre_comp_type == DIAG_PC )
+ if ( control->cm_solver_pre_comp_type == NONE_PC ||
+ control->cm_solver_pre_comp_type == DIAG_PC )
{
/* Givens rotations on the upper-Hessenberg matrix to make it U */
for ( i = 0; i <= j; i++ )
diff --git a/sPuReMD/src/lin_alg.h b/sPuReMD/src/lin_alg.h
index c62a70b6..7085f5a6 100644
--- a/sPuReMD/src/lin_alg.h
+++ b/sPuReMD/src/lin_alg.h
@@ -32,7 +32,7 @@ typedef enum
} TRIANGULARITY;
-void Transpose( const sparse_matrix const *, sparse_matrix const * );
+void Transpose( const sparse_matrix * const, sparse_matrix const * );
void Transpose_I( sparse_matrix * const );
void tri_solve( const sparse_matrix * const, const real * const,
diff --git a/sPuReMD/src/mytypes.h b/sPuReMD/src/mytypes.h
index 8b079ded..3833dff8 100644
--- a/sPuReMD/src/mytypes.h
+++ b/sPuReMD/src/mytypes.h
@@ -206,12 +206,12 @@ enum solver
enum pre_comp
{
- DIAG_PC = 0, ICHOLT_PC = 1, ILU_PAR_PC = 2, ILUT_PAR_PC = 3, ILU_SUPERLU_MT_PC = 4,
+ NONE_PC = 0, DIAG_PC = 1, ICHOLT_PC = 2, ILU_PAR_PC = 3, ILUT_PAR_PC = 4, ILU_SUPERLU_MT_PC = 5,
};
enum pre_app
{
- NONE_PA = 0, TRI_SOLVE_PA = 1, TRI_SOLVE_LEVEL_SCHED_PA = 2, TRI_SOLVE_GC_PA = 3, JACOBI_ITER_PA = 4,
+ TRI_SOLVE_PA = 0, TRI_SOLVE_LEVEL_SCHED_PA = 1, TRI_SOLVE_GC_PA = 2, JACOBI_ITER_PA = 3,
};
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 776bb565..cc1da28a 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -290,6 +290,8 @@ void Trim_Spaces( char *element )
/************ from system_props.c *************/
real Get_Time( )
{
+ struct timeval tim;
+
gettimeofday(&tim, NULL );
return ( tim.tv_sec + (tim.tv_usec / 1000000.0) );
}
@@ -297,6 +299,9 @@ real Get_Time( )
real Get_Timing_Info( real t_start )
{
+ struct timeval tim;
+ real t_end;
+
gettimeofday(&tim, NULL );
t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
return (t_end - t_start);
@@ -305,6 +310,9 @@ real Get_Timing_Info( real t_start )
void Update_Timing_Info( real *t_start, real *timing )
{
+ struct timeval tim;
+ real t_end;
+
gettimeofday(&tim, NULL );
t_end = tim.tv_sec + (tim.tv_usec / 1000000.0);
*timing += (t_end - *t_start);
diff --git a/sPuReMD/src/tool_box.h b/sPuReMD/src/tool_box.h
index 5712152d..aec68895 100644
--- a/sPuReMD/src/tool_box.h
+++ b/sPuReMD/src/tool_box.h
@@ -24,9 +24,6 @@
#include "mytypes.h"
-struct timeval tim;
-real t_end;
-
/* from box.h */
void Transform( rvec, simulation_box*, char, rvec );
void Transform_to_UnitBox( rvec, simulation_box*, char, rvec );
--
GitLab