lin_alg.c

    return i;
}


/* Dual iteration for the Pipelined Preconditioned Conjugate Gradient Method
 * for QEq (2 simaltaneous solves)
 *
 * References:
 * 1) Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm,
 *  P. Ghysels and W. Vanroose, Parallel Computing, 2014.
 * 2) Scalable Non-blocking Preconditioned Conjugate Gradient Methods,
 *  Paul R. Eller and William Gropp, SC '16 Proceedings of the International Conference
 *  for High Performance Computing, Networking, Storage and Analysis, 2016.
 *  */
int dual_PIPECG( reax_system const * const system, control_params const * const control,
        simulation_data * const data,
        storage * const workspace, sparse_matrix * const H, rvec2 * const b,
        real tol, rvec2 * const x, mpi_datatypes * const  mpi_data )
{
    int i, j;
    rvec2 alpha, beta, delta, gamma_old, gamma_new, norm, b_norm;
    real t_start, t_pa, t_spmv, t_vops, t_comm, t_allreduce;
    real timings[5], redux[8];
    MPI_Request req;

    t_pa = 0.0;
    t_spmv = 0.0;
    t_vops = 0.0;
    t_comm = 0.0;
    t_allreduce = 0.0;

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            x, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    dual_Sparse_MatVec_local( H, x, workspace->u2, H->NT );
#else
    dual_Sparse_MatVec_local( H, x, workspace->u2, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->u2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
    //Vector_Sum( workspace->r , 1.0,  b, -1.0, workspace->u, system->n );
    for ( j = 0; j < system->n; ++j )
    {
        workspace->r2[j][0] = b[j][0] - workspace->u2[j][0];
        workspace->r2[j][1] = b[j][1] - workspace->u2[j][1];
    }
    t_vops += MPI_Wtime( ) - t_start;

    /* pre-conditioning */
    if ( control->cm_solver_pre_comp_type == NONE_PC )
    {
        //Vector_Copy( workspace->u, workspace->r, system->n );
        for ( j = 0; j < system->n ; ++j )
        {
            workspace->u2[j][0] = workspace->r2[j][0];
            workspace->u2[j][1] = workspace->r2[j][1];
        }
    }
    else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
    {
        t_start = MPI_Wtime( );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->u2[j][0] = workspace->r2[j][0] * workspace->Hdia_inv[j];
            workspace->u2[j][1] = workspace->r2[j][1] * workspace->Hdia_inv[j];
        }
        t_pa += MPI_Wtime( ) - t_start;
    }
    else if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->r2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );
        
        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->r2, workspace->u2, H->NT );
#else
        dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->r2, workspace->u2, system->n );
#endif
        t_pa += MPI_Wtime( ) - t_start;

        /* no comm part2 because u2 is only local portion */
    }

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            workspace->u2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    dual_Sparse_MatVec_local( H, workspace->u2, workspace->w2, H->NT );
#else
    dual_Sparse_MatVec_local( H, workspace->u2, workspace->w2, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->w2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
    //redux[0] = Dot_local( workspace->w, workspace->u, system->n );
    //redux[1] = Dot_local( workspace->r, workspace->u, system->n );
    //redux[2] = Dot_local( workspace->u, workspace->u, system->n );
    //redux[3] = Dot_local( b, b, system->n );
    for ( j = 0; j < 8; ++j )
    {
        redux[j] = 0.0;
    }
    for( j = 0; j < system->n; ++j )
    {
        redux[0] += workspace->w2[j][0] * workspace->u2[j][0];
        redux[1] += workspace->w2[j][1] * workspace->u2[j][1];

        redux[2] += workspace->r2[j][0] * workspace->u2[j][0];
        redux[3] += workspace->r2[j][1] * workspace->u2[j][1];

        redux[4] += workspace->u2[j][0] * workspace->u2[j][0];
        redux[5] += workspace->u2[j][1] * workspace->u2[j][1];

        redux[6] += b[j][0] * b[j][0];
        redux[7] += b[j][1] * b[j][1];
    }
    t_vops += MPI_Wtime( ) - t_start;

    MPI_Iallreduce( MPI_IN_PLACE, redux, 8, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

    /* pre-conditioning */
    if ( control->cm_solver_pre_comp_type == NONE_PC )
    {
        //Vector_Copy( workspace->m, workspace->w, system->n );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->m2[j][0] = workspace->w2[j][0];
            workspace->m2[j][1] = workspace->w2[j][1];
        }
    }
    else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
    {
        t_start = MPI_Wtime( );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->m2[j][0] = workspace->w2[j][0] * workspace->Hdia_inv[j];
            workspace->m2[j][1] = workspace->w2[j][1] * workspace->Hdia_inv[j];
        }
        t_pa += MPI_Wtime( ) - t_start;
    }
    else if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->w2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );
        
        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->w2, workspace->m2, H->NT );
