From eddce04a89c460bfc4594cac6bea37d8769d5750 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@cse.msu.edu>
Date: Tue, 21 Feb 2017 09:52:57 -0500
Subject: [PATCH] Complete EEM implementation. Fix EEM charge matrix entries.
Other general refactoring.
---
environ/param.gpu.water | 3 +-
sPuReMD/src/allocate.c | 2 +-
sPuReMD/src/charges.c | 33 ++++++++++++++----
sPuReMD/src/forces.c | 74 +++++++++++++++++++++++------------------
sPuReMD/src/geo_tools.c | 6 ++--
sPuReMD/src/init_md.c | 10 +++---
sPuReMD/src/lin_alg.c | 65 +++++++++++++++++++-----------------
7 files changed, 113 insertions(+), 80 deletions(-)
diff --git a/environ/param.gpu.water b/environ/param.gpu.water
index 0b251cf4..bed22006 100644
--- a/environ/param.gpu.water
+++ b/environ/param.gpu.water
@@ -15,7 +15,8 @@ thb_cutoff 0.001 ! cutoff value for three body in
hbond_cutoff 7.50 ! cutoff distance for hydrogen bond interactions (Angstroms)
charge_method 0 ! charge method: 0 = QEq, 1 = EEM, 2 = ACKS2
-cm_solver_type 0 ! iterative linear solver used in charge method
+cm_q_net 0.0 ! net system charge
+cm_solver_type 0 ! iterative linear solver for charge method: 0 = GMRES, 1 = GMRES_H, 2 = CG, 3 = SDM
cm_solver_q_err 1e-6 ! relative residual norm threshold used in solver
cm_domain_sparsity 1.0 ! scalar for scaling cut-off distance, used to sparsify charge matrix (between 0.0 and 1.0)
cm_solver_pre_comp_type 1 ! method used to compute preconditioner, if applicable
diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index c8bdddb0..ed2c07a4 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -280,7 +280,7 @@ void Reallocate( reax_system *system, static_storage *workspace, list **lists,
if ( realloc->Htop > 0 )
{
- Reallocate_Matrix(&(workspace->H), system->N, realloc->Htop * SAFE_ZONE, "H");
+ Reallocate_Matrix(&(workspace->H), system->N_cm, realloc->Htop * SAFE_ZONE, "H");
realloc->Htop = -1;
Deallocate_Matrix( workspace->L );
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index 8ac55b47..d6fc3ff2 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -207,7 +207,8 @@ static void Calculate_Droptol( const sparse_matrix * const A,
#endif
}
- val = A->val[k]; // diagonal entry
+ // diagonal entry
+ val = A->val[k];
#ifdef _OPENMP
droptol_local[tid * A->n + i] += val * val;
#else
@@ -581,10 +582,11 @@ static real diag_pre_comp( const reax_system * const system, real * const Hdia_i
#pragma omp parallel for schedule(static) \
default(none) private(i)
- for ( i = 0; i < system->N_cm; ++i )
+ for ( i = 0; i < system->N; ++i )
{
Hdia_inv[i] = 1.0 / system->reaxprm.sbp[system->atoms[i].type].eta;
}
+ Hdia_inv[system->N_cm - 1] = 1.0;
return Get_Timing_Info( start );
}
@@ -617,7 +619,9 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
memset( U->start, 0, (A->n + 1) * sizeof(unsigned int) );
memset( Utop, 0, A->n * sizeof(unsigned int) );
- //fprintf( stderr, "n: %d\n", A->n );
+// fprintf( stderr, "n: %d\n", A->n );
+// fflush( stderr );
+
for ( i = 0; i < A->n; ++i )
{
L->start[i] = Ltop;
@@ -627,7 +631,9 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
{
j = A->j[pj];
val = A->val[pj];
- //fprintf( stderr, "i: %d, j: %d", i, j );
+
+// fprintf( stderr, " i: %d, j: %d\n", i, j );
+// fflush( stderr );
if ( FABS(val) > droptol[i] )
{
@@ -657,7 +663,9 @@ static real ICHOLT( const sparse_matrix * const A, const real * const droptol,
tmp_val[tmptop] = val;
++tmptop;
}
- //fprintf( stderr, " -- done\n" );
+
+// fprintf( stderr, " -- done\n" );
+// fflush( stderr );
}
// sanity check
@@ -1476,8 +1484,8 @@ static void Init_Charges( const reax_system * const system,
|| workspace->L == NULL))
{
// Print_Linear_System( system, control, workspace, data->step );
- Print_Sparse_Matrix2( workspace->H, "H_0.out" );
- exit(0);
+// Print_Sparse_Matrix2( workspace->H, "H_0.out" );
+// exit(0);
time = Get_Time( );
if ( control->cm_solver_pre_comp_type != DIAG_PC )
@@ -1691,6 +1699,7 @@ static void QEq( reax_system * const system, control_params * const control,
workspace->b_t, control->cm_solver_q_err, workspace->t[0],
out_control->log, FALSE );
break;
+
case GMRES_H_S:
iters = GMRES_HouseHolder( workspace, control, data, workspace->H,
workspace->b_s, control->cm_solver_q_err, workspace->s[0],
@@ -1699,6 +1708,7 @@ static void QEq( reax_system * const system, control_params * const control,
workspace->b_t, control->cm_solver_q_err, workspace->t[0],
out_control->log, 0 );
break;
+
case CG_S:
iters = CG( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
@@ -1711,6 +1721,7 @@ static void QEq( reax_system * const system, control_params * const control,
// workspace->L, workspace->U, workspace->t[0], control->cm_solver_pre_app_type,
// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
break;
+
case SDM_S:
iters = SDM( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
@@ -1723,6 +1734,7 @@ static void QEq( reax_system * const system, control_params * const control,
// workspace->L, workspace->U, workspace->t[0], control->cm_solver_pre_app_type,
// control->cm_solver_pre_app_jacobi_iters, out_control->log ) + 1;
break;
+
default:
fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
exit( INVALID_INPUT );
@@ -1792,20 +1804,24 @@ static void EEM( reax_system * const system, control_params * const control,
out_control->log,
((data->step - data->prev_steps) % control->cm_solver_pre_comp_refactor == 0) ? TRUE : FALSE );
break;
+
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 );
break;
+
case CG_S:
iters = CG( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
break;
+
case SDM_S:
iters = SDM( workspace, workspace->H, workspace->b_s, control->cm_solver_q_err,
workspace->s[0], out_control->log ) + 1;
break;
+
default:
fprintf( stderr, "Unrecognized QEq solver selection. Terminating...\n" );
exit( INVALID_INPUT );
@@ -1834,12 +1850,15 @@ void Compute_Charges( reax_system * const system, control_params * const control
case QEQ_CM:
QEq( system, control, data, workspace, far_nbrs, out_control );
break;
+
case EEM_CM:
EEM( system, control, data, workspace, far_nbrs, out_control );
break;
+
case ACKS2_CM:
//TODO
break;
+
default:
fprintf( stderr, "Invalid charge method. Terminating...\n" );
exit( INVALID_INPUT );
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 91e3430f..5e602253 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -372,7 +372,7 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
/* shielding */
dr3gamij_1 = ( r_ij * r_ij * r_ij +
system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
- dr3gamij_3 = POW( dr3gamij_1 , 0.33333333333333 );
+ dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
break;
@@ -392,33 +392,28 @@ static inline real Init_Charge_Matrix_Entry( reax_system *system,
switch ( pos )
{
case OFF_DIAGONAL:
- Tap = control->Tap7 * r_ij + control->Tap6;
- Tap = Tap * r_ij + control->Tap5;
- Tap = Tap * r_ij + control->Tap4;
- Tap = Tap * r_ij + control->Tap3;
- Tap = Tap * r_ij + control->Tap2;
- Tap = Tap * r_ij + control->Tap1;
- Tap = Tap * r_ij + control->Tap0;
-
-// gamij = system->reaxprm.sbp[system->atoms[i].type].gamma
-// + system->reaxprm.sbp[system->atoms[j].type].gamma;
-// gamij = SQRT( gamij );
-// dr3gamij_1 = ( r_ij * r_ij * r_ij +
-// 1.0 / ( gamij * gamij * gamij ) );
-// dr3gamij_3 = POW( dr3gamij_1 , 0.33333333333333 );
-//
-// ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
-
- /* shielding */
- dr3gamij_1 = ( r_ij * r_ij * r_ij +
- system->reaxprm.tbp[system->atoms[i].type][system->atoms[j].type].gamma );
- dr3gamij_3 = POW( dr3gamij_1 , 0.33333333333333 );
-
- ret = ((i == j) ? 0.5 : 1.0) * Tap * EV_to_KCALpMOL / dr3gamij_3;
+ if ( r_ij < control->r_cut && r_ij > 0.001 )
+ {
+ Tap = control->Tap7 * r_ij + control->Tap6;
+ Tap = Tap * r_ij + control->Tap5;
+ Tap = Tap * r_ij + control->Tap4;
+ Tap = Tap * r_ij + control->Tap3;
+ Tap = Tap * r_ij + control->Tap2;
+ Tap = Tap * r_ij + control->Tap1;
+ Tap = Tap * r_ij + control->Tap0;
+
+ gamij = SQRT( system->reaxprm.sbp[system->atoms[i].type].gamma
+ * system->reaxprm.sbp[system->atoms[j].type].gamma );
+ /* shielding */
+ dr3gamij_1 = POW( r_ij, 3.0 ) + 1.0 / POW( gamij, 3.0 );
+ dr3gamij_3 = POW( dr3gamij_1 , 1.0 / 3.0 );
+
+ ret = Tap * EV_to_KCALpMOL / dr3gamij_3;
+ }
break;
case DIAGONAL:
- ret = 2.0 * system->reaxprm.sbp[system->atoms[i].type].eta;
+ ret = system->reaxprm.sbp[system->atoms[i].type].eta;
break;
default:
@@ -465,19 +460,27 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
break;
case EEM_CM:
- H->start[system->N] = *Htop;
- H_sp->start[system->N] = *H_sp_top;
+ H->start[system->N_cm - 1] = *Htop;
+ H_sp->start[system->N_cm - 1] = *H_sp_top;
- for ( i = 0; i < system->N_cm; ++i )
+ for ( i = 0; i < system->N_cm - 1; ++i )
{
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
- ++(*Htop);
+ *Htop = *Htop + 1;
H_sp->j[*H_sp_top] = i;
H_sp->val[*H_sp_top] = 1.0;
- ++(*H_sp_top);
+ *H_sp_top = *H_sp_top + 1;
}
+
+ H->j[*Htop] = system->N_cm - 1;
+ H->val[*Htop] = 0.0;
+ *Htop = *Htop + 1;
+
+ H_sp->j[*H_sp_top] = system->N_cm - 1;
+ H_sp->val[*H_sp_top] = 0.0;
+ *H_sp_top = *H_sp_top + 1;
break;
case ACKS2_CM:
@@ -538,8 +541,11 @@ void Init_Forces( reax_system *system, control_params *control,
btop_i = End_Index( i, bonds );
sbp_i = &(system->reaxprm.sbp[type_i]);
ihb = ihb_top = -1;
+
if ( control->hb_cut > 0 && (ihb = sbp_i->p_hbond) == 1 )
+ {
ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ }
for ( pj = start_i; pj < end_i; ++pj )
{
@@ -769,14 +775,15 @@ void Init_Forces( reax_system *system, control_params *control,
// i, Start_Index( i, bonds ), End_Index( i, bonds ) );
}
-// printf("Htop = %d\n", Htop);
-// printf("H_sp_top = %d\n", H_sp_top);
+ Init_Charge_Matrix_Remaining_Entries( system, control, H, H_sp, &Htop, &H_sp_top );
// mark the end of j list
- Init_Charge_Matrix_Remaining_Entries( system, control, H, H_sp, &Htop, &H_sp_top );
H->start[system->N_cm] = Htop;
H_sp->start[system->N_cm] = H_sp_top;
+// printf("Htop = %d\n", Htop);
+// printf("H_sp_top = %d\n", H_sp_top);
+
/* validate lists - decide if reallocation is required! */
Validate_Lists( workspace, lists,
data->step, system->N, H->m, Htop, num_bonds, num_hbonds );
@@ -1178,6 +1185,7 @@ void Compute_Forces( reax_system *system, control_params *control,
}
t_elapsed = Get_Timing_Info( t_start );
data->timing.init_forces += t_elapsed;
+
#if defined(DEBUG_FOCUS)
fprintf( stderr, "init_forces - ");
#endif
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 01f404b4..3d7d4ce3 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -381,9 +381,9 @@ char Read_PDB( char* pdb_file, reax_system* system, control_params *control,
// system->my_atoms[top].q, occupancy, temp_factor,
// seg_id, element );
- //fprintf( stderr, "atom( %8.3f %8.3f %8.3f ) --> p%d\n",
- // system->my_atoms[top].x[0], system->my_atoms[top].x[1],
- // system->my_atoms[top].x[2], system->my_rank );
+// fprintf( stderr, "atom( %8.3f %8.3f %8.3f )\n",
+// atom->x[0], atom->x[1],
+// atom->x[2] );
c++;
}
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index aed06740..751df004 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -236,15 +236,14 @@ void Init_Taper( control_params *control )
fprintf( stderr, "Warning: non-zero value for lower Taper-radius cutoff\n" );
}
- if ( swb < 0 )
+ if ( swb < 0.0 )
{
fprintf( stderr, "Negative value for upper Taper-radius cutoff\n" );
exit( INVALID_INPUT );
}
- else if ( swb < 5 )
+ else if ( swb < 5.0 )
{
- fprintf( stderr, "Warning: low value for upper Taper-radius cutoff:%f\n",
- swb );
+ fprintf( stderr, "Warning: low value for upper Taper-radius cutoff:%f\n", swb );
}
d1 = swb - swa;
@@ -352,6 +351,7 @@ void Init_Workspace( reax_system *system, control_params *control,
workspace->b[i + system->N] = -1.0;
}
break;
+
case EEM_CM:
for ( i = 0; i < system->N; ++i )
{
@@ -364,9 +364,11 @@ void Init_Workspace( reax_system *system, control_params *control,
workspace->b_s[system->N] = control->cm_q_net;
workspace->b[system->N] = control->cm_q_net;
break;
+
case ACKS2_CM:
//TODO
break;
+
default:
fprintf( stderr, "Unknown charge method type. Terminating...\n" );
exit( INVALID_INPUT );
diff --git a/sPuReMD/src/lin_alg.c b/sPuReMD/src/lin_alg.c
index 727de37e..2ebf010d 100644
--- a/sPuReMD/src/lin_alg.c
+++ b/sPuReMD/src/lin_alg.c
@@ -225,10 +225,10 @@ void Transpose_I( sparse_matrix * const A )
* Hdia_inv: diagonal inverse preconditioner (constructed using H)
* y: current residual
* x: preconditioned residual
- * N: length of preconditioner and vectors (# rows in H)
+ * N: dimensions of preconditioner and vectors (# rows in H)
*/
static void diag_pre_app( const real * const Hdia_inv, const real * const y,
- real * const x, const int N )
+ real * const x, const int N )
{
unsigned int i;
@@ -245,13 +245,14 @@ static void diag_pre_app( const real * const Hdia_inv, const real * const y,
* LU: lower/upper triangular, stored in CSR
* y: constants in linear system (RHS)
* x: solution
+ * N: dimensions of matrix and vectors
* tri: triangularity of LU (lower/upper)
*
* Assumptions:
* LU has non-zero diagonals
* Each row of LU has at least one non-zero (i.e., no rows with all zeros) */
static void tri_solve( const sparse_matrix * const LU, const real * const y,
- real * const x, const TRIANGULARITY tri )
+ real * const x, const int N, const TRIANGULARITY tri )
{
int i, pj, j, si, ei;
real val;
@@ -260,7 +261,7 @@ static void tri_solve( const sparse_matrix * const LU, const real * const y,
{
if ( tri == LOWER )
{
- for ( i = 0; i < LU->n; ++i )
+ for ( i = 0; i < N; ++i )
{
x[i] = y[i];
si = LU->start[i];
@@ -276,7 +277,7 @@ static void tri_solve( const sparse_matrix * const LU, const real * const y,
}
else
{
- for ( i = LU->n - 1; i >= 0; --i )
+ for ( i = N - 1; i >= 0; --i )
{
x[i] = y[i];
si = LU->start[i];
@@ -299,14 +300,16 @@ static void tri_solve( const sparse_matrix * const LU, const real * const y,
* LU: lower/upper triangular, stored in CSR
* y: constants in linear system (RHS)
* x: solution
+ * N: dimensions of matrix and vectors
* tri: triangularity of LU (lower/upper)
* find_levels: perform level search if positive, otherwise reuse existing levels
*
* Assumptions:
* LU has non-zero diagonals
* Each row of LU has at least one non-zero (i.e., no rows with all zeros) */
-static void tri_solve_level_sched( const sparse_matrix * const LU, const real * const y,
- real * const x, const TRIANGULARITY tri, int find_levels )
+static void tri_solve_level_sched( const sparse_matrix * const LU,
+ const real * const y, real * const x, const int N,
+ const TRIANGULARITY tri, int find_levels )
{
int i, j, pj, local_row, local_level;
@@ -329,9 +332,9 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
if ( row_levels == NULL || level_rows == NULL || level_rows_cnt == NULL )
{
- if ( (row_levels = (unsigned int*) malloc((size_t)LU->n * sizeof(unsigned int))) == NULL
- || (level_rows = (unsigned int*) malloc((size_t)LU->n * sizeof(unsigned int))) == NULL
- || (level_rows_cnt = (unsigned int*) malloc((size_t)(LU->n + 1) * sizeof(unsigned int))) == NULL )
+ if ( (row_levels = (unsigned int*) malloc((size_t)N * sizeof(unsigned int))) == NULL
+ || (level_rows = (unsigned int*) malloc((size_t)N * sizeof(unsigned int))) == NULL
+ || (level_rows_cnt = (unsigned int*) malloc((size_t)(N + 1) * sizeof(unsigned int))) == NULL )
{
fprintf( stderr, "Not enough space for triangular solve via level scheduling. Terminating...\n" );
exit( INSUFFICIENT_MEMORY );
@@ -340,7 +343,7 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
if ( top == NULL )
{
- if ( (top = (unsigned int*) malloc((size_t)(LU->n + 1) * sizeof(unsigned int))) == NULL )
+ if ( (top = (unsigned int*) malloc((size_t)(N + 1) * sizeof(unsigned int))) == NULL )
{
fprintf( stderr, "Not enough space for triangular solve via level scheduling. Terminating...\n" );
exit( INSUFFICIENT_MEMORY );
@@ -350,14 +353,14 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
/* find levels (row dependencies in substitutions) */
if ( find_levels == TRUE )
{
- memset( row_levels, 0, LU->n * sizeof(unsigned int) );
- memset( level_rows_cnt, 0, LU->n * sizeof(unsigned int) );
- memset( top, 0, LU->n * sizeof(unsigned int) );
+ memset( row_levels, 0, N * sizeof(unsigned int) );
+ memset( level_rows_cnt, 0, N * sizeof(unsigned int) );
+ memset( top, 0, N * sizeof(unsigned int) );
levels = 1;
if ( tri == LOWER )
{
- for ( i = 0; i < LU->n; ++i )
+ for ( i = 0; i < N; ++i )
{
local_level = 1;
for ( pj = LU->start[i]; pj < LU->start[i + 1] - 1; ++pj )
@@ -372,12 +375,12 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
//#if defined(DEBUG)
fprintf(stderr, "levels(L): %d\n", levels);
- fprintf(stderr, "NNZ(L): %d\n", LU->start[LU->n]);
+ fprintf(stderr, "NNZ(L): %d\n", LU->start[N]);
//#endif
}
else
{
- for ( i = LU->n - 1; i >= 0; --i )
+ for ( i = N - 1; i >= 0; --i )
{
local_level = 1;
for ( pj = LU->start[i] + 1; pj < LU->start[i + 1]; ++pj )
@@ -392,7 +395,7 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
//#if defined(DEBUG)
fprintf(stderr, "levels(U): %d\n", levels);
- fprintf(stderr, "NNZ(U): %d\n", LU->start[LU->n]);
+ fprintf(stderr, "NNZ(U): %d\n", LU->start[N]);
//#endif
}
@@ -402,7 +405,7 @@ static void tri_solve_level_sched( const sparse_matrix * const LU, const real *
top[i] = level_rows_cnt[i];
}
- for ( i = 0; i < LU->n; ++i )
+ for ( i = 0; i < N; ++i )
{
level_rows[top[row_levels[i] - 1]] = i;
++top[row_levels[i] - 1];
@@ -1043,8 +1046,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
case ICHOLT_PC:
case ILU_PAR_PC:
case ILUT_PAR_PC:
- tri_solve( workspace->L, y, x, LOWER );
- tri_solve( workspace->U, x, x, UPPER );
+ 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" );
@@ -1061,8 +1064,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
case ICHOLT_PC:
case ILU_PAR_PC:
case ILUT_PAR_PC:
- tri_solve_level_sched( workspace->L, y, x, LOWER, fresh_pre );
- tri_solve_level_sched( workspace->U, x, x, UPPER, fresh_pre );
+ 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" );
@@ -1088,8 +1091,8 @@ static void apply_preconditioner( const static_storage * const workspace, const
#pragma omp barrier
permute_vector( y_p, workspace->H->n, FALSE, LOWER );
- tri_solve_level_sched( workspace->L, y_p, x, LOWER, fresh_pre );
- tri_solve_level_sched( workspace->U, x, x, UPPER, fresh_pre );
+ 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:
@@ -1745,8 +1748,8 @@ int PCG( static_storage *workspace, sparse_matrix *A, real *b, real tol,
//Print_Soln( workspace, x, q, b, N );
//fprintf( stderr, "res: %.15e\n", r_norm );
- tri_solve( L, workspace->r, workspace->d, LOWER );
- tri_solve( U, workspace->d, workspace->p, UPPER );
+ tri_solve( L, workspace->r, workspace->d, L->n, LOWER );
+ tri_solve( U, workspace->d, workspace->p, U->n, UPPER );
sig_new = Dot( workspace->r, workspace->p, N );
sig0 = sig_new;
@@ -1764,8 +1767,8 @@ int PCG( static_storage *workspace, sparse_matrix *A, real *b, real tol,
r_norm = Norm(workspace->r, N);
//fprintf( stderr, "res: %.15e\n", r_norm );
- tri_solve( L, workspace->r, workspace->d, LOWER );
- tri_solve( U, workspace->d, workspace->d, UPPER );
+ tri_solve( L, workspace->r, workspace->d, L->n, LOWER );
+ tri_solve( U, workspace->d, workspace->d, U->n, UPPER );
sig_old = sig_new;
sig_new = Dot( workspace->r, workspace->d, N );
beta = sig_new / sig_old;
@@ -1922,8 +1925,8 @@ real condest( const sparse_matrix * const L, const sparse_matrix * const U )
memset( e, 1., N * sizeof(real) );
- tri_solve( L, e, e, LOWER );
- tri_solve( U, e, e, UPPER );
+ tri_solve( L, e, e, L->n, LOWER );
+ tri_solve( U, e, e, U->n, UPPER );
/* compute 1-norm of vector e */
c = FABS(e[0]);
--
GitLab