Newer
Older
Gigg, Martyn Anthony
committed
#ifndef MANTID_KERNEL_MATRIX_H_
#define MANTID_KERNEL_MATRIX_H_
Gigg, Martyn Anthony
committed
#include "MantidKernel/DllConfig.h"
#include <memory>
#include <vector>
Gigg, Martyn Anthony
committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
namespace Mantid {
namespace Kernel {
//-------------------------------------------------------------------------
// Forward declarations
//-------------------------------------------------------------------------
class V3D;
/** Numerical Matrix class. Holds a matrix of variable type and size.
Should work for real and complex objects. Carries out eigenvalue
and inversion if the matrix is square
Copyright © 2008-2011 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
This file is part of Mantid.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>.
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
template <typename T> class DLLExport Matrix {
public:
/// Enable users to retrieve the element type
typedef T value_type;
private:
size_t m_numRowsX; ///< Number of rows (x coordinate)
size_t m_numColumnsY; ///< Number of columns (y coordinate)
/// 1D memory allocation to be wrapped with pointers to the rows
template <class U> using CMemoryArray = std::unique_ptr<U[]>;
/// Pointers to rows in the raw data array to make it appear 2D
template <typename U> using MatrixMemoryPtrs = std::unique_ptr<U*[]>;
/// Pointer to allocated memory
CMemoryArray<T> m_rawDataAlloc;
/// Representation of 2D data
MatrixMemoryPtrs<T> m_rawData;
void lubcmp(int *, int &); ///< starts inversion process
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);
/** Constructor to take two vectors and multiply them to construct a matrix.
* (assuming that we have columns x row vector. */
Matrix(const std::vector<T> &, const std::vector<T> &);
/// Build square matrix from a linear vector. Throw if the vector.size() !=
/// m_numRowsX * m_numRowsX;
/// Build a non-square matrix from vector and dimensions
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 Array accessor. Use, e.g. Matrix[row][col]
const T *operator[](const size_t row) const { return m_rawData[row]; }
/// Array accessor. Use, e.g. Matrix[row][col]
T *operator[](const size_t row) { return m_rawData[row]; };
Matrix<T> &operator+=(const Matrix<T> &); ///< Basic addition operator
Matrix<T> operator+(const Matrix<T> &) const; ///< Basic addition operator
Matrix<T> &operator-=(const Matrix<T> &); ///< Basic subtraction operator
Matrix<T> operator-(const Matrix<T> &) const; ///< Basic subtraction operator
Matrix<T> operator*(const Matrix<T> &)const; ///< Basic matrix multiply
std::vector<T> operator*(const std::vector<T> &)const; ///< Multiply M*Vec
void multiplyPoint(const std::vector<T> &in,
V3D operator*(const V3D &)const; ///< Multiply M*Vec
Matrix<T> operator*(const T &)const; ///< Multiply by constant
Matrix<T> &operator*=(const Matrix<T> &); ///< Basic matrix multipy
Matrix<T> &operator*=(const T &); ///< Multiply by constant
Matrix<T> &operator/=(const T &); ///< Divide by constant
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(const int a, const int b) const {
return m_rawData[a][b]; ///< disallows access
}
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;
}
//
Lamar Moore
committed
void setColumn(const size_t nCol, const std::vector<T> &newCol);
void setRow(const size_t nRow, const std::vector<T> &newRow);
void zeroMatrix(); ///< Set the matrix to zero
void identityMatrix();
Lamar Moore
committed
void setRandom(size_t seed = 0, double rMin = -1,
double rMax = 1); ///< initialize random matrix;
void normVert(); ///< Vertical normalisation
T Trace() const; ///< Trace of the matrix
std::vector<T> Diagonal() const; ///< Returns a vector of the diagonal
Matrix<T>
preMultiplyByDiagonal(const std::vector<T> &) const; ///< pre-multiply D*this
Matrix<T> postMultiplyByDiagonal(
const std::vector<T> &) const; ///< post-multiply this*D
void setMem(const size_t, const size_t);
/// Access matrix sizes
std::pair<size_t, size_t> size() const {
return std::pair<size_t, size_t>(m_numRowsX, m_numColumnsY);
Gigg, Martyn Anthony
committed
}
/// Return the number of rows in the matrix
size_t numRows() const { return m_numRowsX; }
/// Return the number of columns in the matrix
size_t numCols() const { return m_numColumnsY; }
/// Return the smallest matrix size
size_t Ssize() const {
return (m_numRowsX > m_numColumnsY) ? m_numColumnsY : m_numRowsX;
}
void swapRows(const size_t, const size_t); ///< Swap rows (first V index)
void swapCols(const size_t, const size_t); ///< Swap cols (second V index)
T Invert(); ///< LU inversion routine
void averSymmetric(); ///< make Matrix symmetric
int Diagonalise(Matrix<T> &, Matrix<T> &) const; ///< (only Symmetric matrix)
void sortEigen(Matrix<T> &); ///< Sort eigenvectors
Matrix<T> Tprime() const; ///< Transpose the matrix
Matrix<T> &Transpose(); ///< Transpose the matrix
T factor(); ///< Calculate the factor
T determinant() const; ///< Calculate the determinant
void GaussJordan(Matrix<T> &); ///< Create a Gauss-Jordan Inversion
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
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);
};
//-------------------------------------------------------------------------
// Typedefs
//-------------------------------------------------------------------------
/// A matrix of doubles
typedef Mantid::Kernel::Matrix<double> DblMatrix;
/// A matrix of ints
typedef Mantid::Kernel::Matrix<int> IntMatrix;
//-------------------------------------------------------------------------
// Utility methods
//-------------------------------------------------------------------------
template <typename T>
DLLExport std::ostream &operator<<(std::ostream &, const Kernel::Matrix<T> &);
template <typename T>
DLLExport void dumpToStream(std::ostream &, const Kernel::Matrix<T> &,
const char);
template <typename T>
DLLExport std::istream &operator>>(std::istream &, Kernel::Matrix<T> &);
template <typename T>
DLLExport void fillFromStream(std::istream &, Kernel::Matrix<T> &, const char);
} // namespace Kernel
} // namespace Mantid
#endif // MANTID_KERNEL_MATRIX_H_