linalg_ops.hpp

Go to the documentation of this file.

/*
 * Copyright © 2025-2026 Akreon LLC
 * SPDX-License-Identifier: Apache-2.0
 */
 
#pragma once
 
#include "operator/inner_product.hpp"
#include "operator/operator_traits.hpp"
#include "vector/vector.hpp"
 
namespace linalg {
 
template <VectorContainer Vec>
constexpr auto add(const Vec& A, const Vec& B) {
  // (dakre) below is safeguard against non-square matrix acting as vector with
  // duck type T as Vector<Vector<...>>
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  Vector<T, N> C;
  for (std::size_t i = 0; i < N; ++i) C(i) = A(i) + B(i);
  return C;
}
 
template <VectorContainer Vec>
constexpr auto operator+(const Vec& A, const Vec& B) {
  return add(A, B);
}
 
template <MatrixContainer Mat>
constexpr auto add(const Mat& A, const Mat& B) {
  using T = matrix_scalar_t<Mat>;
  constexpr auto Rows = mat_rows_v<Mat>;
  constexpr auto Cols = mat_cols_v<Mat>;
  Matrix<T, Rows, Cols> result;
  for (std::size_t row = 0; row < Rows; ++row)
    for (std::size_t col = 0; col < Cols; ++col)
      result(row, col) = A(row, col) + B(row, col);
  return result;
}
 
template <MatrixContainer Mat>
constexpr auto operator+(const Mat& A, const Mat& B) {
  return add(A, B);
}
 
template <VectorContainer Vec, Field Alpha>
constexpr auto scalar_multiplication(const Alpha& alpha, const Vec& A) {
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  Vector<T, N> C;
  for (std::size_t i = 0; i < N; ++i) C(i) = alpha * A(i);
  return C;
}
 
template <VectorContainer Vec, Field Alpha>
constexpr auto operator*(const Alpha& alpha, const Vec& A) {
  return scalar_multiplication(alpha, A);
}
 
template <MatrixContainer Mat, Field Alpha>
constexpr auto scalar_multiplication(const Alpha& alpha, const Mat& A) {
  using T = matrix_scalar_t<Mat>;
  constexpr auto Rows = mat_rows_v<Mat>;
  constexpr auto Cols = mat_cols_v<Mat>;
  Matrix<T, Rows, Cols> result;
  for (std::size_t row = 0; row < Rows; ++row)
    for (std::size_t col = 0; col < Cols; ++col)
      result(row, col) = alpha * A(row, col);
  return result;
}
 
template <MatrixContainer Mat, Field Alpha>
constexpr auto operator*(const Alpha& alpha, const Mat& A) {
  return scalar_multiplication(alpha, A);
}
 
template <VectorContainer Vec>
constexpr auto operator-(const Vec& A, const Vec& B) {
  return A + vector_scalar_t<Vec>(-1) * B;
}
 
template <MatrixContainer Mat>
constexpr auto operator-(const Mat& A, const Mat& B) {
  return A + matrix_scalar_t<Mat>(-1) * B;
}
 
template <FieldVector Vec>
constexpr auto transpose(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  Matrix<T, N, 1> x_t;
  for (std::size_t i = 0; i < N; ++i) x_t(i, 0) = x(i);
  return x_t;
}
 
template <ComplexVector Vec>
constexpr auto conjugate(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  Vector<T, N> result;
  for (std::size_t i = 0; i < N; ++i) result(i) = std::conj(x(i));
  return result;
}
 
template <FieldVector Vec>
  requires(vector_dim_v<Vec> == 3)
constexpr auto cross(const Vec& x, const Vec& y) {
  return Vec{x(1) * y(2) - x(2) * y(1), -(x(0) * y(2) - x(2) * y(0)),
             x(0) * y(1) - x(1) * y(0)};
}
 
template <FieldVector Vec>
  requires(vector_dim_v<Vec> > 0)
constexpr bool all_zero(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  int cnt = 0;
  for (std::size_t i = 0; i < N; ++i) cnt += x(i) == T(0);
  return cnt == N;
}
 
template <std::size_t idx, FieldVector Vec>
  requires(vector_dim_v<Vec> > 0 && idx < vector_dim_v<Vec>)
constexpr auto splice(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr auto N = vector_dim_v<Vec>;
  Vector<T, idx> a;
  std::size_t i = 0;
  for (; i < idx; ++i) a(i) = x(i);
  Vector<T, N - idx> b;
  for (std::size_t j = 0; i < N; ++i, ++j) b(j) = x(i);
  return std::tuple(a, b);
}
 
template <MatrixContainer MatA, MatrixContainer MatB>
  requires(MatMulCompatible<MatA, MatB>)
constexpr auto outer_product(const MatA& A, const MatB& B) {
  using T = matrix_scalar_t<MatA>;
  constexpr auto Rows = mat_rows_v<MatA>;
  constexpr auto Cols = mat_cols_v<MatB>;
  Matrix<T, Rows, Cols> result;
  for (std::size_t i = 0; i < Rows; ++i)
    for (std::size_t j = 0; j < Cols; ++j) result(i, j) = A(i, 0) * B(0, j);
  return result;
}
 
template <MatrixContainer Mat>
constexpr auto transpose(const Mat& A) {
  using T = matrix_scalar_t<Mat>;
  constexpr auto Cols = mat_cols_v<Mat>;
  constexpr auto Rows = mat_rows_v<Mat>;
  Matrix<T, Cols, Rows> result;
  for (std::size_t i = 0; i < Rows; ++i)
    for (std::size_t j = 0; j < Cols; ++j) result(j, i) = A(i, j);
  return result;
}
 
template <SquareMatrix Mat>
  requires(mat_rows_v<Mat> > 0)
constexpr auto trace(const Mat& A) {
  using T = matrix_scalar_t<Mat>;
  constexpr auto Rows = mat_rows_v<Mat>;
  T sum = 0;
  for (std::size_t i = 0; i < Rows; ++i) sum += A(i, i);
  return sum;
}
 
template <MatrixContainer Mat>
constexpr auto conjugate_transpose(const Mat& A)
  requires(Complex<matrix_scalar_t<Mat>>)
{
  constexpr auto Rows = mat_rows_v<Mat>;
  constexpr auto Cols = mat_cols_v<Mat>;
  auto A_t = transpose(A);
  for (std::size_t i = 0; i < Cols; ++i)
    for (std::size_t j = 0; j < Rows; ++j) A_t(i, j) = std::conj(A_t(i, j));
  return A_t;
}
 
template <MatrixContainer MatA, MatrixContainer MatB>
  requires(MatMulCompatible<MatA, MatB>)
constexpr auto operator*(const MatA& A, const MatB& B) {
  using T = matrix_scalar_t<MatA>;
  constexpr auto Rows = mat_rows_v<MatA>;
  constexpr auto Cols = mat_cols_v<MatB>;
  Matrix<T, Rows, Cols> result;
  constexpr auto Inner = mat_cols_v<MatA>;
  for (std::size_t i = 0; i < Rows; ++i) {
    for (std::size_t j = 0; j < Cols; ++j) {
      T value{};
      for (std::size_t k = 0; k < Inner; ++k) value += A(i, k) * B(k, j);
      result(i, j) = value;
    }
  }
  return result;
}
 
template <MatrixContainer Mat, VectorContainer Vec>
  requires(MatVecCompatible<Mat, Vec>)
constexpr auto operator*(const Mat& A, const Vec& X) {
  using T = matrix_scalar_t<Mat>;
  constexpr auto Rows = mat_rows_v<Mat>;
  Vector<T, Rows> result;
  constexpr auto Cols = mat_cols_v<Mat>;
  for (std::size_t i = 0; i < Rows; ++i) {
    T value{};
    for (std::size_t j = 0; j < Cols; ++j) value += A(i, j) * X(j);
    result(i) = value;
  }
  return result;
}
 
template <VectorContainer Vec>
  requires(vector_dim_v<Vec> > 0)
constexpr auto min(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr std::size_t N = vector_dim_v<Vec>;
  T mn = x(0);
  for (std::size_t i = 1; i < N; ++i) mn = x(i) < mn ? x(i) : mn;
  return mn;
}
 
template <VectorContainer Vec>
  requires(vector_dim_v<Vec> > 0)
constexpr auto max(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr std::size_t N = vector_dim_v<Vec>;
  T mx = x(0);
  for (std::size_t i = 1; i < N; ++i) mx = x(i) > mx ? x(i) : mx;
  return mx;
}
 
template <VectorContainer Vec>
  requires(vector_dim_v<Vec> > 0)
constexpr std::size_t arg_min(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr std::size_t N = vector_dim_v<Vec>;
  T mn = x(0);
  std::size_t idx{0};
  for (std::size_t i = 1; i < N; ++i)
    if (x(i) < mn) mn = x(i), idx = i;
  return idx;
}
 
template <VectorContainer Vec>
  requires(vector_dim_v<Vec> > 0)
constexpr std::size_t arg_max(const Vec& x) {
  using T = vector_element_t<Vec>;
  constexpr std::size_t N = vector_dim_v<Vec>;
  T mx = x(0);
  std::size_t idx{0};
  for (std::size_t i = 1; i < N; ++i)
    if (x(i) > mx) mx = x(i), idx = i;
  return idx;
}
 
} // namespace linalg