MADNESS  0.10.1
operator.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 */
31 
32 #ifndef MADNESS_MRA_OPERATOR_H__INCLUDED
33 #define MADNESS_MRA_OPERATOR_H__INCLUDED
34 
35 /// \file mra/operator.h
36 /// \brief Implements most functionality of separated operators
37 
38 /// \ingroup function
39 
40 #include <type_traits>
41 #include <limits.h>
42 #include <madness/mra/adquad.h>
43 #include <madness/tensor/aligned.h>
44 #include <madness/tensor/tensor_lapack.h>
45 #include <madness/constants.h>
46 
47 #include <madness/mra/simplecache.h>
48 #include <madness/mra/convolution1d.h>
49 #include <madness/mra/displacements.h>
50 #include <madness/mra/function_common_data.h>
51 #include <madness/mra/gfit.h>
52 #include <madness/mra/operatorinfo.h>
53 
54 namespace madness {
55 
56  template<typename T, std::size_t NDIM>
57  class SeparatedConvolution;
58 
59  template<typename T, std::size_t NDIM>
60  class CCPairFunction;
61 
62  template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
63  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
64  apply(const SeparatedConvolution<T,KDIM>& op, const std::vector< Function<R,NDIM> > f);
65 
66  template<typename T, std::size_t NDIM>
67  std::vector<CCPairFunction<T,NDIM>> apply(const SeparatedConvolution<T,NDIM>& op, const std::vector<CCPairFunction<T,NDIM>>& argument);
68 
69  template<typename T, std::size_t NDIM>
70  std::vector<CCPairFunction<T,NDIM>> apply(const SeparatedConvolution<T,NDIM/2>& op, const std::vector<CCPairFunction<T,NDIM>>& argument);
71 
72  template<typename T, std::size_t NDIM>
73  CCPairFunction<T,NDIM> apply(const SeparatedConvolution<T,NDIM>& op, const CCPairFunction<T,NDIM>& argument);
74 
75  template<typename T, std::size_t NDIM>
76  CCPairFunction<T,NDIM> apply(const SeparatedConvolution<T,NDIM/2>& op, const CCPairFunction<T,NDIM>& argument);
77 
78  /// SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement
79  /// Why is this here?? Why don't you just use ConvolutionND in SeparatedConvolutionData??
80  template <typename Q, std::size_t NDIM>
81  struct SeparatedConvolutionInternal {
82  double norm;
83  const ConvolutionData1D<Q>* ops[NDIM];
84  };
85 
86  /// SeparatedConvolutionData keeps data for all terms, all dimensions
87 
88  /// this struct is used to cache the data that are generated by
89  template <typename Q, std::size_t NDIM>
90  struct SeparatedConvolutionData {
91  std::vector< SeparatedConvolutionInternal<Q,NDIM> > muops;
92  double norm;
93 
94  SeparatedConvolutionData(int rank) : muops(rank), norm(0.0) {}
95  SeparatedConvolutionData(const SeparatedConvolutionData<Q,NDIM>& q) {
96  muops = q.muops;
97  norm = q.norm;
98  }
99  };
100 
101 
102  /// Convolutions in separated form (including Gaussian)
103 
104  /* this stuff is very confusing, poorly commented, and extremely poorly named!
105 
106  I think it works like this:
107  We try to apply transition matrices to the compressed form of function coefficients.
108  Most of the code is about caching these transition matrices. They are cached (key of the map is the displacement)
109  in the SimpleCache "data", which is of type SeparatedConvolutionData, which keeps the matrices
110  for all separated terms and dimensions. These SeparatedConvolutionData are constructed using
111  ConvolutionND "ops", which is constructed at the construction of the SeparatedConvolution.
112 
113  SeparatedConvolution (all terms, all dim, all displacements)
114 
115  construction storage
116 
117  SimpleCache<SeparatedConvolutionData>
118  (all terms, all dim) / (all disp)
119  vector<ConvolutionND>
120  (1 term, all dim) / (all terms)
121  vector<SeparatedConvolutionInternal>
122  (1 term, all dim) / (all terms)
123 
124 
125  vector<ConvolutionData1D>
126  (1 term, 1 dim) / (all dim)
127 
128  ConvolutionND and SeparatedConvolutionInternal both point to the same data in ConvolutionData1D.
129  Why we need SeparatedConvolutionInternal in the first place I have no idea. ConvolutionND has the global
130  factor, and SeparatedConvolutionInternal has a norm.
131 
132 
133  */
134 
135  template <typename Q, std::size_t NDIM>
136  class SeparatedConvolution : public WorldObject< SeparatedConvolution<Q,NDIM> > {
137  public:
138 
139  typedef Q opT; ///< The apply function uses this to infer resultT=opT*inputT
140 
141  OperatorInfo info;
142 
143  bool doleaves; ///< If should be applied to leaf coefficients ... false by default
144  bool isperiodicsum;///< If true the operator 1D kernels have been summed over lattice translations
145  ///< and may be non-zero at both ends of the unit cell
146  bool modified_=false; ///< use modified NS form
147  int particle_=1; ///< must only be 1 or 2
148  bool destructive_=false; ///< destroy the argument or restore it (expensive for 6d functions)
149  bool print_timings=false;
150 
151  typedef Key<NDIM> keyT;
152  const static size_t opdim=NDIM;
153  Timer timer_full;
154  Timer timer_low_transf;
155  Timer timer_low_accumulate;
156  Timer timer_stats_accumulate;
157 
158  private:
159 
160 
161  mutable std::vector< ConvolutionND<Q,NDIM> > ops; ///< ConvolutionND keeps data for 1 term, all dimensions, 1 displacement
162  const BoundaryConditions<NDIM> bc;
163  const int k;
164  const FunctionCommonData<Q,NDIM>& cdata;
165  int rank;
166  const std::vector<long> vk;
167  const std::vector<long> v2k;
168  const std::vector<Slice> s0;
169 
170  // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
171  mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, NDIM > data; ///< cache for all terms, dims and displacements
172  mutable SimpleCache< SeparatedConvolutionData<Q,NDIM>, 2*NDIM > mod_data; ///< cache for all terms, dims and displacements
173 
174  public:
175 
176  bool& modified() {return modified_;}
177  const bool& modified() const {return modified_;}
178 
179  int& particle() {return particle_;}
180  const int& particle() const {return particle_;}
181  SeparatedConvolution<Q,NDIM>& set_particle(const int p) {
182  if (p!=1 and p!=2) throw std::runtime_error("particle must be 1 or 2");
183  particle_=p;
184  return *this;
185  }
186 
187  bool& destructive() {return destructive_;}
188  const bool& destructive() const {return destructive_;}
189 
190  const double& gamma() const {return info.mu;}
191  const double& mu() const {return info.mu;}
192 
193  private:
194 
195  /// laziness for calling lists: which terms to apply
196  struct ApplyTerms {
197  ApplyTerms() : r_term(false), t_term(false) {}
198  bool r_term;
199  bool t_term;
200  bool any_terms() const {return r_term or t_term;}
201  };
202 
203  /// too lazy for extended calling lists
204  struct Transformation {
205  long r; // Effective rank of transformation
206  const Q* U; // Ptr to matrix
207  const Q* VT;
208  };
209 
210  static inline std::pair<Tensor<Q>,Tensor<Q>>
211  make_coeff_for_operator(World& world, double mu, double lo, double eps, OpType type,
212  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc()) {
213 
214  OperatorInfo info(mu,lo,eps,type);
215  return make_coeff_for_operator(world, info, bc);
216 // const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
217 // double hi = cell_width.normf(); // Diagonal width of cell
218 // if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
219 //
220 // OperatorInfo info(mu,lo,eps,type);
221 // info.hi=hi;
222 // GFit<Q,NDIM> fit(info);
223 //
224 // Tensor<Q> coeff=fit.coeffs();
225 // Tensor<Q> expnt=fit.exponents();
226 //
227 // if (bc(0,0) == BC_PERIODIC) {
228 // fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
229 // }
230 //
231 // return std::make_pair(coeff,expnt);
232  }
233 
234  static inline std::pair<Tensor<double>,Tensor<double>>
235  make_coeff_for_operator(World& world, OperatorInfo info,
236  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc()) {
237 
238  const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
239  double hi = cell_width.normf(); // Diagonal width of cell
240  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
241 
242  info.hi=hi;
243  GFit<Q,NDIM> fit(info);
244 
245  Tensor<Q> coeff=fit.coeffs();
246  Tensor<Q> expnt=fit.exponents();
247 
248  if (bc(0,0) == BC_PERIODIC) {
249  fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
250  }
251 
252  return std::make_pair(coeff,expnt);
253  }
254 
255 // /// return the right block of the upsampled operator (modified NS only)
256 //
257 // /// unlike the operator matrices on the natural level the upsampled operator
258 // /// matrices are not Toeplitz, so we need more information than just the displacement
259 // ///.@param[in] source the source key
260 // /// @param[in] disp the displacement
261 // /// @param[in] upop the unfiltered operator matrix from scale n-1
262 // /// @return (k,k) patch of the upop(2k,2k) matrix
263 // static Tensor<Q> operator_patch(const Translation& source, const Translation& disp, const Tensor<Q>& upop) {
264 //
265 // // which of the 4 upsampled matrices do we need?
266 // Translation sx=source%2; // source offset
267 // Translation tx=(source+disp)%2; // target offset
268 //
269 // Tensor<Q> rij(k,k);
270 // // those two are equivalent:
271 ///*
272 // if (sx==0 and tx==0) copy_2d_patch(rij.ptr(), k, upop.ptr(), 2*k, k, k);
273 // if (sx==1 and tx==0) copy_2d_patch(rij.ptr() + k, k, upop.ptr(), 2*k, k, k);
274 // if (sx==0 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k, k, upop.ptr(), 2*k, k, k);
275 // if (sx==1 and tx==1) copy_2d_patch(rij.ptr() + 2*k*k + k, k, upop.ptr(), 2*k, k, k);
276 //*/
277 // Slice s0(0,k-1), s1(k,2*k-1);
278 // if (sx==0 and tx==0) rij=Rm(s0,s0);
279 // if (sx==1 and tx==0) rij=Rm(s1,s0);
280 // if (sx==0 and tx==1) rij=Rm(s0,s1);
281 // if (sx==1 and tx==1) rij=Rm(s1,s1);
282 //
283 // return rij;
284 // }
285 
286 
287 
288  /// accumulate into result
289  template <typename T, typename R>
290  void apply_transformation(long dimk,
291  const Transformation trans[NDIM],
292  const Tensor<T>& f,
293  Tensor<R>& work1,
294  Tensor<R>& work2,
295  const Q mufac,
296  Tensor<R>& result) const {
297 
298  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
299  long size = 1;
300  for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
301  long dimi = size/dimk;
302 
303  R* MADNESS_RESTRICT w1=work1.ptr();
304  R* MADNESS_RESTRICT w2=work2.ptr();
305 
306 #ifdef HAVE_IBMBGQ
307  mTxmq_padding(dimi, trans[0].r, dimk, dimk, w1, f.ptr(), trans[0].U);
308 #else
309  mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), trans[0].U, dimk);
310 #endif
311 
312  size = trans[0].r * size / dimk;
313  dimi = size/dimk;
314  for (std::size_t d=1; d<NDIM; ++d) {
315 #ifdef HAVE_IBMBGQ
316  mTxmq_padding(dimi, trans[d].r, dimk, dimk, w2, w1, trans[d].U);
317 #else
318  mTxmq(dimi, trans[d].r, dimk, w2, w1, trans[d].U, dimk);
319 #endif
320  size = trans[d].r * size / dimk;
321  dimi = size/dimk;
322  std::swap(w1,w2);
323  }
324 
325  // If all blocks are full rank we can skip the transposes
326  bool doit = false;
327  for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;
328 
329  if (doit) {
330  for (std::size_t d=0; d<NDIM; ++d) {
331  if (trans[d].VT) {
332  dimi = size/trans[d].r;
333 #ifdef HAVE_IBMBGQ
334  mTxmq_padding(dimi, dimk, trans[d].r, dimk, w2, w1, trans[d].VT);
335 #else
336  mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
337 #endif
338  size = dimk*size/trans[d].r;
339  }
340  else {
341  fast_transpose(dimk, dimi, w1, w2);
342  }
343  std::swap(w1,w2);
344  }
345  }
346  // Assuming here that result is contiguous and aligned
347  aligned_axpy(size, result.ptr(), w1, mufac);
348  }
349 
350 
351  /// accumulate into result
352  template <typename T, typename R>
353  void apply_transformation3(const Tensor<T> trans2[NDIM],
354  const Tensor<T>& f,
355  const Q mufac,
356  Tensor<R>& result) const {
357 
358  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
359 
360  Tensor<R> result2=general_transform(f,trans2);
361  result2.scale(mufac);
362  result+=result2;
363 
364  }
365 
366 
367 
368  /// don't accumulate, since we want to do this at apply()
369  template <typename T, typename R>
370  void apply_transformation2(Level n, long dimk, double tol,
371  const Tensor<T> trans2[NDIM],
372  const GenTensor<T>& f,
373  GenTensor<R>& work1,
374  GenTensor<R>& work2,
375  const Q mufac,
376  GenTensor<R>& result) const {
377 
378  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
379 
380 #if 1
381  result=general_transform(f,trans2);
382  result.scale(mufac);
383 
384 #else
385 
386  long size = 1;
387  for (std::size_t i=0; i<NDIM; ++i) size *= dimk;
388  long dimi = size/dimk;
389 
390  R* MADNESS_RESTRICT w1=work1.ptr();
391  R* MADNESS_RESTRICT w2=work2.ptr();
392 
393  mTxmq(dimi, trans[0].r, dimk, w1, f.ptr(), trans[0].U, dimk);
394  size = trans[0].r * size / dimk;
395  dimi = size/dimk;
396  for (std::size_t d=1; d<NDIM; ++d) {
397  mTxmq(dimi, trans[d].r, dimk, w2, w1, trans[d].U, dimk);
398  size = trans[d].r * size / dimk;
399  dimi = size/dimk;
400  std::swap(w1,w2);
401  }
402 
403  // If all blocks are full rank we can skip the transposes
404  bool doit = false;
405  for (std::size_t d=0; d<NDIM; ++d) doit = doit || trans[d].VT;
406 
407  if (doit) {
408  for (std::size_t d=0; d<NDIM; ++d) {
409  if (trans[d].VT) {
410  dimi = size/trans[d].r;
411  mTxmq(dimi, dimk, trans[d].r, w2, w1, trans[d].VT);
412  size = dimk*size/trans[d].r;
413  }
414  else {
415  fast_transpose(dimk, dimi, w1, w2);
416  }
417  std::swap(w1,w2);
418  }
419  }
420  // Assuming here that result is contiguous and aligned
421  aligned_axpy(size, result.ptr(), w1, mufac);
422  // long one = 1;
423  //daxpy_(&size, &mufac, w1, &one, result.ptr(), &one);
424 #endif
425  }
426 
427 
428  /// Apply one of the separated terms, accumulating into the result
429  template <typename T>
430  void muopxv_fast(ApplyTerms at,
431  const ConvolutionData1D<Q>* const ops_1d[NDIM],
432  const Tensor<T>& f, const Tensor<T>& f0,
433  Tensor<TENSOR_RESULT_TYPE(T,Q)>& result,
434  Tensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
435  double tol,
436  const Q mufac,
437  Tensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
438  Tensor<TENSOR_RESULT_TYPE(T,Q)>& work2) const {
439 
440  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
441  Transformation trans[NDIM];
442  Tensor<T> trans2[NDIM];
443 
444  double Rnorm = 1.0;
445  for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;
446 
447  if (at.r_term and (Rnorm > 1.e-20)) {
448 
449  tol = tol/(Rnorm*NDIM); // Errors are relative within here
450 
451  // Determine rank of SVD to use or if to use the full matrix
452  long twok = 2*k;
453  if (modified()) twok=k;
454 
455  long break_even;
456  if (NDIM==1) break_even = long(0.5*twok);
457  else if (NDIM==2) break_even = long(0.6*twok);
458  else if (NDIM==3) break_even=long(0.65*twok);
459  else break_even=long(0.7*twok);
460  bool rank_is_zero = false;
461  for (std::size_t d=0; d<NDIM; ++d) {
462  long r;
463  for (r=0; r<twok; ++r) {
464  if (ops_1d[d]->Rs[r] < tol) break;
465  }
466  if (r >= break_even) {
467  trans[d].r = twok;
468  trans[d].U = ops_1d[d]->R.ptr();
469  trans[d].VT = 0;
470  }
471  else {
472 
473 #ifdef USE_GENTENSOR
474  r = std::max(2L,r+(r&1L)); // (needed for 6D == when GENTENSOR is on) NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
475 #endif
476  if (r == 0) {
477  rank_is_zero = true;
478  break;
479  }
480  trans[d].r = r;
481  trans[d].U = ops_1d[d]->RU.ptr();
482  trans[d].VT = ops_1d[d]->RVT.ptr();
483  }
484  trans2[d]=ops_1d[d]->R;
485  }
486 
487  if (!rank_is_zero)
488  apply_transformation(twok, trans, f, work1, work2, mufac, result);
489 
490  // apply_transformation2(n, twok, tol, trans2, f, work1, work2, mufac, result);
491 // apply_transformation3(trans2, f, mufac, result);
492  }
493 
494  double Tnorm = 1.0;
495  for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;
496 
497  if (at.t_term and (Tnorm>0.0)) {
498  tol = tol/(Tnorm*NDIM); // Errors are relative within here
499 
500  long break_even;
501  if (NDIM==1) break_even = long(0.5*k);
502  else if (NDIM==2) break_even = long(0.6*k);
503  else if (NDIM==3) break_even=long(0.65*k);
504  else break_even=long(0.7*k);
505  bool rank_is_zero = false;
506  for (std::size_t d=0; d<NDIM; ++d) {
507  long r;
508  for (r=0; r<k; ++r) {
509  if (ops_1d[d]->Ts[r] < tol) break;
510  }
511  if (r >= break_even) {
512  trans[d].r = k;
513  trans[d].U = ops_1d[d]->T.ptr();
514  trans[d].VT = 0;
515  }
516  else {
517 
518 #ifdef USE_GENTENSOR
519  r = std::max(2L,r+(r&1L)); // (needed for 6D == GENTENSOR is USED) NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
520 #endif
521  if (r == 0) {
522  rank_is_zero = true;
523  break;
524  }
525  trans[d].r = r;
526  trans[d].U = ops_1d[d]->TU.ptr();
527  trans[d].VT = ops_1d[d]->TVT.ptr();
528  }
529  trans2[d]=ops_1d[d]->T;
530  }
531  if (!rank_is_zero)
532  apply_transformation(k, trans, f0, work1, work2, -mufac, result0);
533 // apply_transformation2(n, k, tol, trans2, f0, work1, work2, -mufac, result0);
534 // apply_transformation3(trans2, f0, -mufac, result0);
535  }
536  }
537 
538 
539  /// Apply one of the separated terms, accumulating into the result
540  template <typename T>
541  void muopxv_fast2(Level n,
542  const ConvolutionData1D<Q>* const ops_1d[NDIM],
543  const GenTensor<T>& f, const GenTensor<T>& f0,
544  GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result,
545  GenTensor<TENSOR_RESULT_TYPE(T,Q)>& result0,
546  double tol,
547  const Q mufac,
548  GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work1,
549  GenTensor<TENSOR_RESULT_TYPE(T,Q)>& work2) const {
550 
551  PROFILE_MEMBER_FUNC(SeparatedConvolution);
552 // Transformation trans[NDIM];
553  Tensor<T> trans2[NDIM];
554 // MADNESS_EXCEPTION("no muopxv_fast2",1);
555 
556  double Rnorm = 1.0;
557  for (std::size_t d=0; d<NDIM; ++d) Rnorm *= ops_1d[d]->Rnorm;
558  if (Rnorm == 0.0) return;
559 
560  if (Rnorm > 1.e-20) {
561 
562  tol = tol/(Rnorm*NDIM); // Errors are relative within here
563 
564  // Determine rank of SVD to use or if to use the full matrix
565  long twok = 2*k;
566  if (modified()) twok=k;
567 // long break_even;
568 // if (NDIM==1) break_even = long(0.5*twok);
569 // else if (NDIM==2) break_even = long(0.6*twok);
570 // else if (NDIM==3) break_even=long(0.65*twok);
571 // else break_even=long(0.7*twok);
572  for (std::size_t d=0; d<NDIM; ++d) {
573  // long r;
574  // for (r=0; r<twok; ++r) {
575  // if (ops_1d[d]->Rs[r] < tol) break;
576  // }
577 // if (r >= break_even) {
578 // trans[d].r = twok;
579 // trans[d].U = ops_1d[d]->R.ptr();
580 // trans[d].VT = 0;
581 // }
582 // else {
583 // //r += std::max(2L,r&1L); // NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
584 // trans[d].r = r;
585 // trans[d].U = ops_1d[d]->RU.ptr();
586 // trans[d].VT = ops_1d[d]->RVT.ptr();
587 // }
588  trans2[d]=ops_1d[d]->R;
589  }
590  apply_transformation2(n, twok, tol, trans2, f, work1, work2, mufac, result);
591  }
592 
593  double Tnorm = 1.0;
594  for (std::size_t d=0; d<NDIM; ++d) Tnorm *= ops_1d[d]->Tnorm;
595 
596  if (n > 0 and (Tnorm>1.e-20)) {
597 // long break_even;
598 //
599 // if (NDIM==1) break_even = long(0.5*k);
600 // else if (NDIM==2) break_even = long(0.6*k);
601 // else if (NDIM==3) break_even=long(0.65*k);
602 // else break_even=long(0.7*k);
603  for (std::size_t d=0; d<NDIM; ++d) {
604  // long r;
605  // for (r=0; r<k; ++r) {
606  // if (ops_1d[d]->Ts[r] < tol) break;
607  // }
608 // if (r >= break_even) {
609 // trans[d].r = k;
610 // trans[d].U = ops_1d[d]->T.ptr();
611 // trans[d].VT = 0;
612 // }
613 // else {
614 // //r += std::max(2L,r&1L); // NOLONGER NEED TO FORCE OPERATOR RANK TO BE EVEN
615 // trans[d].r = r;
616 // trans[d].U = ops_1d[d]->TU.ptr();
617 // trans[d].VT = ops_1d[d]->TVT.ptr();
618 // }
619  trans2[d]=ops_1d[d]->T;
620  }
621  apply_transformation2(n, k, tol, trans2, f0, work1, work2, -mufac, result0);
622  }
623  }
624 
625 
626  /// Computes the Frobenius norm of one of the separated terms ... WITHOUT FACTOR INCLUDED
627  /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
628  double munorm2(Level n, const ConvolutionData1D<Q>* ops[]) const {
629  if (modified()) return munorm2_modified(n,ops);
630  return munorm2_ns(n,ops);
631  }
632 
633  /// Computes the Frobenius norm of one of the separated terms for the NS form
634  /// ... WITHOUT FACTOR INCLUDED
635  /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
636  double munorm2_ns(Level n, const ConvolutionData1D<Q>* ops[]) const {
637  //PROFILE_MEMBER_FUNC(SeparatedConvolution);
638 
639  double prodR=1.0, prodT=1.0;
640  for (std::size_t d=0; d<NDIM; ++d) {
641  prodR *= ops[d]->Rnormf;
642  prodT *= ops[d]->Tnormf;
643 
644  }
645 // if (n) prodR = sqrt(std::max(prodR*prodR - prodT*prodT,0.0));
646 
647  // this kicks in if the line above has no numerically significant digits.
648 // if (prodR < 1e-8*prodT) {
649  double prod=1.0, sum=0.0;
650  for (std::size_t d=0; d<NDIM; ++d) {
651  double a = ops[d]->NSnormf;
652  double b = ops[d]->Tnormf;
653  double aa = std::min(a,b);
654  double bb = std::max(a,b);
655  prod *= bb;
656  if (bb > 0.0) sum +=(aa/bb);
657  }
658  if (n) prod *= sum;
659  prodR = prod;
660 // }
661 
662  return prodR;
663  }
664 
665 
666  /// Computes the operator norm of one of the separated terms of the modified NS form
667  /// ... WITHOUT FACTOR INCLUDED
668  /// compute for 1 term, all dim, 1 disp, essentially for SeparatedConvolutionInternal
669  double munorm2_modified(Level n, const ConvolutionData1D<Q>* ops_1d[]) const {
670  PROFILE_MEMBER_FUNC(SeparatedConvolution);
671 
672  // follows Eq. (21) ff of Beylkin 2008 (Beylkin Appl. Comput. Harmon. Anal. 24, pp 354)
673 
674  // we have all combinations of difference, upsampled, F terms (d, u, f),
675  // with the constraint that d is in each term exactly once. In the mixed terms (udf)
676  // we just get all possible combinations, in the pure terms (dff, duu) we have
677  // to multiply each term (dff, fdf, ffd) with (NDIM-1)!, to get the right number.
678 
679  double dff = 0.0;
680  double duu = 0.0;
681  double udf = 0.0;
682 
683  // loop over d shifting over the dimensions dxx, xdx, xxd,
684  for (size_t d=0; d<NDIM; ++d) {
685  double dff_tmp = ops_1d[d]->N_diff;
686  double duu_tmp = ops_1d[d]->N_diff;
687  double udf_tmp = ops_1d[d]->N_diff;
688 
689 
690  for (size_t dd=0; dd<NDIM; ++dd) {
691  if (dd!=d) {
692  dff_tmp *= ops_1d[dd]->N_F;
693  duu_tmp *= ops_1d[dd]->N_up;
694 
695  udf_tmp *= ops_1d[dd]->N_F;
696  for (size_t ddd=0; ddd<NDIM; ++ddd) {
697  if (ddd!=dd) udf += udf_tmp * ops_1d[ddd]->N_up;
698  }
699  }
700  }
701 
702  dff+=dff_tmp;
703  duu+=duu_tmp;
704  }
705 
706  // finalize with the factorial
707  double factorial=1.0;
708  for (int i=1; i<static_cast<int>(NDIM)-1; ++i) factorial*=double(i);
709  dff*=factorial;
710  duu*=factorial;
711 
712  // Eq. (23) of Beylkin 2008, for one separated term WITHOUT the factor
713  double norm=(dff + udf + duu) /(factorial * double(NDIM));
714 
715 // // double check
716 // if (NDIM==3) {
717 // Tensor<Q> R_full=outer(ops_1d[0]->R,outer(ops_1d[1]->R,ops_1d[2]->R));
718 // Tensor<Q> T_full=outer(ops_1d[0]->T,outer(ops_1d[1]->T,ops_1d[2]->T));
719 // double n2=(R_full-T_full).normf();
720 //// print("norm estimate, norm",norm, n2, norm<n2);
721 // norm=n2;
722 // }
723 
724  return norm;
725 
726  }
727 
728 
729  /// get the transformation matrices for 1 term and all dimensions and one displacement
730 
731  /// use ConvolutionND, which uses ConvolutionData1D to collect the transformation matrices
732  const SeparatedConvolutionInternal<Q,NDIM> getmuop(int mu, Level n, const Key<NDIM>& disp) const {
733  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
734  SeparatedConvolutionInternal<Q,NDIM> op;
735  for (std::size_t d=0; d<NDIM; ++d) {
736  op.ops[d] = ops[mu].getop(d)->nonstandard(n, disp.translation()[d]);
737  }
738  op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());
739 
740 // double newnorm = munorm2(n, op.ops);
741 // // This rescaling empirically based upon BSH separated expansion
742 // // ... needs more testing. OK also for TDSE.
743 // // All is good except for some 000 blocks which are up to sqrt(k^d) off.
744 // for (int d=0; d<NDIM; ++d) {
745 // if (disp[d] == 0) newnorm *= 0.5;
746 // else if (std::abs(disp[d]) == 1) newnorm *= 0.8;
747 // }
748 // double oldnorm = munorm(n, op.ops);
749 // if (oldnorm > 1e-13 && (newnorm < 0.5*oldnorm || newnorm > 2.0*oldnorm) )
750 // print("munorm", n, disp, mu, newnorm, oldnorm, newnorm/oldnorm);
751 
752  return op;
753  }
754 
755 
756 
757  /// get the transformation matrices for 1 term and all dimensions and one displacement
758 
759  /// use ConvolutionND, which uses ConvolutionData1D to collect the transformation matrices
760  const SeparatedConvolutionInternal<Q,NDIM>
761  getmuop_modified(int mu, Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
762  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
763 
764 
765  // SeparatedConvolutionInternal keeps data for 1 term and all dimensions
766  SeparatedConvolutionInternal<Q,NDIM> op;
767 
768  // in the modified NS form we need not only the displacement, but also the source Translation
769  // for correctly constructing the operator, b/c the operator is not Toeplitz
770 
771  // op.ops is of type ConvolutionData1D (1 term, 1 dim, 1 disp)
772  // ops[mu] is of type ConvolutionND (1 term, all dim, 1 disp)
773  for (std::size_t d=0; d<NDIM; ++d) {
774  Translation sx=source.translation()[d]; // source translation
775  Translation tx=source.translation()[d]+disp.translation()[d]; // target translation
776 
777  Key<2> op_key(n,Vector<Translation,2>{sx,tx});
778  op.ops[d] = ops[mu].getop(d)->mod_nonstandard(op_key);
779  }
780 
781  // works for both modified and not modified NS form
782  op.norm = munorm2(n, op.ops)*std::abs(ops[mu].getfac());
783 // op.norm=1.0;
784  return op;
785  }
786 
787  /// get the data for all terms and all dimensions for one displacement
788  const SeparatedConvolutionData<Q,NDIM>* getop(Level n, const Key<NDIM>& d, const Key<NDIM>& source) const {
789 
790  // in the NS form the operator depends only on the displacement
791  if (not modified()) return getop_ns(n,d);
792  return getop_modified(n, d, source);
793  }
794 
795 
796  /// get the data for all terms and all dimensions for one displacement
797 
798  /// uses SeparatedConvolutionInternal (ConvolutionND, ConvolutionData1D) to construct
799  /// the transformation matrices.
800  /// @param[in] d displacement
801  /// @return pointer to cached operator
802  const SeparatedConvolutionData<Q,NDIM>* getop_ns(Level n, const Key<NDIM>& d) const {
803  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
804  const SeparatedConvolutionData<Q,NDIM>* p = data.getptr(n,d);
805  if (p) return p;
806 
807  // get the data for each term
808  SeparatedConvolutionData<Q,NDIM> op(rank);
809  for (int mu=0; mu<rank; ++mu) {
810  // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
811  // getmuop uses ConvolutionND
812  op.muops[mu] = getmuop(mu, n, d);
813  }
814 
815  double norm = 0.0;
816  for (int mu=0; mu<rank; ++mu) {
817  const double munorm = op.muops[mu].norm;
818  norm += munorm*munorm;
819  }
820  //print("getop", n, d, norm);
821  op.norm = sqrt(norm);
822  data.set(n, d, op);
823  return data.getptr(n,d);
824  }
825 
826 
827 
828  /// get the data for all terms and all dimensions for one displacement (modified NS form)
829 
830  /// remember that the operator in the modified NS form is not Toeplitz, so we need
831  /// information about the displacement and the source key
832  /// @param[in] n level (=scale) (actually redundant, since included in source)
833  /// @param[in] disp displacement key
834  /// @param[in] source source key
835  /// @return pointer to cached operator
836  const SeparatedConvolutionData<Q,NDIM>* getop_modified(Level n, const Key<NDIM>& disp, const Key<NDIM>& source) const {
837  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
838 
839  // in the modified NS form the upsampled part of the operator depends on the modulus of the source
840  Vector<Translation,NDIM> t=source.translation();
841  for (size_t i=0; i<NDIM; ++i) t[i]=t[i]%2;
842  Key<2*NDIM> key=disp.merge_with(Key<NDIM>(source.level(),t));
843 
844  const SeparatedConvolutionData<Q,NDIM>* p = mod_data.getptr(n,key);
845  if (p) return p;
846 
847  // get the data for each term
848  // op.muops is of type SeparatedConvolutionInternal (1 term, all dim, 1 disp)
849  // getmuop uses ConvolutionND
850  SeparatedConvolutionData<Q,NDIM> op(rank);
851  for (int mu=0; mu<rank; ++mu) op.muops[mu] = getmuop_modified(mu, n, disp, source);
852 
853  double norm = 0.0;
854  for (int mu=0; mu<rank; ++mu) {
855  const double munorm = op.muops[mu].norm;
856  norm += munorm*munorm;
857  }
858 
859  op.norm = sqrt(norm);
860  mod_data.set(n, key, op);
861  return mod_data.getptr(n,key);
862  }
863 
864 
865  void check_cubic() {
866  // !!! NB ... cell volume obtained from global defaults
867  const Tensor<double>& cell_width = FunctionDefaults<NDIM>::get_cell_width();
868  // Check that the cell is cubic since currently is assumed
869  for (std::size_t d=1; d<NDIM; ++d) {
870  MADNESS_CHECK(fabs(cell_width(d)-cell_width(0L)) < 1e-14*cell_width(0L));
871  }
872  }
873 
874 
875  /// upsample some of the dimensions of coeff to its child indicated by key
876 
877  /// @param[in] coeff the coeffs of dim 2*NDIM that will be upsampled
878  /// @param[in] key the key indicating the child -- only some dimensions will be "reproductive"
879  /// @param[in] particle if 0: upsample dimensions 0-2
880  /// if 1: upsample dimensions 3-5
881  /// @return a partially upsampled coefficient tensor
882  template<typename T, size_t FDIM>
883  GenTensor<T> partial_upsample(const Key<FDIM>& key, const GenTensor<T>& coeff, const int particle) const {
884 
885  if (coeff.rank()==0) return GenTensor<T>();
886  MADNESS_ASSERT(coeff.dim(0)==k);
887  if (NDIM==coeff.ndim()) {
888  MADNESS_ASSERT(particle==1); // other particle, leave this particle unchanged
889  return coeff;
890  }
891 
892  MADNESS_ASSERT(coeff.ndim()==FDIM);
893  MADNESS_ASSERT(particle==0 or (2*NDIM==FDIM));
894 
895  // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
896  // handle the spectator dimensions with the identity matrix
897  const Tensor<T> h[2] = {cdata.h0, cdata.h1};
898  Tensor<T> identity(k,k);
899  for (int i=0; i<k; ++i) identity(i,i)=1.0;
900  Tensor<T> matrices[2*NDIM];
901 
902  // get the appropriate twoscale coefficients for each dimension
903  if (particle==0) {
904  for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
905  for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=identity;
906  } else if (particle==1) {
907  for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=identity;
908  for (size_t ii=0; ii<NDIM; ++ii) matrices[ii+NDIM]=h[key.translation()[ii+NDIM]%2];
909  } else {
910  MADNESS_EXCEPTION("unknown particle",1);
911  }
912 
913  // transform and accumulate on the result
914  const GenTensor<T> result=general_transform(coeff,matrices);
915  return result;
916  }
917 
918 
919  /// upsample the sum coefficients of level 1 to sum coeffs on level n+1
920 
921  /// specialization of the unfilter method, will transform only the sum coefficients
922  /// @param[in] key key of level n+1
923  /// @param[in] coeff sum coefficients of level n (does NOT belong to key!!)
924  /// @return sum coefficients on level n+1
925  template<typename T, size_t FDIM>
926  GenTensor<T> upsample(const Key<FDIM>& key, const GenTensor<T>& coeff) const {
927 
928  // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
929  // note there are no difference coefficients; if you want that use unfilter
930  const Tensor<T> h[2] = {cdata.h0, cdata.h1};
931  Tensor<T> matrices[FDIM];
932 
933  // get the appropriate twoscale coefficients for each dimension
934  for (size_t ii=0; ii<FDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
935 
936  // transform and accumulate on the result
937  const GenTensor<T> result=general_transform(coeff,matrices);
938  return result;
939  }
940 
941 
942  public:
943 
944  // For separated convolutions with same operator in each direction (isotropic)
945  SeparatedConvolution(World& world,
946  std::vector< std::shared_ptr< Convolution1D<Q> > >& argops,
947  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
948  long k = FunctionDefaults<NDIM>::get_k(),
949  bool doleaves = false)
950  : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
951  , info()
952  , doleaves(doleaves)
953  , isperiodicsum(bc(0,0)==BC_PERIODIC)
954  , modified_(false)
955  , particle_(1)
956  , destructive_(false)
957  , bc(bc)
958  , k(k)
959  , cdata(FunctionCommonData<Q,NDIM>::get(k))
960  , rank(argops.size())
961  , vk(NDIM,k)
962  , v2k(NDIM,2*k)
963  , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
964  {
965  // Presently we must have periodic or non-periodic in all dimensions.
966  for (std::size_t d=1; d<NDIM; ++d) {
967  MADNESS_ASSERT(bc(d,0)==bc(0,0));
968  }
969  check_cubic();
970 
971  for (unsigned int mu=0; mu < argops.size(); ++mu) {
972  this->ops.push_back(ConvolutionND<Q,NDIM>(argops[mu]));
973  }
974 
975  this->process_pending();
976  }
977 
978  // For general convolutions
979  SeparatedConvolution(World& world,
980  std::vector< ConvolutionND<Q,NDIM> >& argops,
981  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
982  long k = FunctionDefaults<NDIM>::get_k(),
983  bool doleaves = false)
984  : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
985  , info()
986  , doleaves(doleaves)
987  , isperiodicsum(bc(0,0)==BC_PERIODIC)
988  , modified_(false)
989  , particle_(1)
990  , destructive_(false)
991  , ops(argops)
992  , bc(bc)
993  , k(k)
994  , cdata(FunctionCommonData<Q,NDIM>::get(k))
995  , rank(argops.size())
996  , vk(NDIM,k)
997  , v2k(NDIM,2*k)
998  , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
999  {
1000  // Presently we must have periodic or non-periodic in all dimensions.
1001  for (std::size_t d=1; d<NDIM; ++d) {
1002  MADNESS_ASSERT(bc(d,0)==bc(0,0));
1003  }
1004  this->process_pending();
1005  }
1006 
1007  /// Constructor for Gaussian Convolutions (mostly for backward compatability)
1008  SeparatedConvolution(World& world, const OperatorInfo info1,
1009  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
1010  int k=FunctionDefaults<NDIM>::get_k(),
1011  bool doleaves = false)
1012  : SeparatedConvolution(world,Tensor<double>(0l),Tensor<double>(0l),info1.lo,info1.thresh,bc,k,doleaves,info1.mu) {
1013  info.type=info1.type;
1014  auto [coeff, expnt] =make_coeff_for_operator(world, info1, bc);
1015  rank=coeff.dim(0);
1016  ops.resize(rank);
1017  initialize(coeff,expnt);
1018  }
1019 
1020  /// Constructor for Gaussian Convolutions (mostly for backward compatability)
1021  SeparatedConvolution(World& world,
1022  const Tensor<Q>& coeff, const Tensor<double>& expnt,
1023  double lo, double thresh,
1024  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
1025  int k=FunctionDefaults<NDIM>::get_k(),
1026  bool doleaves = false,
1027  double mu=0.0)
1028  : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
1029  , info(mu,lo,thresh,OT_UNDEFINED)
1030  , doleaves(doleaves)
1031  , isperiodicsum(bc(0,0)==BC_PERIODIC)
1032  , ops(coeff.dim(0))
1033  , bc(bc)
1034  , k(k)
1035  , cdata(FunctionCommonData<Q,NDIM>::get(k))
1036  , rank(coeff.dim(0))
1037  , vk(NDIM,k)
1038  , v2k(NDIM,2*k)
1039  , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1)) {
1040  initialize(coeff,expnt);
1041  }
1042 
1043  void initialize(const Tensor<Q>& coeff, const Tensor<double>& expnt) {
1044  // Presently we must have periodic or non-periodic in all dimensions.
1045  for (std::size_t d=1; d<NDIM; ++d) {
1046  MADNESS_ASSERT(bc(d,0)==bc(0,0));
1047  }
1048 
1049  const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
1050  const double pi = constants::pi;
1051 
1052  for (int mu=0; mu<rank; ++mu) {
1053  Q c = std::pow(sqrt(expnt(mu)/pi),static_cast<int>(NDIM)); // Normalization coeff
1054 
1055  // We cache the normalized operator so the factor is the value we must multiply
1056  // by to recover the coeff we want.
1057  ops[mu].setfac(coeff(mu)/c);
1058 
1059  for (std::size_t d=0; d<NDIM; ++d) {
1060  ops[mu].setop(d,GaussianConvolution1DCache<Q>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
1061  }
1062  }
1063  }
1064 
1065  /// WSTHORNTON Constructor for Gaussian Convolutions (mostly for backward compatability)
1066  SeparatedConvolution(World& world,
1067  Vector<double,NDIM> args,
1068  const Tensor<Q>& coeff, const Tensor<double>& expnt,
1069  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
1070  int k=FunctionDefaults<NDIM>::get_k(),
1071  bool doleaves=false)
1072  : WorldObject< SeparatedConvolution<Q,NDIM> >(world)
1073  , info(0.0,0.0,0.0,OT_UNDEFINED)
1074  , doleaves(doleaves)
1075  , isperiodicsum(bc(0,0)==BC_PERIODIC)
1076  , modified_(false)
1077  , particle_(1)
1078  , destructive_(false)
1079  , ops(coeff.dim(0))
1080  , bc(bc)
1081  , k(k)
1082  , cdata(FunctionCommonData<Q,NDIM>::get(k))
1083  , rank(coeff.dim(0))
1084  , vk(NDIM,k)
1085  , v2k(NDIM,2*k)
1086  , s0(std::max<std::size_t>(2,NDIM),Slice(0,k-1))
1087  {
1088  // Presently we must have periodic or non-periodic in all dimensions.
1089  for (std::size_t d=1; d<NDIM; ++d) {
1090  MADNESS_ASSERT(bc(d,0)==bc(0,0));
1091  }
1092 
1093  const Tensor<double>& width = FunctionDefaults<NDIM>::get_cell_width();
1094 
1095  for (int mu=0; mu<rank; ++mu) {
1096  double c = std::pow(sqrt(expnt(mu)/madness::constants::pi),static_cast<int>(NDIM)); // Normalization coeff
1097  ops[mu].setfac(coeff(mu)/c);
1098  for (std::size_t d=0; d<NDIM; ++d) {
1099  double c2 = sqrt(expnt[mu]*width[d]*width[d]/madness::constants::pi);
1100  std::shared_ptr<GaussianConvolution1D<double_complex> >
1101  gcptr(new GaussianConvolution1D<double_complex>(k, c2,
1102  expnt(mu)*width[d]*width[d], 0, isperiodicsum, args[d]));
1103  ops[mu].setop(d,gcptr);
1104  }
1105  }
1106  }
1107 
1108  virtual ~SeparatedConvolution() { }
1109 
1110  void print_timer() const {
1111  if (this->get_world().rank()==0) {
1112  timer_full.print("op full tensor ");
1113  timer_low_transf.print("op low rank transform");
1114  timer_low_accumulate.print("op low rank addition ");
1115  }
1116  }
1117 
1118  void reset_timer() const {
1119  if (this->get_world().rank()==0) {
1120  timer_full.reset();
1121  timer_low_transf.reset();
1122  timer_low_accumulate.reset();
1123  }
1124  }
1125 
1126  const BoundaryConditions<NDIM>& get_bc() const {return bc;}
1127 
1128  const std::vector< Key<NDIM> >& get_disp(Level n) const {
1129  return Displacements<NDIM>().get_disp(n, isperiodicsum);
1130  }
1131 
1132  /// return the operator norm for all terms, all dimensions and 1 displacement
1133  double norm(Level n, const Key<NDIM>& d, const Key<NDIM>& source_key) const {
1134  // SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
1135 // return 1.0;
1136  return getop(n, d, source_key)->norm;
1137  }
1138 
1139  /// return that part of a hi-dim key that serves as the base for displacements of this operator
1140 
1141  /// if the function and the operator have the same dimension return key
1142  /// if the function has a higher dimension than the operator (e.g. in the exchange operator)
1143  /// return only that part of key that corresponds to the particle this operator works on
1144  /// @param[in] key hi-dim key
1145  /// @return a lo-dim part of key; typically first or second half
1146  template<size_t FDIM>
1147  typename std::enable_if<FDIM!=NDIM, Key<NDIM> >::type
1148  get_source_key(const Key<FDIM> key) const {
1149  Key<NDIM> source;
1150  Key<FDIM-NDIM> dummykey;
1151  if (particle()==1) key.break_apart(source,dummykey);
1152  if (particle()==2) key.break_apart(dummykey,source);
1153  return source;
1154  }
1155 
1156  /// return that part of a hi-dim key that serves as the base for displacements of this operator
1157 
1158  /// if the function and the operator have the same dimension return key
1159  /// if the function has a higher dimension than the operator (e.g. in the exchange operator)
1160  /// return only that part of key that corresponds to the particle this operator works on
1161  /// @param[in] key hi-dim key
1162  /// @return a lo-dim part of key; typically first or second half
1163  template<size_t FDIM>
1164  typename std::enable_if<FDIM==NDIM, Key<NDIM> >::type
1165  get_source_key(const Key<FDIM> key) const {
1166  return key;
1167  }
1168 
1169  /// apply this operator on a function f
1170 
1171  /// the operator does not need to have the same dimension as the function, e,g,
1172  /// the Poisson kernel for the exchange operator acts only on 1 electron of a
1173  /// given (pair) function.
1174  /// @param[in] f a function of same or different dimension as this operator
1175  /// @return the result function of the same dimensionality as the input function f
1176  template <typename T, size_t FDIM>
1177  Function<TENSOR_RESULT_TYPE(T,Q),FDIM> operator()(const Function<T,FDIM>& f) const {
1178  return madness::apply(*this, f);
1179  }
1180 
1181  /// apply this on a vector of functions
1182  template <typename T, size_t FDIM>
1183  std::vector<Function<TENSOR_RESULT_TYPE(T,Q),FDIM>> operator()(const std::vector<Function<T,FDIM>>& f) const {
1184  return madness::apply(*this, f);
1185  }
1186 
1187  /// apply this operator on a separable function f(1,2) = f(1) f(2)
1188 
1189  /// @param[in] f1 a function of dim LDIM
1190  /// @param[in] f2 a function of dim LDIM
1191  /// @return the result function of dim NDIM=2*LDIM: g(1,2) = G(1,1',2,2') f(1',2')
1192  template <typename T, size_t LDIM>
1193  Function<TENSOR_RESULT_TYPE(T,Q),LDIM+LDIM>
1194  operator()(const Function<T,LDIM>& f1, const Function<Q,LDIM>& f2) const {
1195  return madness::apply(*this, std::vector<Function<Q,LDIM>>({f1}),
1196  std::vector<Function<Q,LDIM>>({f2}));
1197  }
1198 
1199  /// apply this operator on a sum of separable functions f(1,2) = \sum_i f_i(1) f_i(2)
1200 
1201  /// @param[in] f1 a function of dim LDIM
1202  /// @param[in] f2 a function of dim LDIM
1203  /// @return the result function of dim NDIM=2*LDIM: g(1,2) = G(1,1',2,2') f(1',2')
1204  template <typename T, size_t LDIM>
1205  Function<TENSOR_RESULT_TYPE(T,Q),LDIM+LDIM>
1206  operator()(const std::vector<Function<T,LDIM>>& f1, const std::vector<Function<Q,LDIM>>& f2) const {
1207  return madness::apply(*this, f1, f2);
1208  }
1209 
1210  /// apply this onto another suitable argument, returning the same type
1211 
1212  /// argT must implement argT::apply(const SeparatedConvolution& op, const argT& arg)
1213  template<typename argT>
1214  argT operator()(const argT& argument) const {
1215  return madness::apply(*this,argument);
1216  }
1217 
1218 
1219  /// apply this operator on coefficients in full rank form
1220 
1221  /// @param[in] coeff source coeffs in full rank
1222  /// @param[in] source the source key
1223  /// @param[in] shift the displacement, where the source coeffs come from
1224  /// @param[in] tol thresh/#neigh*cnorm
1225  /// @return a tensor of full rank with the result op(coeff)
1226  template <typename T>
1227  Tensor<TENSOR_RESULT_TYPE(T,Q)> apply(const Key<NDIM>& source,
1228  const Key<NDIM>& shift,
1229  const Tensor<T>& coeff,
1230  double tol) const {
1231  //PROFILE_MEMBER_FUNC(SeparatedConvolution); // Too fine grain for routine profiling
1232  MADNESS_ASSERT(coeff.ndim()==NDIM);
1233 
1234  double cpu0=cpu_time();
1235 
1236  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
1237  const Tensor<T>* input = &coeff;
1238  Tensor<T> dummy;
1239 
1240  if (not modified()) {
1241  if (coeff.dim(0) == k) {
1242  // This processes leaf nodes with only scaling
1243  // coefficients ... FuncImpl::apply by default does not
1244  // apply the operator to these since for smoothing operators
1245  // it is not necessary. It is necessary for operators such
1246  // as differentiation and time evolution and will also occur
1247  // if the application of the operator widens the tree.
1248  dummy = Tensor<T>(v2k);
1249  dummy(s0) = coeff;
1250  input = &dummy;
1251  }
1252  else {
1253  MADNESS_ASSERT(coeff.dim(0)==2*k);
1254  }
1255  }
1256 
1257  tol = 0.01*tol/rank; // Error is per separated term
1258  ApplyTerms at;
1259  at.r_term=true;
1260  at.t_term=(source.level()>0);
1261 
1262  /// SeparatedConvolutionData keeps data for all terms and all dimensions and 1 displacement
1263  const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
1264 
1265  //print("sepop",source,shift,op->norm,tol);
1266 
1267  Tensor<resultT> r(v2k), r0(vk);
1268  Tensor<resultT> work1(v2k,false), work2(v2k,false);
1269 
1270  if (modified()) {
1271  r=Tensor<resultT>(vk);
1272  work1=Tensor<resultT>(vk,false);
1273  work2=Tensor<resultT>(vk,false);
1274  }
1275 
1276  const Tensor<T> f0 = copy(coeff(s0));
1277  for (int mu=0; mu<rank; ++mu) {
1278  // SeparatedConvolutionInternal keeps data for 1 term and all dimensions and 1 displacement
1279  const SeparatedConvolutionInternal<Q,NDIM>& muop = op->muops[mu];
1280  if (muop.norm > tol) {
1281  // ops is of ConvolutionND, returns data for 1 term and all dimensions
1282  Q fac = ops[mu].getfac();
1283  muopxv_fast(at, muop.ops, *input, f0, r, r0, tol/std::abs(fac), fac,
1284  work1, work2);
1285  }
1286  }
1287 
1288  r(s0).gaxpy(1.0,r0,1.0);
1289  double cpu1=cpu_time();
1290  timer_full.accumulate(cpu1-cpu0);
1291 
1292  return r;
1293  }
1294 
1295 
1296  /// apply this operator on only 1 particle of the coefficients in low rank form
1297 
1298  /// note the unfortunate mess with NDIM: here NDIM is the operator dimension, and FDIM is the
1299  /// function's dimension, whereas in the function we have OPDIM for the operator and NDIM for
1300  /// the function
1301  /// @tparam T the dimension of the function this operator is applied on. \todo MGR: Make sure info on T is correct. Was previously labeled FDIM.
1302  /// @param[in] coeff source coeffs in SVD (=optimal!) form, in high dimensionality (FDIM)
1303  /// @param[in] source the source key in low dimensionality (NDIM)
1304  /// @param[in] shift the displacement in low dimensionality (NDIM)
1305  /// @param[in] tol thresh/(#neigh*cnorm)
1306  /// @param[in] tol2 thresh/#neigh
1307  /// @return coeff result
1308  template<typename T>
1309  GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2_lowdim(const Key<NDIM>& source,
1310  const Key<NDIM>& shift, const GenTensor<T>& coeff, double tol, double tol2) const {
1311 
1312  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
1313 
1314  // prepare access to the singular vectors
1315  const SVDTensor<T>& svdcoeff=coeff.get_svdtensor();
1316 // std::vector<Slice> s(coeff.config().dim_per_vector()+1,_);
1317  std::vector<Slice> s(svdcoeff.dim_per_vector(particle()-1)+1,_);
1318  // can't use predefined slices and vectors -- they have the wrong dimension
1319  const std::vector<Slice> s00(coeff.ndim(),Slice(0,k-1));
1320 
1321  // some checks
1322  MADNESS_ASSERT(coeff.is_svd_tensor()); // for now
1323  MADNESS_ASSERT(not modified());
1324  MADNESS_ASSERT(not doleaves);
1325  MADNESS_ASSERT(coeff.dim(0)==2*k);
1326  MADNESS_ASSERT(2*NDIM==coeff.ndim());
1327 
1328  double cpu0=cpu_time();
1329  const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
1330 
1331  // some workspace
1332  Tensor<resultT> work1(v2k,false), work2(v2k,false);
1333 
1334  // sliced input and final result
1335  const GenTensor<T> f0 = copy(coeff(s00));
1336  GenTensor<resultT> final=copy(coeff);
1337  GenTensor<resultT> final0=copy(f0);
1338 
1339  tol = tol/rank*0.01; // Error is per separated term
1340  tol2= tol2/rank;
1341 
1342  // the operator norm is missing the identity working on the other particle
1343  // use as (muop.norm*exchange_norm < tol)
1344  // for some reason the screening is not working at all..
1345 // double exchange_norm=std::pow(2.0*k,1.5);
1346 
1347  for (int r=0; r<coeff.rank(); ++r) {
1348 
1349  // get the appropriate singular vector (left or right depends on particle)
1350  // and apply the full tensor muopxv_fast on it, term by term
1351  s[0]=Slice(r,r);
1352  const Tensor<T> chunk=svdcoeff.ref_vector(particle()-1)(s).reshape(v2k);
1353  const Tensor<T> chunk0=f0.get_svdtensor().ref_vector(particle()-1)(s).reshape(vk);
1354 // const double weight=std::abs(coeff.config().weights(r));
1355 
1356  // accumulate all terms of the operator for a specific term of the function
1357  Tensor<resultT> result(v2k), result0(vk);
1358 
1359  ApplyTerms at;
1360  at.r_term=true;
1361  at.t_term=source.level()>0;
1362 
1363  // this loop will return on result and result0 the terms [(P+Q) G (P+Q)]_1,
1364  // and [P G P]_1, respectively
1365  for (int mu=0; mu<rank; ++mu) {
1366  const SeparatedConvolutionInternal<Q,NDIM>& muop = op->muops[mu];
1367  Q fac = ops[mu].getfac();
1368  muopxv_fast(at, muop.ops, chunk, chunk0, result, result0,
1369  tol/std::abs(fac), fac, work1, work2);
1370  }
1371 
1372  // reinsert the transformed terms into result, leaving the other particle unchanged
1373  MADNESS_ASSERT(final.get_svdtensor().has_structure());
1374  final.get_svdtensor().ref_vector(particle()-1)(s)=result;
1375 
1376  if (source.level()>0) {
1377  final0.get_svdtensor().ref_vector(particle()-1)(s)=result0;
1378  } else {
1379  final0.get_svdtensor().ref_vector(0)(s)=0.0;
1380  final0.get_svdtensor().ref_vector(1)(s)=0.0;
1381  }
1382 
1383  }
1384  double cpu1=cpu_time();
1385  timer_low_transf.accumulate(cpu1-cpu0);
1386 
1387  double cpu00=cpu_time();
1388 
1389  final.reduce_rank(tol2*0.5);
1390  final0.reduce_rank(tol2*0.5);
1391  final(s00)+=final0;
1392  final.reduce_rank(tol2);
1393 
1394  double cpu11=cpu_time();
1395  timer_low_accumulate.accumulate(cpu11-cpu00);
1396  return final;
1397  }
1398 
1399  /// apply this operator on coefficients in low rank form
1400 
1401  /// @param[in] coeff source coeffs in SVD (=optimal!) form
1402  /// @param[in] tol thresh/#neigh*cnorm
1403  /// @param[in] tol2 thresh/#neigh
1404  template <typename T>
1405  GenTensor<TENSOR_RESULT_TYPE(T,Q)> apply2(const Key<NDIM>& source,
1406  const Key<NDIM>& shift,
1407  const GenTensor<T>& coeff,
1408  double tol, double tol2) const {
1409  PROFILE_MEMBER_FUNC(SeparatedConvolution);
1410  typedef TENSOR_RESULT_TYPE(T,Q) resultT;
1411 
1412  MADNESS_ASSERT(coeff.ndim()==NDIM);
1413  MADNESS_ASSERT(coeff.is_svd_tensor()); // we use the rank below
1414 // MADNESS_EXCEPTION("no apply2",1);
1415  const TensorType tt=TT_2D;
1416 
1417  const GenTensor<T>* input = &coeff;
1418  GenTensor<T> dummy;
1419 
1420  if (not modified()) {
1421  if (coeff.dim(0) == k) {
1422  // This processes leaf nodes with only scaling
1423  // coefficients ... FuncImpl::apply by default does not
1424  // apply the operator to these since for smoothing operators
1425  // it is not necessary. It is necessary for operators such
1426  // as differentiation and time evolution and will also occur
1427  // if the application of the operator widens the tree.
1428  dummy = GenTensor<T>(v2k,TT_2D);
1429  dummy(s0) += coeff;
1430  input = &dummy;
1431  }
1432  else {
1433  MADNESS_ASSERT(coeff.dim(0)==2*k);
1434  }
1435  }
1436 
1437  tol = tol/rank; // Error is per separated term
1438  tol2= tol2/rank;
1439 
1440  const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
1441 
1442  GenTensor<resultT> r, r0, result, result0;
1443  GenTensor<resultT> work1(v2k,tt), work2(v2k,tt);
1444 
1445  if (modified()) {
1446  r=GenTensor<resultT>(vk,tt);
1447  work1=GenTensor<resultT>(vk,tt);
1448  work2=GenTensor<resultT>(vk,tt);
1449  }
1450 
1451  // collect the results of the individual operator terms
1452  std::list<GenTensor<T> > r_list;
1453  std::list<GenTensor<T> > r0_list;
1454 
1455 // const GenTensor<T> f0 = copy(coeff(s0));
1456  const GenTensor<T> f0 = copy((*input)(s0));
1457  for (int mu=0; mu<rank; ++mu) {
1458  const SeparatedConvolutionInternal<Q,NDIM>& muop = op->muops[mu];
1459  //print("muop",source, shift, mu, muop.norm);
1460 
1461  // delta(g) < delta(T) * || f ||
1462  if (muop.norm > tol) {
1463 
1464  // get maximum rank of coeff to contribute:
1465  // delta(g) < eps < || T || * delta(f)
1466  // delta(coeff) * || T || < tol2
1467  const int r_max=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.get_svdtensor().weights_);
1468  // print("r_max",coeff.config().weights(r_max));
1469 
1470  // note that max_sigma is inclusive!
1471  if (r_max>=0) {
1472  const GenTensor<resultT> chunk=SVDTensor<resultT>(input->get_svdtensor().get_configs(0,r_max));
1473  const GenTensor<resultT> chunk0=SVDTensor<resultT>(f0.get_svdtensor().get_configs(0,r_max));
1474 
1475  double cpu0=cpu_time();
1476 
1477  Q fac = ops[mu].getfac();
1478  muopxv_fast2(source.level(), muop.ops, chunk, chunk0, r, r0,
1479  tol/std::abs(fac), fac, work1, work2);
1480  double cpu1=cpu_time();
1481  timer_low_transf.accumulate(cpu1-cpu0);
1482 
1483  r_list.push_back(r);
1484  r0_list.push_back(r0);
1485  }
1486  }
1487  }
1488 
1489  // finally accumulate all the resultant terms into one tensor
1490  double cpu0=cpu_time();
1491 
1492  result0=reduce(r0_list,tol2*rank);
1493  if (r_list.size()>0) r_list.front()(s0)+=result0;
1494  result=reduce(r_list,tol2*rank);
1495 // result.reduce_rank(tol2*rank);
1496 
1497  double cpu1=cpu_time();
1498  timer_low_accumulate.accumulate(cpu1-cpu0);
1499  timer_stats_accumulate.accumulate(result.rank());
1500  return result;
1501  }
1502 
1503  /// estimate the ratio of cost of full rank versus low rank
1504 
1505  /// @param[in] source source key
1506  /// @param[in] shift displacement
1507  /// @param[in] tol thresh/#neigh/cnorm
1508  /// @param[in] tol2 thresh/#neigh
1509  /// @return cost_ratio r=-1: no terms left
1510  /// 0<r<1: better to do full rank
1511  /// 1<r: better to do low rank
1512  template<typename T>
1513  double estimate_costs(const Key<NDIM>& source,
1514  const Key<NDIM>& shift,
1515  const GenTensor<T>& coeff,
1516  double tol, double tol2) const {
1517 
1518  if (coeff.is_full_tensor()) return 0.5;
1519  if (2*NDIM==coeff.ndim()) return 1.5;
1520  MADNESS_ASSERT(NDIM==coeff.ndim());
1521  MADNESS_ASSERT(coeff.is_svd_tensor());
1522 
1523  const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
1524 
1525  tol = tol/rank; // Error is per separated term
1526  tol2= tol2/rank;
1527 
1528  const double full_operator_cost=pow(coeff.dim(0),NDIM+1);
1529  const double low_operator_cost=pow(coeff.dim(0),NDIM/2+1);
1530  const double low_reduction_cost=pow(coeff.dim(0),NDIM/2);
1531 
1532  double full_cost=0.0;
1533  double low_cost=0.0;
1534 
1535  long initial_rank=0;
1536  long final_rank=sqrt(coeff.size())*0.05; // size=ncol*nrow; final rank is 5% of max rank
1537 
1538  for (int mu=0; mu<rank; ++mu) {
1539  const SeparatedConvolutionInternal<Q,NDIM>& muop = op->muops[mu];
1540 
1541  // delta(g) < delta(T) * || f ||
1542  if (muop.norm > tol) {
1543  // note that max_sigma is inclusive: it returns a slice w(Slice(0,i))
1544  long nterms=SRConf<T>::max_sigma(tol2/muop.norm,coeff.rank(),coeff.get_svdtensor().weights_)+1;
1545 
1546  // take only the first overlap computation of rank reduction into account
1547 // low_cost+=nterms*low_operator_cost + 2.0*nterms*nterms*low_reduction_cost;
1548  initial_rank+=nterms;
1549 
1550  full_cost+=full_operator_cost;
1551  }
1552  }
1553  low_cost=initial_rank*low_operator_cost + initial_rank*final_rank*low_reduction_cost;
1554 
1555  // include random empirical factor of 2
1556  double ratio=-1.0;
1557  if (low_cost>0.0) ratio=full_cost/low_cost;
1558 // print("nterms, full, low, full/low", full_cost, low_cost,shift.distsq(), ratio);
1559  return ratio;
1560 
1561  }
1562 
1563  /// construct the tensortrain representation of the operator
1564 
1565  /// @param[in] source source coefficient box
1566  /// @param[in] shift displacement
1567  /// @param[in] tol threshold for the TT truncation
1568  /// @param[in] do_R compute the R term of the operator (2k^d)
1569  /// @param[in] do_T compute the T term of the operator (k^d), including factor -1
1570  /// Both do_R and do_T may be used simultaneously, then the final
1571  /// operator will have dimensions (2k^d)
1572  TensorTrain<double> make_tt_representation(const Key<NDIM>& source,
1573  const Key<NDIM>& shift, double tol, bool do_R, bool do_T) const {
1574 
1575  if (not (do_R or do_T)) {
1576  print("no operator requested in make_tt_representation??");
1577  MADNESS_EXCEPTION("you're sure you know what you're doing?",1);
1578  }
1579 
1580 
1581  const SeparatedConvolutionData<Q,NDIM>* op = getop(source.level(), shift, source);
1582 
1583  // check for significant ranks since the R/T matrices' construction
1584  // might have been omitted. Tnorm is always smaller than Rnorm
1585  long lo=0,hi=rank;
1586  for (int mu=0; mu<rank; ++mu) {
1587  double Rnorm=1.0;
1588  for (std::size_t d=0; d<NDIM; ++d) Rnorm *= op->muops[mu].ops[d]->Rnorm;
1589  if (Rnorm>1.e-20) hi=mu;
1590  if ((Rnorm<1.e-20) and (mu<hi)) lo=mu;
1591  }
1592  hi++;lo++;
1593 
1594  // think about dimensions
1595  long rank_eff=(hi-lo); // R or T matrices
1596  long step=1;
1597  if (do_R and do_T) { // R and T matrices
1598  rank_eff*=2;
1599  step*=2;
1600  }
1601 
1602  long k2k=k; // T matrices
1603  if (do_R) k2k=2*k; // R matrices
1604 
1605 
1606  // construct empty TT cores and fill them with the significant R/T matrices
1607  std::vector<Tensor<double> > cores(NDIM,Tensor<double>(rank_eff,k2k,k2k,rank_eff));
1608  cores[0]=Tensor<double>(k2k,k2k,rank_eff);
1609  cores[NDIM-1]=Tensor<double>(rank_eff,k2k,k2k);
1610 
1611 
1612  for (int mu=lo, r=0; mu<hi; ++mu, ++r) {
1613  const SeparatedConvolutionInternal<Q,NDIM>& muop = op->muops[mu];
1614  const Q fac = ops[mu].getfac();
1615  const Slice sr0(step*r, step*r, 0);
1616  const Slice sr1(step*r+step-1,step*r+step-1,0);
1617  const Slice s00(0,k-1,1);
1618 
1619  if (do_R) {
1620  cores[0](_, _ ,sr0)=muop.ops[0]->R;
1621  for (std::size_t idim=1; idim<NDIM-1; ++idim) {
1622  cores[idim](sr0,_ ,_ ,sr0)=muop.ops[idim]->R;
1623  }
1624  cores[NDIM-1](sr0,_ ,_ )=muop.ops[NDIM-1]->R*fac;
1625  }
1626 
1627  if (do_T) {
1628  cores[0](s00,s00,sr1)=muop.ops[0]->T;
1629  for (std::size_t idim=1; idim<NDIM-1; ++idim) {
1630  cores[idim](sr1,s00,s00,sr1)=muop.ops[idim]->T;
1631  }
1632  cores[NDIM-1](sr1,s00,s00)=muop.ops[NDIM-1]->T*(-fac);
1633  }
1634  }
1635 
1636  // construct TT representation
1637  TensorTrain<double> tt(cores);
1638 
1639  // need to reshape for the TT truncation
1640  tt.make_tensor();
1641  tt.truncate(tol*GenTensor<double>::fac_reduce());
1642  tt.make_operator();
1643 
1644  return tt;
1645  }
1646 
1647 
1648  static bool can_combine(const SeparatedConvolution<Q,NDIM>& left, const SeparatedConvolution<Q,NDIM>& right) {
1649  return (combine_OT(left,right).type!=OT_UNDEFINED);
1650  }
1651 
1652  /// return operator type and other info of the combined operator (e.g. fg = f(1,2)* g(1,2)
1653  static OperatorInfo combine_OT(const SeparatedConvolution<Q,NDIM>& left, const SeparatedConvolution<Q,NDIM>& right) {
1654  OperatorInfo info=left.info;
1655  if ((left.info.type==OT_F12) and (right.info.type==OT_G12)) {
1656  info.type=OT_FG12;
1657  } else if ((left.info.type==OT_GAUSS) and (right.info.type==OT_GAUSS)) {
1658  info=right.info;
1659  info.type=OT_GAUSS;
1660  info.mu=2.0*right.info.mu;
1661  } else if ((left.info.type==OT_SLATER) and (right.info.type==OT_SLATER)) {
1662  info=right.info;
1663  info.type=OT_SLATER;
1664  info.mu=2.0*right.info.mu;
1665  } else if ((left.info.type==OT_G12) and (right.info.type==OT_F12)) {
1666  info=right.info;
1667  info.type=OT_FG12;
1668  } else if ((left.info.type==OT_G12) and (right.info.type==OT_F212)) {
1669  info=right.info;
1670  info.type=OT_F2G12;
1671  } else if (((left.info.type==OT_F212) and (right.info.type==OT_G12)) or
1672  ((left.info.type==OT_F12) and (right.info.type==OT_FG12)) or
1673  ((left.info.type==OT_FG12) and (right.info.type==OT_F12))) {
1674  info=right.info;
1675  info.type=OT_F2G12;
1676  if (right.info.type!=OT_G12) MADNESS_CHECK(right.info.mu == left.info.mu);
1677  } else if ((left.info.type==OT_F12) and (right.info.type==OT_F12)) {
1678  info.type=OT_F212;
1679  // keep the original gamma
1680  // (f12)^2 = (1- slater12)^2 = 1/(4 gamma) (1 - 2 exp(-gamma) + exp(-2 gamma))
1681  MADNESS_CHECK(right.info.mu == left.info.mu);
1682  } else {
1683  MADNESS_EXCEPTION("unknown combination of SeparatedConvolutions: feel free to extend in operator.h",1);
1684  }
1685  return info;
1686  }
1687 
1688 
1689  /// combine 2 convolution operators to one
1690  static SeparatedConvolution<Q,NDIM> combine(const SeparatedConvolution<Q,NDIM>& left,
1691  const SeparatedConvolution<Q,NDIM>& right) {
1692  MADNESS_CHECK(can_combine(left,right));
1693  MADNESS_CHECK(left.get_world().id()==right.get_world().id());
1694 
1695  auto info=combine_OT(left,right);
1696  return SeparatedConvolution<Q,NDIM>(left.get_world(),info,left.bc,left.k);
1697  }
1698 
1699  /// combine 2 convolution operators to one
1700  friend SeparatedConvolution<Q,NDIM> combine(const std::shared_ptr<SeparatedConvolution<Q,NDIM>> left,
1701  const std::shared_ptr<SeparatedConvolution<Q,NDIM>> right) {
1702  SeparatedConvolution<Q,NDIM> result;
1703  if (left and right) {
1704  return combine(*left, *right);
1705  } else if (left) {
1706  return *left;
1707  } else if (right) {
1708  return *right;
1709  } else {
1710  MADNESS_EXCEPTION("can't combine empty SeparatedConvolutions",1);
1711  }
1712  return result;
1713  }
1714  };
1715 
1716 
1717 
1718  /// Factory function generating separated kernel for convolution with 1/r in 3D.
1719  static
1720  inline
1721  SeparatedConvolution<double_complex,3> PeriodicHFExchangeOperator(World& world,
1722  Vector<double,3> args,
1723  double lo,
1724  double eps,
1725  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1726  int k=FunctionDefaults<3>::get_k())
1727  {
1728  const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
1729  double hi = cell_width.normf(); // Diagonal width of cell
1730  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
1731 
1732  GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
1733  Tensor<double> coeff=fit.coeffs();
1734  Tensor<double> expnt=fit.exponents();
1735 
1736  if (bc(0,0) == BC_PERIODIC) {
1737  fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), true);
1738  }
1739 
1740  return SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k, false);
1741 // return SeparatedConvolution<double_complex,3>(world, coeff, expnt, bc, k);
1742 
1743  }
1744 
1745  /// Factory function generating separated kernel for convolution with 1/r in 3D.
1746  static
1747  inline
1748  SeparatedConvolution<double,3> CoulombOperator(World& world,
1749  double lo,
1750  double eps,
1751  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1752 
1753  int k=FunctionDefaults<3>::get_k())
1754  {
1755  return SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k);
1756  }
1757 
1758 
1759  /// Factory function generating separated kernel for convolution with 1/r in 3D.
1760  static
1761  inline
1762  SeparatedConvolution<double,3>* CoulombOperatorPtr(World& world,
1763  double lo,
1764  double eps,
1765  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1766  int k=FunctionDefaults<3>::get_k())
1767  {
1768  return new SeparatedConvolution<double,3>(world,OperatorInfo(0.0,lo,eps,OT_G12),bc,k);
1769  }
1770 
1771 
1772  /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
1773  template <std::size_t NDIM>
1774  static inline
1775  SeparatedConvolution<double,NDIM>
1776  BSHOperator(World& world, double mu, double lo, double eps,
1777  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
1778  int k=FunctionDefaults<NDIM>::get_k()) {
1779  if (eps>1.e-4) {
1780  if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps);
1781  MADNESS_EXCEPTION("0",1);
1782  }
1783  return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
1784  }
1785 
1786  /// Factory function generating separated kernel for convolution with BSH kernel in general NDIM
1787  template <std::size_t NDIM>
1788  static inline
1789  SeparatedConvolution<double,NDIM>*
1790  BSHOperatorPtr(World& world, double mu, double lo, double eps,
1791  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
1792  int k=FunctionDefaults<NDIM>::get_k()) {
1793  if (eps>1.e-4) {
1794  if (world.rank()==0) print("the accuracy in BSHOperator is too small, tighten the threshold",eps);
1795  MADNESS_EXCEPTION("0",1);
1796  }
1797  return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
1798  }
1799 
1800 
1801  /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
1802  static inline SeparatedConvolution<double,3>
1803  BSHOperator3D(World& world, double mu, double lo, double eps,
1804  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1805  int k=FunctionDefaults<3>::get_k()) {
1806  return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_BSH),bc,k);
1807  }
1808 
1809  /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
1810  static
1811  inline
1812  SeparatedConvolution<double_complex,3> PeriodicBSHOperator3D(World& world,
1813  Vector<double,3> args,
1814  double mu,
1815  double lo,
1816  double eps,
1817  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1818  int k=FunctionDefaults<3>::get_k())
1819 
1820  {
1821  const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
1822  double hi = cell_width.normf(); // Diagonal width of cell
1823  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
1824 
1825  GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
1826  Tensor<double> coeff=fit.coeffs();
1827  Tensor<double> expnt=fit.exponents();
1828 
1829  if (bc(0,0) == BC_PERIODIC) {
1830  fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
1831  }
1832  return SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k);
1833  }
1834 
1835  /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
1836  static inline
1837  SeparatedConvolution<double_complex,3>* PeriodicBSHOperatorPtr3D(World& world,
1838  Vector<double,3> args,
1839  double mu,
1840  double lo,
1841  double eps,
1842  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1843  int k=FunctionDefaults<3>::get_k())
1844 
1845  {
1846  const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
1847  double hi = cell_width.normf(); // Diagonal width of cell
1848  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
1849 
1850  GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
1851  Tensor<double> coeff=fit.coeffs();
1852  Tensor<double> expnt=fit.exponents();
1853 
1854  if (bc(0,0) == BC_PERIODIC) {
1855  fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
1856  }
1857  return new SeparatedConvolution<double_complex,3>(world, args, coeff, expnt, bc, k);
1858  }
1859 
1860 
1861  static inline SeparatedConvolution<double,3>
1862  SlaterF12Operator(World& world, double mu, double lo, double eps,
1863  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(), int k=FunctionDefaults<3>::get_k()) {
1864  return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k);
1865  }
1866 
1867  static inline SeparatedConvolution<double,3> SlaterF12sqOperator(World& world,
1868  double mu, double lo, double eps,
1869  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1870  int k=FunctionDefaults<3>::get_k()) {
1871  return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k);
1872  }
1873 
1874  static inline SeparatedConvolution<double,3>* SlaterF12sqOperatorPtr(World& world,
1875  double mu, double lo, double eps,
1876  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1877  int k=FunctionDefaults<3>::get_k()) {
1878  return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F212),bc,k);
1879  }
1880 
1881  /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
1882  template<std::size_t NDIM=3>
1883  static inline SeparatedConvolution<double,NDIM> SlaterOperator(World& world,
1884  double mu, double lo, double eps,
1885  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
1886  int k=FunctionDefaults<NDIM>::get_k()) {
1887  return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
1888  }
1889 
1890  /// Factory function generating separated kernel for convolution with exp(-mu*r*r)
1891 
1892  /// lo and eps are not used here
1893  template<std::size_t NDIM>
1894  static inline SeparatedConvolution<double,NDIM> GaussOperator(World& world,
1895  double mu, double lo=0.0, double eps=0.0,
1896  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
1897  int k=FunctionDefaults<NDIM>::get_k()) {
1898  return SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k);
1899  }
1900 
1901  /// Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D
1902 
1903  /// lo and eps are not used here
1904  template<std::size_t NDIM>
1905  static inline SeparatedConvolution<double, NDIM>* GaussOperatorPtr(World& world,
1906  double mu, double lo=0.0, double eps=0.0,
1907  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
1908  int k = FunctionDefaults<NDIM>::get_k()) {
1909  return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_GAUSS),bc,k);
1910  }
1911 
1912 
1913  /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
1914  /// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function
1915  template<std::size_t NDIM>
1916  static inline SeparatedConvolution<double, NDIM>* SlaterOperatorPtr_ND(World& world,
1917  double mu, double lo, double eps,
1918  const BoundaryConditions<NDIM>& bc = FunctionDefaults<NDIM>::get_bc(),
1919  int k = FunctionDefaults<NDIM>::get_k()) {
1920  return new SeparatedConvolution<double,NDIM>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
1921  }
1922 
1923  /// Factory function generating separated kernel for convolution with exp(-mu*r) in 3D
1924  /// Note that the 1/(2mu) factor of SlaterF12Operator is not included, this is just the exponential function
1925  static inline SeparatedConvolution<double, 3>* SlaterOperatorPtr(World& world,
1926  double mu, double lo, double eps,
1927  const BoundaryConditions<3>& bc = FunctionDefaults<3>::get_bc(),
1928  int k = FunctionDefaults<3>::get_k()) {
1929  return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_SLATER),bc,k);
1930  }
1931 
1932  /// Factory function generating separated kernel for convolution with (1 - exp(-mu*r))/(2 mu) in 3D
1933 
1934  /// includes the factor 1/(2 mu)
1935  static inline SeparatedConvolution<double,3>* SlaterF12OperatorPtr(World& world,
1936  double mu, double lo, double eps,
1937  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1938  int k=FunctionDefaults<3>::get_k()) {
1939  return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F12),bc,k);
1940  }
1941 
1942 
1943  /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D
1944 
1945  /// fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh
1946  /// includes the factor 1/(2 mu)
1947  static inline SeparatedConvolution<double,3>
1948  FGOperator(World& world, double mu, double lo, double eps,
1949  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1950  int k=FunctionDefaults<3>::get_k()) {
1951  return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k);
1952  }
1953 
1954  /// Factory function generating separated kernel for convolution with 1/(2 mu)*(1 - exp(-mu*r))/r in 3D
1955 
1956  /// fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh
1957  /// includes the factor 1/(2 mu)
1958  static inline SeparatedConvolution<double,3>*
1959  FGOperatorPtr(World& world, double mu, double lo, double eps,
1960  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1961  int k=FunctionDefaults<3>::get_k()) {
1962  return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_FG12),bc,k);
1963  }
1964 
1965  /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
1966 
1967  /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2 / r12
1968  /// = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ]
1969  /// includes the factor 1/(2 mu)^2
1970  static inline SeparatedConvolution<double,3>*
1971  F2GOperatorPtr(World& world, double mu, double lo, double eps,
1972  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1973  int k=FunctionDefaults<3>::get_k()) {
1974  return new SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k);
1975  }
1976 
1977  /// Factory function generating separated kernel for convolution with (1/(2 mu)*(1 - exp(-mu*r)))^2/r in 3D
1978 
1979  /// f2g = (1/(2 gamma) (1 - exp(-gamma r12)))^2 / r12
1980  /// = 1/(4 gamma) * [ 1/r12 - 2 exp(-gamma r12)/r12 + exp(-2 gamma r12)/r12 ]
1981  /// includes the factor 1/(2 mu)^2
1982  static inline SeparatedConvolution<double,3>
1983  F2GOperator(World& world, double mu, double lo, double eps,
1984  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1985  int k=FunctionDefaults<3>::get_k()) {
1986  return SeparatedConvolution<double,3>(world,OperatorInfo(mu,lo,eps,OT_F2G12),bc,k);
1987  }
1988 
1989 
1990  /// Factory function generating separated kernel for convolution a normalized
1991  /// Gaussian (aka a widened delta function)
1992  static inline SeparatedConvolution<double,3> SmoothingOperator3D(World& world,
1993  double eps,
1994  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
1995  int k=FunctionDefaults<3>::get_k()) {
1996 
1997  double exponent = 1.0/(2.0*eps);
1998  Tensor<double> coeffs(1), exponents(1);
1999  exponents(0L) = exponent;
2000  coeffs(0L)=pow(exponent/M_PI,0.5*3.0); // norm of the gaussian
2001  return SeparatedConvolution<double,3>(world, coeffs, exponents, 1.e-8, eps);
2002 
2003  }
2004 
2005  /// Factory function generating separated kernel for convolution a normalized
2006  /// Gaussian (aka a widened delta function)
2007  template<std::size_t NDIM>
2008  static inline SeparatedConvolution<double,NDIM> SmoothingOperator(World& world,
2009  double eps,
2010  const BoundaryConditions<NDIM>& bc=FunctionDefaults<NDIM>::get_bc(),
2011  int k=FunctionDefaults<NDIM>::get_k()) {
2012 
2013  double exponent = 1.0/(2.0*eps);
2014  Tensor<double> coeffs(1), exponents(1);
2015  exponents(0L) = exponent;
2016  coeffs(0L)=pow(exponent/M_PI,0.5*NDIM); // norm of the gaussian
2017  return SeparatedConvolution<double,NDIM>(world, coeffs, exponents, 1.e-8, eps);
2018  }
2019 
2020  /// Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D
2021  static
2022  inline
2023  SeparatedConvolution<double,3>* BSHOperatorPtr3D(World& world,
2024  double mu,
2025  double lo,
2026  double eps,
2027  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
2028  int k=FunctionDefaults<3>::get_k())
2029  {
2030  const Tensor<double>& cell_width = FunctionDefaults<3>::get_cell_width();
2031  double hi = cell_width.normf(); // Diagonal width of cell
2032  if (bc(0,0) == BC_PERIODIC) hi *= 100; // Extend range for periodic summation
2033 
2034  GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
2035  Tensor<double> coeff=fit.coeffs();
2036  Tensor<double> expnt=fit.exponents();
2037 
2038  if (bc(0,0) == BC_PERIODIC) {
2039  fit.truncate_periodic_expansion(coeff, expnt, cell_width.max(), false);
2040  }
2041  return new SeparatedConvolution<double,3>(world, coeff, expnt, lo, eps, bc, k);
2042  }
2043 
2044 
2045  /// Factory function generating operator for convolution with grad(1/r) in 3D
2046 
2047  /// Returns a 3-vector containing the convolution operator for the
2048  /// x, y, and z components of grad(1/r)
2049  static
2050  inline
2051  std::vector< std::shared_ptr< SeparatedConvolution<double,3> > >
2052  GradCoulombOperator(World& world,
2053  double lo,
2054  double eps,
2055  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
2056  int k=FunctionDefaults<3>::get_k())
2057  {
2058  typedef SeparatedConvolution<double,3> real_convolution_3d;
2059  typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
2060  const double pi = constants::pi;
2061  const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
2062  double hi = width.normf(); // Diagonal width of cell
2063  const bool isperiodicsum = (bc(0,0)==BC_PERIODIC);
2064  if (isperiodicsum) hi *= 100; // Extend range for periodic summation
2065 
2066  GFit<double,3> fit=GFit<double,3>::CoulombFit(lo,hi,eps,false);
2067  Tensor<double> coeff=fit.coeffs();
2068  Tensor<double> expnt=fit.exponents();
2069 
2070  if (bc(0,0) == BC_PERIODIC) {
2071  fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
2072  }
2073 
2074  int rank = coeff.dim(0);
2075 
2076  std::vector<real_convolution_3d_ptr> gradG(3);
2077 
2078  for (int dir=0; dir<3; dir++) {
2079  std::vector< ConvolutionND<double,3> > ops(rank);
2080  for (int mu=0; mu<rank; mu++) {
2081  // We cache the normalized operator so the factor is the value we must multiply
2082  // by to recover the coeff we want.
2083  double c = std::pow(sqrt(expnt(mu)/pi),3); // Normalization coeff
2084  ops[mu].setfac(coeff(mu)/c/width[dir]);
2085 
2086  for (int d=0; d<3; d++) {
2087  if (d != dir)
2088  ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
2089  }
2090  ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum));
2091  }
2092  gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
2093  }
2094 
2095  return gradG;
2096  }
2097 
2098  /// Factory function generating operator for convolution with grad(bsh) in 3D
2099 
2100  /// Returns a 3-vector containing the convolution operator for the
2101  /// x, y, and z components of grad(bsh)
2102  static
2103  inline
2104  std::vector< std::shared_ptr< SeparatedConvolution<double,3> > >
2105  GradBSHOperator(World& world,
2106  double mu,
2107  double lo,
2108  double eps,
2109  const BoundaryConditions<3>& bc=FunctionDefaults<3>::get_bc(),
2110  int k=FunctionDefaults<3>::get_k())
2111  {
2112  typedef SeparatedConvolution<double,3> real_convolution_3d;
2113  typedef std::shared_ptr<real_convolution_3d> real_convolution_3d_ptr;
2114  const double pi = constants::pi;
2115  const Tensor<double> width = FunctionDefaults<3>::get_cell_width();
2116  double hi = width.normf(); // Diagonal width of cell
2117  const bool isperiodicsum = (bc(0,0)==BC_PERIODIC);
2118  if (isperiodicsum) hi *= 100; // Extend range for periodic summation
2119 
2120  GFit<double,3> fit=GFit<double,3>::BSHFit(mu,lo,hi,eps,false);
2121  Tensor<double> coeff=fit.coeffs();
2122  Tensor<double> expnt=fit.exponents();
2123 
2124  if (bc(0,0) == BC_PERIODIC) {
2125  fit.truncate_periodic_expansion(coeff, expnt, width.max(), true);
2126  }
2127 
2128  int rank = coeff.dim(0);
2129 
2130  std::vector<real_convolution_3d_ptr> gradG(3);
2131 
2132  for (int dir=0; dir<3; dir++) {
2133  std::vector< ConvolutionND<double,3> > ops(rank);
2134  for (int mu=0; mu<rank; mu++) {
2135  // We cache the normalized operator so the factor is the value we must multiply
2136  // by to recover the coeff we want.
2137  double c = std::pow(sqrt(expnt(mu)/pi),3); // Normalization coeff
2138  ops[mu].setfac(coeff(mu)/c/width[dir]);
2139 
2140  for (int d=0; d<3; d++) {
2141  if (d != dir)
2142  ops[mu].setop(d,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[d]*width[d], 0, isperiodicsum));
2143  }
2144  ops[mu].setop(dir,GaussianConvolution1DCache<double>::get(k, expnt(mu)*width[dir]*width[dir], 1, isperiodicsum));
2145  }
2146  gradG[dir] = real_convolution_3d_ptr(new SeparatedConvolution<double,3>(world, ops));
2147  }
2148 
2149  return gradG;
2150  }
2151 
2152 
2153 
2154  namespace archive {
2155  template <class Archive, class T, std::size_t NDIM>
2156  struct ArchiveLoadImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
2157  static inline void load(const Archive& ar, const SeparatedConvolution<T,NDIM>*& ptr) {
2158  WorldObject< SeparatedConvolution<T,NDIM> >* p = nullptr;
2159  ar & p;
2160  ptr = static_cast< const SeparatedConvolution<T,NDIM>* >(p);
2161  }
2162  };
2163 
2164  template <class Archive, class T, std::size_t NDIM>
2165  struct ArchiveStoreImpl<Archive,const SeparatedConvolution<T,NDIM>*> {
2166  static inline void store(const Archive& ar, const SeparatedConvolution<T,NDIM>*const& ptr) {
2167  ar & static_cast< const WorldObject< SeparatedConvolution<T,NDIM> >* >(ptr);
2168  }
2169  };
2170  }
2171 
2172 }
2173 
2174 
2175 
2176 
2177 #endif // MADNESS_MRA_OPERATOR_H__INCLUDED
q
double q(double t)
Definition: DKops.h:18
aligned.h
Provides routines for internal use optimized for aligned data.
madness::BaseTensor::dim
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
madness::BaseTensor::ndim
long ndim() const
Returns the number of dimensions in the tensor.
Definition: basetensor.h:144
madness::BoundaryConditions
This class is used to specify boundary conditions for all operators.
Definition: funcdefaults.h:101
madness::Convolution1D
Provides the common functionality/interface of all 1D convolutions.
Definition: convolution1d.h:257
madness::ConvolutionND
Array of 1D convolutions (one / dimension)
Definition: convolution1d.h:543
madness::Displacements
Holds displacements for applying operators to avoid replicating for all operators.
Definition: displacements.h:39
madness::Displacements::get_disp
const std::vector< Key< NDIM > > & get_disp(Level n, bool isperiodicsum)
Definition: displacements.h:191
madness::FunctionCommonData::h0
Tensor< double > h0
Definition: function_common_data.h:105
madness::FunctionCommonData::h1
Tensor< double > h1
Definition: function_common_data.h:105
madness::FunctionDefaults
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
madness::FunctionDefaults::get_cell_width
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
madness::Function
A multiresolution adaptive numerical function.
Definition: mra.h:122
madness::GFit
Definition: gfit.h:57
madness::GFit::BSHFit
static GFit BSHFit(double mu, double lo, double hi, double eps, bool prnt=false)
return a fit for the bound-state Helmholtz function
Definition: gfit.h:117
madness::GFit::CoulombFit
static GFit CoulombFit(double lo, double hi, double eps, bool prnt=false)
return a fit for the Coulomb function
Definition: gfit.h:102
madness::GaussianConvolution1D
1D convolution with (derivative) Gaussian; coeff and expnt given in simulation coordinates [0,...
Definition: convolution1d.h:683
madness::GenTensor
Definition: lowranktensor.h:59
madness::GenTensor::dim
long dim(const int i) const
return the number of entries in dimension i
Definition: lowranktensor.h:391
madness::GenTensor::ndim
long ndim() const
Definition: lowranktensor.h:386
madness::GenTensor::is_full_tensor
constexpr bool is_full_tensor() const
Definition: gentensor.h:224
madness::GenTensor::scale
IsSupported< TensorTypeData< Q >, GenTensor< T > & >::type scale(Q fac)
Inplace multiplication by scalar of supported type (legacy name)
Definition: lowranktensor.h:426
madness::GenTensor::ptr
const BaseTensor * ptr() const
might return a NULL pointer!
Definition: lowranktensor.h:709
madness::GenTensor::reduce_rank
void reduce_rank(const double &eps)
Definition: gentensor.h:217
madness::GenTensor::rank
long rank() const
Definition: gentensor.h:212
madness::GenTensor::size
long size() const
Definition: lowranktensor.h:482
madness::GenTensor::get_svdtensor
SVDTensor< T > & get_svdtensor()
Definition: gentensor.h:228
madness::GenTensor::is_svd_tensor
constexpr bool is_svd_tensor() const
Definition: gentensor.h:222
madness::Key
Key is the index for a node of the 2^NDIM-tree.
Definition: key.h:66
madness::Key::translation
const Vector< Translation, NDIM > & translation() const
Definition: key.h:164
madness::Key::merge_with
Key< NDIM+LDIM > merge_with(const Key< LDIM > &rhs) const
merge with other key (ie concatenate), use level of rhs, not of this
Definition: key.h:313
madness::Key::break_apart
void break_apart(Key< LDIM > &key1, Key< KDIM > &key2) const
break key into two low-dimensional keys
Definition: key.h:269
madness::SRConf::dim_per_vector
int dim_per_vector(int idim) const
return the number of physical dimensions
Definition: srconf.h:665
madness::SRConf::max_sigma
static int max_sigma(const double &thresh, const long &rank, const Tensor< double > &w)
Definition: srconf.h:109
madness::SRConf::ref_vector
Tensor< T > & ref_vector(const unsigned int &idim)
return reference to one of the vectors F
Definition: srconf.h:530
madness::SVDTensor
Definition: SVDTensor.h:42
madness::SeparatedConvolution
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
madness::SeparatedConvolution::SeparatedConvolution
SeparatedConvolution(World &world, const OperatorInfo info1, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k(), bool doleaves=false)
Constructor for Gaussian Convolutions (mostly for backward compatability)
Definition: operator.h:1008
madness::SeparatedConvolution::timer_low_transf
Timer timer_low_transf
Definition: operator.h:154
madness::SeparatedConvolution::particle
const int & particle() const
Definition: operator.h:180
madness::SeparatedConvolution::destructive_
bool destructive_
destroy the argument or restore it (expensive for 6d functions)
Definition: operator.h:148
madness::SeparatedConvolution::particle_
int particle_
must only be 1 or 2
Definition: operator.h:147
madness::SeparatedConvolution::isperiodicsum
bool isperiodicsum
Definition: operator.h:144
madness::SeparatedConvolution::mu
const double & mu() const
Definition: operator.h:191
madness::SeparatedConvolution::muopxv_fast2
void muopxv_fast2(Level n, const ConvolutionData1D< Q > *const ops_1d[NDIM], const GenTensor< T > &f, const GenTensor< T > &f0, GenTensor< TENSOR_RESULT_TYPE(T, Q)> &result, GenTensor< TENSOR_RESULT_TYPE(T, Q)> &result0, double tol, const Q mufac, GenTensor< TENSOR_RESULT_TYPE(T, Q)> &work1, GenTensor< TENSOR_RESULT_TYPE(T, Q)> &work2) const
Apply one of the separated terms, accumulating into the result.
Definition: operator.h:541
madness::SeparatedConvolution::get_source_key
std::enable_if< FDIM==NDIM, Key< NDIM > >::type get_source_key(const Key< FDIM > key) const
return that part of a hi-dim key that serves as the base for displacements of this operator
Definition: operator.h:1165
madness::SeparatedConvolution::SeparatedConvolution
SeparatedConvolution(World &world, std::vector< ConvolutionND< Q, NDIM > > &argops, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), long k=FunctionDefaults< NDIM >::get_k(), bool doleaves=false)
Definition: operator.h:979
madness::SeparatedConvolution::ops
std::vector< ConvolutionND< Q, NDIM > > ops
ConvolutionND keeps data for 1 term, all dimensions, 1 displacement.
Definition: operator.h:161
madness::SeparatedConvolution::make_coeff_for_operator
static std::pair< Tensor< Q >, Tensor< Q > > make_coeff_for_operator(World &world, double mu, double lo, double eps, OpType type, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc())
Definition: operator.h:211
madness::SeparatedConvolution::combine
static SeparatedConvolution< Q, NDIM > combine(const SeparatedConvolution< Q, NDIM > &left, const SeparatedConvolution< Q, NDIM > &right)
combine 2 convolution operators to one
Definition: operator.h:1690
madness::SeparatedConvolution::upsample
GenTensor< T > upsample(const Key< FDIM > &key, const GenTensor< T > &coeff) const
upsample the sum coefficients of level 1 to sum coeffs on level n+1
Definition: operator.h:926
madness::SeparatedConvolution::muopxv_fast
void muopxv_fast(ApplyTerms at, const ConvolutionData1D< Q > *const ops_1d[NDIM], const Tensor< T > &f, const Tensor< T > &f0, Tensor< TENSOR_RESULT_TYPE(T, Q)> &result, Tensor< TENSOR_RESULT_TYPE(T, Q)> &result0, double tol, const Q mufac, Tensor< TENSOR_RESULT_TYPE(T, Q)> &work1, Tensor< TENSOR_RESULT_TYPE(T, Q)> &work2) const
Apply one of the separated terms, accumulating into the result.
Definition: operator.h:430
madness::SeparatedConvolution::combine
friend SeparatedConvolution< Q, NDIM > combine(const std::shared_ptr< SeparatedConvolution< Q, NDIM >> left, const std::shared_ptr< SeparatedConvolution< Q, NDIM >> right)
combine 2 convolution operators to one
Definition: operator.h:1700
madness::SeparatedConvolution::get_bc
const BoundaryConditions< NDIM > & get_bc() const
Definition: operator.h:1126
madness::SeparatedConvolution::munorm2
double munorm2(Level n, const ConvolutionData1D< Q > *ops[]) const
Definition: operator.h:628
madness::SeparatedConvolution::apply
Tensor< TENSOR_RESULT_TYPE(T, Q)> apply(const Key< NDIM > &source, const Key< NDIM > &shift, const Tensor< T > &coeff, double tol) const
apply this operator on coefficients in full rank form
Definition: operator.h:1227
madness::SeparatedConvolution::timer_stats_accumulate
Timer timer_stats_accumulate
Definition: operator.h:156
madness::SeparatedConvolution::data
SimpleCache< SeparatedConvolutionData< Q, NDIM >, NDIM > data
cache for all terms, dims and displacements
Definition: operator.h:171
madness::SeparatedConvolution::munorm2_ns
double munorm2_ns(Level n, const ConvolutionData1D< Q > *ops[]) const
Definition: operator.h:636
madness::SeparatedConvolution::~SeparatedConvolution
virtual ~SeparatedConvolution()
Definition: operator.h:1108
madness::SeparatedConvolution::mod_data
SimpleCache< SeparatedConvolutionData< Q, NDIM >, 2 *NDIM > mod_data
cache for all terms, dims and displacements
Definition: operator.h:172
madness::SeparatedConvolution::operator()
argT operator()(const argT &argument) const
apply this onto another suitable argument, returning the same type
Definition: operator.h:1214
madness::SeparatedConvolution::getop
const SeparatedConvolutionData< Q, NDIM > * getop(Level n, const Key< NDIM > &d, const Key< NDIM > &source) const
get the data for all terms and all dimensions for one displacement
Definition: operator.h:788
madness::SeparatedConvolution::timer_full
Timer timer_full
Definition: operator.h:153
madness::SeparatedConvolution::apply2
GenTensor< TENSOR_RESULT_TYPE(T, Q)> apply2(const Key< NDIM > &source, const Key< NDIM > &shift, const GenTensor< T > &coeff, double tol, double tol2) const
apply this operator on coefficients in low rank form
Definition: operator.h:1405
madness::SeparatedConvolution::opT
Q opT
The apply function uses this to infer resultT=opT*inputT.
Definition: operator.h:139
madness::SeparatedConvolution::operator()
Function< TENSOR_RESULT_TYPE(T, Q), LDIM+LDIM > operator()(const std::vector< Function< T, LDIM >> &f1, const std::vector< Function< Q, LDIM >> &f2) const
apply this operator on a sum of separable functions f(1,2) = \sum_i f_i(1) f_i(2)
Definition: operator.h:1206
madness::SeparatedConvolution::v2k
const std::vector< long > v2k
Definition: operator.h:167
madness::SeparatedConvolution::doleaves
bool doleaves
If should be applied to leaf coefficients ... false by default.
Definition: operator.h:143
madness::SeparatedConvolution::operator()
Function< TENSOR_RESULT_TYPE(T, Q), FDIM > operator()(const Function< T, FDIM > &f) const
apply this operator on a function f
Definition: operator.h:1177
madness::SeparatedConvolution::print_timer
void print_timer() const
Definition: operator.h:1110
madness::SeparatedConvolution::getop_modified
const SeparatedConvolutionData< Q, NDIM > * getop_modified(Level n, const Key< NDIM > &disp, const Key< NDIM > &source) const
get the data for all terms and all dimensions for one displacement (modified NS form)
Definition: operator.h:836
madness::SeparatedConvolution::SeparatedConvolution
SeparatedConvolution(World &world, std::vector< std::shared_ptr< Convolution1D< Q > > > &argops, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), long k=FunctionDefaults< NDIM >::get_k(), bool doleaves=false)
Definition: operator.h:945
madness::SeparatedConvolution::apply_transformation3
void apply_transformation3(const Tensor< T > trans2[NDIM], const Tensor< T > &f, const Q mufac, Tensor< R > &result) const
accumulate into result
Definition: operator.h:353
madness::SeparatedConvolution::modified
bool & modified()
Definition: operator.h:176
madness::SeparatedConvolution::getmuop_modified
const SeparatedConvolutionInternal< Q, NDIM > getmuop_modified(int mu, Level n, const Key< NDIM > &disp, const Key< NDIM > &source) const
get the transformation matrices for 1 term and all dimensions and one displacement
Definition: operator.h:761
madness::SeparatedConvolution::vk
const std::vector< long > vk
Definition: operator.h:166
madness::SeparatedConvolution::norm
double norm(Level n, const Key< NDIM > &d, const Key< NDIM > &source_key) const
return the operator norm for all terms, all dimensions and 1 displacement
Definition: operator.h:1133
madness::SeparatedConvolution::make_coeff_for_operator
static std::pair< Tensor< double >, Tensor< double > > make_coeff_for_operator(World &world, OperatorInfo info, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc())
Definition: operator.h:235
madness::SeparatedConvolution::destructive
bool & destructive()
Definition: operator.h:187
madness::SeparatedConvolution::operator()
std::vector< Function< TENSOR_RESULT_TYPE(T, Q), FDIM > > operator()(const std::vector< Function< T, FDIM >> &f) const
apply this on a vector of functions
Definition: operator.h:1183
madness::SeparatedConvolution::apply2_lowdim
GenTensor< TENSOR_RESULT_TYPE(T, Q)> apply2_lowdim(const Key< NDIM > &source, const Key< NDIM > &shift, const GenTensor< T > &coeff, double tol, double tol2) const
apply this operator on only 1 particle of the coefficients in low rank form
Definition: operator.h:1309
madness::SeparatedConvolution::reset_timer
void reset_timer() const
Definition: operator.h:1118
madness::SeparatedConvolution::combine_OT
static OperatorInfo combine_OT(const SeparatedConvolution< Q, NDIM > &left, const SeparatedConvolution< Q, NDIM > &right)
return operator type and other info of the combined operator (e.g. fg = f(1,2)* g(1,...
Definition: operator.h:1653
madness::SeparatedConvolution::SeparatedConvolution
SeparatedConvolution(World &world, const Tensor< Q > &coeff, const Tensor< double > &expnt, double lo, double thresh, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k(), bool doleaves=false, double mu=0.0)
Constructor for Gaussian Convolutions (mostly for backward compatability)
Definition: operator.h:1021
madness::SeparatedConvolution::gamma
const double & gamma() const
Definition: operator.h:190
madness::SeparatedConvolution::particle
int & particle()
Definition: operator.h:179
madness::SeparatedConvolution::print_timings
bool print_timings
Definition: operator.h:149
madness::SeparatedConvolution::timer_low_accumulate
Timer timer_low_accumulate
Definition: operator.h:155
madness::SeparatedConvolution::estimate_costs
double estimate_costs(const Key< NDIM > &source, const Key< NDIM > &shift, const GenTensor< T > &coeff, double tol, double tol2) const
estimate the ratio of cost of full rank versus low rank
Definition: operator.h:1513
madness::SeparatedConvolution::s0
const std::vector< Slice > s0
Definition: operator.h:168
madness::SeparatedConvolution::apply_transformation2
void apply_transformation2(Level n, long dimk, double tol, const Tensor< T > trans2[NDIM], const GenTensor< T > &f, GenTensor< R > &work1, GenTensor< R > &work2, const Q mufac, GenTensor< R > &result) const
don't accumulate, since we want to do this at apply()
Definition: operator.h:370
madness::SeparatedConvolution::destructive
const bool & destructive() const
Definition: operator.h:188
madness::SeparatedConvolution::partial_upsample
GenTensor< T > partial_upsample(const Key< FDIM > &key, const GenTensor< T > &coeff, const int particle) const
upsample some of the dimensions of coeff to its child indicated by key
Definition: operator.h:883
madness::SeparatedConvolution::cdata
const FunctionCommonData< Q, NDIM > & cdata
Definition: operator.h:164
madness::SeparatedConvolution::k
const int k
Definition: operator.h:163
madness::SeparatedConvolution::operator()
Function< TENSOR_RESULT_TYPE(T, Q), LDIM+LDIM > operator()(const Function< T, LDIM > &f1, const Function< Q, LDIM > &f2) const
apply this operator on a separable function f(1,2) = f(1) f(2)
Definition: operator.h:1194
madness::SeparatedConvolution::info
OperatorInfo info
Definition: operator.h:141
madness::SeparatedConvolution::set_particle
SeparatedConvolution< Q, NDIM > & set_particle(const int p)
Definition: operator.h:181
madness::SeparatedConvolution::make_tt_representation
TensorTrain< double > make_tt_representation(const Key< NDIM > &source, const Key< NDIM > &shift, double tol, bool do_R, bool do_T) const
construct the tensortrain representation of the operator
Definition: operator.h:1572
madness::SeparatedConvolution::getop_ns
const SeparatedConvolutionData< Q, NDIM > * getop_ns(Level n, const Key< NDIM > &d) const
get the data for all terms and all dimensions for one displacement
Definition: operator.h:802
madness::SeparatedConvolution::keyT
Key< NDIM > keyT
Definition: operator.h:151
madness::SeparatedConvolution::opdim
static const size_t opdim
Definition: operator.h:152
madness::SeparatedConvolution::apply_transformation
void apply_transformation(long dimk, const Transformation trans[NDIM], const Tensor< T > &f, Tensor< R > &work1, Tensor< R > &work2, const Q mufac, Tensor< R > &result) const
Definition: operator.h:290
madness::SeparatedConvolution::check_cubic
void check_cubic()
Definition: operator.h:865
madness::SeparatedConvolution::initialize
void initialize(const Tensor< Q > &coeff, const Tensor< double > &expnt)
Definition: operator.h:1043
madness::SeparatedConvolution::get_disp
const std::vector< Key< NDIM > > & get_disp(Level n) const
Definition: operator.h:1128
madness::SeparatedConvolution::modified_
bool modified_
use modified NS form
Definition: operator.h:146
madness::SeparatedConvolution::rank
int rank
Definition: operator.h:165
madness::SeparatedConvolution::bc
const BoundaryConditions< NDIM > bc
Definition: operator.h:162
madness::SeparatedConvolution::munorm2_modified
double munorm2_modified(Level n, const ConvolutionData1D< Q > *ops_1d[]) const
Definition: operator.h:669
madness::SeparatedConvolution::SeparatedConvolution
SeparatedConvolution(World &world, Vector< double, NDIM > args, const Tensor< Q > &coeff, const Tensor< double > &expnt, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k(), bool doleaves=false)
WSTHORNTON Constructor for Gaussian Convolutions (mostly for backward compatability)
Definition: operator.h:1066
madness::SeparatedConvolution::getmuop
const SeparatedConvolutionInternal< Q, NDIM > getmuop(int mu, Level n, const Key< NDIM > &disp) const
get the transformation matrices for 1 term and all dimensions and one displacement
Definition: operator.h:732
madness::SeparatedConvolution::modified
const bool & modified() const
Definition: operator.h:177
madness::SeparatedConvolution::can_combine
static bool can_combine(const SeparatedConvolution< Q, NDIM > &left, const SeparatedConvolution< Q, NDIM > &right)
Definition: operator.h:1648
madness::SeparatedConvolution::get_source_key
std::enable_if< FDIM!=NDIM, Key< NDIM > >::type get_source_key(const Key< FDIM > key) const
return that part of a hi-dim key that serves as the base for displacements of this operator
Definition: operator.h:1148
madness::SimpleCache
Simplified interface around hash_map to cache stuff for 1D.
Definition: simplecache.h:46
madness::Slice
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
madness::TensorTrain
Definition: tensortrain.h:123
madness::TensorTrain::make_operator
TensorTrain< T > & make_operator()
convert this into an operator representation (r,k',k,r)
Definition: tensortrain.h:1188
madness::TensorTrain::make_tensor
TensorTrain< T > & make_tensor()
convert this into a tensor representation (r,k,r)
Definition: tensortrain.h:1176
madness::TensorTrain::truncate
std::enable_if<!std::is_arithmetic< R >::value, void >::type truncate(double eps)
recompress and truncate this TT representation
Definition: tensortrain.h:883
madness::Tensor::ptr
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
madness::Tensor::normf
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
madness::Tensor::max
T max(long *ind=0) const
Return the maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1703
madness::Tensor::scale
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
madness::Timer
Definition: function_common_data.h:169
madness::Timer::print
void print(std::string line="") const
print timer
Definition: function_common_data.h:216
madness::Timer::accumulate
void accumulate(const double time) const
accumulate timer
Definition: function_common_data.h:183
madness::Timer::reset
void reset() const
Definition: function_common_data.h:210
madness::Vector
A simple, fixed dimension vector.
Definition: vector.h:64
madness::WorldObject
Implements most parts of a globally addressable object (via unique ID).
Definition: world_object.h:364
madness::WorldObject< SeparatedConvolution< Q, NDIM > >::world
World & world
The World this object belongs to. (Think globally, act locally).
Definition: world_object.h:381
madness::WorldObject< SeparatedConvolution< Q, NDIM > >::process_pending
void process_pending()
To be called from derived constructor to process pending messages.
Definition: world_object.h:656
madness::WorldObject< SeparatedConvolution< Q, NDIM > >::get_world
World & get_world() const
Returns a reference to the world.
Definition: world_object.h:717
madness::World
A parallel world class.
Definition: world.h:132
madness::World::rank
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
madness::World::id
unsigned long id() const
Definition: world.h:313
constants.h
Defines common mathematical and physical constants.
convolution1d.h
Compuates most matrix elements over 1D operators (including Gaussians)
R
static const double R
Definition: csqrt.cc:46
f1
double(* f1)(const coord_3d &)
Definition: derivatives.cc:55
f2
double(* f2)(const coord_3d &)
Definition: derivatives.cc:56
p
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition: derivatives.cc:72
lo
static double lo
Definition: dirac-hatom.cc:23
shift
static double shift
Definition: dirac-hatom.cc:19
function_common_data.h
gfit.h
fit isotropic functions to a set of Gaussians with controlled precision
T
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
op
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
pow
static double pow(const double *a, const double *b)
Definition: lda.h:74
max
#define max(a, b)
Definition: lda.h:51
MADNESS_CHECK
#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
MADNESS_EXCEPTION
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
MADNESS_ASSERT
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
madness::constants::pi
const double pi
Mathematical constant .
Definition: constants.h:48
madness
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
madness::BC_PERIODIC
@ BC_PERIODIC
Definition: funcdefaults.h:56
madness::FGOperatorPtr
static SeparatedConvolution< double, 3 > * FGOperatorPtr(World &world, double mu, 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/(2 mu)*(1 - exp(-mu*r))/r in 3D.
Definition: operator.h:1959
madness::PeriodicBSHOperatorPtr3D
static SeparatedConvolution< double_complex, 3 > * PeriodicBSHOperatorPtr3D(World &world, Vector< double, 3 > args, double mu, 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 exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:1837
madness::real_convolution_3d_ptr
std::shared_ptr< real_convolution_3d > real_convolution_3d_ptr
Definition: functypedefs.h:135
madness::cpu_time
static double cpu_time()
Returns the cpu time in seconds relative to an arbitrary origin.
Definition: timers.h:127
madness::BSHOperatorPtr
static SeparatedConvolution< double, NDIM > * BSHOperatorPtr(World &world, double mu, double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Factory function generating separated kernel for convolution with BSH kernel in general NDIM.
Definition: operator.h:1790
madness::real_convolution_3d
SeparatedConvolution< double, 3 > real_convolution_3d
Definition: functypedefs.h:121
madness::SlaterF12Operator
static SeparatedConvolution< double, 3 > SlaterF12Operator(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Definition: operator.h:1862
madness::F2GOperatorPtr
static SeparatedConvolution< double, 3 > * F2GOperatorPtr(World &world, double mu, 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/(2 mu)*(1 - exp(-mu*r)))^2/r in ...
Definition: operator.h:1971
madness::general_transform
GenTensor< TENSOR_RESULT_TYPE(R, Q)> general_transform(const GenTensor< R > &t, const Tensor< Q > c[])
Definition: gentensor.h:274
madness::GradCoulombOperator
static std::vector< std::shared_ptr< SeparatedConvolution< double, 3 > > > GradCoulombOperator(World &world, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating operator for convolution with grad(1/r) in 3D.
Definition: operator.h:2052
madness::apply
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
madness::fast_transpose
void fast_transpose(long n, long m, const T *a, T *MADNESS_RESTRICT b)
a(n,m) --> b(m,n) ... optimized for smallish matrices
Definition: convolution1d.h:69
madness::GradBSHOperator
static std::vector< std::shared_ptr< SeparatedConvolution< double, 3 > > > GradBSHOperator(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating operator for convolution with grad(bsh) in 3D.
Definition: operator.h:2105
madness::BSHOperator3D
static SeparatedConvolution< double, 3 > BSHOperator3D(World &world, double mu, 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 exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:1803
madness::SlaterOperatorPtr_ND
static SeparatedConvolution< double, NDIM > * SlaterOperatorPtr_ND(World &world, double mu, double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Definition: operator.h:1916
madness::FGOperator
static SeparatedConvolution< double, 3 > FGOperator(World &world, double mu, 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/(2 mu)*(1 - exp(-mu*r))/r in 3D.
Definition: operator.h:1948
madness::copy
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
madness::Translation
int64_t Translation
Definition: key.h:54
madness::SmoothingOperator3D
static SeparatedConvolution< double, 3 > SmoothingOperator3D(World &world, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Definition: operator.h:1992
madness::mTxmq_padding
void mTxmq_padding(long dimi, long dimj, long dimk, long ext_b, cT *c, const aT *a, const bT *b)
Definition: mtxmq.h:96
madness::PeriodicBSHOperator3D
static SeparatedConvolution< double_complex, 3 > PeriodicBSHOperator3D(World &world, Vector< double, 3 > args, double mu, 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 exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:1812
madness::PeriodicHFExchangeOperator
static SeparatedConvolution< double_complex, 3 > PeriodicHFExchangeOperator(World &world, Vector< double, 3 > args, 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:1721
madness::GaussOperator
static SeparatedConvolution< double, NDIM > GaussOperator(World &world, double mu, double lo=0.0, double eps=0.0, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Factory function generating separated kernel for convolution with exp(-mu*r*r)
Definition: operator.h:1894
madness::_
static const Slice _(0,-1, 1)
madness::CoulombOperator
static SeparatedConvolution< double, 3 > CoulombOperator(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:1748
madness::Level
int Level
Definition: key.h:55
madness::OpType
OpType
operator types
Definition: operatorinfo.h:11
madness::OT_FG12
@ OT_FG12
1-exp(-r)
Definition: operatorinfo.h:18
madness::OT_SLATER
@ OT_SLATER
1/r
Definition: operatorinfo.h:15
madness::OT_GAUSS
@ OT_GAUSS
exp(-r)
Definition: operatorinfo.h:16
madness::OT_BSH
@ OT_BSH
(1-exp(-r))^2/r = 1/r + exp(-2r)/r - 2 exp(-r)/r
Definition: operatorinfo.h:21
madness::OT_F12
@ OT_F12
exp(-r2)
Definition: operatorinfo.h:17
madness::OT_F212
@ OT_F212
(1-exp(-r))/r
Definition: operatorinfo.h:19
madness::OT_UNDEFINED
@ OT_UNDEFINED
Definition: operatorinfo.h:12
madness::OT_G12
@ OT_G12
indicates the identity
Definition: operatorinfo.h:14
madness::OT_F2G12
@ OT_F2G12
(1-exp(-r))^2
Definition: operatorinfo.h:20
madness::SlaterF12sqOperator
static SeparatedConvolution< double, 3 > SlaterF12sqOperator(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Definition: operator.h:1867
madness::GaussOperatorPtr
static SeparatedConvolution< double, NDIM > * GaussOperatorPtr(World &world, double mu, double lo=0.0, double eps=0.0, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Factory function generating separated kernel for convolution with exp(-mu*r*r) in 3D.
Definition: operator.h:1905
madness::TENSOR_RESULT_TYPE
TENSOR_RESULT_TYPE(T, R) inner(const Function< T
Computes the scalar/inner product between two functions.
Definition: vmra.h:1059
madness::SlaterOperatorPtr
static SeparatedConvolution< double, 3 > * SlaterOperatorPtr(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Definition: operator.h:1925
madness::SmoothingOperator
static SeparatedConvolution< double, NDIM > SmoothingOperator(World &world, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Definition: operator.h:2008
madness::print
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
madness::TensorType
TensorType
low rank representations of tensors (see gentensor.h)
Definition: gentensor.h:120
madness::TT_2D
@ TT_2D
Definition: gentensor.h:120
madness::SlaterF12sqOperatorPtr
static SeparatedConvolution< double, 3 > * SlaterF12sqOperatorPtr(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Definition: operator.h:1874
madness::sum
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
madness::f
NDIM & f
Definition: mra.h:2416
madness::reduce
GenTensor< T > reduce(std::list< GenTensor< T > > &addends, double eps, bool are_optimal=false)
add all the GenTensors of a given list
Definition: gentensor.h:246
madness::F2GOperator
static SeparatedConvolution< double, 3 > F2GOperator(World &world, double mu, 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/(2 mu)*(1 - exp(-mu*r)))^2/r in ...
Definition: operator.h:1983
madness::CoulombOperatorPtr
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
madness::type
std::string type(const PairType &n)
Definition: PNOParameters.h:18
madness::SlaterOperator
static SeparatedConvolution< double, NDIM > SlaterOperator(World &world, double mu, double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Factory function generating separated kernel for convolution with exp(-mu*r) in 3D.
Definition: operator.h:1883
madness::BSHOperatorPtr3D
static SeparatedConvolution< double, 3 > * BSHOperatorPtr3D(World &world, double mu, 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 exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:2023
madness::BSHOperator
static SeparatedConvolution< double, NDIM > BSHOperator(World &world, double mu, double lo, double eps, const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), int k=FunctionDefaults< NDIM >::get_k())
Factory function generating separated kernel for convolution with BSH kernel in general NDIM.
Definition: operator.h:1776
madness::swap
void swap(Vector< T, N > &l, Vector< T, N > &r)
Swap the contents of two Vectors.
Definition: vector.h:497
madness::SlaterF12OperatorPtr
static SeparatedConvolution< double, 3 > * SlaterF12OperatorPtr(World &world, double mu, 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 - exp(-mu*r))/(2 mu) in 3D.
Definition: operator.h:1935
madness::mTxmq
void mTxmq(long dimi, long dimj, long dimk, T *MADNESS_RESTRICT c, const T *a, const T *b, long ldb=-1)
Matrix = Matrix transpose * matrix ... MKL interface version.
Definition: mxm.h:257
madness::dir
static const string dir
Definition: corepotential.cc:249
madness::aligned_axpy
static void aligned_axpy(long n, T *MADNESS_RESTRICT a, const T *MADNESS_RESTRICT b, Q s)
Definition: aligned.h:75
std
Definition: mraimpl.h:50
std::abs
static long abs(long a)
Definition: tensor.h:218
mu
const double mu
Definition: navstokes_cosines.cc:95
b
static const double b
Definition: nonlinschro.cc:119
a
static const double a
Definition: nonlinschro.cc:118
Q
double Q(double a)
Definition: relops.cc:20
c
static const double c
Definition: relops.cc:10
L
static const double L
Definition: rk.cc:46
thresh
static const double thresh
Definition: rk.cc:45
k
static const long k
Definition: rk.cc:44
data
Definition: test_ccpairfunction.cc:22
madness::ConvolutionData1D
!!! Note that if Rnormf is zero then ALL of the tensors are empty
Definition: convolution1d.h:161
madness::ConvolutionData1D::N_up
double N_up
Definition: convolution1d.h:172
madness::ConvolutionData1D::N_F
double N_F
the norms according to Beylkin 2008, Eq. (21) ff
Definition: convolution1d.h:172
madness::ConvolutionData1D::N_diff
double N_diff
Definition: convolution1d.h:172
madness::GaussianConvolution1DCache
Definition: convolution1d.h:849
madness::OperatorInfo
Definition: operatorinfo.h:58
madness::OperatorInfo::hi
double hi
Definition: operatorinfo.h:65
madness::OperatorInfo::type
OpType type
introspection
Definition: operatorinfo.h:64
madness::OperatorInfo::mu
double mu
some introspection
Definition: operatorinfo.h:61
madness::SeparatedConvolutionData
SeparatedConvolutionData keeps data for all terms, all dimensions.
Definition: operator.h:90
madness::SeparatedConvolutionData::muops
std::vector< SeparatedConvolutionInternal< Q, NDIM > > muops
Definition: operator.h:91
madness::SeparatedConvolutionData::SeparatedConvolutionData
SeparatedConvolutionData(int rank)
Definition: operator.h:94
madness::SeparatedConvolutionData::norm
double norm
Definition: operator.h:92
madness::SeparatedConvolutionData::SeparatedConvolutionData
SeparatedConvolutionData(const SeparatedConvolutionData< Q, NDIM > &q)
Definition: operator.h:95
madness::SeparatedConvolutionInternal
Definition: operator.h:81
madness::SeparatedConvolutionInternal::norm
double norm
Definition: operator.h:82
madness::SeparatedConvolutionInternal::ops
const ConvolutionData1D< Q > * ops[NDIM]
Definition: operator.h:83
madness::SeparatedConvolution::ApplyTerms
laziness for calling lists: which terms to apply
Definition: operator.h:196
madness::SeparatedConvolution::ApplyTerms::t_term
bool t_term
Definition: operator.h:199
madness::SeparatedConvolution::ApplyTerms::r_term
bool r_term
Definition: operator.h:198
madness::SeparatedConvolution::ApplyTerms::any_terms
bool any_terms() const
Definition: operator.h:200
madness::SeparatedConvolution::ApplyTerms::ApplyTerms
ApplyTerms()
Definition: operator.h:197
madness::SeparatedConvolution::Transformation
too lazy for extended calling lists
Definition: operator.h:204
madness::SeparatedConvolution::Transformation::VT
const Q * VT
Definition: operator.h:207
madness::SeparatedConvolution::Transformation::U
const Q * U
Definition: operator.h:206
madness::SeparatedConvolution::Transformation::r
long r
Definition: operator.h:205
madness::archive::ArchiveLoadImpl< Archive, const SeparatedConvolution< T, NDIM > * >::load
static void load(const Archive &ar, const SeparatedConvolution< T, NDIM > *&ptr)
Definition: operator.h:2157
madness::archive::ArchiveLoadImpl
Default load of an object via serialize(ar, t).
Definition: archive.h:666
madness::archive::ArchiveStoreImpl< Archive, const SeparatedConvolution< T, NDIM > * >::store
static void store(const Archive &ar, const SeparatedConvolution< T, NDIM > *const &ptr)
Definition: operator.h:2166
madness::archive::ArchiveStoreImpl
Default store of an object via serialize(ar, t).
Definition: archive.h:611
madness::particle
Definition: lowrankfunction.h:332
doit
void doit(World &world)
Definition: tdse.cc:921
tensor_lapack.h
Prototypes for a partial interface from Tensor to LAPACK.
factorial
int factorial(int n)
Definition: test_BSHApply.cc:14
d
void d()
Definition: test_sig.cc:79
e
void e()
Definition: test_sig.cc:75
aa
double aa
Definition: testbsh.cc:68
pi
static const double pi
Definition: testcosine.cc:6
fit
std::vector< double > fit(size_t m, size_t n, const std::vector< double > N, const std::vector< double > &f)
Definition: testfuns.cc:36
h
double h(const coord_1d &r)
Definition: testgconv.cc:68
NDIM
static const std::size_t NDIM
Definition: testpdiff.cc:42
source
double source(const coordT &r)
Definition: testperiodic.cc:48
PROFILE_MEMBER_FUNC
#define PROFILE_MEMBER_FUNC(classname)
Definition: worldprofile.h:210
Generated on Sat Oct 26 2024 16:26:21 for MADNESS by doxygen 1.9.1