Matrix.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2008 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 "MantidKernel/DllConfig.h"
#include <cfloat>
#include <iosfwd>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
 
// Our aim with Matrix is to provide the definition for
// predefined types and avoid multiple copies of the same
// definition in multiple shared libraries.
// On MSVC template declarations cannot have
// both extern and __declspec(dllexport) but gcc
// requires the __attribute__((visibility("default")))
// attribute on the class definition.
#if defined(_MSC_VER)
#define KERNEL_MATRIX_CLASS_DLL
#else
#define KERNEL_MATRIX_CLASS_DLL MANTID_KERNEL_DLL
#endif
 
namespace Mantid {
 
namespace Kernel {
//-------------------------------------------------------------------------
// Forward declarations
//-------------------------------------------------------------------------
class V3D;
 
template <typename T> class KERNEL_MATRIX_CLASS_DLL Matrix {
public:
  using value_type = T;
 
private:
  size_t m_numRows;    
  size_t m_numColumns; 
 
  template <class U> using CMemoryArray = std::unique_ptr<U[]>;
  template <typename U> using MatrixMemoryPtrs = std::unique_ptr<U *[]>;
 
  CMemoryArray<T> m_rawDataAlloc;
  MatrixMemoryPtrs<T> m_rawData;
 
  void lubcmp(int *, int &); 
  void lubksb(int const *, double *);
  void rotate(double const, double const, int const, int const, int const, int const);
 
public:
  Matrix(const size_t nrow = 0, const size_t ncol = 0, bool const makeIdentity = false);
  Matrix(const std::vector<T> &, const std::vector<T> &);
  Matrix(const std::vector<T> &);
  Matrix(const std::vector<T> &, const size_t nrow, const size_t ncol);
 
  Matrix(const Matrix<T> &, const size_t nrow, const size_t ncol);
 
  Matrix(const Matrix<T> &);
  Matrix<T> &operator=(const Matrix<T> &);
  Matrix(Matrix<T> &&) noexcept;
  Matrix<T> &operator=(Matrix<T> &&) noexcept;
 
  const T *operator[](const size_t row) const { return m_rawData[row]; }
 
  T *operator[](const size_t row) { return m_rawData[row]; }
 
  Matrix<T> &operator+=(const Matrix<T> &);     
  Matrix<T> operator+(const Matrix<T> &) const; 
 
  Matrix<T> &operator-=(const Matrix<T> &);     
  Matrix<T> operator-(const Matrix<T> &) const; 
 
  Matrix<T> operator*(const Matrix<T> &) const;           
  std::vector<T> operator*(const std::vector<T> &) const; 
  void multiplyPoint(const std::vector<T> &in,
                     std::vector<T> &out) const; 
  V3D operator*(const V3D &) const;              
  Matrix<T> operator*(const T &) const;          
 
  Matrix<T> &operator*=(const Matrix<T> &); 
  Matrix<T> &operator*=(const T &);         
  Matrix<T> &operator/=(const T &);         
 
  bool operator<(const Matrix<T> &) const;
  bool operator>=(const Matrix<T> &) const;
  bool operator!=(const Matrix<T> &) const;
  bool operator==(const Matrix<T> &) const;
  bool equals(const Matrix<T> &A, const double Tolerance = FLT_EPSILON) const;
  T item(size_t row, size_t col) const { return m_rawData[row][col]; }
 
  void print() const;
  void write(std::ostream &, int const = 0) const;
  std::string str() const;
 
  // returns this matrix in 1D vector representation
  std::vector<T> getVector() const;
  // explicit conversion into the vector
  operator std::vector<T>() const {
    std::vector<T> tmp = this->getVector();
    return tmp;
  }
  //
  void setColumn(const size_t nCol, const std::vector<T> &newCol);
  void setRow(const size_t nRow, const std::vector<T> &newRow);
  void zeroMatrix(); 
  void identityMatrix();
  void setRandom(size_t seed = 0, double rMin = -1,
                 double rMax = 1); 
  void normVert();                 
  T Trace() const;                 
 
  std::vector<T> Diagonal() const;                                
  Matrix<T> preMultiplyByDiagonal(const std::vector<T> &) const;  
  Matrix<T> postMultiplyByDiagonal(const std::vector<T> &) const; 
 
  void setMem(const size_t, const size_t);
 
  std::pair<size_t, size_t> size() const { return std::pair<size_t, size_t>(m_numRows, m_numColumns); }
 
  size_t numRows() const { return m_numRows; }
 
  size_t numCols() const { return m_numColumns; }
 
  size_t Ssize() const { return (m_numRows > m_numColumns) ? m_numColumns : m_numRows; }
 
  void swapRows(const size_t, const size_t); 
  void swapCols(const size_t, const size_t); 
 
  T Invert(); 
  // Optimized inversion routine for tridiagonal matrices
  void invertTridiagonal();
  void averSymmetric();                            
  int Diagonalise(Matrix<T> &, Matrix<T> &) const; 
  void sortEigen(Matrix<T> &);                     
  Matrix<T> Tprime() const;                        
  Matrix<T> &Transpose();                          
 
  T factor();            
  T determinant() const; 
 
  void GaussJordan(Matrix<T> &); 
  T compSum() const;
 
  // Check if a rotation matrix
  bool isRotation() const;
  // Check if orthogonal
  bool isOrthogonal() const;
  // Transform to a rotation matrix
  std::vector<T> toRotation();
 
private:
  template <typename TYPE> friend void dumpToStream(std::ostream &, const Kernel::Matrix<TYPE> &, const char);
  template <typename TYPE> friend void fillFromStream(std::istream &, Kernel::Matrix<TYPE> &, const char);
};
 
template <>
inline void Matrix<int>::GaussJordan(Kernel::Matrix<int> &)
{
  throw std::invalid_argument("Gauss-Jordan inversion not valid for integer matrix");
}
 
// Explicit declarations required by Visual C++. Symbols provided by matching
// explicit instantiations in source file
#if defined(_MSC_VER)
extern template class Matrix<double>;
extern template class Matrix<int>;
extern template class Matrix<float>;
#endif
 
// clean up
#undef KERNEL_MATRIX_CLASS_DLL
 
//-------------------------------------------------------------------------
// Typedefs
//-------------------------------------------------------------------------
using DblMatrix = Mantid::Kernel::Matrix<double>;
using IntMatrix = Mantid::Kernel::Matrix<int>;
using FloatMatrix = Mantid::Kernel::Matrix<float>;
 
//-------------------------------------------------------------------------
// Utility methods
//-------------------------------------------------------------------------
template <typename T> std::ostream &operator<<(std::ostream &, const Kernel::Matrix<T> &);
template <typename T> void dumpToStream(std::ostream &, const Kernel::Matrix<T> &, const char);
 
template <typename T> std::istream &operator>>(std::istream &, Kernel::Matrix<T> &);
template <typename T> void fillFromStream(std::istream &, Kernel::Matrix<T> &, const char);
} // namespace Kernel
} // namespace Mantid