#pragma once

#include <cstddef>
#include <stdexcept>
#include <vector>

namespace qc {

class Matrix {
public:
  Matrix() = default;
  Matrix(std::size_t rows, std::size_t cols, double fill = 0.0);

  std::size_t rows() const { return rows_; }
  std::size_t cols() const { return cols_; }
  double* data() { return data_.data(); }
  const double* data() const { return data_.data(); }

  double& operator()(std::size_t r, std::size_t c) {
    return data_[r * cols_ + c];
  }
  double operator()(std::size_t r, std::size_t c) const {
    return data_[r * cols_ + c];
  }

  void resize(std::size_t rows, std::size_t cols, double fill = 0.0);

  Matrix transposed() const;
  Matrix symmetrize_lower() const;

  static Matrix identity(std::size_t n);

private:
  std::size_t rows_{0};
  std::size_t cols_{0};
  std::vector<double> data_;
};

class Vector {
public:
  Vector() = default;
  explicit Vector(std::size_t n, double fill = 0.0);

  std::size_t size() const { return data_.size(); }
  double* data() { return data_.data(); }
  const double* data() const { return data_.data(); }

  double& operator[](std::size_t i) { return data_[i]; }
  double operator[](std::size_t i) const { return data_[i]; }

  void resize(std::size_t n, double fill = 0.0);

private:
  std::vector<double> data_;
};

void gemm(double alpha, const Matrix& A, const Matrix& B, double beta, Matrix& C);

void axpy(double a, const Vector& x, Vector& y);
double dot(const Vector& x, const Vector& y);
double norm2(const Vector& x);

struct CholeskyResult {
  Matrix L;
  bool ok{false};
};

CholeskyResult cholesky_lower(const Matrix& A);

void cholesky_solve_lower(const Matrix& L, Vector& x);

struct SymmetricEigenDecomp {
  Vector eigenvalues;
  Matrix eigenvectors;
};

SymmetricEigenDecomp symmetric_eigen_jacobi(Matrix A, double tol = 1e-14,
                                             int max_sweeps = 50);

void generalized_symmetric_eigen(const Matrix& F, const Matrix& S,
                                 Vector& evals, Matrix& evecs);

Matrix matrix_add(const Matrix& A, const Matrix& B);
void matrix_scale(double s, Matrix& A);
Matrix matrix_copy(const Matrix& A);

}  // namespace qc