#else
        dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->w2, workspace->m2, system->n );
#endif
        t_pa += MPI_Wtime( ) - t_start;

        /* no comm part2 because m2 is only local portion */
    }

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            workspace->m2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    dual_Sparse_MatVec_local( H, workspace->m2, workspace->n2, H->NT );
#else
    dual_Sparse_MatVec_local( H, workspace->m2, workspace->n2, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->n2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

    t_start = MPI_Wtime( );
    MPI_Wait( &req, MPI_STATUS_IGNORE );
    t_allreduce += MPI_Wtime( ) - t_start;
    delta[0] = redux[0];
    delta[1] = redux[1];
    gamma_new[0] = redux[2];
    gamma_new[1] = redux[3];
    norm[0] = SQRT( redux[4] );
    norm[1] = SQRT( redux[5] );
    b_norm[0] = SQRT( redux[6] );
    b_norm[1] = SQRT( redux[7] );

    for ( i = 0; i < control->cm_solver_max_iters; ++i )
    {
        if ( norm[0] / b_norm[0] <= tol || norm[1] / b_norm[1] <= tol )
        {
            break;
        }
        if ( i > 0 )
        {
            beta[0] = gamma_new[0] / gamma_old[0];
            beta[1] = gamma_new[1] / gamma_old[1];
            alpha[0] = gamma_new[0] / (delta[0] - beta[0] / alpha[0] * gamma_new[0]);
            alpha[1] = gamma_new[1] / (delta[1] - beta[1] / alpha[1] * gamma_new[1]);
        }
        else
        {
            beta[0] = 0.0;
            beta[1] = 0.0;
            alpha[0] = gamma_new[0] / delta[0];
            alpha[1] = gamma_new[1] / delta[1];
        }

        t_start = MPI_Wtime( );
        //Vector_Sum( workspace->z, 1.0, workspace->n, beta, workspace->z, system->n );
        //Vector_Sum( workspace->q, 1.0, workspace->m, beta, workspace->q, system->n );
        //Vector_Sum( workspace->p, 1.0, workspace->u, beta, workspace->p, system->n );
        //Vector_Sum( workspace->d, 1.0, workspace->w, beta, workspace->d, system->n );
        //Vector_Sum( x, 1.0, x, alpha, workspace->p, system->n );
        //Vector_Sum( workspace->u, 1.0, workspace->u, -alpha, workspace->q, system->n );
        //Vector_Sum( workspace->w, 1.0, workspace->w, -alpha, workspace->z, system->n );
        //Vector_Sum( workspace->r, 1.0, workspace->r, -alpha, workspace->d, system->n );
        //redux[0] = Dot_local( workspace->w, workspace->u, system->n );
        //redux[1] = Dot_local( workspace->r, workspace->u, system->n );
        //redux[2] = Dot_local( workspace->u, workspace->u, system->n );
        for ( j = 0; j < 6; ++j )
        {
            redux[j] = 0.0;
        }
        for ( j = 0; j < system->n; ++j )
        {
            workspace->z2[j][0] = workspace->n2[j][0] + beta[0] * workspace->z2[j][0];
            workspace->z2[j][1] = workspace->n2[j][1] + beta[1] * workspace->z2[j][1];

            workspace->q2[j][0] = workspace->m2[j][0] + beta[0] * workspace->q2[j][0];
            workspace->q2[j][1] = workspace->m2[j][1] + beta[1] * workspace->q2[j][1];

            workspace->p2[j][0] = workspace->u2[j][0] + beta[0] * workspace->p2[j][0];
            workspace->p2[j][1] = workspace->u2[j][1] + beta[1] * workspace->p2[j][1];

            workspace->d2[j][0] = workspace->w2[j][0] + beta[0] * workspace->d2[j][0];
            workspace->d2[j][1] = workspace->w2[j][1] + beta[1] * workspace->d2[j][1];

            x[j][0] += alpha[0] * workspace->p2[j][0];
            x[j][1] += alpha[1] * workspace->p2[j][1];

            workspace->u2[j][0] -= alpha[0] * workspace->q2[j][0];
            workspace->u2[j][1] -= alpha[1] * workspace->q2[j][1];

            workspace->w2[j][0] -= alpha[0] * workspace->z2[j][0];
            workspace->w2[j][1] -= alpha[1] * workspace->z2[j][1];

            workspace->r2[j][0] -= alpha[0] * workspace->d2[j][0];
            workspace->r2[j][1] -= alpha[1] * workspace->d2[j][1];

            redux[0] += workspace->w2[j][0] * workspace->u2[j][0];
            redux[1] += workspace->w2[j][1] * workspace->u2[j][1];
            
            redux[2] += workspace->r2[j][0] * workspace->u2[j][0];
            redux[3] += workspace->r2[j][1] * workspace->u2[j][1];
            
            redux[4] += workspace->u2[j][0] * workspace->u2[j][0];
            redux[5] += workspace->u2[j][1] * workspace->u2[j][1];

        }
        t_vops += MPI_Wtime( ) - t_start;

        MPI_Iallreduce( MPI_IN_PLACE, redux, 6, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

        /* pre-conditioning */
        if ( control->cm_solver_pre_comp_type == NONE_PC )
        {
            //Vector_Copy( workspace->m, workspace->w, system->n );
            for ( j = 0; j < system->n; ++j )
            {
                workspace->m2[j][0] = workspace->w2[j][0];
                workspace->m2[j][1] = workspace->w2[j][1];
            }
        }
        else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
        {
            t_start = MPI_Wtime( );
            for ( j = 0; j < system->n; ++j )
            {
                workspace->m2[j][0] = workspace->w2[j][0] * workspace->Hdia_inv[j];
                workspace->m2[j][1] = workspace->w2[j][1] * workspace->Hdia_inv[j];
            }
            t_pa += MPI_Wtime( ) - t_start;
        }
        else if ( control->cm_solver_pre_comp_type == SAI_PC )
        {
            t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                    workspace->w2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );
            
            t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
            dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->w2, workspace->m2, H->NT );
#else
            dual_Sparse_MatVec_local( &workspace->H_app_inv, workspace->w2, workspace->m2, system->n );
#endif
            t_pa += MPI_Wtime( ) - t_start;

            /* no comm part2 because m2 is only local portion */
        }

        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->m2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        dual_Sparse_MatVec_local( H, workspace->m2, workspace->n2, H->NT );
#else
        dual_Sparse_MatVec_local( H, workspace->m2, workspace->n2, system->N );
#endif
        t_spmv += MPI_Wtime( ) - t_start;

        t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
                H->format, workspace->n2, RVEC2_PTR_TYPE, mpi_data->mpi_rvec2 );

        gamma_old[0] = gamma_new[0];
        gamma_old[1] = gamma_new[1];

        t_start = MPI_Wtime( );
        MPI_Wait( &req, MPI_STATUS_IGNORE );
        t_allreduce += MPI_Wtime( ) - t_start;
        delta[0] = redux[0];
        delta[1] = redux[1];
        gamma_new[0] = redux[2];
        gamma_new[1] = redux[3];
        norm[0] = SQRT( redux[4] );
        norm[1] = SQRT( redux[5] );
    }

    timings[0] = t_pa;
    timings[1] = t_spmv;
    timings[2] = t_vops;
    timings[3] = t_comm;
    timings[4] = t_allreduce;

    if ( system->my_rank == MASTER_NODE )
    {
        MPI_Reduce( MPI_IN_PLACE, timings, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );

        data->timing.cm_solver_pre_app += timings[0] / control->nprocs;
        data->timing.cm_solver_spmv += timings[1] / control->nprocs;
        data->timing.cm_solver_vector_ops += timings[2] / control->nprocs;
        data->timing.cm_solver_comm += timings[3] / control->nprocs;
        data->timing.cm_solver_allreduce += timings[4] / control->nprocs;
    }
    else
    {
        MPI_Reduce( timings, NULL, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );
    }

    /* continue to solve the system that has not converged yet */
    if ( norm[0] / b_norm[0] > tol )
    {
        for ( j = 0; j < system->n; ++j )
        {
            workspace->s[j] = workspace->x[j][0];
        }

        i += PIPECG( system, control, data, workspace,
                H, workspace->b_s, tol, workspace->s, mpi_data );

        for ( j = 0; j < system->n; ++j )
        {
            workspace->x[j][0] = workspace->s[j];
        }
    }
    else if ( norm[1] / b_norm[1] > tol )
    {
        for ( j = 0; j < system->n; ++j )
        {
            workspace->t[j] = workspace->x[j][1];
        }

        i += PIPECG( system, control, data, workspace,
                H, workspace->b_t, tol, workspace->t, mpi_data );

        for ( j = 0; j < system->n; ++j )
        {
            workspace->x[j][1] = workspace->t[j];
        }
    }

    if ( i >= control->cm_solver_max_iters && system->my_rank == MASTER_NODE )
    {
        fprintf( stderr, "[WARNING] PIPECG convergence failed!\n" );
        return i;
    }

    return i;
}


