#pragma once

#include "spectra_reml/types.hpp"

#include <cstddef>
#include <string>
#include <vector>

namespace spectra::reml {

// C = A' B.  A is n x k and B is n x m; C is k x m.
[[nodiscard]] ColMajorMatrix cross_product(const ColMajorMatrix& a,
                                           const ColMajorMatrix& b);

// C = A' B, specialized for an orthogonal n x n matrix A.  This is kept as a
// distinct API because it is the dominant level-3 BLAS operation in a batch.
[[nodiscard]] ColMajorMatrix rotate_to_eigenspace(
    const ColMajorMatrix& eigenvectors, const ColMajorMatrix& variables);

// Overwrites a symmetric matrix with its eigenvectors (columns) and returns
// eigenvalues in ascending order.  Only the lower triangle is inspected.
[[nodiscard]] std::vector<double> symmetric_eigen_decomposition(
    ColMajorMatrix& symmetric_matrix);

// Cholesky helpers.  The factor is lower triangular and stored in the lower
// triangle of a full column-major matrix.
[[nodiscard]] bool cholesky_factor_in_place(std::vector<double>& matrix,
                                            std::size_t order,
                                            std::string* error = nullptr);
[[nodiscard]] bool cholesky_solve_in_place(
    const std::vector<double>& factor, std::size_t order, double* rhs,
    std::size_t rhs_columns, std::string* error = nullptr);
[[nodiscard]] bool cholesky_inverse(
    const std::vector<double>& factor, std::size_t order,
    std::vector<double>& inverse, std::string* error = nullptr);

[[nodiscard]] double dot(const std::vector<double>& a,
                         const std::vector<double>& b);
[[nodiscard]] double infinity_norm(const std::vector<double>& x);

}  // namespace spectra::reml