MADNESS  0.10.1
exchangeoperator.h
Go to the documentation of this file.
1 #ifndef SRC_APPS_CHEM_EXCHANGEOPERATOR_H_
2 #define SRC_APPS_CHEM_EXCHANGEOPERATOR_H_
3 
4 #include<madness.h>
5 #include<madness/world/cloud.h>
8 
9 namespace madness {
10 
11 // forward declaration
12 class SCF;
13 class Nemo;
14 
15 
16 template<typename T, std::size_t NDIM>
19  typedef std::vector<functionT> vecfuncT;
20 
21  static inline std::atomic<long> apply_timer;
22  static inline std::atomic<long> mul2_timer;
23  static inline std::atomic<long> mul1_timer; ///< timing
24 
25  static void reset_timer() {
26  mul1_timer = 0l;
27  mul2_timer = 0l;
28  apply_timer = 0l;
29  }
30 
31  static void print_timer(World& world) {
32  double t1 = double(mul1_timer) * 0.001;
33  double t2 = double(apply_timer) * 0.001;
34  double t3 = double(mul2_timer) * 0.001;
35  world.gop.sum(t1);
36  world.gop.sum(t2);
37  world.gop.sum(t3);
38  if (world.rank() == 0) {
39  printf(" cpu time spent in multiply1 %8.2fs\n", t1);
40  printf(" cpu time spent in apply %8.2fs\n", t2);
41  printf(" cpu time spent in multiply2 %8.2fs\n", t3);
42  }
43  }
44 
45 public:
46 
49 
50  /// default ctor
51  ExchangeImpl(World& world, const double lo, const double thresh) : world(world), lo(lo), thresh(thresh) {}
52 
53  /// ctor with a conventional calculation
54  ExchangeImpl(World& world, const SCF *calc, const int ispin) ;
55 
56  /// ctor with a nemo calculation
57  ExchangeImpl(World& world, const Nemo *nemo, const int ispin);
58 
59  /// set the bra and ket orbital spaces, and the occupation
60 
61  /// @param[in] bra bra space, must be provided as complex conjugate
62  /// @param[in] ket ket space
63  void set_bra_and_ket(const vecfuncT& bra, const vecfuncT& ket) {
64  mo_bra = copy(world, bra);
65  mo_ket = copy(world, ket);
66  }
67 
68  std::string info() const {return "K";}
69 
70  static auto set_poisson(World& world, const double lo, const double econv = FunctionDefaults<3>::get_thresh()) {
71  return std::shared_ptr<real_convolution_3d>(CoulombOperatorPtr(world, lo, econv));
72  }
73 
74  /// apply the exchange operator on a vector of functions
75 
76  /// note that only one spin is used (either alpha or beta orbitals)
77  /// @param[in] vket the orbitals |i> that the operator is applied on
78  /// @return a vector of orbitals K| i>
79  vecfuncT operator()(const vecfuncT& vket) const;
80 
81  bool is_symmetric() const { return symmetric_; }
82 
83  ExchangeImpl& set_taskq(std::shared_ptr<MacroTaskQ> taskq1) {
84  this->taskq=taskq1;
85  return *this;
86  }
87 
88  ExchangeImpl& symmetric(const bool flag) {
89  symmetric_ = flag;
90  return *this;
91  }
92 
94  algorithm_ = alg;
95  return *this;
96  }
97 
98  ExchangeImpl& set_printlevel(const long& level) {
99  printlevel=level;
100  return *this;
101  }
102 
103 private:
104 
105  /// exchange using macrotasks, i.e. apply K on a function in individual worlds
106  vecfuncT K_macrotask_efficient(const vecfuncT& vket, const double mul_tol = 0.0) const;
107 
108  /// computing the full square of the double sum (over vket and the K orbitals)
109  vecfuncT K_small_memory(const vecfuncT& vket, const double mul_tol = 0.0) const;
110 
111  /// computing the upper triangle of the double sum (over vket and the K orbitals)
112  vecfuncT K_large_memory(const vecfuncT& vket, const double mul_tol = 0.0) const;
113 
114  /// computing the upper triangle of the double sum (over vket and the K orbitals)
115  static vecfuncT compute_K_tile(World& world, const vecfuncT& mo_bra, const vecfuncT& mo_ket,
116  const vecfuncT& vket, std::shared_ptr<real_convolution_3d> poisson,
117  const bool symmetric, const double mul_tol = 0.0);
118 
119  inline bool do_print_timings() const { return (world.rank() == 0) and (printlevel >= 3); }
120 
121  inline bool printdebug() const {return printlevel >= 10; }
122 
124  std::shared_ptr<MacroTaskQ> taskq;
125  bool symmetric_ = false; /// is the exchange matrix symmetric? K phi_i = \sum_k \phi_k \int \phi_k \phi_i
126  vecfuncT mo_bra, mo_ket; ///< MOs for bra and ket
127  double lo = 1.e-4;
129  long printlevel = 0;
130  double mul_tol = 0.0;
131 
133 
134  long nresult;
135  double lo = 1.e-4;
136  double mul_tol = 1.e-7;
137  bool symmetric = false;
138 
139  /// custom partitioning for the exchange operator in exchangeoperator.h
140 
141  /// arguments are: result[i] += sum_k vket[k] \int 1/r vbra[k] f[i]
142  /// with f and vbra being batched, result and vket being passed on as a whole
144  public:
145  MacroTaskPartitionerExchange(const bool symmetric) : symmetric(symmetric) {
146  max_batch_size=30;
147  }
148 
149  bool symmetric = false;
150 
151  partitionT do_partitioning(const std::size_t& vsize1, const std::size_t& vsize2,
152  const std::string policy) const override {
153 
154  partitionT partition1 = do_1d_partition(vsize1, policy);
155  partitionT partition2 = do_1d_partition(vsize2, policy);
156  partitionT result;
157  for (auto i = partition1.begin(); i != partition1.end(); ++i) {
158  if (symmetric) {
159  for (auto j = i; j != partition1.end(); ++j) {
160  Batch batch(i->first.input[0], j->first.input[0], _);
161  double priority=compute_priority(batch);
162  result.push_back(std::make_pair(batch,priority));
163  }
164  } else {
165  for (auto j = partition2.begin(); j != partition2.end(); ++j) {
166  Batch batch(i->first.input[0], j->first.input[0], _);
167  double priority=compute_priority(batch);
168  result.push_back(std::make_pair(batch,priority));
169  }
170  }
171  }
172  return result;
173  }
174 
175  /// compute the priority of this task for non-dumb scheduling
176 
177  /// \return the priority as double number (no limits)
178  double compute_priority(const Batch& batch) const override {
179  MADNESS_CHECK(batch.input.size() == 2); // must be quadratic batches
180  long nrow = batch.input[0].size();
181  long ncol = batch.input[1].size();
182  return double(nrow * ncol);
183  }
184  };
185 
186  public:
187  MacroTaskExchangeSimple(const long nresult, const double lo, const double mul_tol, const bool symmetric)
188  : nresult(nresult), lo(lo), mul_tol(mul_tol), symmetric(symmetric) {
189  partitioner.reset(new MacroTaskPartitionerExchange(symmetric));
190  }
191 
192 
193  // you need to define the exact argument(s) of operator() as tuple
194  typedef std::tuple<const std::vector<Function<T, NDIM>>&,
195  const std::vector<Function<T, NDIM>>&,
196  const std::vector<Function<T, NDIM>>&> argtupleT;
197 
198  using resultT = std::vector<Function<T, NDIM>>;
199 
200  // you need to define an empty constructor for the result
201  // resultT must implement operator+=(const resultT&)
202  resultT allocator(World& world, const argtupleT& argtuple) const {
203  std::size_t n = std::get<0>(argtuple).size();
204  resultT result = zero_functions_compressed<T, NDIM>(world, n);
205  return result;
206  }
207 
208  std::vector<Function<T, NDIM>>
209  operator()(const std::vector<Function<T, NDIM>>& vf_batch, // will be batched (column)
210  const std::vector<Function<T, NDIM>>& bra_batch, // will be batched (row)
211  const std::vector<Function<T, NDIM>>& vket) { // will not be batched
212 
213  World& world = vf_batch.front().world();
214  resultT Kf = zero_functions_compressed<T, NDIM>(world, nresult);
215 
216  bool diagonal_block = batch.input[0] == batch.input[1];
217  auto& bra_range = batch.input[1]; // corresponds to vbra
218  auto& vf_range = batch.input[0]; // corresponds to vf_batch
219 
220  if (vf_range.is_full_size()) vf_range.end = vf_batch.size();
221  if (bra_range.is_full_size()) bra_range.end = bra_batch.size();
222 
223  MADNESS_CHECK(vf_range.end <= nresult);
224  if (symmetric) MADNESS_CHECK(bra_range.end <= nresult);
225 
226  if (symmetric and diagonal_block) {
227  auto ket_batch = bra_range.copy_batch(vket);
228  vecfuncT resultcolumn = compute_diagonal_batch_in_symmetric_matrix(world, ket_batch, bra_batch,
229  vf_batch);
230 
231  for (int i = vf_range.begin; i < vf_range.end; ++i)
232  Kf[i] += resultcolumn[i - vf_range.begin];
233 
234  } else if (symmetric and not diagonal_block) {
235  auto[resultcolumn, resultrow]=compute_offdiagonal_batch_in_symmetric_matrix(world, vket, bra_batch,
236  vf_batch);
237 
238  for (int i = bra_range.begin; i < bra_range.end; ++i)
239  Kf[i] += resultcolumn[i - bra_range.begin];
240  for (int i = vf_range.begin; i < vf_range.end; ++i)
241  Kf[i] += resultrow[i - vf_range.begin];
242  } else {
243  auto ket_batch = bra_range.copy_batch(vket);
244  vecfuncT resultcolumn = compute_batch_in_asymmetric_matrix(world, ket_batch, bra_batch, vf_batch);
245  for (int i = vf_range.begin; i < vf_range.end; ++i)
246  Kf[i] += resultcolumn[i - vf_range.begin];
247  }
248  return Kf;
249  }
250 
251  /// compute a batch of the exchange matrix, with identical ranges, exploiting the matrix symmetry
252 
253  /// \param subworld the world we're computing in
254  /// \param cloud where to store the results
255  /// \param bra_batch the bra batch of orbitals (including the nuclear correlation factor square)
256  /// \param ket_batch the ket batch of orbitals, i.e. the orbitals to premultiply with
257  /// \param vf_batch the argument of the exchange operator
259  const vecfuncT& ket_batch, // is batched
260  const vecfuncT& bra_batch, // is batched
261  const vecfuncT& vf_batch // is batched
262  ) const {
263  double mul_tol = 0.0;
264  double symmetric = true;
265  auto poisson = Exchange<double, 3>::ExchangeImpl::set_poisson(subworld, lo);
266  return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, poisson, symmetric,
267  mul_tol);
268  }
269 
270  /// compute a batch of the exchange matrix, with non-identical ranges
271 
272  /// \param subworld the world we're computing in
273  /// \param cloud where to store the results
274  /// \param bra_batch the bra batch of orbitals (including the nuclear correlation factor square)
275  /// \param ket_batch the ket batch of orbitals, i.e. the orbitals to premultiply with
276  /// \param vf_batch the argument of the exchange operator
278  const vecfuncT& ket_batch,
279  const vecfuncT& bra_batch,
280  const vecfuncT& vf_batch) const {
281  double mul_tol = 0.0;
282  double symmetric = false;
283  auto poisson = Exchange<double, 3>::ExchangeImpl::set_poisson(subworld, lo);
284  return Exchange<T, NDIM>::ExchangeImpl::compute_K_tile(subworld, bra_batch, ket_batch, vf_batch, poisson, symmetric,
285  mul_tol);
286  }
287 
288  /// compute a batch of the exchange matrix, with non-identical ranges
289 
290  /// \param subworld the world we're computing in
291  /// \param cloud where to store the results
292  /// \param bra_batch the bra batch of orbitals (including the nuclear correlation factor square)
293  /// \param ket_batch the ket batch of orbitals, i.e. the orbitals to premultiply with
294  /// \param vf_batch the argument of the exchange operator
295  std::pair<vecfuncT, vecfuncT> compute_offdiagonal_batch_in_symmetric_matrix(World& subworld,
296  const vecfuncT& ket, // not batched
297  const vecfuncT& bra_batch, // batched
298  const vecfuncT& vf_batch) const; // batched
299 
300  };
301 
302 };
303 
304 } /* namespace madness */
305 
306 #endif /* SRC_APPS_CHEM_EXCHANGEOPERATOR_H_ */
Operators for the molecular HF and DFT code.
a batch consists of a 2D-input batch and a 1D-output batch: K-batch <- (I-batch, J-batch)
Definition: macrotaskpartitioner.h:132
std::vector< Batch_1D > input
Definition: macrotaskpartitioner.h:135
custom partitioning for the exchange operator in exchangeoperator.h
Definition: exchangeoperator.h:143
double compute_priority(const Batch &batch) const override
compute the priority of this task for non-dumb scheduling
Definition: exchangeoperator.h:178
MacroTaskPartitionerExchange(const bool symmetric)
Definition: exchangeoperator.h:145
partitionT do_partitioning(const std::size_t &vsize1, const std::size_t &vsize2, const std::string policy) const override
override this if you want your own partitioning
Definition: exchangeoperator.h:151
vecfuncT compute_diagonal_batch_in_symmetric_matrix(World &subworld, const vecfuncT &ket_batch, const vecfuncT &bra_batch, const vecfuncT &vf_batch) const
compute a batch of the exchange matrix, with identical ranges, exploiting the matrix symmetry
Definition: exchangeoperator.h:258
MacroTaskExchangeSimple(const long nresult, const double lo, const double mul_tol, const bool symmetric)
Definition: exchangeoperator.h:187
long nresult
Definition: exchangeoperator.h:134
std::vector< Function< T, NDIM > > resultT
Definition: exchangeoperator.h:198
vecfuncT compute_batch_in_asymmetric_matrix(World &subworld, const vecfuncT &ket_batch, const vecfuncT &bra_batch, const vecfuncT &vf_batch) const
compute a batch of the exchange matrix, with non-identical ranges
Definition: exchangeoperator.h:277
std::tuple< const std::vector< Function< T, NDIM > > &, const std::vector< Function< T, NDIM > > &, const std::vector< Function< T, NDIM > > & > argtupleT
Definition: exchangeoperator.h:196
resultT allocator(World &world, const argtupleT &argtuple) const
Definition: exchangeoperator.h:202
std::vector< Function< T, NDIM > > operator()(const std::vector< Function< T, NDIM >> &vf_batch, const std::vector< Function< T, NDIM >> &bra_batch, const std::vector< Function< T, NDIM >> &vket)
Definition: exchangeoperator.h:209
Definition: exchangeoperator.h:17
static std::atomic< long > mul1_timer
timing
Definition: exchangeoperator.h:23
static void print_timer(World &world)
Definition: exchangeoperator.h:31
Exchange< T, NDIM >::Algorithm Algorithm
Definition: exchangeoperator.h:47
ExchangeImpl & set_taskq(std::shared_ptr< MacroTaskQ > taskq1)
Definition: exchangeoperator.h:83
ExchangeImpl & set_algorithm(const Algorithm &alg)
Definition: exchangeoperator.h:93
static void reset_timer()
Definition: exchangeoperator.h:25
vecfuncT mo_bra
is the exchange matrix symmetric? K phi_i = \sum_k \phi_k \int \phi_k \phi_i
Definition: exchangeoperator.h:126
World & world
Definition: exchangeoperator.h:123
std::shared_ptr< MacroTaskQ > taskq
Definition: exchangeoperator.h:124
ExchangeImpl & symmetric(const bool flag)
Definition: exchangeoperator.h:88
Function< T, NDIM > functionT
Definition: exchangeoperator.h:18
bool is_symmetric() const
Definition: exchangeoperator.h:81
ExchangeImpl & set_printlevel(const long &level)
Definition: exchangeoperator.h:98
std::vector< functionT > vecfuncT
Definition: exchangeoperator.h:19
static vecfuncT compute_K_tile(World &world, const vecfuncT &mo_bra, const vecfuncT &mo_ket, const vecfuncT &vket, std::shared_ptr< real_convolution_3d > poisson, const bool symmetric, const double mul_tol=0.0)
computing the upper triangle of the double sum (over vket and the K orbitals)
Definition: exchangeoperator.cc:167
static std::atomic< long > apply_timer
Definition: exchangeoperator.h:21
ExchangeImpl(World &world, const double lo, const double thresh)
default ctor
Definition: exchangeoperator.h:51
std::string info() const
Definition: exchangeoperator.h:68
bool do_print_timings() const
Definition: exchangeoperator.h:119
bool printdebug() const
Definition: exchangeoperator.h:121
static auto set_poisson(World &world, const double lo, const double econv=FunctionDefaults< 3 >::get_thresh())
Definition: exchangeoperator.h:70
static std::atomic< long > mul2_timer
Definition: exchangeoperator.h:22
void set_bra_and_ket(const vecfuncT &bra, const vecfuncT &ket)
set the bra and ket orbital spaces, and the occupation
Definition: exchangeoperator.h:63
Definition: SCFOperators.h:104
Function< T, NDIM > operator()(const Function< T, NDIM > &ket) const
Definition: SCFOperators.h:150
Algorithm
Definition: SCFOperators.h:116
@ multiworld_efficient
Definition: SCFOperators.h:117
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
A multiresolution adaptive numerical function.
Definition: mra.h:122
Definition: macrotaskq.h:716
partition one (two) vectors into 1D (2D) batches.
Definition: macrotaskpartitioner.h:190
std::list< std::pair< Batch, double > > partitionT
Definition: macrotaskpartitioner.h:194
The Nemo class.
Definition: nemo.h:326
std::vector< functionT > vecfuncT
Definition: SCFOperators.h:62
std::shared_ptr< MacroTaskQ > taskq
Definition: SCFOperators.h:70
Definition: SCF.h:187
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:870
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
WorldGopInterface & gop
Global operations.
Definition: world.h:205
Declares the Cloud class for storing data and transfering them between worlds.
static double lo
Definition: dirac-hatom.cc:23
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
Declares the macrotaskq and MacroTaskBase classes.
General header file for using MADNESS.
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:190
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
static const Slice _(0,-1, 1)
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition: operator.h:1762
static const double thresh
Definition: rk.cc:45
static const std::size_t NDIM
Definition: testpdiff.cc:42