/* Pipelined Preconditioned Conjugate Gradient Method
 *
 * References:
 * 1) Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm,
 *  P. Ghysels and W. Vanroose, Parallel Computing, 2014.
 * 2) Scalable Non-blocking Preconditioned Conjugate Gradient Methods,
 *  Paul R. Eller and William Gropp, SC '16 Proceedings of the International Conference
 *  for High Performance Computing, Networking, Storage and Analysis, 2016.
 *  */
int PIPECG( reax_system const * const system, control_params const * const control,
        simulation_data * const data,
        storage * const workspace, sparse_matrix * const H, real * const b,
        real tol, real * const x, mpi_datatypes * const  mpi_data )
{
    int i, j;
    real alpha, beta, delta, gamma_old, gamma_new, norm, b_norm;
    real t_start, t_pa, t_spmv, t_vops, t_comm, t_allreduce;
    real timings[5], redux[4];
    MPI_Request req;

    t_pa = 0.0;
    t_spmv = 0.0;
    t_vops = 0.0;
    t_comm = 0.0;
    t_allreduce = 0.0;

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            x, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    Sparse_MatVec_local( H, x, workspace->u, H->NT );
#else
    Sparse_MatVec_local( H, x, workspace->u, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->u, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
    Vector_Sum( workspace->r , 1.0,  b, -1.0, workspace->u, system->n );
    t_vops += MPI_Wtime( ) - t_start;

    /* pre-conditioning */
    if ( control->cm_solver_pre_comp_type == NONE_PC )
    {
        Vector_Copy( workspace->u, workspace->r, system->n );
    }
    else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
    {
        t_start = MPI_Wtime( );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->u[j] = workspace->r[j] * workspace->Hdia_inv[j];
        }
        t_pa += MPI_Wtime( ) - t_start;
    }
    else if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->r, REAL_PTR_TYPE, MPI_DOUBLE );
        
        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->r, workspace->u, H->NT );
