MPQC  3.0.0-alpha
ordered_shells.h
1 //
2 // ordered_shells.h
3 //
4 // Copyright (C) 2014 David Hollman
5 //
6 // Author: David Hollman
7 // Maintainer: DSH
8 // Created: May 5, 2014
9 //
10 // This file is part of the SC Toolkit.
11 //
12 // The SC Toolkit is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Library General Public License as published by
14 // the Free Software Foundation; either version 2, or (at your option)
15 // any later version.
16 //
17 // The SC Toolkit is distributed in the hope that it will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Library General Public License for more details.
21 //
22 // You should have received a copy of the GNU Library General Public License
23 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
24 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
25 //
26 // The U.S. Government is granted a limited license as per AL 91-7.
27 //
28 
29 #ifndef _chemistry_qc_scf_cadf_ordered_shells_h
30 #define _chemistry_qc_scf_cadf_ordered_shells_h
31 
32 #include <array>
33 #include <mutex>
34 
35 #include <Eigen/Dense>
36 
37 #include "iters.h"
38 
39 namespace sc {
40 
42 
43  public:
44 
45  typedef std::vector<ShellIndexWithValue> index_list;
46  typedef std::unordered_set<
50  > index_set;
51 
52  typedef index_list::const_iterator index_iterator;
55 
56  private:
57 
58  // TODO Make this a shared mutex
59  mutable std::mutex insert_mtx_;
60  mutable std::mutex aux_vector_mtx_;
61 
62  index_list indices_;
63  index_set idx_set_;
64  bool sorted_ = false;
65  bool sort_by_value_ = true;
66 
67  // An auxiliary value and vector to tag along with the list
68  //Eigen::VectorXd aux_vector_;
69  //bool aux_vector_initialized_ = false;
70  double aux_value_ = 0.0;
71 
72  public:
73 
74  GaussianBasisSet* basis_ = 0;
75  GaussianBasisSet* dfbasis_ = 0;
76 
77  int nbf = 0;
78 
79  explicit OrderedShellList(bool sort_by_value = true)
80  : sort_by_value_(sort_by_value)
81  {
82  aux_value_ = 0.0;
83  //aux_vector_initialized_ = false;
84  }
85 
86  template<typename Iterable>
88  const Iterable& indices_in,
89  GaussianBasisSet* basis,
90  GaussianBasisSet* dfbasis = 0
91  )
92  : sort_by_value_(false),
93  basis_(basis), dfbasis_(dfbasis)
94  {
95  std::copy(indices_in.begin(), indices_in.end(), std::back_inserter(indices_));
96  sort(false);
97  aux_value_ = 0.0;
98  //aux_vector_initialized_ = false;
99  }
100 
102  : indices_(other.indices_),
103  idx_set_(other.idx_set_),
104  sorted_(other.sorted_),
105  sort_by_value_(other.sort_by_value_),
106  nbf(other.nbf)
107  {
108  //aux_value_ = other.aux_value_;
109  //aux_vector_initialized_ = other.aux_vector_initialized_;
110  }
111 
112  void add_to_aux_value(const double add_val) {
113  std::lock_guard<std::mutex> lg(aux_vector_mtx_);
114  aux_value_ += add_val;
115  }
116 
117  //template <typename Derived>
118  //void add_to_aux_value_vector(const double add_val, const Eigen::MatrixBase<Derived>& to_add) {
119  // std::lock_guard<std::mutex> lg(aux_vector_mtx_);
120  // if(not aux_vector_initialized_) {
121  // aux_vector_.resize(to_add.rows());
122  // aux_vector_ = decltype(aux_vector_)::Zero(to_add.rows());
123  // aux_vector_initialized_ = true;
124  // }
125  // aux_value_ += add_val;
126  // aux_vector_.noalias() += to_add;
127  //}
128 
129  //template <typename Derived>
130  //void add_to_aux_vector(const Eigen::MatrixBase<Derived>& to_add) {
131  // std::lock_guard<std::mutex> lg(aux_vector_mtx_);
132  // if(not aux_vector_initialized_) {
133  // aux_vector_.resize(to_add.rows());
134  // aux_vector_ = decltype(aux_vector_)::Zero(to_add.rows());
135  // aux_vector_initialized_ = true;
136  // }
137  // aux_vector_.noalias() += to_add;
138  //}
139 
140  //void set_aux_value(const double add_val) {
141  // aux_value_ = add_val;
142  //}
143 
144  double get_aux_value() const {
145  return aux_value_;
146  }
147 
148  //const decltype(aux_vector_)& get_aux_vector() const {
149  // MPQC_ASSERT(aux_vector_initialized_);
150  // return aux_vector_;
151  //}
152 
153  void set_basis(GaussianBasisSet* basis, GaussianBasisSet* dfbasis = 0)
154  {
155  basis_ = basis;
156  if(dfbasis) dfbasis_ = dfbasis;
157  }
158 
159  void insert(const ShellData& ish, double value = 0, int nbf_in = 0) {
160  //----------------------------------------//
161  assert(basis_ == 0 || ish.basis == basis_);
162  if(basis_ == 0) basis_ = ish.basis;
163  assert(dfbasis_ == 0 || ish.dfbasis == dfbasis_);
164  if(dfbasis_ == 0) dfbasis_ = ish.dfbasis;
165  //----------------------------------------//
166  std::lock_guard<std::mutex> lg(insert_mtx_);
167  sorted_ = false;
168  ShellIndexWithValue insert_val(ish, value);
169  const auto& found = idx_set_.find(insert_val);
170  if(found != idx_set_.end()) {
171  const double old_val = found->value;
172  if(value > old_val) found->value = value;
173  }
174  else {
175  idx_set_.insert(insert_val);
176  nbf += nbf_in;
177  }
178  }
179 
180  size_t size() const {
181  if(sorted_) return indices_.size();
182  else return idx_set_.size();
183  }
184 
185  void set_sort_by_value(bool new_srt=true) {
186  if(new_srt != sort_by_value_) {
187  sort_by_value_ = new_srt;
188  sorted_ = false;
189  }
190  }
191 
192  double value_for_index(int index) {
193  if(idx_set_.size() < indices_.size()) {
194  idx_set_.clear();
195  for(auto&& idx : indices_) {
196  idx_set_.insert(idx);
197  }
198  }
199  std::lock_guard<std::mutex> lg(insert_mtx_);
200  ShellIndexWithValue find_val(index);
201  const auto& found = idx_set_.find(find_val);
202  if(found != idx_set_.end()) {
203  return found->value;
204  }
205  else {
206  throw AlgorithmException("Index not found in list", __FILE__, __LINE__);
207  }
208  }
209 
210  void sort(bool transfer_idx_set = true) {
211  std::lock_guard<std::mutex> lg(insert_mtx_);
212  if(transfer_idx_set) {
213  assert(idx_set_.size() >= indices_.size());
214  indices_.clear();
215  std::copy(idx_set_.begin(), idx_set_.end(), std::back_inserter(indices_));
216  }
217  if(sort_by_value_) {
218  std::sort(indices_.begin(), indices_.end(),
219  [](const ShellIndexWithValue& a, const ShellIndexWithValue& b){
220  return a.value > b.value or (a.value == b.value and a.index < b.index);
221  }
222  );
223  }
224  else {
225  std::sort(indices_.begin(), indices_.end(),
226  [](const ShellIndexWithValue& a, const ShellIndexWithValue& b){
227  return a.index < b.index;
228  }
229  );
230  }
231  sorted_ = true;
232  }
233 
234  void acquire_and_sort(
235  ShellIndexWithValue* new_data,
236  size_t n_new_data
237  )
238  {
239  assert(size() == 0);
240  std::copy(new_data, new_data + n_new_data, std::back_inserter(indices_));
241  sort(false);
242  }
243 
244  template <typename ValueIterable>
245  void acquire_and_sort(
246  const ValueIterable& val_iter,
247  double cutoff = 0.0
248  )
249  {
250  // There is probably a faster way to do this, particularly for Eigen data structures
251  int index = 0;
252  indices_.clear();
253  idx_set_.clear();
254  sorted_ = false;
255  for(auto&& val : val_iter) {
256  if(val > cutoff) {
257  indices_.emplace_back(index++, val);
258  }
259  }
260  sort(false);
261  }
262 
263  void acquire_and_sort(
264  const double* vals,
265  size_t n_vals,
266  double cutoff = 0.0,
267  bool sort_by_val = true
268  )
269  {
270  indices_.clear();
271  idx_set_.clear();
272  sorted_ = false;
273  for(size_t idx = 0; idx < n_vals; ++idx) {
274  if(vals[idx] > cutoff) {
275  indices_.emplace_back(idx, vals[idx]);
276  }
277  }
278  set_sort_by_value(sort_by_val);
279  sort(false);
280  }
281 
282  template<typename IndexIterable>
283  void acquire_and_sort(
284  const IndexIterable& idxs_in,
285  const double* vals_in,
286  double cutoff = 0.0,
287  bool sort_by_val = true
288  )
289  {
290  indices_.clear();
291  idx_set_.clear();
292  typename IndexIterable::iterator idxiter = idxs_in.begin();
293  int val_off = 0;
294  for(; idxiter != idxs_in.end(); ++idxiter, ++val_off) {
295  if(vals_in[val_off] > cutoff) {
296  indices_.emplace_back(*idxiter, vals_in[val_off]);
297  }
298  }
299  set_sort_by_value(sort_by_val);
300  sort(false);
301  }
302 
303  template <typename IndexIterable>
305  intersection_with(const IndexIterable& other_indices) {
306  if(!sorted_) {
307  sort();
308  }
309  OrderedShellList rv(false);
310  std::copy(indices_.begin(), indices_.end(), std::back_inserter(rv.indices_));
311  if(sort_by_value_) {
312  rv.sort(false);
313  }
314  IndexIterable other_copy;
315  std::copy(other_indices.begin(), other_indices.end(), std::back_inserter(other_copy));
316  std::sort(other_copy.begin(), other_copy.end(), std::less<int>());
317 
318  // do the intersection
319  std::vector<ShellIndexWithValue> intersect;
320  std::set_intersection(
321  rv.indices_.begin(), rv.indices_.end(),
322  other_copy.begin(), other_copy.end(),
323  std::back_inserter(intersect), std::less<int>()
324  );
325 
326  rv.indices_.clear();
327  std::move(intersect.begin(), intersect.end(), std::back_inserter(rv.indices_));
328 
329  if(sort_by_value_) {
330  rv.sort_by_value_ = true;
331  rv.sorted_ = false;
332  rv.sort(false);
333  }
334 
335  rv.basis_ = basis_;
336  rv.dfbasis_ = dfbasis_;
337 
338  return rv;
339 
340  }
341 
342  iterator begin()
343  {
344  assert(sorted_);
345  return iterator(
346  basis_, dfbasis_,
347  index_begin()
348  );
349  }
350 
351  const_iterator begin() const
352  {
353  assert(sorted_);
354  return const_iterator(
355  basis_, dfbasis_,
356  index_begin()
357  );
358  }
359 
360  index_iterator index_begin() const
361  {
362  assert(sorted_);
363  return indices_.cbegin();
364  }
365 
366  iterator end()
367  {
368  assert(sorted_);
369  return iterator(
370  basis_, dfbasis_,
371  index_end()
372  );
373  }
374 
375  const_iterator end() const
376  {
377  assert(sorted_);
378  return const_iterator(
379  basis_, dfbasis_,
380  index_end()
381  );
382  }
383 
384  index_iterator index_end() const
385  {
386  assert(sorted_);
387  return indices_.cend();
388  }
389 
390  const index_list& unsorted_indices() {
391  std::lock_guard<std::mutex> lg(insert_mtx_);
392  assert(!sorted_);
393  if(indices_.size() != idx_set_.size()) {
394  indices_.clear();
395  std::copy(idx_set_.begin(), idx_set_.end(), std::back_inserter(indices_));
396  }
397  return indices_;
398  }
399 
400  void clear()
401  {
402  indices_.clear();
403  idx_set_.clear();
404  }
405 
406  friend std::ostream&
407  operator <<(std::ostream& out, const OrderedShellList& list);
408 
409 };
410 
412 
413 typedef enum {
414  L3 = 0,
415  L3_star = 1,
416  LB = 2,
417  Ld_over = 3,
418  Ld_under = 4
419 } LinKListName;
420 
422 {
423  std::array<OrderedShellList, 5> lists_;
424 
425  public:
426 
427  OrderedShellList& operator[](LinKListName name) { return lists_[name]; }
428 
429 };
430 
432 
433 inline range_of_shell_blocks<typename OrderedShellList::index_iterator>
435  const OrderedShellList& shlist,
436  int requirements=SameCenter,
437  int target_size=DEFAULT_TARGET_BLOCK_SIZE
438 )
439 {
440  return boost::make_iterator_range(
442  shlist.index_begin(), shlist.index_end(),
443  shlist.basis_, shlist.dfbasis_, requirements, target_size
444  ),
446  shlist.index_end(), shlist.index_end(),
447  shlist.basis_, shlist.dfbasis_, requirements, target_size
448  )
449  );
450 }
451 
453 
454 
455 } // end namespace sc
456 
457 #endif /* _chemistry_qc_scf_cadf_ordered_shells_h */
sc::basis_element_with_value_iterator
Definition: iters.h:394
sc::detail::index_equal_
Definition: iters.h:1266
sc::ShellData
Definition: iters.h:135
sc::AlgorithmException
This exception is thrown whenever a problem with an algorithm is encountered.
Definition: scexception.h:342
shell_block_iterable
Definition: cadf_attic.h:406
sc::detail::hash_< ShellIndexWithValue >
Definition: iters.h:1258
sc::OrderedShellList
Definition: ordered_shells.h:41
sc::other
SpinCase1 other(SpinCase1 S)
given 1-spin return the other 1-spin
sc::LinKListGroup
Definition: ordered_shells.h:421
sc::intersect
std::pair< I, I > intersect(const std::pair< I, I > &range1, const std::pair< I, I > &range2)
return intersect of two ranges defined as pair<start,fence>, i.e. [start, fence)
Definition: mtensor.h:50
sc::shell_block_iterator
Definition: iters.h:840
sc::GaussianBasisSet
The GaussianBasisSet class is used describe a basis set composed of atomic gaussian orbitals.
Definition: gaussbas.h:141
sc
Contains all MPQC code up to version 3.
Definition: mpqcin.h:14
sc::ShellIndexWithValue
binds an integer index + real annotation, e.g. Shell index + associated operator norm
Definition: iters.h:1229

Generated at Sun Jan 26 2020 23:23:59 for MPQC 3.0.0-alpha using the documentation package Doxygen 1.8.16.