#[cfg(test)]
mod cg {

    use rstest::*;
    use mtxlib::{*, st_ops::*};

    use krysys::cg::*;

    struct SymM<'a>(SymmetricMatrixCompact<'a, f64>);
    struct GenM<'a>(GenericMatrix<'a, f64>);

    macro_rules! impl_linear_map {
        ($t:ty) => {
            impl<H> LinearMap<f64, H> for $t {
                #[inline(always)]
                fn apply(&self, result: &mut [f64], input: &[f64], _helper: &mut H) {
                    result.reset();
                    let Self(core) = self;
                    core.mv_mul(result, input);
                }
            }
        };
    }

    impl_linear_map!(SymM<'_>);
    impl_linear_map!(GenM<'_>);

    const TOLERANCE: f64 = 1e-12;
    const REDUCED_TOL: f64 = 1e-8;

    #[rstest]
    #[case(0.0, -3.0)]
    #[case(-3.0, -3.0)]
    #[case(-3.0, 5.0)]
    fn identity(#[case] x_vals: f64, #[case] b_vals: f64) {
        const N: usize = 10;

        let rows: Vec<usize> = (0..N).collect();
        let cols: Vec<usize> = (0..N).collect();
        let vals = [1.0; N];

        let b = [b_vals; N];

        let (mut x, mut error) = ([x_vals; N], [0.0; N]);

        let mut plugin = DoNothing;

        let matrix = GenM(GenericMatrix::new(&rows, &cols, &vals));
        let opts = Settings { abs_tolerance: TOLERANCE, max_iterations: N, ..Default::default() };

        let _report =
            matrix.cg(&mut x, &b, &opts, &mut plugin, Buckets::new(&mut [0.0; N], &mut [0.0; N], &mut [0.0; N]));

        // Compute true residual
        matrix.apply(&mut error, &x, &mut plugin);
        let true_residual = error.scale_add(-1.0, 1.0, &b).norm_2();

        // Compute solution error
        let solution_error = error.linear_comb(1.0, &x, -1.0, &b).norm_2();

        assert!(true_residual <= TOLERANCE);
        assert!(solution_error <= REDUCED_TOL);
    }

    #[rstest]
    fn test_1by1() {
        let maxiters: usize = 100;
        let mut plugin = DoNothing;
        const SIZE: usize = 1;

        let rows: [usize; 1] = [0];
        let cols: [usize; 1] = [0];
        let vals: [f64; 1] = [2.0; 1];

        let matrix = GenM(GenericMatrix::new(&rows, &cols, &vals));
        let opts = Settings { abs_tolerance: TOLERANCE, max_iterations: maxiters, ..Default::default() };

        let sol: [f64; SIZE] = [1.0];

        let mut b: [f64; SIZE] = [0.0; SIZE];
        matrix.apply(&mut b, &sol, &mut plugin);

        let (mut x, mut error) = (vec![0.0; SIZE], vec![0.0; SIZE]);

        let _report = matrix.cg(
            &mut x,
            &b,
            &opts,
            &mut plugin,
            Buckets::new(&mut [0.0; SIZE], &mut [0.0; SIZE], &mut [0.0; SIZE]),
        );

        // Compute true residual
        matrix.apply(&mut error, &x, &mut plugin);
        let true_residual = error.scale_add(-1.0, 1.0, &b).norm_2();

        // Compute solution error
        let solution_error = error.linear_comb(1.0, &x, -1.0, &sol).norm_2();

        assert!(true_residual <= TOLERANCE);
        assert!(solution_error <= REDUCED_TOL);
    }

    #[rstest]
    fn test_2by2() {
        let maxiters: usize = 100;
        let mut plugin = DoNothing;
        const SIZE: usize = 2;

        let rows: [usize; 3] = [0, 1, 0];
        let cols: [usize; 3] = [0, 0, 1];
        let vals: [f64; 3] = [1.0, 1.0, 1.0];

        let matrix = GenM(GenericMatrix::new(&rows, &cols, &vals));
        let opts = Settings { abs_tolerance: TOLERANCE, max_iterations: maxiters, ..Default::default() };

        let sol: [f64; SIZE] = [1.0, 1.0];

        let mut rhs: [f64; SIZE] = [0.0; SIZE];
        matrix.apply(&mut rhs, &sol, &mut plugin);

        let (mut x, mut error) = ([0.0; SIZE], [0.0; SIZE]);
        let (mut a, mut b, mut c) = (x, x, x);
        let buckets = Buckets::new(&mut a, &mut b, &mut c);
        let _report = matrix.cg(&mut x, &rhs, &opts, &mut plugin, buckets);

        // Compute true residual
        matrix.apply(&mut error, &x, &mut plugin);
        let true_residual = error.scale_add(-1.0, 1.0, &rhs).norm_2();

        // Compute solution error
        let solution_error = error.linear_comb(1.0, &x, -1.0, &sol).norm_2();

        assert!(true_residual <= TOLERANCE);
        assert!(solution_error <= REDUCED_TOL);
    }

    #[rstest]
    fn test_4by4() {
        const SIZE: usize = 4;
        let maxiters: usize = 2 * SIZE;
        let mut plugin = DoNothing;

        let rows: [usize; 6] = [0, 2, 1, 3, 0, 1];
        let cols: [usize; 6] = [0, 0, 1, 1, 2, 3];
        let vals: [f64; 6] = [3.0, 1.0, 1.0, 1.0, 1.0, 1.0];

        let matrix = GenM(GenericMatrix::new(&rows, &cols, &vals));
        let opts = Settings { abs_tolerance: TOLERANCE, max_iterations: maxiters, ..Default::default() };

        let sol: [f64; SIZE] = [1.0; SIZE];

        let mut rhs: [f64; SIZE] = [0.0; SIZE];
        matrix.apply(&mut rhs, &sol, &mut plugin);

        let (mut x, mut error) = ([0.0; SIZE], [0.0; SIZE]);

        let (mut a, mut b, mut c) = (x, x, x);
        let buckets = Buckets::new(&mut a, &mut b, &mut c);
        let _report = matrix.cg(&mut x, &rhs, &opts, &mut plugin, buckets);

        // Compute true residual
        matrix.apply(&mut error, &x, &mut plugin);
        let true_residual = error.scale_add(-1.0, 1.0, &rhs).norm_2();

        // Compute solution error
        let solution_error = error.linear_comb(1.0, &x, -1.0, &sol).norm_2();

        assert!(true_residual <= TOLERANCE);
        assert!(solution_error <= REDUCED_TOL);
    }

    #[test]
    fn singular_matrix_2by2() {
        let rows = [0, 0, 1, 1];
        let cols = [0, 1, 0, 1];
        let vals = [4.0, -2.0, -2.0, 1.0 + 1e-7];

        let matrix = GenM(GenericMatrix::new(&rows, &cols, &vals));

        const SIZE: usize = 2;

        let rhs = [-2.0, 1.0];
        let mut x = [0.0; SIZE];

        let (mut a, mut b, mut c) = (x, x, x);
        let mut plugin = DoNothing;

        let opts = Settings { abs_tolerance: TOLERANCE, max_iterations: 100, ..Default::default() };

        let buckets = Buckets::new(&mut a, &mut b, &mut c);
        let _report = matrix.cg(&mut x, &rhs, &opts, &mut plugin, buckets);

        let mut error = [0.0; SIZE];

        // Compute true residual
        matrix.apply(&mut error, &x, &mut plugin);
        let true_residual = error.scale_add(-1.0, 1.0, &rhs).norm_2();

        assert!((true_residual <= TOLERANCE) | (true_residual <= rhs.norm_2() * opts.rel_tolerance));
    }

    #[rstest]
    #[case("./tests/data/positive-definite.well-behaved/Trefethen_20.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/Trefethen_20b.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/nos4.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/mesh1em6.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/mesh1em1.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/mesh1e1.mtx")]
    #[case("./tests/data/positive-definite.well-behaved/bcsstm02.mtx")]
    #[case("./tests/data/full-rank.well-behaved/bfwb62.mtx")]
    #[case("./tests/data/full-rank.well-behaved/Stranke94.mtx")]
    #[case("./tests/data/full-rank.well-behaved/nemeth26.mtx")]
    #[case("./tests/data/singular.well-behaved/bcsstm01.mtx")]
    #[case("./tests/data/singular.well-behaved/GD99_b.mtx")]
    #[case("./tests/data/singular.well-behaved/Sandi_authors.mtx")]
    #[case("./tests/data/positive-definite.stiff/bcsstk01.mtx")]
    #[case("./tests/data/positive-definite.stiff/LF10.mtx")]
    #[case("./tests/data/positive-definite.stiff/LFAT5.mtx")]
    #[case("./tests/data/full-rank.stiff/spaceStation_1.mtx")]
    #[case("./tests/data/full-rank.stiff/meg4.mtx")]
    #[case("./tests/data/singular.stiff/GD97_b.mtx")]
    #[case("./tests/data/singular.stiff/lesmis.mtx")]
    fn matrix_market(#[case] file_path: &str) {
        let file_content = std::fs::read_to_string(file_path).unwrap();
        let owned_mat: SparseMatrix<f64> = file_content.parse().unwrap();

        if owned_mat.pattern.symmetry != SymmetryInfo::Symmetric {
            panic!("Provided matrix is not symmetric");
        }

        let (_, size) = owned_mat.dims();
        let matrix = SymM(SymmetricMatrixCompact::new(owned_mat.coords(), &owned_mat.vals));

        let mut plugin = DoNothing;
        let mut x = vec![0.0; size];
        let sol = {
            let mut aux = vec![1.0; size];

            matrix.apply(&mut aux, &vec![1.0; size], &mut plugin);

            let norm = aux.norm_2();
            aux.iter_mut().for_each(|elm| *elm /= norm);
            aux
        };

        let mut rhs: Vec<f64> = vec![0.0; size];
        matrix.apply(&mut rhs, &sol, &mut ());

        let opts = Settings { abs_tolerance: TOLERANCE, rel_tolerance: TOLERANCE, max_iterations: 9999 };
        let (mut a, mut b, mut c) = (x.clone(), x.clone(), x.clone());
        let buckets = Buckets::new(&mut a, &mut b, &mut c);
        let report = matrix.cg(&mut x, &rhs, &opts, &mut plugin, buckets);

        let mut error = vec![0.0; size];

        // // Compute true residual
        matrix.apply(&mut error, &x, &mut ());
        let true_residual = error.scale_add(-1.0, 1.0, &rhs).norm_2();

        // // Compute solution error
        let solution_error = error.linear_comb(1.0, &x, -1.0, &sol).norm_2();

        let scaled_reltol = rhs.norm_2() * opts.rel_tolerance;

        println!();
        println!("Problem : {file_path}");
        println!("Total Iterations: {}", report.iterations);
        println!("Size: {}", size);
        println!(
            "Problem is Stiff: {}",
            (report.iterations > 2 * size)
                | !((solution_error <= opts.abs_tolerance) | (solution_error <= scaled_reltol))
        );
        println!("Estimated Residual: {:.6e}", report.residual);
        println!("True Residual: {:.6e}", true_residual);
        println!("Solution error: {:.6e}", solution_error);
        println!();
        println!();

        assert!((true_residual <= opts.abs_tolerance) | (true_residual <= scaled_reltol));
    }
}