#else
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->r, workspace->u, system->n );
#endif
        t_pa += MPI_Wtime( ) - t_start;

        /* no comm part2 because u is only local portion */
    }

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            workspace->u, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    Sparse_MatVec_local( H, workspace->u, workspace->w, H->NT );
#else
    Sparse_MatVec_local( H, workspace->u, workspace->w, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->w, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
    redux[0] = Dot_local( workspace->w, workspace->u, system->n );
    redux[1] = Dot_local( workspace->r, workspace->u, system->n );
    redux[2] = Dot_local( workspace->u, workspace->u, system->n );
    redux[3] = Dot_local( b, b, system->n );
    t_vops += MPI_Wtime( ) - t_start;

    MPI_Iallreduce( MPI_IN_PLACE, redux, 4, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

    /* pre-conditioning */
    if ( control->cm_solver_pre_comp_type == NONE_PC )
    {
        Vector_Copy( workspace->m, workspace->w, system->n );
    }
    else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
    {
        t_start = MPI_Wtime( );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->m[j] = workspace->w[j] * workspace->Hdia_inv[j];
        }
        t_pa += MPI_Wtime( ) - t_start;
    }
    else if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->w, REAL_PTR_TYPE, MPI_DOUBLE );
        
        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, H->NT );
#else
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, system->n );
#endif
        t_pa += MPI_Wtime( ) - t_start;

        /* no comm part2 because m is only local portion */
    }

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            workspace->m, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    Sparse_MatVec_local( H, workspace->m, workspace->n, H->NT );
#else
    Sparse_MatVec_local( H, workspace->m, workspace->n, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->n, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
    MPI_Wait( &req, MPI_STATUS_IGNORE );
    t_allreduce += MPI_Wtime( ) - t_start;
    delta = redux[0];
    gamma_new = redux[1];
    norm = SQRT( redux[2] );
    b_norm = SQRT( redux[3] );

    for ( i = 0; i < control->cm_solver_max_iters && norm / b_norm > tol; ++i )
    {
        if ( i > 0 )
        {
            beta = gamma_new / gamma_old;
            alpha = gamma_new / (delta - beta / alpha * gamma_new);
        }
        else
        {
            beta = 0.0;
            alpha = gamma_new / delta;
        }

        t_start = MPI_Wtime( );
        Vector_Sum( workspace->z, 1.0, workspace->n, beta, workspace->z, system->n );
        Vector_Sum( workspace->q, 1.0, workspace->m, beta, workspace->q, system->n );
        Vector_Sum( workspace->p, 1.0, workspace->u, beta, workspace->p, system->n );
        Vector_Sum( workspace->d, 1.0, workspace->w, beta, workspace->d, system->n );
        Vector_Sum( x, 1.0, x, alpha, workspace->p, system->n );
        Vector_Sum( workspace->u, 1.0, workspace->u, -alpha, workspace->q, system->n );
        Vector_Sum( workspace->w, 1.0, workspace->w, -alpha, workspace->z, system->n );
        Vector_Sum( workspace->r, 1.0, workspace->r, -alpha, workspace->d, system->n );
        redux[0] = Dot_local( workspace->w, workspace->u, system->n );
        redux[1] = Dot_local( workspace->r, workspace->u, system->n );
        redux[2] = Dot_local( workspace->u, workspace->u, system->n );
        t_vops += MPI_Wtime( ) - t_start;

        MPI_Iallreduce( MPI_IN_PLACE, redux, 3, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

        /* pre-conditioning */
        if ( control->cm_solver_pre_comp_type == NONE_PC )
        {
            Vector_Copy( workspace->m, workspace->w, system->n );
        }
        else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
        {
            t_start = MPI_Wtime( );
            for ( j = 0; j < system->n; ++j )
            {
                workspace->m[j] = workspace->w[j] * workspace->Hdia_inv[j];
            }
            t_pa += MPI_Wtime( ) - t_start;
        }
        else if ( control->cm_solver_pre_comp_type == SAI_PC )
        {
            t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                    workspace->w, REAL_PTR_TYPE, MPI_DOUBLE );
            
            t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
            Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, H->NT );
#else
            Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, system->n );
#endif
            t_pa += MPI_Wtime( ) - t_start;

            /* no comm part2 because m is only local portion */
        }

        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->m, REAL_PTR_TYPE, MPI_DOUBLE );

        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        Sparse_MatVec_local( H, workspace->m, workspace->n, H->NT );
