EigenComplexMatrix.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/EigenComplexMatrix.h"
 
namespace Mantid::CurveFitting {
 
ComplexMatrix::ComplexMatrix() {}
ComplexMatrix::ComplexMatrix(const size_t nx, const size_t ny) : m_matrix(Eigen::MatrixXcd(nx, ny)) { zero(); }
 
ComplexMatrix::ComplexMatrix(const ComplexMatrix &M) : m_matrix(M.eigen()) {}
 
ComplexMatrix::ComplexMatrix(const ComplexMatrix &M, const size_t row, const size_t col, const size_t nRows,
                             const size_t nCols) {
  if (row + nRows > M.size1() || col + nCols > M.size2()) {
    throw std::runtime_error("Submatrix exceeds matrix size.");
  }
 
  m_matrix = complex_matrix_map_type(M.eigen().data() + (col * M.size1()) + row, nRows, nCols,
                                     dynamic_stride(M.eigen().outerStride(), M.eigen().innerStride()));
}
 
ComplexMatrix::ComplexMatrix(ComplexMatrix &&m) noexcept : m_matrix(std::move(m.m_matrix)) {}
 
ComplexMatrix::ComplexMatrix(Eigen::MatrixXcd &&m) noexcept : m_matrix(std::move(m)) {}
 
ComplexMatrix &ComplexMatrix::operator=(const ComplexMatrix &M) {
  m_matrix = M.eigen();
  return *this;
}
 
ComplexMatrix &ComplexMatrix::operator=(ComplexMatrix &&M) {
  m_matrix = std::move(M.m_matrix);
  return *this;
}
 
ComplexMatrix &ComplexMatrix::operator=(const Eigen::MatrixXcd eigenMatrix) {
  m_matrix = eigenMatrix;
  return *this;
}
 
bool ComplexMatrix::isEmpty() const { return m_matrix.size() == 0; }
 
void ComplexMatrix::resize(const size_t nx, const size_t ny) {
  if (nx == size1() && ny == size2()) {
    return;
  }
  m_matrix.resize(nx, ny);
  zero();
}
 
size_t ComplexMatrix::size1() const { return m_matrix.rows(); }
 
size_t ComplexMatrix::size2() const { return m_matrix.cols(); }
 
void ComplexMatrix::set(const size_t i, const size_t j, const ComplexType value) {
  if (i < size1() && j < size2()) {
    m_matrix(i, j) = value;
  } else {
    throw std::out_of_range("ComplexMatrix indices are out of range.");
  }
}
 
ComplexType ComplexMatrix::get(const size_t i, const size_t j) const {
  if (i < size1() && j < size2()) {
    return m_matrix(i, j);
  }
  throw std::out_of_range("ComplexMatrix indices are out of range.");
}
 
ComplexType ComplexMatrix::operator()(const size_t i, const size_t j) const {
  return const_cast<ComplexMatrix &>(*this)(i, j);
}
 
ComplexType &ComplexMatrix::operator()(const size_t i, const size_t j) { return m_matrix(i, j); }
 
void ComplexMatrix::identity() { m_matrix.setIdentity(); }
 
void ComplexMatrix::zero() { m_matrix.setZero(); }
 
void ComplexMatrix::diag(const ComplexVector &d) {
  const auto n = d.size();
  resize(n, n);
  zero();
  for (size_t i = 0; i < n; ++i) {
    set(i, i, d.get(i));
  }
}
 
ComplexMatrix &ComplexMatrix::operator+=(const ComplexMatrix &M) {
  m_matrix += M.eigen();
  return *this;
}
ComplexMatrix &ComplexMatrix::operator+=(const ComplexType &d) {
  m_matrix.array() += d;
  return *this;
}
ComplexMatrix &ComplexMatrix::operator-=(const ComplexMatrix &M) {
  m_matrix -= M.eigen();
  return *this;
}
ComplexMatrix &ComplexMatrix::operator*=(const ComplexType &d) {
  m_matrix *= d;
  return *this;
}
 
ComplexMatrix ComplexMatrix::operator*(const EigenMatrix &m) const {
  if (m.size1() != size2()) {
    throw std::invalid_argument("Matrix by matrix multiplication: matricies are of incompatible sizes.");
  }
 
  ComplexMatrix res(m.size1(), size2());
  res.eigen() = eigen() * m.inspector();
  return res;
}
 
ComplexMatrix ComplexMatrix::operator*(const ComplexMatrix &m) const {
  if (m.size1() != size2()) {
    throw std::invalid_argument("Matrix by matrix multiplication: matricies are of incompatible sizes.");
  }
 
  ComplexMatrix res(m.size1(), size2());
  res.eigen() = eigen() * m.eigen();
  return res;
}
 
void ComplexMatrix::solve(const ComplexVector &rhs, ComplexVector &x) {
  if (size1() != size2()) {
    throw std::invalid_argument("System of linear equations: the matrix must be square.");
  }
  size_t n = size1();
  if (rhs.size() != n) {
    throw std::invalid_argument("System of linear equations: right-hand side vector has wrong size.");
  }
  if (det() == ComplexType(0, 0)) {
    throw std::invalid_argument("Matrix A is singular.");
  }
 
  Eigen::MatrixXcd b = rhs.eigen();
  Eigen::ColPivHouseholderQR<Eigen::MatrixXcd> dec(eigen());
  Eigen::VectorXcd res = dec.solve(b);
  x = res;
 
  if (!rhs.eigen().isApprox(eigen() * x.eigen())) {
    throw std::runtime_error("Matrix Solution Error: solution does not exist.");
  }
}
 
void ComplexMatrix::invert() {
  if (size1() != size2()) {
    throw std::runtime_error("Matrix inverse: the matrix must be square.");
  }
  m_matrix = m_matrix.inverse();
}
 
ComplexType ComplexMatrix::det() {
  if (size1() != size2()) {
    throw std::runtime_error("Matrix inverse: the matrix must be square.");
  }
 
  return m_matrix.determinant();
}
 
void ComplexMatrix::eigenSystemHermitian(EigenVector &eigenValues, ComplexMatrix &eigenVectors) {
  size_t n = size1();
  if (n != size2()) {
    throw std::runtime_error("Matrix eigenSystem: the matrix must be square.");
  }
  eigenValues.resize(n);
 
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXcd> solver;
  solver.compute(m_matrix);
 
  eigenValues = solver.eigenvalues();
  eigenVectors = solver.eigenvectors();
}
 
ComplexVector ComplexMatrix::copyRow(const size_t i) const {
  if (i >= size1()) {
    throw std::out_of_range("EigenMatrix row index is out of range.");
  }
 
  Eigen::VectorXcd row = complex_vector_map_type(eigen().data() + i, size2(), dynamic_stride(0, eigen().outerStride()));
  return ComplexVector(row);
}
 
ComplexVector ComplexMatrix::copyColumn(const size_t i) const {
  if (i >= size2()) {
    throw std::out_of_range("ComplexMatrix column index is out of range.");
  }
 
  Eigen::VectorXcd col =
      complex_vector_map_type(eigen().data() + (i * size1()), size1(), dynamic_stride(0, eigen().innerStride()));
  return ComplexVector(col);
}
 
void ComplexMatrix::sortColumns(const std::vector<size_t> &indices) {
  Eigen::MatrixXcd matrix = Eigen::MatrixXcd(size1(), size2());
  for (size_t col = 0; col < size2(); ++col) {
    auto col1 = indices[col];
    for (size_t row = 0; row < size1(); ++row) {
      matrix(row, col) = m_matrix(row, col1);
    }
  }
  m_matrix = matrix;
}
 
std::vector<double> ComplexMatrix::packToStdVector() const {
  const size_t n1 = size1();
  const size_t n2 = size2();
 
  std::vector<double> packed(2 * n1 * n2);
  for (size_t i = 0; i < n1; ++i) {
    for (size_t j = 0; j < n2; ++j) {
      auto k = 2 * (i * n2 + j);
      ComplexType value = get(i, j);
      packed[k] = value.real();
      packed[k + 1] = value.imag();
    }
  }
  return packed;
}
 
void ComplexMatrix::unpackFromStdVector(const std::vector<double> &packed) {
  const size_t n1 = size1();
  const size_t n2 = size2();
  if (2 * n1 * n2 != packed.size()) {
    throw std::runtime_error("Cannot unpack vector into ComplexMatrix: size mismatch.");
  }
  for (size_t i = 0; i < n1; ++i) {
    for (size_t j = 0; j < n2; ++j) {
      auto k = 2 * (i * n2 + j);
      ComplexType value(packed[k], packed[k + 1]);
      set(i, j, value);
    }
  }
}
 
ComplexMatrix ComplexMatrix::tr() const {
  ComplexMatrix res = *this;
  res.eigen() = m_matrix.transpose();
  return res;
}
 
ComplexMatrix ComplexMatrix::ctr() const {
  ComplexMatrix res = *this;
  res.eigen() = m_matrix.adjoint();
  return res;
}
 
} // namespace Mantid::CurveFitting