#![allow(non_snake_case)]

use krysys::minres::*;
use mtxlib::{*, st_ops::*};

use std::io::Write;

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

impl<H> LinearMap<f64, H> for SymM<'_> {
    #[inline(always)]
    fn apply(&self, result: &mut [f64], input: &[f64], _helper: &mut H) {
        result.reset();
        let Self(core) = self;
        core.mv_mul(result, input);
    }
}

fn main() {
    let file_path: &str = "./tests/data/positive-definite.well-behaved/nos4.mtx";

    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 (_, sz) = owned_mat.dims();
    let matrix = SymM(SymmetricMatrixCompact::new(owned_mat.coords(), &owned_mat.vals));

    // Initialize dense solution matrix
    let mut inverse_matrix: Vec<Vec<f64>> = vec![vec![0.0; sz]; sz];

    // Initialize right hand side
    let mut b: Vec<f64> = vec![0.0; sz];

    // Specify tolerance
    let tol: f64 = 1e-8;

    // Specify max iterations allowed
    let max_iters = sz.pow(2);

    let options = Settings { abs_tolerance: tol, rel_tolerance: tol, max_iterations: max_iters, ..Default::default() };

    // Select an appropriate plugin.
    let mut do_nothing = DoNothing;
    let mut plugin = StopWatch::extend(&mut do_nothing);
    // let mut console_logger = ConsoleLogger::extend(&mut do_nothing);
    // let mut plugin = StopWatch::extend(&mut console_logger);

    let (mut aux1, mut aux2, mut aux3, mut aux4, mut aux5) = (b.clone(), b.clone(), b.clone(), b.clone(), b.clone());

    let mut average_solver_time: f64 = 0.0;

    for kk in 0..sz {
        b.reset();
        b[kk] = 1.0;

        let report = matrix.minres(
            &mut inverse_matrix[kk],
            &b,
            &options,
            &mut plugin,
            Buckets::new(&mut aux1, &mut aux2, &mut aux3, &mut aux4, &mut aux5),
        );

        if report.status == ExitStatus::MaxIterations {
            println!("Failed to obtain solution at iteration {kk}");
        }

        average_solver_time =
            (plugin.elapsed_time_in_seconds() + average_solver_time * (kk as f64 + 1.0)) / (kk as f64 + 2.0);
    }

    // ----
    // Write the first (top left) 10x10 block of the matrix inverse to a csv file.

    let f = std::fs::File::options()
        .read(true)
        .write(true)
        .create(true)
        .truncate(true)
        .open("matrix_inverse.csv")
        .ok()
        .unwrap();
    let mut f = std::io::BufWriter::new(f);

    println!(" ");
    println!("Top left 9x9 block of the matrix inverse:");
    println!(" ");
    for row in inverse_matrix.iter().take(9) {
        for val in row.iter().take(9) {
            write!(&mut f, "{val:<+.10} ").ok();
            // print!("{val: <+14.7e} ");
            print!("{val:<+.10} ");
        }
        writeln!(&mut f).ok();
        println!();
    }

    println!("Total solver executions: {sz}.");
    println!("Average solver execution time: {:.7e} s.", average_solver_time);
}