#else
        Sparse_MatVec_local( H, workspace->m, workspace->n, system->N );
#endif
        t_spmv += MPI_Wtime( ) - t_start;

        t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
                H->format, workspace->n, REAL_PTR_TYPE, MPI_DOUBLE );

        gamma_old = gamma_new;

        t_start = MPI_Wtime( );
        MPI_Wait( &req, MPI_STATUS_IGNORE );
        t_allreduce += MPI_Wtime( ) - t_start;
        delta = redux[0];
        gamma_new = redux[1];
        norm = SQRT( redux[2] );
    }

    timings[0] = t_pa;
    timings[1] = t_spmv;
    timings[2] = t_vops;
    timings[3] = t_comm;
    timings[4] = t_allreduce;

    if ( system->my_rank == MASTER_NODE )
    {
        MPI_Reduce( MPI_IN_PLACE, timings, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );

        data->timing.cm_solver_pre_app += timings[0] / control->nprocs;
        data->timing.cm_solver_spmv += timings[1] / control->nprocs;
        data->timing.cm_solver_vector_ops += timings[2] / control->nprocs;
        data->timing.cm_solver_comm += timings[3] / control->nprocs;
        data->timing.cm_solver_allreduce += timings[4] / control->nprocs;
    }
    else
    {
        MPI_Reduce( timings, NULL, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );
    }

    if ( i >= control->cm_solver_max_iters && system->my_rank == MASTER_NODE )
    {
        fprintf( stderr, "[WARNING] PIPECG convergence failed!\n" );
        return i;
    }

    return i;
}


