ampsci
c++ program for high-precision atomic structure calculations of single-valence systems
Loading...
Searching...
No Matches
Vector.ipp
1#pragma once
2
3namespace LinAlg {
4//==============================================================================
5
6template <typename T>
7Vector<T> Vector<T>::conj() const {
8 static_assert(is_complex_v<T>, "conj() only available for complex Vector");
9 std::vector<T> conj_data;
10 conj_data.reserve(this->size());
11 for (std::size_t i = 0; i < this->size(); ++i) {
12 conj_data.push_back(std::conj(this->data()[i]));
13 }
14 return Vector<T>{std::move(conj_data)};
15}
16
17template <typename T>
18auto Vector<T>::real() const {
19 static_assert(is_complex_v<T>, "real() only available for complex Vector");
20 std::vector<typename T::value_type> real_data;
21 real_data.reserve(this->size());
22 for (std::size_t i = 0; i < this->size(); ++i) {
23 real_data.push_back(std::real(this->data()[i]));
24 }
25 return Vector<typename T::value_type>{std::move(real_data)};
26}
27template <typename T>
28auto Vector<T>::imag() const {
29 static_assert(is_complex_v<T>, "imag() only available for complex Vector");
30 std::vector<typename T::value_type> imag_data;
31 imag_data.reserve(this->size());
32 for (std::size_t i = 0; i < this->size(); ++i) {
33 imag_data.push_back(std::imag(this->data()[i]));
34 }
35 return Vector<typename T::value_type>{std::move(imag_data)};
36}
37template <typename T>
38auto Vector<T>::complex() const {
39 static_assert(!is_complex_v<T>, "complex() only available for real Vector");
40 // use move constructor to avoid default Matrix construction:
41 std::vector<std::complex<T>> new_data;
42 new_data.reserve(this->m_data.size());
43 for (std::size_t i = 0; i < this->m_data.size(); ++i) {
44 new_data.push_back(this->m_data[i]);
45 }
46 return Vector<std::complex<T>>{std::move(new_data)};
47}
48
49//==============================================================================
50template <typename T>
51Vector<T> &Vector<T>::operator+=(const Vector<T> &rhs) {
52 assert(this->rows() == rhs.rows() && this->cols() == rhs.cols());
53 using namespace qip::overloads;
54 this->m_data += rhs.m_data;
55 return *this;
56}
57template <typename T>
58Vector<T> &Vector<T>::operator-=(const Vector<T> rhs) {
59 assert(this->rows() == rhs.rows() && this->cols() == rhs.cols());
60 using namespace qip::overloads;
61 this->m_data -= rhs.m_data;
62 return *this;
63}
64template <typename T>
65Vector<T> &Vector<T>::operator*=(const T x) {
66 using namespace qip::overloads;
67 this->m_data *= x;
68 return *this;
69}
70template <typename T>
71Vector<T> &Vector<T>::operator/=(const T x) {
72 using namespace qip::overloads;
73 this->m_data /= x;
74 return *this;
75}
76
77//==============================================================================
78template <typename T>
80 const auto size = std::max(this->rows(), this->cols());
81 if constexpr (std::is_same_v<T, double>) {
82 return gsl_vector_view_array(this->data(), size);
83 } else if constexpr (std::is_same_v<T, float>) {
84 return gsl_vector_float_view_array(this->data(), size);
85 } else if constexpr (std::is_same_v<T, std::complex<double>>) {
86 return gsl_vector_complex_view_array(
87 reinterpret_cast<double *>(this->data()), size);
88 } else if constexpr (std::is_same_v<T, std::complex<float>>) {
89 return gsl_vector_complex_float_view_array(
90 reinterpret_cast<float *>(this->data()), size);
91 } else {
92 assert(false &&
93 "as_gsl_view only for double/float (or complex double/float)");
94 }
95}
96
97template <typename T>
99 const auto size = std::max(this->rows(), this->cols());
100 if constexpr (std::is_same_v<T, double>) {
101 return gsl_vector_const_view_array(this->data(), size);
102 } else if constexpr (std::is_same_v<T, float>) {
103 return gsl_vector_float_const_view_array(this->data(), size);
104 } else if constexpr (std::is_same_v<T, std::complex<double>>) {
105 return gsl_vector_complex_const_view_array(
106 reinterpret_cast<const double *>(this->data()), size);
107 } else if constexpr (std::is_same_v<T, std::complex<float>>) {
108 return gsl_vector_complex_float_const_view_array(
109 reinterpret_cast<const float *>(this->data()), size);
110 } else {
111 assert(false &&
112 "as_gsl_view only for double/float (or complex double/float)");
113 }
114}
115
116} // namespace LinAlg
auto as_gsl_view()
Returns gsl_vector_view (or _float_view, _complex_view, _complex_float_view). Call ....
Definition Vector.ipp:79
Defines Matrix, Vector classes, and linear some algebra functions.
Definition Matrix.hpp:26
namespace qip::overloads provides operator overloads for std::vector
Definition Vector.hpp:450