1 // 2 // BAGEL - Brilliantly Advanced General Electronic Structure Library 3 // Filename: vectorb.h 4 // Copyright (C) 2014 Toru Shiozaki 5 // 6 // Author: Toru Shiozaki <shiozaki@northwestern.edu> 7 // Maintainer: Shiozaki group 8 // 9 // This file is part of the BAGEL package. 10 // 11 // This program is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // This program is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU General Public License 22 // along with this program. If not, see <http://www.gnu.org/licenses/>. 23 // 24 25 26 #ifndef __SRC_MATH_BVECTOR_H 27 #define __SRC_MATH_BVECTOR_H 28 29 #include <src/util/math/algo.h> 30 #include <src/util/math/btas_interface.h> 31 #include <src/util/parallel/mpi_interface.h> 32 33 namespace bagel { 34 35 template<typename DataType> 36 class Vector_; 37 38 template<typename DataType> 39 class VecView_ : public btas::TensorView1<DataType> { 40 public: 41 using data_type = DataType; 42 using btas::TensorView1<DataType>::begin; 43 using btas::TensorView1<DataType>::cbegin; 44 using btas::TensorView1<DataType>::end; 45 using btas::TensorView1<DataType>::cend; 46 47 public: VecView_(VecView_<DataType> & o)48 VecView_(VecView_<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_(VecView_<DataType> && o)49 VecView_(VecView_<DataType>&& o) : btas::TensorView1<DataType>(std::move(o)) { } VecView_(const VecView_<DataType> & o)50 VecView_(const VecView_<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_(btas::TensorView1<DataType> & o)51 VecView_(btas::TensorView1<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_(btas::TensorView1<DataType> && o)52 VecView_(btas::TensorView1<DataType>&& o) : btas::TensorView1<DataType>(std::move(o)) { } VecView_(const btas::TensorView1<DataType> & o)53 VecView_(const btas::TensorView1<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_(Vector_<DataType> & o)54 VecView_(Vector_<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_(const Vector_<DataType> & o)55 VecView_(const Vector_<DataType>& o) : btas::TensorView1<DataType>(o) { } VecView_()56 VecView_() { } ~VecView_()57 virtual ~VecView_() { } 58 size()59 size_t size() const { return this->extent(0); } rms()60 double rms() const { return std::sqrt(detail::real(btas::dotc(*this, *this))/size()); } 61 data()62 DataType* data() { /*assert(contiguous());*/ return &*begin(); } data()63 const DataType* data() const { /*assert(contiguous());*/ return &*cbegin(); } 64 operator()65 DataType& operator()(const int i) { return *(data()+i); } operator()66 const DataType& operator()(const int i) const { return *(data()+i); } 67 contiguous()68 bool contiguous() const { return this->range().ordinal().contiguous(); } 69 70 template<typename U = DataType, class = typename std::enable_if<std::is_same<std::complex<double>,U>::value>::type> get_real_part()71 std::shared_ptr<Vector_<double>> get_real_part() const { 72 assert(contiguous()); 73 auto out = std::make_shared<Vector_<double>>(size()); 74 auto iter = out->begin(); 75 for (auto& i : *this) *iter++ = std::real(i); 76 return out; 77 } 78 79 template<typename U = DataType, class = typename std::enable_if<std::is_same<std::complex<double>,U>::value>::type> get_imag_part()80 std::shared_ptr<Vector_<double>> get_imag_part() const { 81 assert(contiguous()); 82 auto out = std::make_shared<Vector_<double>>(size()); 83 auto iter = out->begin(); 84 for (auto& i : *this) *iter++ = std::imag(i); 85 return out; 86 } 87 }; 88 89 using VecView = VecView_<double>; 90 using ZVecView = VecView_<std::complex<double>>; 91 92 template<typename DataType> 93 class Vector_ : public btas::Tensor1<DataType> { 94 public: 95 using data_type = DataType; 96 using btas::Tensor1<DataType>::data; 97 using btas::Tensor1<DataType>::begin; 98 using btas::Tensor1<DataType>::cbegin; 99 using btas::Tensor1<DataType>::end; 100 using btas::Tensor1<DataType>::cend; 101 102 private: 103 // serialization 104 friend class boost::serialization::access; 105 template<class Archive> serialize(Archive & ar,const unsigned int file_version)106 void serialize(Archive& ar, const unsigned int file_version) { 107 ar & boost::serialization::base_object<btas::Tensor1<DataType>>(*this); 108 } 109 110 public: Vector_(const size_t n)111 Vector_(const size_t n) : btas::Tensor1<DataType>(n) { this->fill(0.0); } Vector_(const Vector_<DataType> & o)112 Vector_(const Vector_<DataType>& o) : btas::Tensor1<DataType>(o) { } Vector_(const btas::TensorView1<DataType> & o)113 Vector_(const btas::TensorView1<DataType>& o) : btas::Tensor1<DataType>(o) { } Vector_(btas::TensorView1<DataType> && o)114 Vector_(btas::TensorView1<DataType>&& o) : btas::Tensor1<DataType>(std::move(o)) { } Vector_(Vector_ && o)115 Vector_(Vector_&& o) : btas::Tensor1<DataType>(std::move(o)) { } Vector_()116 Vector_() { } ~Vector_()117 virtual ~Vector_() { } 118 119 template<typename U = DataType, class = typename std::enable_if<std::is_same<std::complex<double>,U>::value>::type> Vector_(const Vector_<double> & r,const Vector_<double> & i)120 Vector_(const Vector_<double>& r, const Vector_<double>& i) : btas::Tensor1<std::complex<double>>(r.size()) { 121 assert(r.size() == i.size()); 122 auto riter = r.begin(); auto iiter = i.begin(); 123 for (auto& i : *this) 124 i = *riter++ + *iiter++ * (std::complex<double>(0.0,1.0)); 125 } 126 clone()127 std::shared_ptr<Vector_<DataType>> clone() const { return std::make_shared<Vector_<DataType>>(size()); } copy()128 std::shared_ptr<Vector_<DataType>> copy() const { return std::make_shared<Vector_<DataType>>(*this); } 129 operator()130 DataType& operator()(size_t i) { return *(data()+i); } operator()131 const DataType& operator()(size_t i) const { return *(data()+i); } 132 slice(const int mstart,const int mend)133 VecView_<DataType> slice(const int mstart, const int mend) { 134 auto low = {mstart}; 135 auto up = {mend}; 136 assert(mstart >= 0 && mend <= size()); 137 return VecView_<DataType>(btas::make_rwview(this->range().slice(low, up), this->storage())); 138 } 139 slice(const int mstart,const int mend)140 const VecView_<DataType> slice(const int mstart, const int mend) const { 141 auto low = {mstart}; 142 auto up = {mend}; 143 assert(mstart >= 0 && mend <= size()); 144 return VecView_<DataType>(btas::make_rwview(this->range().slice(low, up), this->storage())); 145 } 146 size()147 size_t size() const { return this->extent(0); } norm()148 double norm() const { return std::sqrt(detail::real(btas::dotc(*this, *this))); } rms()149 double rms() const { return std::sqrt(detail::real(btas::dotc(*this, *this))/size()); } dot_product(std::shared_ptr<const Vector_<DataType>> o)150 DataType dot_product(std::shared_ptr<const Vector_<DataType>> o) const { return blas::dot_product(data(), size(), o->data()); } dot_product(const Vector_<DataType> & o)151 DataType dot_product(const Vector_<DataType>& o) const { return blas::dot_product(data(), size(), o.data()); } zero()152 void zero() { this->fill(0.0); } scale(const DataType & a)153 void scale(const DataType& a) { blas::scale_n(a, data(), size()); } ax_plus_y(const DataType a,std::shared_ptr<const Vector_<DataType>> o)154 void ax_plus_y(const DataType a, std::shared_ptr<const Vector_<DataType>> o) { btas::axpy(a, *o, *this); } ax_plus_y(const DataType a,const Vector_<DataType> & o)155 void ax_plus_y(const DataType a, const Vector_<DataType>& o) { btas::axpy(a, o, *this); } 156 157 Vector_<DataType>& operator=(const Vector_<DataType>& o) { btas::Tensor1<DataType>::operator=(o); return *this; } 158 Vector_<DataType>& operator=(Vector_<DataType>&& o) { btas::Tensor1<DataType>::operator=(std::move(o)); return *this; } 159 160 template<typename U = DataType, class = typename std::enable_if<std::is_same<std::complex<double>,U>::value>::type> get_real_part()161 std::shared_ptr<Vector_<double>> get_real_part() const { 162 auto out = std::make_shared<Vector_<double>>(size()); 163 auto iter = out->begin(); 164 for (auto& i : *this) *iter++ = std::real(i); 165 return out; 166 } 167 168 template<typename U = DataType, class = typename std::enable_if<std::is_same<std::complex<double>,U>::value>::type> get_imag_part()169 std::shared_ptr<Vector_<double>> get_imag_part() const { 170 auto out = std::make_shared<Vector_<double>>(size()); 171 auto iter = out->begin(); 172 for (auto& i : *this) *iter++ = std::imag(i); 173 return out; 174 } 175 176 void print(const std::string tag = "") const { 177 if (tag != "") 178 std::cout << std::endl << " ++ " << tag << " ++" << std::endl << std::endl; 179 180 for (int i = 0; i != size(); ++i) 181 std::cout << std::setw(6) << i << std::setw(20) << std::setprecision(10) << *(data()+i) << std::endl; 182 } 183 allreduce()184 void allreduce() { 185 mpi__->allreduce(data(), size()); 186 } 187 188 }; 189 190 using VectorB = Vector_<double>; 191 using ZVectorB = Vector_<std::complex<double>>; 192 193 } 194 195 #endif 196