/* Pipelined Preconditioned Conjugate Residual Method
 *
 * References:
 * 1) Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm,
 *  P. Ghysels and W. Vanroose, Parallel Computing, 2014.
 *  */
int PIPECR( reax_system const * const system, control_params const * const control,
        simulation_data * const data,
        storage * const workspace, sparse_matrix * const H, real * const b,
        real tol, real * const x, mpi_datatypes * const  mpi_data )
{
    int i, j;
    real alpha, beta, delta, gamma_old, gamma_new, norm, b_norm;
    real t_start, t_pa, t_spmv, t_vops, t_comm, t_allreduce;
    real timings[5], redux[3];
    MPI_Request req;

    t_pa = 0.0;
    t_spmv = 0.0;
    t_vops = 0.0;
    t_comm = 0.0;
    t_allreduce = 0.0;

    t_start = MPI_Wtime( );
    redux[0] = Dot_local( b, b, system->n );
    t_vops += MPI_Wtime( ) - t_start;

    MPI_Iallreduce( MPI_IN_PLACE, redux, 1, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            x, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    Sparse_MatVec_local( H, x, workspace->u, H->NT );
#else
    Sparse_MatVec_local( H, x, workspace->u, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->u, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
    Vector_Sum( workspace->r , 1.0,  b, -1.0, workspace->u, system->n );
    t_vops += MPI_Wtime( ) - t_start;

    /* pre-conditioning */
    if ( control->cm_solver_pre_comp_type == NONE_PC )
    {
        Vector_Copy( workspace->u, workspace->r, system->n );
    }
    else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
    {
        t_start = MPI_Wtime( );
        for ( j = 0; j < system->n; ++j )
        {
            workspace->u[j] = workspace->r[j] * workspace->Hdia_inv[j];
        }
        t_pa += MPI_Wtime( ) - t_start;
    }
    else if ( control->cm_solver_pre_comp_type == SAI_PC )
    {
        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->r, REAL_PTR_TYPE, MPI_DOUBLE );
        
        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->r, workspace->u, H->NT );
#else
        Sparse_MatVec_local( &workspace->H_app_inv, workspace->r, workspace->u, system->n );
#endif
        t_pa += MPI_Wtime( ) - t_start;

        /* no comm part2 because u is only local portion */
    }

    t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
            workspace->u, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
    Sparse_MatVec_local( H, workspace->u, workspace->w, H->NT );
#else
    Sparse_MatVec_local( H, workspace->u, workspace->w, system->N );
#endif
    t_spmv += MPI_Wtime( ) - t_start;

    t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
            H->format, workspace->w, REAL_PTR_TYPE, MPI_DOUBLE );

    t_start = MPI_Wtime( );
    MPI_Wait( &req, MPI_STATUS_IGNORE );
    t_allreduce += MPI_Wtime( ) - t_start;
    b_norm = SQRT( redux[0] );

    t_start =  MPI_Wtime( );
    norm = Parallel_Norm( workspace->u, system->n, mpi_data->world );
    t_allreduce += MPI_Wtime( ) - t_start;

    for ( i = 0; i < control->cm_solver_max_iters && norm / b_norm > tol; ++i )
    {
        /* pre-conditioning */
        if ( control->cm_solver_pre_comp_type == NONE_PC )
        {
            Vector_Copy( workspace->m, workspace->w, system->n );
        }
        else if ( control->cm_solver_pre_comp_type == JACOBI_PC )
        {
            t_start = MPI_Wtime( );
            for ( j = 0; j < system->n; ++j )
            {
                workspace->m[j] = workspace->w[j] * workspace->Hdia_inv[j];
            }
            t_pa += MPI_Wtime( ) - t_start;
        }
        else if ( control->cm_solver_pre_comp_type == SAI_PC )
        {
            t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                    workspace->w, REAL_PTR_TYPE, MPI_DOUBLE );
            
            t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
            Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, H->NT );
