ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
Matrix.hpp
1#pragma once
2#include "qip/Vector.hpp" // for std::vector overloads
3#include <array>
4#include <cassert>
5#include <complex>
6#include <gsl/gsl_blas.h>
7#include <gsl/gsl_complex.h>
8#include <gsl/gsl_complex_math.h>
9#include <gsl/gsl_eigen.h>
10#include <gsl/gsl_linalg.h>
11#include <gsl/gsl_math.h>
12#include <iostream>
13#include <type_traits>
14#include <utility>
15#include <vector>
16
17template <typename T>
18struct is_complex : std::false_type {};
19template <typename T>
20struct is_complex<std::complex<T>> : std::true_type {};
21template <typename T>
22constexpr bool is_complex_v = is_complex<T>::value;
23
24//==============================================================================
26namespace LinAlg {
27
29template <typename T>
30class View;
31
32//==============================================================================
34template <typename T = double>
35class Matrix {
36
37protected:
38 std::size_t m_rows;
39 std::size_t m_cols;
40 std::vector<T> m_data{};
41
42public:
44 Matrix() : m_rows(0), m_cols(0) {}
45
47 Matrix(std::size_t rows, std::size_t cols)
48 : m_rows(rows), m_cols(cols), m_data(rows * cols) {}
49
51 // excplicit, since don't alow flot->int converions
52 explicit Matrix(std::size_t dimension)
53 : m_rows(dimension), m_cols(dimension), m_data(dimension * dimension) {}
54
57 Matrix(std::initializer_list<std::initializer_list<T>> ll)
58 : m_rows(ll.size()), m_cols(ll.begin()->size()), m_data{} {
59 // way to avoid copy?
60 m_data.reserve(m_rows * m_cols);
61 for (auto &l : ll) {
62 m_data.insert(m_data.end(), l.begin(), l.end());
63 }
64 }
65
67 Matrix(std::size_t rows, std::size_t cols, std::initializer_list<T> l)
68 : m_rows(rows), m_cols(cols), m_data{l} {
69 assert(m_data.size() == rows * cols &&
70 "initializer_list must be rows*cols");
71 }
72
74 Matrix(std::size_t rows, std::size_t cols, std::vector<T> &&v)
75 : m_rows(rows), m_cols(cols), m_data{std::forward<std::vector<T>>(v)} {
76 assert(m_data.size() == rows * cols &&
77 "initializer_list must be rows*cols");
78 }
80 Matrix(std::size_t rows, std::size_t cols, const std::vector<T> &v)
81 : m_rows(rows), m_cols(cols), m_data{v} {
82 assert(m_data.size() == rows * cols &&
83 "initializer_list must be rows*cols");
84 }
85
86 //============================================================================
87
89 std::size_t rows() const { return m_rows; }
91 std::size_t cols() const { return m_cols; }
93 std::size_t size() const { return m_data.size(); }
94
97 T *data() { return m_data.data(); }
99 const T *data() const { return m_data.data(); }
100
101 //============================================================================
102
104 const T *operator[](std::size_t i) const { return &(m_data[i * m_cols]); }
106 T *operator[](std::size_t i) { return &(m_data[i * m_cols]); }
107
109 T &at(std::size_t row_i, std::size_t col_j) {
110 assert(row_i < m_rows && col_j < m_cols);
111 return m_data[row_i * m_cols + col_j];
112 }
114 T at(std::size_t row_i, std::size_t col_j) const {
115 assert(row_i < m_rows && col_j < m_cols);
116 return m_data[row_i * m_cols + col_j];
117 }
119 T &operator()(std::size_t i, std::size_t j) { return at(i, j); }
121 T operator()(std::size_t i, std::size_t j) const { return at(i, j); }
122
123 //============================================================================
124
126 auto begin() { return m_data.begin(); }
127 auto cbegin() const { return m_data.cbegin(); }
128 auto end() { return m_data.end(); }
129 auto cend() const { return m_data.cend(); }
130
132 // [[deprecated]]
133 const T *row(std::size_t row) const {
134 return m_data.data() + long(row * m_cols);
135 }
136
138 [[nodiscard]] View<T> row_view(std::size_t row) {
139 return View<T>(this->data(), row * m_cols, m_cols, 1ul);
140 }
142 [[nodiscard]] View<const T> row_view(std::size_t row) const {
143 return View<const T>(this->data(), row * m_cols, m_cols, 1ul);
144 }
146 [[nodiscard]] View<T> column_view(std::size_t col) {
147 return View<T>(this->data(), col, m_rows, m_rows);
148 }
150 [[nodiscard]] View<const T> column_view(std::size_t col) const {
151 return View<const T>(this->data(), col, m_rows, m_rows);
152 }
153
154 //============================================================================
159 [[nodiscard]] auto as_gsl_view();
160
162 [[nodiscard]] auto as_gsl_view() const;
163
164 //============================================================================
165 // Basic matrix operations:
166 //============================================================================
167 //============================================================================
168
171 [[nodiscard]] T determinant() const;
172
176
179 [[nodiscard]] Matrix<T> inverse() const;
181 [[nodiscard]] Matrix<T> transpose() const;
182
183 //============================================================================
184
188 Matrix<T> &zero();
189 // //! M -> M + aI, for I=identity (adds a to diag elements), in place
190 // Matrix<T> &plusIdent(T a = T(1));
191 //============================================================================
192
194 [[nodiscard]] Matrix<T> conj() const;
196 [[nodiscard]] auto real() const;
198 [[nodiscard]] auto imag() const;
200 [[nodiscard]] auto complex() const;
201
204
205 //============================================================================
208
210 [[nodiscard]] friend Matrix<T> mult_elements(Matrix<T> a,
211 const Matrix<T> &b) {
212 return a.mult_elements_by(b);
213 }
214
215 //============================================================================
216 // Operator overloads: +,-, scalar */
218 Matrix<T> &operator+=(const Matrix<T> &rhs);
219 Matrix<T> &operator-=(const Matrix<T> &rhs);
220 Matrix<T> &operator*=(const T x);
221 Matrix<T> &operator/=(const T x);
222
223 //============================================================================
224 // nb: these are defined inline here to avoid ambiguous overload?
225 [[nodiscard]] friend Matrix<T> operator+(Matrix<T> lhs,
226 const Matrix<T> &rhs) {
227 return (lhs += rhs);
228 }
229 [[nodiscard]] friend Matrix<T> operator-(Matrix<T> lhs,
230 const Matrix<T> &rhs) {
231 return (lhs -= rhs);
232 }
233 [[nodiscard]] friend Matrix<T> operator*(const T x, Matrix<T> rhs) {
234 return (rhs *= x);
235 }
236 [[nodiscard]] friend Matrix<T> operator*(Matrix<T> lhs, const T x) {
237 return (lhs *= x);
238 }
239 [[nodiscard]] friend Matrix<T> operator/(Matrix<T> lhs, const T x) {
240 return (lhs /= x);
241 }
242
243 //============================================================================
245 Matrix<T> &operator+=(T aI);
247 Matrix<T> &operator-=(T aI);
248
249 [[nodiscard]] friend Matrix<T> operator+(Matrix<T> M, T aI) {
250 return (M += aI);
251 }
252 [[nodiscard]] friend Matrix<T> operator-(Matrix<T> M, T aI) {
253 return (M -= aI);
254 }
255
256 //============================================================================
258 template <typename U>
259 friend Matrix<U> operator*(const Matrix<U> &a, const Matrix<U> &b);
260
261 //============================================================================
262 template <typename U>
263 friend std::ostream &operator<<(std::ostream &os, const Matrix<U> &a);
264};
265
266//==============================================================================
267//==============================================================================
268//==============================================================================
269
270//==============================================================================
271//==============================================================================
272//==============================================================================
273
274template <typename T>
275constexpr auto myEps();
278template <typename T>
279bool equal(const Matrix<T> &lhs, const Matrix<T> &rhs, T eps = myEps<T>());
280
281} // namespace LinAlg
282
283#include "Matrix.ipp"
Matrix class; row-major.
Definition Matrix.hpp:35
T at(std::size_t row_i, std::size_t col_j) const
As above, but const.
Definition Matrix.hpp:114
const T * row(std::size_t row) const
Returns raw c pointer to start of a row.
Definition Matrix.hpp:133
Matrix< T > & make_identity()
Constructs a diagonal unit matrix (identity), in place; only for square.
Definition Matrix.ipp:143
Matrix< T > & invert_in_place()
Inverts the matrix, in place. Uses GSL; via LU decomposition. Only works for double/complex<double>.
Definition Matrix.ipp:77
Matrix(std::size_t rows, std::size_t cols, const std::vector< T > &v)
Initialise a matrix from single initialiser list. {...}.
Definition Matrix.hpp:80
std::size_t size() const
Return rows*columns [total array size].
Definition Matrix.hpp:93
const T * data() const
As above, but const.
Definition Matrix.hpp:99
T * data()
Returns pointer to first element. Note: for std::complex<T>, this is a pointer to complex<T>,...
Definition Matrix.hpp:97
std::size_t rows() const
Return rows [major index size].
Definition Matrix.hpp:89
friend Matrix< U > operator*(const Matrix< U > &a, const Matrix< U > &b)
Matrix multiplication: C_ij = sum_k A_ik*B_kj.
Matrix< T > conj() const
Returns conjugate of matrix.
Definition Matrix.ipp:163
View< T > row_view(std::size_t row)
Returns a mutable 'View' of a row.
Definition Matrix.hpp:138
Matrix< T > & operator+=(const Matrix< T > &rhs)
Overload standard operators: do what expected.
Definition Matrix.ipp:218
auto begin()
iterators for underlying std::vector (entire data)
Definition Matrix.hpp:126
T & at(std::size_t row_i, std::size_t col_j)
() index access (with range checking). (i,j) returns ith row, jth col
Definition Matrix.hpp:109
auto complex() const
Converts a real to complex matrix (changes type; returns a complex matrix)
Definition Matrix.ipp:205
auto imag() const
Returns imag part of complex matrix (changes type; returns a real matrix)
Definition Matrix.ipp:194
auto as_gsl_view()
Returns gsl_matrix_view (or _float_view, _complex_view, _complex_float_view). Call ....
Definition Matrix.ipp:310
Matrix< T > & conj_in_place()
Conjugates matrix, in place.
Definition Matrix.ipp:174
T * operator[](std::size_t i)
As above, but const.
Definition Matrix.hpp:106
Matrix(std::size_t dimension)
Initialise a blank square matrix dimension*dimension, filled with 0.
Definition Matrix.hpp:52
Matrix()
Default initialiser.
Definition Matrix.hpp:44
Matrix(std::size_t rows, std::size_t cols, std::initializer_list< T > l)
Initialise a matrix from single initialiser list. {...}.
Definition Matrix.hpp:67
T & operator()(std::size_t i, std::size_t j)
() index access (with range checking). (i,j) returns ith row, jth col
Definition Matrix.hpp:119
Matrix< T > transpose() const
Returns transpose of matrix.
Definition Matrix.ipp:111
Matrix< T > & zero()
Sets all elements to zero, in place.
Definition Matrix.ipp:154
Matrix< T > inverse() const
Returns inverse of the matrix. Leaves original matrix intact. Uses GSL; via LU decomposition....
Definition Matrix.ipp:104
View< const T > row_view(std::size_t row) const
Returns an immutable 'View' of a row.
Definition Matrix.hpp:142
std::size_t cols() const
Return columns [minor index size].
Definition Matrix.hpp:91
View< T > column_view(std::size_t col)
Returns a mutable 'View' of a column.
Definition Matrix.hpp:146
T determinant() const
Returns the determinant. Uses GSL; via LU decomposition. Only works for double/complex<double>
Definition Matrix.ipp:49
Matrix< T > & mult_elements_by(const Matrix< T > &a)
Muplitplies all the elements by those of matrix a, in place: M_ij *= a_ij.
Definition Matrix.ipp:270
const T * operator[](std::size_t i) const
[] index access (with no range checking). [i][j] returns ith row, jth col
Definition Matrix.hpp:104
auto real() const
Returns real part of complex matrix (changes type; returns a real matrix)
Definition Matrix.ipp:183
Matrix(std::initializer_list< std::initializer_list< T > > ll)
Initialise a matrix from initialiser list. {{},{},{}}. Each row must be same length.
Definition Matrix.hpp:57
Matrix(std::size_t rows, std::size_t cols, std::vector< T > &&v)
Initialise a matrix from single initialiser list. {...}.
Definition Matrix.hpp:74
T operator()(std::size_t i, std::size_t j) const
As above, but const.
Definition Matrix.hpp:121
Matrix(std::size_t rows, std::size_t cols)
Initialise a blank matrix rows*cols, filled with 0.
Definition Matrix.hpp:47
friend Matrix< T > mult_elements(Matrix< T > a, const Matrix< T > &b)
Returns new matrix, C_ij = A_ij*B_ij.
Definition Matrix.hpp:210
View< const T > column_view(std::size_t col) const
Returns an immutable 'View' of a column.
Definition Matrix.hpp:150
Proved a "view" onto an array.
Definition Matrix.ipp:7
constexpr bool is_complex_v
User-defined type-trait: Checks whether T is a std::complex type.
Definition AdamsMoulton.hpp:109
Defines Matrix, Vector classes, and linear some algebra functions.
Definition Matrix.hpp:26
bool equal(const Matrix< T > &lhs, const Matrix< T > &rhs, T eps=myEps< T >())
Compares two matrices; returns true iff all elements compare relatively to better than eps.
Definition Matrix.ipp:381