EigenFortranMatrix.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2016 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 +
#pragma once
 
#include <cstddef>
#include <stdexcept>
 
namespace Mantid {
namespace CurveFitting {
 
template <class MatrixClass> class FortranMatrix : public MatrixClass {
  int m_base1;
  int m_base2;
  using ElementConstType = decltype(std::declval<const MatrixClass>().operator()(0, 0));
  using ElementRefType = decltype(std::declval<MatrixClass>().operator()(0, 0));
 
public:
  FortranMatrix();
  FortranMatrix(const int nx, const int ny);
  FortranMatrix(const int iFrom, const int iTo, const int jFrom, const int jTo);
  void allocate(const int iFrom, const int iTo, const int jFrom, const int jTo);
  void allocate(const int nx, const int ny);
  int len1() const;
  int len2() const;
  ElementConstType operator()(const int i, const int j) const;
  ElementRefType operator()(const int i, const int j);
  FortranMatrix<MatrixClass> &operator=(const MatrixClass &m);
 
  MatrixClass moveToBaseMatrix();
  FortranMatrix<MatrixClass> transpose() const;
 
private:
  static size_t makeSize(const int firstIndex, const int lastIndex);
};
 
//-----------------  Method implementations -------------------------//
 
template <class MatrixClass> size_t FortranMatrix<MatrixClass>::makeSize(const int firstIndex, const int lastIndex) {
  if (lastIndex < firstIndex) {
    throw std::invalid_argument("Matrix defined with invalid index range.");
  }
  return static_cast<size_t>(lastIndex - firstIndex + 1);
}
 
template <class MatrixClass>
FortranMatrix<MatrixClass>::FortranMatrix()
    : MatrixClass(this->makeSize(1, 1), this->makeSize(1, 1)), m_base1(1), m_base2(1) {}
 
template <class MatrixClass>
FortranMatrix<MatrixClass>::FortranMatrix(const int nx, const int ny)
    : MatrixClass(this->makeSize(1, nx), this->makeSize(1, ny)), m_base1(1), m_base2(1) {}
 
template <class MatrixClass>
FortranMatrix<MatrixClass>::FortranMatrix(const int iFirst, const int iLast, const int jFirst, const int jLast)
    : MatrixClass(this->makeSize(iFirst, iLast), this->makeSize(jFirst, jLast)), m_base1(iFirst), m_base2(jFirst) {}
 
template <class MatrixClass>
void FortranMatrix<MatrixClass>::allocate(const int iFirst, const int iLast, const int jFirst, const int jLast) {
  m_base1 = iFirst;
  m_base2 = jFirst;
  this->resize(this->makeSize(iFirst, iLast), this->makeSize(jFirst, jLast));
}
 
template <class MatrixClass> void FortranMatrix<MatrixClass>::allocate(const int nx, const int ny) {
  m_base1 = 1;
  m_base2 = 1;
  this->resize(this->makeSize(1, nx), this->makeSize(1, ny));
}
 
template <class MatrixClass>
typename FortranMatrix<MatrixClass>::ElementConstType FortranMatrix<MatrixClass>::operator()(const int i,
                                                                                             const int j) const {
  return this->MatrixClass::operator()(static_cast<size_t>(i - m_base1), static_cast<size_t>(j - m_base2));
}
 
template <class MatrixClass>
typename FortranMatrix<MatrixClass>::ElementRefType FortranMatrix<MatrixClass>::operator()(const int i, const int j) {
  return this->MatrixClass::operator()(static_cast<size_t>(i - m_base1), static_cast<size_t>(j - m_base2));
}
 
template <class MatrixClass> MatrixClass FortranMatrix<MatrixClass>::moveToBaseMatrix() {
  return MatrixClass(std::move(*this));
}
 
template <class MatrixClass> int FortranMatrix<MatrixClass>::len1() const { return static_cast<int>(this->size1()); }
 
template <class MatrixClass> int FortranMatrix<MatrixClass>::len2() const { return static_cast<int>(this->size2()); }
 
template <class MatrixClass> FortranMatrix<MatrixClass> FortranMatrix<MatrixClass>::transpose() const {
  FortranMatrix<MatrixClass> res;
  res = this->tr();
  return res;
}
 
template <class MatrixClass> FortranMatrix<MatrixClass> &FortranMatrix<MatrixClass>::operator=(const MatrixClass &m) {
  this->resize(m.size1(), m.size2());
  for (size_t i = 0; i < this->size1(); i++) {
    for (size_t j = 0; j < this->size2(); j++) {
      this->operator()((int)i + 1, (int)j + 1) = m.get(i, j);
    }
  }
  return *this;
}
 
} // namespace CurveFitting
} // namespace Mantid