#else
            Sparse_MatVec_local( &workspace->H_app_inv, workspace->w, workspace->m, system->n );
#endif
            t_pa += MPI_Wtime( ) - t_start;

            /* no comm part2 because m is only local portion */
        }

        t_start = MPI_Wtime( );
        redux[0] = Dot_local( workspace->w, workspace->u, system->n );
        redux[1] = Dot_local( workspace->m, workspace->w, system->n );
        redux[2] = Dot_local( workspace->u, workspace->u, system->n );
        t_vops += MPI_Wtime( ) - t_start;

        MPI_Iallreduce( MPI_IN_PLACE, redux, 3, MPI_DOUBLE, MPI_SUM, mpi_data->world, &req );

        t_comm += Sparse_MatVec_Comm_Part1( system, control, mpi_data,
                workspace->m, REAL_PTR_TYPE, MPI_DOUBLE );

        t_start = MPI_Wtime( );
#if defined(NEUTRAL_TERRITORY)
        Sparse_MatVec_local( H, workspace->m, workspace->n, H->NT );
#else
        Sparse_MatVec_local( H, workspace->m, workspace->n, system->N );
#endif
        t_spmv += MPI_Wtime( ) - t_start;

        t_comm += Sparse_MatVec_Comm_Part2( system, control, mpi_data,
                H->format, workspace->n, REAL_PTR_TYPE, MPI_DOUBLE );

        t_start = MPI_Wtime( );
        MPI_Wait( &req, MPI_STATUS_IGNORE );
        t_allreduce += MPI_Wtime( ) - t_start;
        gamma_new = redux[0];
        delta = redux[1];
        norm = SQRT( redux[2] );

        if ( i > 0 )
        {
            beta = gamma_new / gamma_old;
            alpha = gamma_new / (delta - beta / alpha * gamma_new);
        }
        else
        {
            beta = 0.0;
            alpha = gamma_new / delta;
        }

        t_start = MPI_Wtime( );
        Vector_Sum( workspace->z, 1.0, workspace->n, beta, workspace->z, system->n );
        Vector_Sum( workspace->q, 1.0, workspace->m, beta, workspace->q, system->n );
        Vector_Sum( workspace->p, 1.0, workspace->u, beta, workspace->p, system->n );
        Vector_Sum( workspace->d, 1.0, workspace->w, beta, workspace->d, system->n );
        Vector_Sum( x, 1.0, x, alpha, workspace->p, system->n );
        Vector_Sum( workspace->u, 1.0, workspace->u, -alpha, workspace->q, system->n );
        Vector_Sum( workspace->w, 1.0, workspace->w, -alpha, workspace->z, system->n );
        Vector_Sum( workspace->r, 1.0, workspace->r, -alpha, workspace->d, system->n );
        t_vops += MPI_Wtime( ) - t_start;

        gamma_old = gamma_new;
    }

    timings[0] = t_pa;
    timings[1] = t_spmv;
    timings[2] = t_vops;
    timings[3] = t_comm;
    timings[4] = t_allreduce;

    if ( system->my_rank == MASTER_NODE )
    {
        MPI_Reduce( MPI_IN_PLACE, timings, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );

        data->timing.cm_solver_pre_app += timings[0] / control->nprocs;
        data->timing.cm_solver_spmv += timings[1] / control->nprocs;
        data->timing.cm_solver_vector_ops += timings[2] / control->nprocs;
        data->timing.cm_solver_comm += timings[3] / control->nprocs;
        data->timing.cm_solver_allreduce += timings[4] / control->nprocs;
    }
    else
    {
        MPI_Reduce( timings, NULL, 5, MPI_DOUBLE, MPI_SUM, MASTER_NODE, mpi_data->world );
    }

    if ( i >= control->cm_solver_max_iters && system->my_rank == MASTER_NODE )
    {
        fprintf( stderr, "[WARNING] PIPECR convergence failed!\n" );
        return i;
    }

    return i;
}