math.h

Go to the documentation of this file.

#ifndef COTILA_MATRIX_MATH_H_
#define COTILA_MATRIX_MATH_H_

#include<algorithm>

#include <cotila/scalar/math.h>
#include <cotila/vector/vector.h>
#include <cotila/vector/math.h>
#include <cotila/matrix/matrix.h>
#include <cotila/matrix/utility.h>
#include <cotila/detail/assert.h>

namespace cotila {

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<T, M, N> conj(const matrix<T, M, N> &m) {
  return elementwise(cotila::conj<T>, m);
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<detail::remove_complex_t<T>, M, N>
real(const matrix<T, M, N> &m) {
  return elementwise([](auto i) { return std::real(i); }, m);
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<detail::remove_complex_t<T>, M, N>
imag(const matrix<T, M, N> &m) {
  return elementwise([](auto i) { return std::imag(i); }, m);
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<T, N, M> transpose(const matrix<T, M, N> &m) {
  return generate<N, M>([&m](auto i, auto j){return m[j][i];});
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<T, N, M> hermitian(const matrix<T, M, N> &m) {
  return transpose(conj(m));
}

template <typename T, std::size_t M, std::size_t N, std::size_t P>
constexpr matrix<T, M, P> matmul(const matrix<T, M, N> &a,
                                 const matrix<T, N, P> &b) {
  return generate<M, P>([&a, &b](auto i, auto j){return cotila::sum(a.row(i)*b.column(j));});
}

template <typename T, std::size_t M, std::size_t N,
                      std::size_t P, std::size_t Q>
constexpr matrix<T, M * P, N * Q> kron(const matrix<T, M, N> &a,
                                       const matrix<T, P, Q> &b) {
  return generate<M * P, N * Q>([&a, &b](auto i, auto j) { return a[i / P][j / Q] * b[i % P][j % Q]; });
}

template <typename T, std::size_t M, std::size_t N>
constexpr T macs(const matrix<T, M, N> &m) {
  return max(
      generate<N>([&m](std::size_t i) { return sum(abs(m.column(i))); }));
}

template <typename T, std::size_t M, std::size_t N>
constexpr T mars(const matrix<T, M, N> &m) {
  return max(generate<M>([&m](std::size_t i) { return sum(abs(m.row(i))); }));
}

template <typename T, std::size_t M, std::size_t N>
constexpr std::tuple<matrix<T, M, N>, std::size_t, T>
gauss_jordan_impl(matrix<T, M, N> m, T tolerance) {
  COTILA_DETAIL_ASSERT_FLOATING_POINT(T)
  COTILA_DETAIL_ASSERT_REAL(T)

  // Define function for determining if an element is negligible
  auto negligible = [&tolerance](const T &v) { return abs(v) < tolerance; };

  T det = 1;
  std::size_t rank = 0;
  std::size_t i = 0, j = 0;
  while (i < M && j < N) {
    // Choose largest magnitude as pivot to avoid adding different magnitudes
    for (std::size_t ip = i + 1; ip < M; ++ip) {
      if (abs(m[ip][j]) > abs(m[i][j])) {
        for (std::size_t jp = 0; jp < N; ++jp) {
          auto tmp = m[ip][jp];
          m[ip][jp] = m[i][jp];
          m[i][jp] = tmp;
        }
        det *= -1;
        break;
      }
    }

    // If m_ij is still 0, continue to the next column
    if (!negligible(m[i][j])) {
      // Scale m_ij to 1
      auto s = m[i][j];
      for (std::size_t jp = 0; jp < N; ++jp)
        m[i][jp] /= s;
      det /= s;

      // Eliminate other values in the column
      for (std::size_t ip = 0; ip < M; ++ip) {
        if (ip == i)
          continue;
        if (!negligible(m[ip][j])) {
          auto s = m[ip][j];
          [&]() { // wrap this in a lambda to get around a gcc bug
            for (std::size_t jp = 0; jp < N; ++jp)
              m[ip][jp] -= s * m[i][jp];
          }();
        }
      }

      // Increment rank
      ++rank;

      // Select next row
      ++i;
    }
    ++j;
  }
  det = (rank == M) ? det : 0;
  return {m, rank, det};
}

template <typename T, std::size_t M, std::size_t N>
constexpr std::tuple<matrix<T, M, N>, std::size_t, T>
gauss_jordan_impl(const matrix<T, M, N> &m) {
  T tol = std::max(N, M) * std::numeric_limits<T>::epsilon() * mars(m);
  return gauss_jordan_impl(m, tol);
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<T, M, N> rref(const matrix<T, M, N> &m) {
  return std::get<0>(gauss_jordan_impl(m));
}

template <typename T, std::size_t M, std::size_t N>
constexpr matrix<T, M, N> rref(const matrix<T, M, N> &m, T tolerance) {
  return std::get<0>(gauss_jordan_impl(m), tolerance);
}

template <typename T, std::size_t M, std::size_t N>
constexpr std::size_t rank(const matrix<T, M, N> &m) {
  return std::get<1>(gauss_jordan_impl(m));
}

template <typename T, std::size_t M>
constexpr T det(const matrix<T, M, M> &m) {
  return std::get<2>(gauss_jordan_impl(m));
}

template <typename T, std::size_t M>
constexpr matrix<T, M, M> inverse(const matrix<T, M, M> &m) {
  if (rank(m) < M)
    throw "matrix is not invertible";
  return submat<M, M>(rref(horzcat(m, identity<T, M>)), 0, M);
}

template <typename T, std::size_t M>
constexpr T trace(const matrix<T, M, M> &m) {
  return sum(generate<M>([&m](std::size_t i){ return m[i][i]; }));
}

} // namespace cotila

#endif // COTILA_MATRIX_MATH_H_