ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
Vector.hpp
1#pragma once
2#include "Matrix.hpp"
3
4namespace LinAlg {
5
6//==============================================================================
8template <typename T = double>
9class Vector : public Matrix<T> {
10public:
12 Vector() : Matrix<T>() {}
13
15 Vector(std::size_t dimension) : Matrix<T>(dimension, 1) {}
16
19 Vector(std::initializer_list<T> l) : Matrix<T>(l.size(), 1, l) {}
20
22 Vector(std::vector<T> &&v)
23 : Matrix<T>(v.size(), 1, std::forward<std::vector<T>>(v)) {}
24
26 Vector(const std::vector<T> &v) : Matrix<T>(v.size(), 1, v) {}
27
29 Vector(const Matrix<T> &&m) : Matrix<T>(std::move(m)) {
30 assert(m.cols() == 1 && "Can only convert Matrix to Vector if matrix has 1 "
31 "column. Traspose first?");
32 }
33
35 Vector(const Matrix<T> &m) : Matrix<T>(m) {
36 assert(m.cols() == 1 && "Can only convert Matrix to Vector if matrix has 1 "
37 "column. Traspose first?");
38 }
39
40 //============================================================================
41
43 T operator[](std::size_t i) const { return this->data()[i]; }
45 T &operator[](std::size_t i) { return this->data()[i]; }
47 T &at(std::size_t i) {
48 assert(i < this->size());
49 return this->data()[i];
50 }
52 T at(std::size_t i) const {
53 assert(i < this->size());
54 return this->data()[i];
55 }
57 T &operator()(std::size_t i) { return at(i); }
59 T operator()(std::size_t i) const { return at(i); }
60
61 //============================================================================
62 [[nodiscard]] Vector<T> conj() const;
63 [[nodiscard]] auto real() const;
64 [[nodiscard]] auto imag() const;
65 [[nodiscard]] auto complex() const;
66
67 Vector<T> transpose() const = delete;
68
69 //============================================================================
70 Vector<T> &operator+=(const Vector<T> &rhs);
71 Vector<T> &operator-=(const Vector<T> rhs);
72 Vector<T> &operator*=(const T x);
73 Vector<T> &operator/=(const T x);
74
75 [[nodiscard]] friend Vector<T> operator+(Vector<T> lhs,
76 const Matrix<T> &rhs) {
77 return (lhs += rhs);
78 }
79 [[nodiscard]] friend Vector<T> operator-(Vector<T> lhs,
80 const Matrix<T> &rhs) {
81 return (lhs -= rhs);
82 }
83 [[nodiscard]] friend Vector<T> operator*(const T x, Vector<T> rhs) {
84 return (rhs *= x);
85 }
86 [[nodiscard]] friend Vector<T> operator*(Vector<T> lhs, const T x) {
87 return (lhs *= x);
88 }
89 [[nodiscard]] friend Vector<T> operator/(Vector<T> lhs, const T x) {
90 return (lhs /= x);
91 }
92
93 //============================================================================
95 [[nodiscard]] friend Vector<T> operator*(const Matrix<T> &a,
96 const Vector<T> &b) {
97 // https://www.gnu.org/software/gsl/doc/html/blas.html
98 assert(a.cols() == b.rows());
99 Vector<T> product(b.rows());
100 const auto a_gsl = a.as_gsl_view();
101 const auto b_gsl = b.as_gsl_view();
102 auto product_gsl = product.as_gsl_view();
103 if constexpr (std::is_same_v<T, double>) {
104 gsl_blas_dgemv(CblasNoTrans, 1.0, &a_gsl.matrix, &b_gsl.vector, 0.0,
105 &product_gsl.vector);
106 } else if constexpr (std::is_same_v<T, float>) {
107 gsl_blas_sgemv(CblasNoTrans, 1.0f, &a_gsl.matrix, &b_gsl.vector, 0.0f,
108 &product_gsl.vector);
109 } else if constexpr (std::is_same_v<T, std::complex<double>>) {
110 gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, &a_gsl.matrix,
111 &b_gsl.vector, GSL_COMPLEX_ZERO, &product_gsl.vector);
112 } else if constexpr (std::is_same_v<T, std::complex<float>>) {
113 const gsl_complex_float one{1.0f, 0.0f};
114 const gsl_complex_float zero{0.0f, 0.0f};
115 gsl_blas_cgemv(CblasNoTrans, one, &a_gsl.matrix, &b_gsl.vector, zero,
116 &product_gsl.vector);
117 }
118
119 return product;
120 }
121
123 [[nodiscard]] friend T operator*(const Vector<T> &a, const Vector<T> &b) {
124 // https://www.gnu.org/software/gsl/doc/html/blas.html
125 assert(a.rows() == b.rows());
126 const auto a_gsl = a.as_gsl_view();
127 const auto b_gsl = b.as_gsl_view();
128 T product = T(0);
129 if constexpr (std::is_same_v<T, double>) {
130 gsl_blas_ddot(&a_gsl.vector, &b_gsl.vector, &product);
131 } else if constexpr (std::is_same_v<T, float>) {
132 gsl_blas_sdot(&a_gsl.vector, &b_gsl.vector, &product);
133 } else if constexpr (std::is_same_v<T, std::complex<double>>) {
134 gsl_blas_zdotu(&a_gsl.vector, &b_gsl.vector,
135 reinterpret_cast<gsl_complex *>(&product));
136 } else if constexpr (std::is_same_v<T, std::complex<float>>) {
137 gsl_blas_cdotu(&a_gsl.vector, &b_gsl.vector,
138 reinterpret_cast<gsl_complex_float *>(&product));
139 }
140
141 return product;
142 }
143
145 [[nodiscard]] friend Matrix<T> outer_product(const Vector<T> &a,
146 const Vector<T> &b) {
147 Matrix<T> op(a.rows(), b.rows());
148 for (std::size_t i = 0; i < op.rows(); ++i) {
149 for (std::size_t j = 0; j < op.cols(); ++j) {
150 op[i][j] = a[i] * b[j];
151 }
152 }
153 return op;
154 }
155
156 //============================================================================
161 [[nodiscard]] auto as_gsl_view();
162
164 [[nodiscard]] auto as_gsl_view() const;
165};
166
167} // namespace LinAlg
168
169#include "Vector.ipp"
Matrix class; row-major.
Definition Matrix.hpp:35
std::size_t size() const
Return rows*columns [total array size].
Definition Matrix.hpp:93
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
auto as_gsl_view()
Returns gsl_matrix_view (or _float_view, _complex_view, _complex_float_view). Call ....
Definition Matrix.ipp:310
Matrix< T > & zero()
Sets all elements to zero, in place.
Definition Matrix.ipp:154
std::size_t cols() const
Return columns [minor index size].
Definition Matrix.hpp:91
Vector class (inherits from Matrix)
Definition Vector.hpp:9
Vector(std::size_t dimension)
Initialise a blank square matrix dimension*dimension, filled with 0.
Definition Vector.hpp:15
T & at(std::size_t i)
() index access (with range checking). (i,j) returns ith row, jth col
Definition Vector.hpp:47
Vector()
Default construct.
Definition Vector.hpp:12
auto as_gsl_view()
Returns gsl_vector_view (or _float_view, _complex_view, _complex_float_view). Call ....
Definition Vector.ipp:79
Vector(const Matrix< T > &&m)
Initialise from Matrix by move.
Definition Vector.hpp:29
Vector(std::vector< T > &&v)
Initialise from std::vector using move()
Definition Vector.hpp:22
friend Matrix< T > outer_product(const Vector< T > &a, const Vector< T > &b)
Outer product:M_ij = a_i*b_j.
Definition Vector.hpp:145
friend T operator*(const Vector< T > &a, const Vector< T > &b)
Inner product: = sum_i a_i*b_i.
Definition Vector.hpp:123
Vector(std::initializer_list< T > l)
Initialise a matrix from initialiser list. {{},{},{}}. Each row must be same length.
Definition Vector.hpp:19
Vector(const std::vector< T > &v)
Initialise from std::vector by copy.
Definition Vector.hpp:26
T operator()(std::size_t i) const
As above, but const.
Definition Vector.hpp:59
T at(std::size_t i) const
As above, but const.
Definition Vector.hpp:52
Vector(const Matrix< T > &m)
Initialise from Matrix by copy.
Definition Vector.hpp:35
friend Vector< T > operator*(const Matrix< T > &a, const Vector< T > &b)
Matrix*Vector multiplication: v_i = sum_j A_ij*B_j.
Definition Vector.hpp:95
T operator[](std::size_t i) const
[] index access (with no range checking). [i][j] returns ith row, jth col
Definition Vector.hpp:43
T & operator()(std::size_t i)
() index access (with range checking). (i,j) returns ith row, jth col
Definition Vector.hpp:57
T & operator[](std::size_t i)
As above, but const.
Definition Vector.hpp:45
Defines Matrix, Vector classes, and linear some algebra functions.
Definition Matrix.hpp:26