#include <fesa/solvers/linear/pardiso_linear_solver.hpp>

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <span>
#include <string_view>
#include <type_traits>
#include <vector>

#include <gtest/gtest.h>

namespace {

static_assert(
    !std::is_copy_constructible_v<fesa::PardisoLinearSolver>);
static_assert(
    !std::is_copy_assignable_v<fesa::PardisoLinearSolver>);

fesa::SymmetricCsr spd_matrix() {
    return {
        3,
        {0, 2, 4, 5},
        {0, 1, 1, 2, 2},
        {4.0, 1.0, 3.0, 1.0, 2.0},
    };
}

double independent_relative_residual(
    const fesa::SymmetricCsr& matrix,
    const std::span<const double> rhs,
    const std::span<const double> solution) {
    std::vector<double> residual(rhs.begin(), rhs.end());
    for (double& value : residual) {
        value = -value;
    }

    for (std::size_t row = 0; row < matrix.order; ++row) {
        const std::size_t begin =
            static_cast<std::size_t>(matrix.row_offsets[row]);
        const std::size_t end =
            static_cast<std::size_t>(matrix.row_offsets[row + 1]);
        for (std::size_t entry = begin; entry < end; ++entry) {
            const std::size_t column =
                static_cast<std::size_t>(matrix.column_indices[entry]);
            const double value = matrix.values[entry];
            residual[row] += value * solution[column];
            if (column != row) {
                residual[column] += value * solution[row];
            }
        }
    }

    double residual_squared = 0.0;
    double rhs_squared = 0.0;
    for (std::size_t index = 0; index < rhs.size(); ++index) {
        residual_squared += residual[index] * residual[index];
        rhs_squared += rhs[index] * rhs[index];
    }
    const double residual_norm = std::sqrt(residual_squared);
    const double rhs_norm = std::sqrt(rhs_squared);
    return rhs_norm == 0.0 ? residual_norm : residual_norm / rhs_norm;
}

bool has_diagnostic(
    const fesa::LinearSolveResult& result,
    const std::string_view code) {
    return std::ranges::find(
               result.diagnostics, code, &fesa::Diagnostic::code) !=
           result.diagnostics.end();
}

TEST(PardisoLinearSolver, SolvesThreeByThreeSpdSystem) {
    fesa::PardisoLinearSolver pardiso;
    fesa::LinearSolver& solver = pardiso;
    const fesa::SymmetricCsr matrix = spd_matrix();
    const std::vector<double> rhs{6.0, 10.0, 8.0};

    const fesa::LinearSolveResult result = solver.solve(matrix, rhs);

    ASSERT_TRUE(result.diagnostics.empty());
    ASSERT_EQ(result.solution.size(), 3);
    EXPECT_NEAR(result.solution[0], 1.0, 1.0e-12);
    EXPECT_NEAR(result.solution[1], 2.0, 1.0e-12);
    EXPECT_NEAR(result.solution[2], 3.0, 1.0e-12);
    const double independent =
        independent_relative_residual(matrix, rhs, result.solution);
    EXPECT_NEAR(result.relative_residual, independent, 1.0e-15);
    EXPECT_LT(result.relative_residual, 1.0e-12);
}

TEST(PardisoLinearSolver, SupportsRepeatedSolve) {
    fesa::PardisoLinearSolver solver;
    const fesa::SymmetricCsr matrix = spd_matrix();

    const fesa::LinearSolveResult first =
        solver.solve(matrix, std::vector<double>{6.0, 10.0, 8.0});
    const fesa::LinearSolveResult second =
        solver.solve(matrix, std::vector<double>{-3.5, 2.5, 4.5});

    ASSERT_TRUE(first.diagnostics.empty());
    ASSERT_TRUE(second.diagnostics.empty());
    ASSERT_EQ(second.solution.size(), 3);
    EXPECT_NEAR(second.solution[0], -1.0, 1.0e-12);
    EXPECT_NEAR(second.solution[1], 0.5, 1.0e-12);
    EXPECT_NEAR(second.solution[2], 2.0, 1.0e-12);
}

TEST(LinearSolver, RejectsInvalidUpperTriangleCsr) {
    fesa::PardisoLinearSolver solver;
    fesa::SymmetricCsr invalid = spd_matrix();
    invalid.column_indices[2] = 0;

    const fesa::LinearSolveResult result =
        solver.solve(invalid, std::vector<double>{6.0, 10.0, 8.0});

    EXPECT_TRUE(result.solution.empty());
    EXPECT_TRUE(has_diagnostic(result, "solver.invalid_csr"));
}

TEST(LinearSolver, RejectsRhsDimensionMismatch) {
    fesa::PardisoLinearSolver solver;

    const fesa::LinearSolveResult result =
        solver.solve(spd_matrix(), std::vector<double>{1.0, 2.0});

    EXPECT_TRUE(result.solution.empty());
    EXPECT_TRUE(has_diagnostic(result, "solver.dimension_mismatch"));
}

TEST(PardisoLinearSolver, ReportsSingularMatrix) {
    fesa::PardisoLinearSolver solver;
    const fesa::SymmetricCsr singular{
        2,
        {0, 2, 3},
        {0, 1, 1},
        {1.0, 1.0, 1.0},
    };

    const fesa::LinearSolveResult result =
        solver.solve(singular, std::vector<double>{2.0, 2.0});

    EXPECT_TRUE(result.solution.empty());
    EXPECT_FALSE(result.diagnostics.empty());
    EXPECT_TRUE(std::ranges::all_of(
        result.diagnostics,
        [](const fesa::Diagnostic& diagnostic) {
            return diagnostic.stage == fesa::DiagnosticStage::solver &&
                   diagnostic.severity == fesa::Severity::error;
        }));
}

}  // namespace