MADNESS  0.10.1
mra.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_MRA_H__INCLUDED
33 #define MADNESS_MRA_MRA_H__INCLUDED
34 
35 /*!
36  \file mra/mra.h
37  \brief Main include file for MADNESS and defines \c Function interface
38 
39  \addtogroup mra
40 
41 */
42 
43 
44 #include <madness/world/MADworld.h>
45 #include <madness/misc/misc.h>
46 #include <madness/tensor/tensor.h>
47 
48 #define FUNCTION_INSTANTIATE_1
49 #define FUNCTION_INSTANTIATE_2
50 #define FUNCTION_INSTANTIATE_3
51 #if !defined(HAVE_IBMBGP) || !defined(HAVE_IBMBGQ)
52 #define FUNCTION_INSTANTIATE_4
53 #define FUNCTION_INSTANTIATE_5
54 #define FUNCTION_INSTANTIATE_6
55 #endif
56 
57 static const bool VERIFY_TREE = false; //true
58 
59 
60 namespace madness {
61  void startup(World& world, int argc, char** argv, bool doprint=false, bool make_stdcout_nice_to_reals = true);
62  std::string get_mra_data_dir();
63 }
64 
65 #include <madness/mra/key.h>
66 #include <madness/mra/twoscale.h>
67 #include <madness/mra/legendre.h>
68 #include <madness/mra/indexit.h>
70 #include <madness/world/worlddc.h>
73 #include <madness/mra/lbdeux.h>
74 #include <madness/mra/funcimpl.h>
75 
76 // some forward declarations
77 namespace madness {
78 
79  template<typename T, std::size_t NDIM>
80  class FunctionImpl;
81 
82  template<typename T, std::size_t NDIM>
83  class Function;
84 
85  template<typename T, std::size_t NDIM>
86  class FunctionNode;
87 
88  template<typename T, std::size_t NDIM>
89  class FunctionFactory;
90 
91  template<typename T, std::size_t NDIM>
92  class FunctionFunctorInterface;
93 
94  template<typename T, std::size_t NDIM>
95  struct leaf_op;
96 
97  template<typename T, std::size_t NDIM>
98  struct mul_leaf_op;
99 
100  template<typename T, std::size_t NDIM>
101  struct hartree_leaf_op;
102 
103  template<typename T, std::size_t NDIM, std::size_t LDIM, typename opT>
105 
106  template<typename T, std::size_t NDIM, typename opT>
107  struct op_leaf_op;
108 
109  template<typename T, std::size_t NDIM>
111 
112 }
113 
114 
115 namespace madness {
116 
117  /// \ingroup mra
118  /// \addtogroup function
119 
120  /// A multiresolution adaptive numerical function
121  template <typename T, std::size_t NDIM>
123  // We make all of the content of Function and FunctionImpl
124  // public with the intent of avoiding the cumbersome forward
125  // and friend declarations. However, this open access should
126  // not be abused.
127 
128  private:
129  std::shared_ptr< FunctionImpl<T,NDIM> > impl;
130 
131  public:
132  bool impl_initialized()const{
133  if(impl==NULL) return false;
134  else return true;
135  }
139  typedef Vector<double,NDIM> coordT; ///< Type of vector holding coordinates
140 
141  /// Asserts that the function is initialized
142  inline void verify() const {
144  }
145 
146  /// Returns true if the function is initialized
147  bool is_initialized() const {
148  return impl.get();
149  }
150 
151  /// Default constructor makes uninitialized function. No communication.
152 
153  /// An uninitialized function can only be assigned to. Any other operation will throw.
154  Function() : impl() {}
155 
156 
157  /// Constructor from FunctionFactory provides named parameter idiom. Possible non-blocking communication.
158  Function(const factoryT& factory)
159  : impl(new FunctionImpl<T,NDIM>(factory)) {
161  }
162 
163 
164  /// Copy constructor is \em shallow. No communication, works in either basis.
166  : impl(f.impl) {
167  }
168 
169 
170  /// Assignment is \em shallow. No communication, works in either basis.
173  if (this != &f) impl = f.impl;
174  return *this;
175  }
176 
177  /// Destruction of any underlying implementation is deferred to next global fence.
179 
180  /// Evaluates the function at a point in user coordinates. Possible non-blocking comm.
181 
182  /// Only the invoking process will receive the result via the future
183  /// though other processes may be involved in the evaluation.
184  ///
185  /// Throws if function is not initialized.
186  Future<T> eval(const coordT& xuser) const {
188  const double eps=1e-15;
189  verify();
191  coordT xsim;
192  user_to_sim(xuser,xsim);
193  // If on the boundary, move the point just inside the
194  // volume so that the evaluation logic does not fail
195  for (std::size_t d=0; d<NDIM; ++d) {
196  if (xsim[d] < -eps) {
197  MADNESS_EXCEPTION("eval: coordinate lower-bound error in dimension", d);
198  }
199  else if (xsim[d] < eps) {
200  xsim[d] = eps;
201  }
202 
203  if (xsim[d] > 1.0+eps) {
204  MADNESS_EXCEPTION("eval: coordinate upper-bound error in dimension", d);
205  }
206  else if (xsim[d] > 1.0-eps) {
207  xsim[d] = 1.0-eps;
208  }
209  }
210 
211  Future<T> result;
212  impl->eval(xsim, impl->key0(), result.remote_ref(impl->world));
213  return result;
214  }
215 
216  /// Evaluate function only if point is local returning (true,value); otherwise return (false,0.0)
217 
218  /// maxlevel is the maximum depth to search down to --- the max local depth can be
219  /// computed with max_local_depth();
220  std::pair<bool,T> eval_local_only(const Vector<double,NDIM>& xuser, Level maxlevel) const {
221  const double eps=1e-15;
222  verify();
224  coordT xsim;
225  user_to_sim(xuser,xsim);
226  // If on the boundary, move the point just inside the
227  // volume so that the evaluation logic does not fail
228  for (std::size_t d=0; d<NDIM; ++d) {
229  if (xsim[d] < -eps) {
230  MADNESS_EXCEPTION("eval: coordinate lower-bound error in dimension", d);
231  }
232  else if (xsim[d] < eps) {
233  xsim[d] = eps;
234  }
235 
236  if (xsim[d] > 1.0+eps) {
237  MADNESS_EXCEPTION("eval: coordinate upper-bound error in dimension", d);
238  }
239  else if (xsim[d] > 1.0-eps) {
240  xsim[d] = 1.0-eps;
241  }
242  }
243  return impl->eval_local_only(xsim,maxlevel);
244  }
245 
246  /// Only the invoking process will receive the result via the future
247  /// though other processes may be involved in the evaluation.
248  ///
249  /// Throws if function is not initialized.
250  ///
251  /// This function is a minimally-modified version of eval()
252  Future<Level> evaldepthpt(const coordT& xuser) const {
254  const double eps=1e-15;
255  verify();
257  coordT xsim;
258  user_to_sim(xuser,xsim);
259  // If on the boundary, move the point just inside the
260  // volume so that the evaluation logic does not fail
261  for (std::size_t d=0; d<NDIM; ++d) {
262  if (xsim[d] < -eps) {
263  MADNESS_EXCEPTION("eval: coordinate lower-bound error in dimension", d);
264  }
265  else if (xsim[d] < eps) {
266  xsim[d] = eps;
267  }
268 
269  if (xsim[d] > 1.0+eps) {
270  MADNESS_EXCEPTION("eval: coordinate upper-bound error in dimension", d);
271  }
272  else if (xsim[d] > 1.0-eps) {
273  xsim[d] = 1.0-eps;
274  }
275  }
276 
277  Future<Level> result;
278  impl->evaldepthpt(xsim, impl->key0(), result.remote_ref(impl->world));
279  return result;
280  }
281 
282 
283  /// Evaluates the function rank at a point in user coordinates. Possible non-blocking comm.
284 
285  /// Only the invoking process will receive the result via the future
286  /// though other processes may be involved in the evaluation.
287  ///
288  /// Throws if function is not initialized.
289  Future<long> evalR(const coordT& xuser) const {
291  const double eps=1e-15;
292  verify();
294  coordT xsim;
295  user_to_sim(xuser,xsim);
296  // If on the boundary, move the point just inside the
297  // volume so that the evaluation logic does not fail
298  for (std::size_t d=0; d<NDIM; ++d) {
299  if (xsim[d] < -eps) {
300  MADNESS_EXCEPTION("eval: coordinate lower-bound error in dimension", d);
301  }
302  else if (xsim[d] < eps) {
303  xsim[d] = eps;
304  }
305 
306  if (xsim[d] > 1.0+eps) {
307  MADNESS_EXCEPTION("eval: coordinate upper-bound error in dimension", d);
308  }
309  else if (xsim[d] > 1.0-eps) {
310  xsim[d] = 1.0-eps;
311  }
312  }
313 
314  Future<long> result;
315  impl->evalR(xsim, impl->key0(), result.remote_ref(impl->world));
316  return result;
317  }
318 
319  /// Evaluates a cube/slice of points (probably for plotting) ... collective but no fence necessary
320 
321  /// All processes recieve the entire result (which is a rather severe limit
322  /// on the size of the cube that is possible).
323 
324  /// Set eval_refine=true to return the refinment levels of
325  /// the given function.
326 
327  /// @param[in] cell A Tensor describe the cube where the function to be evaluated in
328  /// @param[in] npt How many points to evaluate in each dimension
329  /// @param[in] eval_refine Wether to return the refinment levels of the given function
331  const std::vector<long>& npt,
332  bool eval_refine = false) const {
333  MADNESS_ASSERT(static_cast<std::size_t>(cell.dim(0))>=NDIM && cell.dim(1)==2 && npt.size()>=NDIM);
335  const double eps=1e-14;
336  verify();
337  reconstruct();
338  coordT simlo, simhi;
339  for (std::size_t d=0; d<NDIM; ++d) {
340  simlo[d] = cell(d,0);
341  simhi[d] = cell(d,1);
342  }
343  user_to_sim(simlo, simlo);
344  user_to_sim(simhi, simhi);
345 
346  // Move the bounding box infintesimally inside dyadic
347  // points so that the evaluation logic does not fail
348  for (std::size_t d=0; d<NDIM; ++d) {
349  MADNESS_ASSERT(simhi[d] >= simlo[d]);
350  MADNESS_ASSERT(simlo[d] >= 0.0);
351  MADNESS_ASSERT(simhi[d] <= 1.0);
352 
353  double delta = eps*(simhi[d]-simlo[d]);
354  simlo[d] += delta;
355  simhi[d] -= 2*delta; // deliberate asymmetry
356  }
357  return impl->eval_plot_cube(simlo, simhi, npt, eval_refine);
358  }
359 
360 
361  /// Evaluates the function at a point in user coordinates. Collective operation.
362 
363  /// Throws if function is not initialized.
364  ///
365  /// This function calls eval, blocks until the result is
366  /// available and then broadcasts the result to everyone.
367  /// Therefore, if you are evaluating many points in parallel
368  /// it is \em vastly less efficient than calling eval
369  /// directly, saving the futures, and then forcing all of the
370  /// results.
371  T operator()(const coordT& xuser) const {
373  verify();
374  if (!is_reconstructed()) reconstruct();
375  T result;
376  if (impl->world.rank() == 0) result = eval(xuser).get();
377  impl->world.gop.broadcast(result);
378  //impl->world.gop.fence();
379  return result;
380  }
381 
382  /// Evaluates the function at a point in user coordinates. Collective operation.
383 
384  /// See "operator()(const coordT& xuser)" for more info
385  T operator()(double x, double y=0, double z=0, double xx=0, double yy=0, double zz=0) const {
386  coordT r;
387  r[0] = x;
388  if (NDIM>=2) r[1] = y;
389  if (NDIM>=3) r[2] = z;
390  if (NDIM>=4) r[3] = xx;
391  if (NDIM>=5) r[4] = yy;
392  if (NDIM>=6) r[5] = zz;
393  return (*this)(r);
394  }
395 
396  /// Throws if function is not initialized.
397  ///
398  /// This function mimics operator() by going through the
399  /// tree looking for the depth of the tree at the point.
400  /// It blocks until the result is
401  /// available and then broadcasts the result to everyone.
402  /// Therefore, if you are evaluating many points in parallel
403  /// it is \em vastly less efficient than calling evaldepthpt
404  /// directly, saving the futures, and then forcing all of the
405  /// results.
406  Level depthpt(const coordT& xuser) const {
408  verify();
409  if (!is_reconstructed()) reconstruct();
410  Level result;
411  if (impl->world.rank() == 0) result = evaldepthpt(xuser).get();
412  impl->world.gop.broadcast(result);
413  //impl->world.gop.fence();
414  return result;
415  }
416 
417  /// Returns an estimate of the difference ||this-func||^2 from local data
418 
419  /// No communication is performed. If the function is not
420  /// reconstructed, it throws an exception. To get the global
421  /// value either do a global sum of the local values or call
422  /// errsq
423  /// @param[in] func Templated interface to the a user specified function
424  template <typename funcT>
425  double errsq_local(const funcT& func) const {
427  verify();
428  if (!is_reconstructed()) MADNESS_EXCEPTION("Function:errsq_local:not reconstructed",0);
429  return impl->errsq_local(func);
430  }
431 
432 
433  /// Returns an estimate of the difference ||this-func|| ... global sum performed
434 
435  /// If the function is compressed, it is reconstructed first. For efficient use
436  /// especially with many functions, reconstruct them all first, and use errsq_local
437  /// instead so you can perform a global sum on all at the same time.
438  /// @param[in] func Templated interface to the a user specified function
439  template <typename funcT>
440  double err(const funcT& func) const {
442  verify();
443  if (VERIFY_TREE) verify_tree();
444  if (!is_reconstructed()) reconstruct();
445  if (VERIFY_TREE) verify_tree();
446  double local = impl->errsq_local(func);
447  impl->world.gop.sum(local);
448  impl->world.gop.fence();
449  return sqrt(local);
450  }
451 
452  /// Verifies the tree data structure ... global sync implied
453  void verify_tree() const {
455  if (impl) impl->verify_tree();
456  }
457 
458 
459  /// Returns true if compressed, false otherwise. No communication.
460 
461  /// If the function is not initialized, returns false.
462  bool is_compressed() const {
464  if (impl)
465  return impl->is_compressed();
466  else
467  return false;
468  }
469 
470  /// Returns true if reconstructed, false otherwise. No communication.
471 
472  /// If the function is not initialized, returns false.
473  bool is_reconstructed() const {
475  if (impl)
476  return impl->is_reconstructed();
477  else
478  return false;
479  }
480 
481  /// Returns true if nonstandard-compressed, false otherwise. No communication.
482 
483  /// If the function is not initialized, returns false.
484  bool is_nonstandard() const {
486  return impl ? impl->is_nonstandard() : false;
487  }
488 
489 
490  /// Returns the number of nodes in the function tree ... collective global sum
491  std::size_t tree_size() const {
493  if (!impl) return 0;
494  return impl->tree_size();
495  }
496 
497  /// print some info about this
498  void print_size(const std::string name) const {
499  if (!impl) {
500  print("function",name,"not assigned yet");
501  } else {
502  impl->print_size(name);
503  }
504  }
505 
506  /// Returns the maximum depth of the function tree ... collective global sum
507  std::size_t max_depth() const {
509  if (!impl) return 0;
510  return impl->max_depth();
511  }
512 
513 
514  /// Returns the maximum local depth of the function tree ... no communications
515  std::size_t max_local_depth() const {
517  if (!impl) return 0;
518  return impl->max_local_depth();
519  }
520 
521 
522  /// Returns the max number of nodes on a processor
523  std::size_t max_nodes() const {
525  if (!impl) return 0;
526  return impl->max_nodes();
527  }
528 
529  /// Returns the min number of nodes on a processor
530  std::size_t min_nodes() const {
532  if (!impl) return 0;
533  return impl->min_nodes();
534  }
535 
536 
537  /// Returns the number of coefficients in the function ... collective global sum
538  std::size_t size() const {
540  if (!impl) return 0;
541  return impl->size();
542  }
543 
544  /// Retunrs
545 
546 
547  /// Returns value of autorefine flag. No communication.
548  bool autorefine() const {
550  if (!impl) return true;
551  return impl->get_autorefine();
552  }
553 
554 
555  /// Sets the value of the autorefine flag. Optional global fence.
556 
557  /// A fence is required to ensure consistent global state.
558  void set_autorefine(bool value, bool fence = true) {
560  verify();
561  impl->set_autorefine(value);
562  if (fence) impl->world.gop.fence();
563  }
564 
565 
566  /// Returns value of truncation threshold. No communication.
567  double thresh() const {
569  if (!impl) return 0.0;
570  return impl->get_thresh();
571  }
572 
573 
574  /// Sets the value of the truncation threshold. Optional global fence.
575 
576  /// A fence is required to ensure consistent global state.
577  void set_thresh(double value, bool fence = true) {
579  verify();
580  impl->set_thresh(value);
581  if (fence) impl->world.gop.fence();
582  }
583 
584 
585  /// Returns the number of multiwavelets (k). No communication.
586  int k() const {
588  verify();
589  return impl->get_k();
590  }
591 
592 
593  /// Truncate the function with optional fence. Compresses with fence if not compressed.
594 
595  /// If the truncation threshold is less than or equal to zero the default value
596  /// specified when the function was created is used.
597  /// If the function is not initialized, it just returns.
598  ///
599  /// Returns this for chaining.
600  /// @param[in] tol Tolerance for truncating the coefficients. Default 0.0 means use the implementation's member value \c thresh instead.
601  /// @param[in] fence Do fence
602  Function<T,NDIM>& truncate(double tol = 0.0, bool fence = true) {
604  if (!impl) return *this;
605  verify();
606 // if (!is_compressed()) compress();
607  impl->truncate(tol,fence);
608  if (VERIFY_TREE) verify_tree();
609  return *this;
610  }
611 
612 
613  /// Returns a shared-pointer to the implementation
614  const std::shared_ptr< FunctionImpl<T,NDIM> >& get_impl() const {
616  verify();
617  return impl;
618  }
619 
620  /// Replace current FunctionImpl with provided new one
621  void set_impl(const std::shared_ptr< FunctionImpl<T,NDIM> >& impl) {
623  this->impl = impl;
624  }
625 
626 
627  /// Replace the current functor with the provided new one
628 
629  /// presumably the new functor will be a CompositeFunctor, which will
630  /// change the behavior of the function: multiply the functor with the function
631  void set_functor(const std::shared_ptr<FunctionFunctorInterface<T, NDIM> > functor) {
632  this->impl->set_functor(functor);
633  print("set functor in mra.h");
634  }
635 
636  bool is_on_demand() const {return this->impl->is_on_demand();}
637 
638  /// Replace current FunctionImpl with a new one using the same parameters & map as f
639 
640  /// If zero is true the function is initialized to zero, otherwise it is empty
641  template <typename R>
642  void set_impl(const Function<R,NDIM>& f, bool zero = true) {
643  impl = std::shared_ptr<implT>(new implT(*f.get_impl(), f.get_pmap(), zero));
644  if (zero) world().gop.fence();
645  }
646 
647  /// Returns the world
648  World& world() const {
650  verify();
651  return impl->world;
652  }
653 
654 
655  /// Returns a shared pointer to the process map
656  const std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& get_pmap() const {
658  verify();
659  return impl->get_pmap();
660  }
661 
662  /// replicate this function, generating a unique pmap
663  void replicate(bool fence=true) const {
664  verify();
665  impl->replicate(fence);
666  }
667 
668  /// distribute this function according to newmap
669  void distribute(std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > > newmap) const {
670  verify();
671  impl->distribute(newmap);
672  }
673 
674 
675  /// Returns the square of the norm of the local function ... no communication
676 
677  /// Works in either basis
678  double norm2sq_local() const {
680  verify();
681  return impl->norm2sq_local();
682  }
683 
684 
685  /// Returns the 2-norm of the function ... global sum ... works in either basis
686 
687  /// See comments for err() w.r.t. applying to many functions.
688  double norm2() const {
690  verify();
691  if (VERIFY_TREE) verify_tree();
692  double local = impl->norm2sq_local();
693 
694  impl->world.gop.sum(local);
695  impl->world.gop.fence();
696  return sqrt(local);
697  }
698 
699 
700  /// Initializes information about the function norm at all length scales
701  void norm_tree(bool fence = true) const {
703  verify();
704  if (VERIFY_TREE) verify_tree();
705  if (!is_reconstructed()) reconstruct();
706  const_cast<Function<T,NDIM>*>(this)->impl->norm_tree(fence);
707  }
708 
709 
710  /// Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
711 
712  /// By default fence=true meaning that this operation completes before returning,
713  /// otherwise if fence=false it returns without fencing and the user must invoke
714  /// world.gop.fence() to assure global completion before using the function
715  /// for other purposes.
716  ///
717  /// Noop if already compressed or if not initialized.
718  ///
719  /// Since reconstruction/compression do not discard information we define them
720  /// as const ... "logical constness" not "bitwise constness".
721  const Function<T,NDIM>& compress(bool fence = true) const {
723  }
724 
725 
726  /// Compresses the function retaining scaling function coeffs. Possible non-blocking comm.
727 
728  /// By default fence=true meaning that this operation completes before returning,
729  /// otherwise if fence=false it returns without fencing and the user must invoke
730  /// world.gop.fence() to assure global completion before using the function
731  /// for other purposes.
732  ///
733  /// Noop if already compressed or if not initialized.
734  void make_nonstandard(bool keepleaves, bool fence=true) const {
736  if (keepleaves) newstate=nonstandard_with_leaves;
737  change_tree_state(newstate,fence);
738  }
739 
740  /// Converts the function standard compressed form. Possible non-blocking comm.
741 
742  /// By default fence=true meaning that this operation completes before returning,
743  /// otherwise if fence=false it returns without fencing and the user must invoke
744  /// world.gop.fence() to assure global completion before using the function
745  /// for other purposes.
746  ///
747  /// Must be already compressed.
748  void standard(bool fence = true) {
750  }
751 
752  /// Converts the function to redundant form, i.e. sum coefficients on all levels
753 
754  /// By default fence=true meaning that this operation completes before returning,
755  /// otherwise if fence=false it returns without fencing and the user must invoke
756  /// world.gop.fence() to assure global completion before using the function
757  /// for other purposes.
758  ///
759  /// Must be already compressed.
760  void make_redundant(bool fence = true) {
762  }
763 
764  /// Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
765 
766  /// By default fence=true meaning that this operation completes before returning,
767  /// otherwise if fence=false it returns without fencing and the user must invoke
768  /// world.gop.fence() to assure global completion before using the function
769  /// for other purposes.
770  ///
771  /// Noop if already reconstructed or if not initialized.
772  ///
773  /// Since reconstruction/compression do not discard information we define them
774  /// as const ... "logical constness" not "bitwise constness".
775  const Function<T,NDIM>& reconstruct(bool fence = true) const {
777  }
778 
779  /// changes tree state to given state
780 
781  /// Since reconstruction/compression do not discard information we define them
782  /// as const ... "logical constness" not "bitwise constness".
783  /// @param[in] finalstate The final state of the tree
784  /// @param[in] fence Fence after the operation (might not be respected!!!)
785  const Function<T,NDIM>& change_tree_state(const TreeState finalstate, bool fence = true) const {
787  if (not impl) return *this;
788  TreeState current_state = impl->get_tree_state();
789  if (finalstate == current_state) return *this;
790  MADNESS_CHECK_THROW(current_state != TreeState::unknown, "unknown tree state");
791 
792  impl->change_tree_state(finalstate, fence);
793  if (fence && VERIFY_TREE) verify_tree();
794  return *this;
795  }
796 
797  /// Sums scaling coeffs down tree restoring state with coeffs only at leaves. Optional fence. Possible non-blocking comm.
798  void sum_down(bool fence = true) const {
800  verify();
801  MADNESS_CHECK_THROW(impl->get_tree_state()==redundant_after_merge, "sum_down requires a redundant_after_merge state");
802  const_cast<Function<T,NDIM>*>(this)->impl->sum_down(fence);
803  const_cast<Function<T,NDIM>*>(this)->impl->set_tree_state(reconstructed);
804 
805  if (fence && VERIFY_TREE) verify_tree(); // Must be after in case nonstandard
806  }
807 
808 
809  /// Inplace autorefines the function. Optional fence. Possible non-blocking comm.
810  template <typename opT>
811  void refine_general(const opT& op, bool fence = true) const {
813  verify();
814  if (!is_reconstructed()) reconstruct();
815  impl->refine(op, fence);
816  }
817 
818 
820  bool operator()(implT* impl, const Key<NDIM>& key, const nodeT& t) const {
821  return impl->autorefine_square_test(key, t);
822  }
823 
824  template <typename Archive> void serialize (Archive& ar) {}
825  };
826 
827  /// Inplace autorefines the function using same test as for squaring.
828 
829  /// return this for chaining
830  const Function<T,NDIM>& refine(bool fence = true) const {
832  return *this;
833  }
834 
835  /// Inplace broadens support in scaling function basis
837  bool fence = true) const {
838  verify();
839  reconstruct();
840  impl->broaden(bc.is_periodic(), fence);
841  }
842 
843 
844  /// Clears the function as if constructed uninitialized. Optional fence.
845 
846  /// Any underlying data will not be freed until the next global fence.
847  void clear(bool fence = true) {
849  if (impl) {
850  World& world = impl->world;
851  impl.reset();
852  if (fence) world.gop.fence();
853  }
854  }
855 
856  /// Process 0 prints a summary of all nodes in the tree (collective)
857  void print_tree(std::ostream& os = std::cout) const {
859  if (impl) impl->print_tree(os);
860  }
861 
862  /// same as print_tree() but produces JSON-formatted string
863  /// @warning enclose the result in braces to make it a valid JSON object
864  void print_tree_json(std::ostream& os = std::cout) const {
866  if (impl) impl->print_tree_json(os);
867  }
868 
869  /// Process 0 prints a graphviz-formatted output of all nodes in the tree (collective)
870  void print_tree_graphviz(std::ostream& os = std::cout) const {
872  os << "digraph G {" << std::endl;
873  if (impl) impl->print_tree_graphviz(os);
874  os << "}" << std::endl;
875  }
876 
877  /// Print a summary of the load balancing info
878 
879  /// This is serial and VERY expensive
880  void print_info() const {
882  if (impl) impl->print_info();
883  }
884 
886  T (*f)(T);
888  void operator()(const Key<NDIM>& key, Tensor<T>& t) const {
889  UNARY_OPTIMIZED_ITERATOR(T, t, *_p0 = f(*_p0));
890  }
891  template <typename Archive> void serialize(Archive& ar) {}
892  };
893 
894  /// Inplace unary operation on function values
895  void unaryop(T (*f)(T)) {
896  // Must fence here due to temporary object on stack
897  // stopping us returning before complete
899  }
900 
901 
902  /// Inplace unary operation on function values
903  template <typename opT>
904  void unaryop(const opT& op, bool fence=true) {
906  verify();
907  reconstruct();
908  impl->unary_op_value_inplace(op, fence);
909  }
910 
911 
912  /// Unary operation applied inplace to the coefficients
913  template <typename opT>
914  void unaryop_coeff(const opT& op,
915  bool fence = true) {
917  verify();
918  impl->unary_op_coeff_inplace(op, fence);
919  }
920 
921 
922  /// Unary operation applied inplace to the nodes
923  template <typename opT>
924  void unaryop_node(const opT& op,
925  bool fence = true) {
927  verify();
928  impl->unary_op_node_inplace(op, fence);
929  }
930 
931 
932 
933 
934  static void doconj(const Key<NDIM>, Tensor<T>& t) {
936  t.conj();
937  }
938 
939  /// Inplace complex conjugate. No communication except for optional fence.
940 
941  /// Returns this for chaining. Works in either basis.
942  Function<T,NDIM> conj(bool fence = true) {
945  return *this;
946  }
947 
948 
949  /// Inplace, scale the function by a constant. No communication except for optional fence.
950 
951  /// Works in either basis. Returns reference to this for chaining.
952  template <typename Q>
953  Function<T,NDIM>& scale(const Q q, bool fence=true) {
955  verify();
956  if (VERIFY_TREE) verify_tree();
957  impl->scale_inplace(q,fence);
958  return *this;
959  }
960 
961 
962  /// Inplace add scalar. No communication except for optional fence.
965  verify();
966  if (VERIFY_TREE) verify_tree();
967  impl->add_scalar_inplace(t,fence);
968  return *this;
969  }
970 
971 
972  /// Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
973 
974  /// If the functions are not in the wavelet basis an exception is thrown since this routine
975  /// is intended to be fast and unexpected compression is assumed to be a performance bug.
976  ///
977  /// Returns this for chaining.
978  ///
979  /// this <-- this*alpha + other*beta
980  template <typename Q, typename R>
982  const Function<Q,NDIM>& other, const R& beta, bool fence=true) {
984  verify();
985  other.verify();
986  MADNESS_CHECK_THROW(impl->get_tree_state() == other.get_impl()->get_tree_state(),
987  "gaxpy requires both functions to be in the same tree state");
988  bool same_world=this->world().id()==other.world().id();
989  MADNESS_CHECK(same_world or is_compressed());
990 
991  if (not same_world) {
992  impl->gaxpy_inplace(alpha,*other.get_impl(),beta,fence);
993  } else {
994  if (is_compressed()) impl->gaxpy_inplace(alpha, *other.get_impl(), beta, fence);
995  if (is_reconstructed()) impl->gaxpy_inplace_reconstructed(alpha,*other.get_impl(),beta,fence);
996  }
997  return *this;
998  }
999 
1000 
1001  /// Inplace addition of functions in the wavelet basis
1002 
1003  /// Using operator notation forces a global fence after every operation.
1004  /// Functions don't need to be compressed, it's the caller's responsibility
1005  /// to choose an appropriate state with performance, usually compressed for 3d,
1006  /// reconstructed for 6d)
1007  template <typename Q>
1010  if (NDIM<=3) {
1011  compress();
1012  other.compress();
1013  } else {
1014  reconstruct();
1015  other.reconstruct();
1016  }
1017  MADNESS_ASSERT(impl->get_tree_state() == other.get_impl()->get_tree_state());
1018  if (VERIFY_TREE) verify_tree();
1019  if (VERIFY_TREE) other.verify_tree();
1020  return gaxpy(T(1.0), other, Q(1.0), true);
1021  }
1022 
1023 
1024  /// Inplace subtraction of functions in the wavelet basis
1025 
1026  /// Using operator notation forces a global fence after every operation
1027  template <typename Q>
1030  if (NDIM<=3) {
1031  compress();
1032  other.compress();
1033  } else {
1034  reconstruct();
1035  other.reconstruct();
1036  }
1037  MADNESS_ASSERT(impl->get_tree_state() == other.get_impl()->get_tree_state());
1038  if (VERIFY_TREE) verify_tree();
1039  if (VERIFY_TREE) other.verify_tree();
1040  return gaxpy(T(1.0), other, Q(-1.0), true);
1041  }
1042 
1043 
1044  /// Inplace scaling by a constant
1045 
1046  /// Using operator notation forces a global fence after every operation
1047  template <typename Q>
1049  operator*=(const Q q) {
1051  scale(q,true);
1052  return *this;
1053  }
1054 
1055 
1056  /// Inplace squaring of function ... global comm only if not reconstructed
1057 
1058  /// Returns *this for chaining.
1061  if (!is_reconstructed()) reconstruct();
1062  if (VERIFY_TREE) verify_tree();
1063  impl->square_inplace(fence);
1064  return *this;
1065  }
1066 
1067  /// Returns *this for chaining.
1068  Function<T,NDIM>& abs(bool fence = true) {
1070  if (!is_reconstructed()) reconstruct();
1071  if (VERIFY_TREE) verify_tree();
1072  impl->abs_inplace(fence);
1073  return *this;
1074  }
1075 
1076  /// Returns *this for chaining.
1079  if (!is_reconstructed()) reconstruct();
1080  if (VERIFY_TREE) verify_tree();
1081  impl->abs_square_inplace(fence);
1082  return *this;
1083  }
1084 
1085  /// Returns local contribution to \c int(f(x),x) ... no communication
1086 
1087  /// In the wavelet basis this is just the coefficient of the first scaling
1088  /// function which is a constant. In the scaling function basis we
1089  /// must add up contributions from each box.
1090  T trace_local() const {
1092  if (!impl) return 0.0;
1093  if (VERIFY_TREE) verify_tree();
1094  return impl->trace_local();
1095  }
1096 
1097 
1098  /// Returns global value of \c int(f(x),x) ... global comm required
1099  T trace() const {
1101  if (!impl) return 0.0;
1102  T sum = impl->trace_local();
1103  impl->world.gop.sum(sum);
1104  impl->world.gop.fence();
1105  return sum;
1106  }
1107 
1108 
1109  /// Returns local part of inner product ... throws if both not compressed
1110  template <typename R>
1111  TENSOR_RESULT_TYPE(T,R) inner_local(const Function<R,NDIM>& g) const {
1114  MADNESS_ASSERT(g.is_compressed());
1116  if (VERIFY_TREE) g.verify_tree();
1117  return impl->inner_local(*(g.get_impl()));
1118  }
1119 
1120 
1121  /// With this being an on-demand function, fill the MRA tree according to different criteria
1122 
1123  /// @param[in] g the function after which the MRA structure is modeled (any basis works)
1124  template<typename R>
1126  MADNESS_ASSERT(g.is_initialized());
1128 
1129  // clear what we have
1130  impl->get_coeffs().clear();
1131 
1132  //leaf_op<T,NDIM> gnode_is_leaf(g.get_impl().get());
1133  Leaf_op_other<T,NDIM> gnode_is_leaf(g.get_impl().get());
1134  impl->make_Vphi(gnode_is_leaf,fence);
1135  return *this;
1136 
1137  }
1138 
1139  /// With this being an on-demand function, fill the MRA tree according to different criteria
1140 
1141  /// @param[in] op the convolution operator for screening
1142  template<typename opT>
1143  Function<T,NDIM>& fill_tree(const opT& op, bool fence=true) {
1145  // clear what we have
1146  impl->get_coeffs().clear();
1147  Specialbox_op<T,NDIM> sbox;
1149  impl ->make_Vphi(leaf_op,fence);
1150  return *this;
1151  }
1152 
1153  /// With this being an on-demand function, fill the MRA tree according to different criteria
1156  // clear what we have
1157  impl->get_coeffs().clear();
1159  impl->make_Vphi(leaf_op,fence);
1160  return *this;
1161  }
1162 
1163  /// Special refinement on 6D boxes where the electrons come close (meet)
1164  /// @param[in] op the convolution operator for screening
1165  template<typename opT>
1166  Function<T,NDIM>& fill_cuspy_tree(const opT& op,const bool fence=true){
1168  // clear what we have
1169  impl->get_coeffs().clear();
1171 
1173  impl ->make_Vphi(leaf_op,fence);
1174 
1175  return *this;
1176  }
1177 
1178  /// Special refinement on 6D boxes where the electrons come close (meet)
1181  // clear what we have
1182  impl->get_coeffs().clear();
1184 
1186  impl ->make_Vphi(leaf_op,fence);
1187 
1188  return *this;
1189  }
1190 
1191  /// Special refinement on 6D boxes for the nuclear potentials (regularized with cusp, non-regularized with singularity)
1192  /// @param[in] op the convolution operator for screening
1193  template<typename opT>
1194  Function<T,NDIM>& fill_nuclear_cuspy_tree(const opT& op,const size_t particle,const bool fence=true){
1196  // clear what we have
1197  impl->get_coeffs().clear();
1199 
1201  impl ->make_Vphi(leaf_op,fence);
1202 
1203  return *this;
1204  }
1205 
1206  /// Special refinement on 6D boxes for the nuclear potentials (regularized with cusp, non-regularized with singularity)
1207  Function<T,NDIM>& fill_nuclear_cuspy_tree(const size_t particle,const bool fence=true){
1209  // clear what we have
1210  impl->get_coeffs().clear();
1212 
1214  impl ->make_Vphi(leaf_op,fence);
1215 
1216  return *this;
1217  }
1218 
1219  /// perform the hartree product of f*g, invoked by result
1220  template<size_t LDIM, size_t KDIM, typename opT>
1221  void do_hartree_product(const std::vector<std::shared_ptr<FunctionImpl<T,LDIM>>> left,
1222  const std::vector<std::shared_ptr<FunctionImpl<T,KDIM>>> right,
1223  const opT* op) {
1224 
1225  // get the right leaf operator
1227  impl->hartree_product(left,right,leaf_op,true);
1228  impl->finalize_sum();
1229 // this->truncate();
1230 
1231  }
1232 
1233  /// perform the hartree product of f*g, invoked by result
1234  template<size_t LDIM, size_t KDIM>
1235  void do_hartree_product(const std::vector<std::shared_ptr<FunctionImpl<T,LDIM>>> left,
1236  const std::vector<std::shared_ptr<FunctionImpl<T,KDIM>>> right) {
1237 
1238 // hartree_leaf_op<T,KDIM+LDIM> leaf_op(impl.get(),cdata.s0);
1240  impl->hartree_product(left,right,leaf_op,true);
1241  impl->finalize_sum();
1242 // this->truncate();
1243 
1244  }
1245 
1246  /// Returns the inner product
1247 
1248  /// Not efficient for computing multiple inner products
1249  /// @param[in] g Function, optionally on-demand
1250  template <typename R>
1253 
1254  // fast return if possible
1255  if (not this->is_initialized()) return 0.0;
1256  if (not g.is_initialized()) return 0.0;
1257 
1258  // if this and g are the same, use norm2()
1259  if (this->get_impl()==g.get_impl()) {
1260  TreeState state=this->get_impl()->get_tree_state();
1262  double norm=this->norm2();
1263  return norm*norm;
1264  }
1265 
1266  // do it case-by-case
1267  if (this->is_on_demand()) return g.inner_on_demand(*this);
1268  if (g.is_on_demand()) return this->inner_on_demand(g);
1269 
1271  if (VERIFY_TREE) g.verify_tree();
1272 
1273  // compression is more efficient for 3D
1274  TreeState state=this->get_impl()->get_tree_state();
1275  TreeState gstate=g.get_impl()->get_tree_state();
1276  if (NDIM<=3) {
1278  g.change_tree_state(compressed,false);
1279  impl->world.gop.fence();
1280  }
1281 
1282  if (this->is_compressed() and g.is_compressed()) {
1283  } else {
1285  g.change_tree_state(redundant,false);
1286  impl->world.gop.fence();
1287  }
1288 
1289 
1290  TENSOR_RESULT_TYPE(T,R) local = impl->inner_local(*g.get_impl());
1291  impl->world.gop.sum(local);
1292  impl->world.gop.fence();
1293 
1294  // restore state
1296  g.change_tree_state(gstate,false);
1297  impl->world.gop.fence();
1298 
1299  return local;
1300  }
1301 
1302  /// Return the local part of inner product with external function ... no communication.
1303  /// If you are going to be doing a bunch of inner_ext calls, set
1304  /// keep_redundant to true and then manually undo_redundant when you
1305  /// are finished.
1306  /// @param[in] f Pointer to function of type T that take coordT arguments. This is the externally provided function
1307  /// @param[in] leaf_refine boolean switch to turn on/off refinement past leaf nodes
1308  /// @param[in] keep_redundant boolean switch to turn on/off undo_redundant
1309  /// @return Returns local part of the inner product, i.e. over the domain of all function nodes on this compute node.
1310  T inner_ext_local(const std::shared_ptr< FunctionFunctorInterface<T,NDIM> > f, const bool leaf_refine=true, const bool keep_redundant=false) const {
1313  T local = impl->inner_ext_local(f, leaf_refine);
1314  if (not keep_redundant) change_tree_state(reconstructed);
1315  return local;
1316  }
1317 
1318  /// Return the inner product with external function ... requires communication.
1319  /// If you are going to be doing a bunch of inner_ext calls, set
1320  /// keep_redundant to true and then manually undo_redundant when you
1321  /// are finished.
1322  /// @param[in] f Reference to FunctionFunctorInterface. This is the externally provided function
1323  /// @param[in] leaf_refine boolean switch to turn on/off refinement past leaf nodes
1324  /// @param[in] keep_redundant boolean switch to turn on/off undo_redundant
1325  /// @return Returns the inner product
1326  T inner_ext(const std::shared_ptr< FunctionFunctorInterface<T,NDIM> > f, const bool leaf_refine=true, const bool keep_redundant=false) const {
1329  T local = impl->inner_ext_local(f, leaf_refine);
1330  impl->world.gop.sum(local);
1331  impl->world.gop.fence();
1332  if (not keep_redundant) change_tree_state(reconstructed);
1333  return local;
1334  }
1335 
1336  /// Return the inner product with external function ... requires communication.
1337  /// If you are going to be doing a bunch of inner_ext calls, set
1338  /// keep_redundant to true and then manually undo_redundant when you
1339  /// are finished.
1340  /// @param[in] f Reference to FunctionFunctorInterface. This is the externally provided function
1341  /// @param[in] leaf_refine boolean switch to turn on/off refinement past leaf nodes
1342  /// @return Returns the inner product
1344  const bool leaf_refine=true) const {
1346  reconstruct();
1347  T local = impl->inner_adaptive_local(f, leaf_refine);
1348  impl->world.gop.sum(local);
1349  impl->world.gop.fence();
1350  return local;
1351  }
1352 
1353  /// Return the local part of gaxpy with external function, this*alpha + f*beta ... no communication.
1354  /// @param[in] alpha prefactor for this Function
1355  /// @param[in] f Pointer to function of type T that take coordT arguments. This is the externally provided function
1356  /// @param[in] beta prefactor for f
1357  template <typename L>
1358  void gaxpy_ext(const Function<L,NDIM>& left, T (*f)(const coordT&), T alpha, T beta, double tol, bool fence=true) const {
1360  if (!left.is_reconstructed()) left.reconstruct();
1361  impl->gaxpy_ext(left.get_impl().get(), f, alpha, beta, tol, fence);
1362  }
1363 
1364  /// Returns the inner product for one on-demand function
1365 
1366  /// It does work, but it might not give you the precision you expect.
1367  /// The assumption is that the function g returns proper sum
1368  /// coefficients on the MRA tree of this. This might not be the case if
1369  /// g is constructed with an implicit multiplication, e.g.
1370  /// result = <this|g>, with g = 1/r12 | gg>
1371  /// @param[in] g on-demand function
1372  template<typename R>
1373  TENSOR_RESULT_TYPE(T, R) inner_on_demand(const Function<R, NDIM>& g) const {
1374  MADNESS_ASSERT(g.is_on_demand() and (not this->is_on_demand()));
1375 
1376  constexpr std::size_t LDIM=std::max(NDIM/2,std::size_t(1));
1377  auto func=dynamic_cast<CompositeFunctorInterface<T,NDIM,LDIM>* >(g.get_impl()->get_functor().get());
1379  func->make_redundant(true);
1380  func->replicate_low_dim_functions(true);
1381  this->reconstruct(); // if this == &g we don't need g to be redundant
1382 
1384 
1385  TENSOR_RESULT_TYPE(T, R) local = impl->inner_local_on_demand(*g.get_impl());
1386  impl->world.gop.sum(local);
1387  impl->world.gop.fence();
1388 
1389  return local;
1390  }
1391 
1392  /// project this on the low-dim function g: h(x) = <f(x,y) | g(y)>
1393 
1394  /// @param[in] g low-dim function
1395  /// @param[in] dim over which dimensions to be integrated: 0..LDIM-1 or LDIM..NDIM-1
1396  /// @return new function of dimension NDIM-LDIM
1397  template <typename R, size_t LDIM>
1399  if (NDIM<=LDIM) MADNESS_EXCEPTION("confused dimensions in project_out?",1);
1400  MADNESS_CHECK_THROW(dim==0 or dim==1,"dim must be 0 or 1 in project_out");
1401  verify();
1402  typedef TENSOR_RESULT_TYPE(T,R) resultT;
1403  static const size_t KDIM=NDIM-LDIM;
1404 
1406  .k(g.k()).thresh(g.thresh());
1407  Function<resultT,KDIM> result=factory; // no empty() here!
1408 
1410  g.change_tree_state(redundant,false);
1411  world().gop.fence();
1412  this->get_impl()->project_out(result.get_impl().get(),g.get_impl().get(),dim,true);
1413 // result.get_impl()->project_out2(this->get_impl().get(),gimpl,dim);
1414  result.world().gop.fence();
1415  g.change_tree_state(reconstructed,false);
1416  result.get_impl()->trickle_down(false);
1417  result.get_impl()->set_tree_state(reconstructed);
1418  result.world().gop.fence();
1419  return result;
1420  }
1421 
1422  Function<T,NDIM/2> dirac_convolution(const bool fence=true) const {
1423  constexpr std::size_t LDIM=NDIM/2;
1424  MADNESS_CHECK_THROW(NDIM==2*LDIM,"NDIM must be even");
1425 // // this will be the result function
1427  Function<T,LDIM> f = factory;
1428  if(!is_reconstructed()) this->reconstruct();
1429  this->get_impl()->do_dirac_convolution(f.get_impl().get(),fence);
1430  return f;
1431  }
1432 
1433  /// Replaces this function with one loaded from an archive using the default processor map
1434 
1435  /// Archive can be sequential or parallel.
1436  ///
1437  /// The & operator for serializing will only work with parallel archives.
1438  template <typename Archive>
1439  void load(World& world, Archive& ar) {
1441  // Type checking since we are probably circumventing the archive's own type checking
1442  long magic = 0l, id = 0l, ndim = 0l, k = 0l;
1443  ar & magic & id & ndim & k;
1444  MADNESS_ASSERT(magic == 7776768); // Mellow Mushroom Pizza tel.# in Knoxville
1446  MADNESS_ASSERT(ndim == NDIM);
1447 
1448  impl.reset(new implT(FunctionFactory<T,NDIM>(world).k(k).empty()));
1449 
1450  impl->load(ar);
1451  }
1452 
1453 
1454  /// Stores the function to an archive
1455 
1456  /// Archive can be sequential or parallel.
1457  ///
1458  /// The & operator for serializing will only work with parallel archives.
1459  template <typename Archive>
1460  void store(Archive& ar) const {
1462  verify();
1463  // For type checking, etc.
1464  ar & long(7776768) & long(TensorTypeData<T>::id) & long(NDIM) & long(k());
1465 
1466  impl->store(ar);
1467  }
1468 
1469  /// change the tensor type of the coefficients in the FunctionNode
1470 
1471  /// @param[in] targs target tensor arguments (threshold and full/low rank)
1472  void change_tensor_type(const TensorArgs& targs, bool fence=true) {
1473  if (not impl) return;
1474  impl->change_tensor_type1(targs,fence);
1475  }
1476 
1477 
1478  /// This is replaced with left*right ... private
1479  template <typename Q, typename opT>
1481  const opT& op, bool fence) {
1483  func.verify();
1484  MADNESS_ASSERT(func.is_reconstructed());
1485  if (VERIFY_TREE) func.verify_tree();
1486  impl.reset(new implT(*func.get_impl(), func.get_pmap(), false));
1487  impl->unaryXX(func.get_impl().get(), op, fence);
1488  return *this;
1489  }
1490 
1491  /// Returns vector of FunctionImpl pointers corresponding to vector of functions
1492  template <typename Q, std::size_t D>
1493  static std::vector< std::shared_ptr< FunctionImpl<Q,D> > > vimpl(const std::vector< Function<Q,D> >& v) {
1495  std::vector< std::shared_ptr< FunctionImpl<Q,D> > > r(v.size());
1496  for (unsigned int i=0; i<v.size(); ++i) r[i] = v[i].get_impl();
1497  return r;
1498  }
1499 
1500  /// This is replaced with op(vector of functions) ... private
1501  template <typename opT>
1502  Function<T,NDIM>& multiop_values(const opT& op, const std::vector< Function<T,NDIM> >& vf) {
1503  std::vector<implT*> v(vf.size(),NULL);
1504  for (unsigned int i=0; i<v.size(); ++i) {
1505  if (vf[i].is_initialized()) v[i] = vf[i].get_impl().get();
1506  }
1507  impl->multiop_values(op, v);
1508  world().gop.fence();
1509  if (VERIFY_TREE) verify_tree();
1510 
1511  return *this;
1512  }
1513 
1514  /// apply op on the input vector yielding an output vector of functions
1515 
1516  /// (*this) is just a dummy Function to be able to call internal methods in FuncImpl
1517  /// @param[in] op the operator working on vin
1518  /// @param[in] vin vector of input Functions
1519  /// @param[out] vout vector of output Functions vout = op(vin)
1520  template <typename opT>
1521  void multi_to_multi_op_values(const opT& op,
1522  const std::vector< Function<T,NDIM> >& vin,
1523  std::vector< Function<T,NDIM> >& vout,
1524  const bool fence=true) {
1525  std::vector<implT*> vimplin(vin.size(),NULL);
1526  for (unsigned int i=0; i<vin.size(); ++i) {
1527  if (vin[i].is_initialized()) vimplin[i] = vin[i].get_impl().get();
1528  }
1529  std::vector<implT*> vimplout(vout.size(),NULL);
1530  for (unsigned int i=0; i<vout.size(); ++i) {
1531  if (vout[i].is_initialized()) vimplout[i] = vout[i].get_impl().get();
1532  }
1533 
1534  impl->multi_to_multi_op_values(op, vimplin, vimplout, fence);
1535  if (VERIFY_TREE) verify_tree();
1536 
1537  }
1538 
1539 
1540  /// Multiplication of function * vector of functions using recursive algorithm of mulxx
1541  template <typename L, typename R>
1542  void vmulXX(const Function<L,NDIM>& left,
1543  const std::vector< Function<R,NDIM> >& right,
1544  std::vector< Function<T,NDIM> >& result,
1545  double tol,
1546  bool fence) {
1548 
1549  std::vector<FunctionImpl<T,NDIM>*> vresult(right.size());
1550  std::vector<const FunctionImpl<R,NDIM>*> vright(right.size());
1551  for (unsigned int i=0; i<right.size(); ++i) {
1552  result[i].set_impl(left,false);
1553  vresult[i] = result[i].impl.get();
1554  vright[i] = right[i].get_impl().get();
1555  }
1556 
1557  left.world().gop.fence(); // Is this still essential? Yes.
1558  vresult[0]->mulXXvec(left.get_impl().get(), vright, vresult, tol, fence);
1559  }
1560 
1561  /// Same as \c operator* but with optional fence and no automatic reconstruction
1562 
1563  /// f or g are on-demand functions
1564  template<typename L, typename R>
1565  void mul_on_demand(const Function<L,NDIM>& f, const Function<R,NDIM>& g, bool fence=true) {
1566  const FunctionImpl<L,NDIM>* fimpl=f.get_impl().get();
1567  const FunctionImpl<R,NDIM>* gimpl=g.get_impl().get();
1568  if (fimpl->is_on_demand() and gimpl->is_on_demand()) {
1569  MADNESS_EXCEPTION("can't multiply two on-demand functions",1);
1570  }
1571 
1572  if (fimpl->is_on_demand()) {
1573  leaf_op<R,NDIM> leaf_op1(gimpl);
1574  impl->multiply(leaf_op1,gimpl,fimpl,fence);
1575  } else {
1576  leaf_op<L,NDIM> leaf_op1(fimpl);
1577  impl->multiply(leaf_op1,fimpl,gimpl,fence);
1578  }
1579  }
1580 
1581  /// sparse transformation of a vector of functions ... private
1582  template <typename R, typename Q>
1583  void vtransform(const std::vector< Function<R,NDIM> >& v,
1584  const Tensor<Q>& c,
1585  std::vector< Function<T,NDIM> >& vresult,
1586  double tol,
1587  bool fence=true) {
1589  vresult[0].impl->vtransform(vimpl(v), c, vimpl(vresult), tol, fence);
1590  }
1591 
1592  /// This is replaced with alpha*left + beta*right ... private
1593  template <typename L, typename R>
1595  T beta, const Function<R,NDIM>& right, bool fence) {
1597  left.verify();
1598  right.verify();
1599  MADNESS_ASSERT(left.is_compressed() && right.is_compressed());
1600  if (VERIFY_TREE) left.verify_tree();
1601  if (VERIFY_TREE) right.verify_tree();
1602  impl.reset(new implT(*left.get_impl(), left.get_pmap(), false));
1603  impl->gaxpy(alpha,*left.get_impl(),beta,*right.get_impl(),fence);
1604  return *this;
1605  }
1606 
1607  /// This is replaced with mapdim(f) ... private
1608  Function<T,NDIM>& mapdim(const Function<T,NDIM>& f, const std::vector<long>& map, bool fence) {
1610  f.verify();
1611  if (VERIFY_TREE) f.verify_tree();
1612  for (std::size_t i=0; i<NDIM; ++i) MADNESS_ASSERT(map[i]>=0 && static_cast<std::size_t>(map[i])<NDIM);
1613  impl.reset(new implT(*f.impl, f.get_pmap(), false));
1614  impl->mapdim(*f.impl,map,fence);
1615  return *this;
1616  }
1617 
1618  /// This is replaced with mirror(f) ... private
1619 
1620  /// similar to mapdim, but maps from x to -x, y to -y, and so on
1621  /// Example: mirror a 3d function on the xy plane: mirror={1,1,-1}
1622  /// @param[in] mirror array of -1 and 1, corresponding to mirror or not
1623  Function<T,NDIM>& mirror(const Function<T,NDIM>& f, const std::vector<long>& mirrormap, bool fence) {
1625  f.verify();
1626  if (VERIFY_TREE) f.verify_tree();
1627  for (std::size_t i=0; i<NDIM; ++i) MADNESS_ASSERT((mirrormap[i]==1) or (mirrormap[i]==-1));
1628  impl.reset(new implT(*f.impl, f.get_pmap(), false));
1629  impl->mirror(*f.impl,mirrormap,fence);
1630  return *this;
1631  }
1632 
1633  /// This is replaced with mirror(map(f)) ... private
1634 
1635  /// first map then mirror!
1636  /// mirror is similar to mapdim, but maps from x to -x, y to -y, and so on
1637  /// Example: mirror a 3d function on the xy plane: mirror={1,1,-1}
1638  /// Example: c4 rotation of a 3d function around the z axis:
1639  /// x->y, y->-x, z->z: map(1,0,2); mirror(-1,1,1)
1640  /// @param[in] map array holding dimensions
1641  /// @param[in] mirror array of -1 and 1, corresponding to mirror or not
1643  const std::vector<long>& map, const std::vector<long>& mirror,
1644  bool fence) {
1646  f.verify();
1647  if (VERIFY_TREE) f.verify_tree();
1648  for (std::size_t i=0; i<mirror.size(); ++i) MADNESS_ASSERT((mirror[i]==1) or (mirror[i]==-1));
1649  for (std::size_t i=0; i<map.size(); ++i) MADNESS_ASSERT(map[i]>=0 && static_cast<std::size_t>(map[i])<NDIM);
1650 
1651  impl.reset(new implT(*f.impl, f.get_pmap(), false));
1652  impl->map_and_mirror(*f.impl,map,mirror,fence);
1653  return *this;
1654  }
1655 
1656 
1657  /// check symmetry of a function by computing the 2nd derivative
1658  double check_symmetry() const {
1659 
1661  if (VERIFY_TREE) verify_tree();
1662  double local = impl->check_symmetry_local();
1663  impl->world.gop.sum(local);
1664  impl->world.gop.fence();
1665  double asy=sqrt(local);
1666  if (this->world().rank()==0) print("asymmetry wrt particle",asy);
1668  return asy;
1669  }
1670 
1671  /// reduce the rank of the coefficient tensors
1672  Function<T,NDIM>& reduce_rank(const double thresh=0.0, const bool fence=true) {
1673  verify();
1674  double thresh1= (thresh==0.0) ? impl->get_tensor_args().thresh : thresh;
1675  impl->reduce_rank(thresh1,fence);
1676  return *this;
1677  }
1678 
1679  /// remove all nodes with level higher than n
1680  Function<T,NDIM>& chop_at_level(const int n, const bool fence=true) {
1681  verify();
1683  impl->chop_at_level(n,true);
1685  return *this;
1686  }
1687  };
1688 
1689 // template <typename T, typename opT, std::size_t NDIM>
1690  template <typename T, typename opT, std::size_t NDIM>
1691  Function<T,NDIM> multiop_values(const opT& op, const std::vector< Function<T,NDIM> >& vf) {
1692  Function<T,NDIM> r;
1693  r.set_impl(vf[0], false);
1694  r.multiop_values(op, vf);
1695  return r;
1696  }
1697 
1698  /// Returns new function equal to alpha*f(x) with optional fence
1699  template <typename Q, typename T, std::size_t NDIM>
1700  Function<TENSOR_RESULT_TYPE(Q,T),NDIM>
1701  mul(const Q alpha, const Function<T,NDIM>& f, bool fence=true) {
1702  PROFILE_FUNC;
1703  f.verify();
1704  if (VERIFY_TREE) f.verify_tree();
1705  Function<TENSOR_RESULT_TYPE(Q,T),NDIM> result;
1706  result.set_impl(f, false);
1707  result.get_impl()->scale_oop(alpha,*f.get_impl(),fence);
1708  return result;
1709  }
1710 
1711 
1712  /// Returns new function equal to f(x)*alpha with optional fence
1713  template <typename Q, typename T, std::size_t NDIM>
1714  Function<TENSOR_RESULT_TYPE(Q,T),NDIM>
1715  mul(const Function<T,NDIM>& f, const Q alpha, bool fence=true) {
1716  PROFILE_FUNC;
1717  return mul(alpha,f,fence);
1718  }
1719 
1720 
1721  /// Returns new function equal to f(x)*alpha
1722 
1723  /// Using operator notation forces a global fence after each operation
1724  template <typename Q, typename T, std::size_t NDIM>
1725  Function<TENSOR_RESULT_TYPE(Q,T),NDIM>
1727  return mul(alpha, f, true);
1728  }
1729 
1730  /// Returns new function equal to alpha*f(x)
1731 
1732  /// Using operator notation forces a global fence after each operation
1733  template <typename Q, typename T, std::size_t NDIM>
1734  Function<TENSOR_RESULT_TYPE(Q,T),NDIM>
1736  return mul(alpha, f, true);
1737  }
1738 
1739  /// Sparse multiplication --- left and right must be reconstructed and if tol!=0 have tree of norms already created
1740  template <typename L, typename R,std::size_t NDIM>
1741  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1742  mul_sparse(const Function<L,NDIM>& left, const Function<R,NDIM>& right, double tol, bool fence=true) {
1743  PROFILE_FUNC;
1744  left.verify();
1745  right.verify();
1746  MADNESS_ASSERT(left.is_reconstructed() and right.is_reconstructed());
1747  if (VERIFY_TREE) left.verify_tree();
1748  if (VERIFY_TREE) right.verify_tree();
1749 
1750  Function<TENSOR_RESULT_TYPE(L,R),NDIM> result;
1751  result.set_impl(left, false);
1752  result.get_impl()->mulXX(left.get_impl().get(), right.get_impl().get(), tol, fence);
1753  return result;
1754  }
1755 
1756  /// Same as \c operator* but with optional fence and no automatic reconstruction
1757  template <typename L, typename R,std::size_t NDIM>
1758  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1759  mul(const Function<L,NDIM>& left, const Function<R,NDIM>& right, bool fence=true) {
1760  return mul_sparse(left,right,0.0,fence);
1761  }
1762 
1763  /// Generate new function = op(left,right) where op acts on the function values
1764  template <typename L, typename R, typename opT, std::size_t NDIM>
1765  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1766  binary_op(const Function<L,NDIM>& left, const Function<R,NDIM>& right, const opT& op, bool fence=true) {
1767  PROFILE_FUNC;
1768  if (!left.is_reconstructed()) left.reconstruct();
1769  if (!right.is_reconstructed()) right.reconstruct();
1770 
1771  Function<TENSOR_RESULT_TYPE(L,R),NDIM> result;
1772  result.set_impl(left, false);
1773  result.get_impl()->binaryXX(left.get_impl().get(), right.get_impl().get(), op, fence);
1774  return result;
1775  }
1776 
1777  /// Out of place application of unary operation to function values with optional fence
1778  template <typename Q, typename opT, std::size_t NDIM>
1779  Function<typename opT::resultT, NDIM>
1780  unary_op(const Function<Q,NDIM>& func, const opT& op, bool fence=true) {
1781  if (!func.is_reconstructed()) func.reconstruct();
1783  if (VERIFY_TREE) func.verify_tree();
1784  result.set_impl(func, false);
1785  result.get_impl()->unaryXXvalues(func.get_impl().get(), op, fence);
1786  return result;
1787  }
1788 
1789 
1790  /// Out of place application of unary operation to scaling function coefficients with optional fence
1791  template <typename Q, typename opT, std::size_t NDIM>
1792  Function<typename opT::resultT, NDIM>
1793  unary_op_coeffs(const Function<Q,NDIM>& func, const opT& op, bool fence=true) {
1794  if (!func.is_reconstructed()) func.reconstruct();
1796  return result.unary_op_coeffs(func,op,fence);
1797  }
1798 
1799  /// Use the vmra/mul(...) interface instead
1800 
1801  /// This so that we don't have to have friend functions in a different header.
1802  ///
1803  /// If using sparsity (tol != 0) you must have created the tree of norms
1804  /// already for both left and right.
1805  template <typename L, typename R, std::size_t D>
1806  std::vector< Function<TENSOR_RESULT_TYPE(L,R),D> >
1807  vmulXX(const Function<L,D>& left, const std::vector< Function<R,D> >& vright, double tol, bool fence=true) {
1808  if (vright.size() == 0) return std::vector< Function<TENSOR_RESULT_TYPE(L,R),D> >();
1809  std::vector< Function<TENSOR_RESULT_TYPE(L,R),D> > vresult(vright.size());
1810  vresult[0].vmulXX(left, vright, vresult, tol, fence);
1811  return vresult;
1812  }
1813 
1814  /// Multiplies two functions with the new result being of type TensorResultType<L,R>
1815 
1816  /// Using operator notation forces a global fence after each operation but also
1817  /// enables us to automatically reconstruct the input functions as required.
1818  template <typename L, typename R, std::size_t NDIM>
1819  Function<TENSOR_RESULT_TYPE(L,R), NDIM>
1820  operator*(const Function<L,NDIM>& left, const Function<R,NDIM>& right) {
1821  if (!left.is_reconstructed()) left.reconstruct();
1822  if (!right.is_reconstructed()) right.reconstruct();
1823  MADNESS_ASSERT(not (left.is_on_demand() or right.is_on_demand()));
1824  return mul(left,right,true);
1825  }
1826 
1827  /// Performs a Hartree/outer product on the two given low-dimensional function vectors
1828 
1829  /// @return result(x,y) = \sum_i f_i(x) g_i(y)
1830  template<typename T, std::size_t KDIM, std::size_t LDIM>
1831  Function<T,KDIM+LDIM>
1832  hartree_product(const std::vector<Function<T,KDIM>>& left, const std::vector<Function<T,LDIM>>& right) {
1833 
1834  MADNESS_CHECK_THROW(left.size()==right.size(), "hartree_product: left and right must have same size");
1835  if (left.size()==0) return Function<T,KDIM+LDIM>();
1836 
1838 
1840  .k(left.front().k()).thresh(thresh);
1841  Function<T,KDIM+LDIM> result=factory.empty();
1842 
1843  // some prep work
1846  std::vector<std::shared_ptr<FunctionImpl<T,KDIM>>> vleft=get_impl(left);
1847  std::vector<std::shared_ptr<FunctionImpl<T,LDIM>>> vright=get_impl(right);
1848 
1849  result.do_hartree_product(vleft,vright);
1850 
1851  return result;
1852 
1853  }
1854 
1855  /// Performs a Hartree product on the two given low-dimensional functions
1856  template<typename T, std::size_t KDIM, std::size_t LDIM>
1857  Function<T,KDIM+LDIM>
1858  hartree_product(const Function<T,KDIM>& left2, const Function<T,LDIM>& right2) {
1859  typedef std::vector<Function<T,KDIM>> vector;
1860  return hartree_product(vector({left2}),vector({right2}));
1861  }
1862 
1863  /// Performs a Hartree product on the two given low-dimensional functions
1864  template<typename T, std::size_t KDIM, std::size_t LDIM, typename opT>
1865  Function<T,KDIM+LDIM>
1867  const opT& op) {
1868 
1869  // we need both sum and difference coeffs for error estimation
1870  Function<T,KDIM>& left = const_cast< Function<T,KDIM>& >(left2);
1871  Function<T,LDIM>& right = const_cast< Function<T,LDIM>& >(right2);
1872 
1874 
1876  .k(left.k()).thresh(thresh);
1877  Function<T,KDIM+LDIM> result=factory.empty();
1878 
1879  if (result.world().rank()==0) {
1880  print("incomplete FunctionFactory in Function::hartree_product");
1881  print("thresh: ", thresh);
1882  }
1883  bool same=(left2.get_impl()==right2.get_impl());
1884 
1885  // some prep work
1886  left.make_nonstandard(true, true);
1887  right.make_nonstandard(true, true);
1888 
1889  std::vector<std::shared_ptr<FunctionImpl<T,KDIM>>> vleft;
1890  std::vector<std::shared_ptr<FunctionImpl<T,LDIM>>> vright;
1891  vleft.push_back(left.get_impl());
1892  vright.push_back(right.get_impl());
1893  result.do_hartree_product(vleft,right,&op);
1894 
1895  left.standard(false);
1896  if (not same) right.standard(false);
1897  left2.world().gop.fence();
1898 
1899  return result;
1900  }
1901 
1902  /// adds beta*right only left: alpha*left + beta*right optional fence and no automatic compression
1903 
1904  /// left and right might live in different worlds, the accumulation is non-blocking
1905  template <typename L, typename R,std::size_t NDIM>
1906  void
1908  TENSOR_RESULT_TYPE(L,R) beta, const Function<R,NDIM>& right, bool fence=true) {
1909  PROFILE_FUNC;
1910  Function<TENSOR_RESULT_TYPE(L,R),NDIM> result;
1911  left.gaxpy(alpha, right, beta, fence);
1912  }
1913 
1914  /// Returns new function alpha*left + beta*right optional fence and no automatic compression
1915  template <typename L, typename R,std::size_t NDIM>
1916  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1918  TENSOR_RESULT_TYPE(L,R) beta, const Function<R,NDIM>& right, bool fence=true) {
1919  PROFILE_FUNC;
1920  Function<TENSOR_RESULT_TYPE(L,R),NDIM> result;
1921  return result.gaxpy_oop(alpha, left, beta, right, fence);
1922  }
1923 
1924  /// Same as \c operator+ but with optional fence and no automatic compression
1925  template <typename L, typename R,std::size_t NDIM>
1926  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1927  add(const Function<L,NDIM>& left, const Function<R,NDIM>& right, bool fence=true) {
1928  return gaxpy_oop(TENSOR_RESULT_TYPE(L,R)(1.0), left,
1929  TENSOR_RESULT_TYPE(L,R)(1.0), right, fence);
1930  }
1931 
1932 
1933  /// Returns new function alpha*left + beta*right optional fence, having both addends reconstructed
1934  template<typename T, std::size_t NDIM>
1936  const double beta, const Function<T,NDIM>& right, const bool fence=true) {
1937  Function<T,NDIM> result;
1938  result.set_impl(right,false);
1939 
1942  result.get_impl()->gaxpy_oop_reconstructed(alpha,*left.get_impl(),beta,*right.get_impl(),fence);
1943  return result;
1944 
1945  }
1946 
1947  /// Adds two functions with the new result being of type TensorResultType<L,R>
1948 
1949  /// Using operator notation forces a global fence after each operation
1950  template <typename L, typename R, std::size_t NDIM>
1951  Function<TENSOR_RESULT_TYPE(L,R), NDIM>
1952  operator+(const Function<L,NDIM>& left, const Function<R,NDIM>& right) {
1953  if (VERIFY_TREE) left.verify_tree();
1954  if (VERIFY_TREE) right.verify_tree();
1955 
1956  // no compression for high-dimensional functions
1957  if (NDIM==6) {
1958  left.reconstruct();
1959  right.reconstruct();
1960  return gaxpy_oop_reconstructed(1.0,left,1.0,right,true);
1961  } else {
1962  if (!left.is_compressed()) left.compress();
1963  if (!right.is_compressed()) right.compress();
1964  return add(left,right,true);
1965  }
1966  }
1967 
1968  /// Same as \c operator- but with optional fence and no automatic compression
1969  template <typename L, typename R,std::size_t NDIM>
1970  Function<TENSOR_RESULT_TYPE(L,R),NDIM>
1971  sub(const Function<L,NDIM>& left, const Function<R,NDIM>& right, bool fence=true) {
1972  return gaxpy_oop(TENSOR_RESULT_TYPE(L,R)(1.0), left,
1973  TENSOR_RESULT_TYPE(L,R)(-1.0), right, fence);
1974  }
1975 
1976 
1977  /// Subtracts two functions with the new result being of type TensorResultType<L,R>
1978 
1979  /// Using operator notation forces a global fence after each operation
1980  template <typename L, typename R, std::size_t NDIM>
1981  Function<TENSOR_RESULT_TYPE(L,R), NDIM>
1982  operator-(const Function<L,NDIM>& left, const Function<R,NDIM>& right) {
1983  PROFILE_FUNC;
1984  // no compression for high-dimensional functions
1985  if (NDIM==6) {
1986  left.reconstruct();
1987  right.reconstruct();
1988  return gaxpy_oop_reconstructed(1.0,left,-1.0,right,true);
1989  } else {
1990  if (!left.is_compressed()) left.compress();
1991  if (!right.is_compressed()) right.compress();
1992  return sub(left,right,true);
1993  }
1994  }
1995 
1996  /// Create a new copy of the function with different distribution and optional fence
1997 
1998  /// Works in either basis. Different distributions imply
1999  /// asynchronous communication and the optional fence is
2000  /// collective.
2001  template <typename T, std::size_t NDIM>
2003  const std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& pmap,
2004  bool fence = true) {
2005  PROFILE_FUNC;
2006  f.verify();
2007  Function<T,NDIM> result;
2008  typedef FunctionImpl<T,NDIM> implT;
2009  result.set_impl(std::shared_ptr<implT>(new implT(*f.get_impl(), pmap, false)));
2010  result.get_impl()->copy_coeffs(*f.get_impl(), fence);
2011  if (VERIFY_TREE) result.verify_tree();
2012  return result;
2013  }
2014 
2015  /// Create a new copy of the function with the same distribution and optional fence
2016  template <typename T, std::size_t NDIM>
2017  Function<T,NDIM> copy(const Function<T,NDIM>& f, bool fence = true) {
2018  PROFILE_FUNC;
2019  return copy(f, f.get_pmap(), fence);
2020  }
2021 
2022  /// Type conversion implies a deep copy. No communication except for optional fence.
2023 
2024  /// Works in either basis but any loss of precision may result in different errors
2025  /// in applied in a different basis.
2026  ///
2027  /// The new function is formed with the options from the default constructor.
2028  ///
2029  /// There is no automatic type conversion since this is generally a rather dangerous
2030  /// thing and because there would be no way to make the fence optional.
2031  template <typename T, typename Q, std::size_t NDIM>
2032  Function<Q,NDIM> convert(const Function<T,NDIM>& f, bool fence = true) {
2033  PROFILE_FUNC;
2034  f.verify();
2035  Function<Q,NDIM> result;
2036  result.set_impl(f, false);
2037  result.get_impl()->copy_coeffs(*f.get_impl(), fence);
2038  return result;
2039  }
2040 
2041 
2042  /// Return the complex conjugate of the input function with the same distribution and optional fence
2043 
2044  /// !!! The fence is actually not optional in the current implementation !!!
2045  template <typename T, std::size_t NDIM>
2046  Function<T,NDIM> conj(const Function<T,NDIM>& f, bool fence = true) {
2047  PROFILE_FUNC;
2048  Function<T,NDIM> result = copy(f,true);
2049  return result.conj(fence);
2050  }
2051 
2052  /// Apply operator on a hartree product of two low-dimensional functions
2053 
2054  /// Supposed to be something like result= G( f(1)*f(2))
2055  /// the hartree product is never constructed explicitly, but its coeffs are
2056  /// constructed on the fly and processed immediately.
2057  /// @param[in] op the operator
2058  /// @param[in] f1 function of particle 1
2059  /// @param[in] f2 function of particle 2
2060  /// @param[in] fence if we shall fence
2061  /// @return a function of dimension NDIM=LDIM+LDIM
2062  template <typename opT, typename T, std::size_t LDIM>
2063  Function<TENSOR_RESULT_TYPE(typename opT::opT,T), LDIM+LDIM>
2064  apply(const opT& op, const std::vector<Function<T,LDIM>>& f1, const std::vector<Function<T,LDIM>>& f2, bool fence=true) {
2065 
2066  World& world=f1.front().world();
2067 
2068  typedef TENSOR_RESULT_TYPE(T,typename opT::opT) resultT;
2069  typedef std::vector<Function<T,LDIM>> vecfuncL;
2070 
2071  vecfuncL& ff1 = const_cast< vecfuncL& >(f1);
2072  vecfuncL& ff2 = const_cast< vecfuncL& >(f2);
2073 
2074  bool same=(ff1[0].get_impl()==ff2[0].get_impl());
2075 
2076  reconstruct(world,f1,false);
2077  reconstruct(world,f2,false);
2078  world.gop.fence();
2079  // keep the leaves! They are assumed to be there later
2080  // even for modified op we need NS form for the hartree_leaf_op
2081  for (auto& f : f1) f.make_nonstandard(true,false);
2082  for (auto& f : f2) f.make_nonstandard(true,false);
2083  world.gop.fence();
2084 
2085 
2088  Function<resultT,LDIM+LDIM> result=factory.empty().fence();
2089 
2090  result.get_impl()->reset_timer();
2091  op.reset_timer();
2092 
2093  // will fence here
2094  for (size_t i=0; i<f1.size(); ++i)
2095  result.get_impl()->recursive_apply(op, f1[i].get_impl().get(),f2[i].get_impl().get(),false);
2096  world.gop.fence();
2097 
2098  if (op.print_timings) {
2099  result.get_impl()->print_timer();
2100  op.print_timer();
2101  }
2102 
2103  result.get_impl()->finalize_apply(); // need fence before reconstruct
2104 
2105  if (op.modified()) {
2106  result.get_impl()->trickle_down(true);
2107  } else {
2108  result.get_impl()->reconstruct(true);
2109  }
2110  standard(world,ff1,false);
2111  if (not same) standard(world,ff2,false);
2112 
2113  return result;
2114  }
2115 
2116 
2117  /// Apply operator ONLY in non-standard form - required other steps missing !!
2118  template <typename opT, typename R, std::size_t NDIM>
2119  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM>
2120  apply_only(const opT& op, const Function<R,NDIM>& f, bool fence=true) {
2121  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> result;
2122 
2123  constexpr std::size_t OPDIM=opT::opdim;
2124  constexpr bool low_dim=(OPDIM*2==NDIM); // apply on some dimensions only
2125 
2126  // specialized version for 3D
2127  if (NDIM <= 3 and (not low_dim)) {
2128  result.set_impl(f, false);
2129  result.get_impl()->apply(op, *f.get_impl(), fence);
2130 
2131  } else { // general version for higher dimension
2132  //bool print_timings=false;
2133  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> r1;
2134 
2135  result.set_impl(f, false);
2136  r1.set_impl(f, false);
2137 
2138  result.get_impl()->reset_timer();
2139  op.reset_timer();
2140 
2141  result.get_impl()->apply_source_driven(op, *f.get_impl(), fence);
2142 
2143  // recursive_apply is about 20% faster than apply_source_driven
2144  //result.get_impl()->recursive_apply(op, f.get_impl().get(),
2145  // r1.get_impl().get(),true); // will fence here
2146 
2147  }
2148 
2149  return result;
2150  }
2151 
2152  /// Apply operator in non-standard form
2153 
2154  /// Returns a new function with the same distribution
2155  ///
2156  /// !!! For the moment does NOT respect fence option ... always fences
2157  /// if the operator acts on one particle only the result will be sorted as
2158  /// g.particle=1: g(f) = \int g(x,x') f(x',y) dx' = result(x,y)
2159  /// g.particle=2: g(f) = \int g(y,y') f(x,y') dy' = result(x,y)
2160  /// for the second case it will notably *not* be as it is implemented in the partial inner product!
2161  /// g.particle=2 g(f) = result(x,y)
2162  /// inner(g(y,y'),f(x,y'),1,1) = result(y,x)
2163  /// also note the confusion with the counting of the particles/integration variables
2164  template <typename opT, typename R, std::size_t NDIM>
2165  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM>
2166  apply(const opT& op, const Function<R,NDIM>& f, bool fence=true) {
2167 
2168  typedef TENSOR_RESULT_TYPE(typename opT::opT,R) resultT;
2169  Function<R,NDIM>& ff = const_cast< Function<R,NDIM>& >(f);
2170  Function<resultT, NDIM> result;
2171 
2172  MADNESS_ASSERT(not f.is_on_demand());
2173  bool print_timings=op.print_timings;
2174 
2175  if (VERIFY_TREE) ff.verify_tree();
2176  ff.reconstruct();
2177  if (print_timings) ff.print_size("ff in apply after reconstruct");
2178 
2179  if (op.modified()) {
2180 
2182 // ff.get_impl()->make_redundant(true);
2183  result = apply_only(op, ff, fence);
2184  ff.get_impl()->undo_redundant(false);
2185  result.get_impl()->trickle_down(true);
2186 
2187  } else {
2188 
2189  // saves the standard() step, which is very expensive in 6D
2190 // Function<R,NDIM> fff=copy(ff);
2191  Function<R,NDIM> fff=(ff);
2192  fff.make_nonstandard(op.doleaves, true);
2193  if (print_timings) fff.print_size("ff in apply after make_nonstandard");
2194  if ((print_timings) and (f.world().rank()==0)) {
2195  fff.get_impl()->timer_filter.print("filter");
2196  fff.get_impl()->timer_compress_svd.print("compress_svd");
2197  }
2198  result = apply_only(op, fff, fence);
2199  result.get_impl()->set_tree_state(nonstandard_after_apply);
2200  ff.world().gop.fence();
2201  if (print_timings) result.print_size("result after apply_only");
2202 
2203  // svd-tensors need some post-processing
2204  if (result.get_impl()->get_tensor_type()==TT_2D) {
2205  double elapsed=result.get_impl()->finalize_apply();
2206  if (print_timings) printf("time in finalize_apply %8.2f\n",elapsed);
2207  }
2208  if (print_timings) {
2209  result.get_impl()->print_timer();
2210  op.print_timer();
2211  }
2212 
2213  result.get_impl()->reconstruct(true);
2214 
2215 // fff.clear();
2216  if (op.destructive()) {
2217  ff.world().gop.fence();
2218  ff.clear();
2219  } else {
2220  ff.standard();
2221  }
2222 
2223  }
2224  if (print_timings) result.print_size("result after reconstruction");
2225  return result;
2226  }
2227 
2228 
2229  template <typename opT, typename R, std::size_t NDIM>
2230  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM>
2231  apply_1d_realspace_push(const opT& op, const Function<R,NDIM>& f, int axis, bool fence=true) {
2232  PROFILE_FUNC;
2233  Function<R,NDIM>& ff = const_cast< Function<R,NDIM>& >(f);
2234  if (VERIFY_TREE) ff.verify_tree();
2235  ff.reconstruct();
2236 
2237  Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> result;
2238 
2239  result.set_impl(ff, false);
2240  result.get_impl()->apply_1d_realspace_push(op, ff.get_impl().get(), axis, fence);
2241  result.get_impl()->set_tree_state(redundant_after_merge);
2242  return result;
2243  }
2244 
2245 
2246  /// Generate a new function by reordering dimensions ... optional fence
2247 
2248  /// You provide an array of dimension NDIM that maps old to new dimensions
2249  /// according to
2250  /// \code
2251  /// newdim = mapdim[olddim]
2252  /// \endcode
2253  /// Works in either scaling function or wavelet basis.
2254  ///
2255  /// Would be easy to modify this to also change the procmap here
2256  /// if desired but presently it uses the same procmap as f.
2257  template <typename T, std::size_t NDIM>
2258  Function<T,NDIM>
2259  mapdim(const Function<T,NDIM>& f, const std::vector<long>& map, bool fence=true) {
2260  PROFILE_FUNC;
2261  Function<T,NDIM> result;
2262  return result.mapdim(f,map,fence);
2263  }
2264 
2265  /// Generate a new function by mirroring within the dimensions .. optional fence
2266 
2267  /// similar to mapdim
2268  /// @param[in] mirror array with -1 and 1, corresponding to mirror this dimension or not
2269  template <typename T, std::size_t NDIM>
2270  Function<T,NDIM>
2271  mirror(const Function<T,NDIM>& f, const std::vector<long>& mirrormap, bool fence=true) {
2272  PROFILE_FUNC;
2273  Function<T,NDIM> result;
2274  return result.mirror(f,mirrormap,fence);
2275  }
2276 
2277  /// This is replaced with mirror(map(f)), optional fence
2278 
2279  /// first map then mirror!
2280  /// mirror is similar to mapdim, but maps from x to -x, y to -y, and so on
2281  /// Example: mirror a 3d function on the xy plane: mirror={1,1,-1}
2282  /// Example: c4 rotation of a 3d function around the z axis:
2283  /// x->y, y->-x, z->z: map(1,0,2); mirror(-1,1,1)
2284  /// @param[in] map array holding dimensions
2285  /// @param[in] mirror array of -1 and 1, corresponding to mirror or not
2286  template <typename T, std::size_t NDIM>
2287  Function<T,NDIM>
2288  map_and_mirror(const Function<T,NDIM>& f, const std::vector<long>& map,
2289  const std::vector<long>& mirror, bool fence=true) {
2290  PROFILE_FUNC;
2291  Function<T,NDIM> result;
2292  return result.map_and_mirror(f,map,mirror,fence);
2293  }
2294 
2295 
2296  /// swap particles 1 and 2
2297 
2298  /// param[in] f a function of 2 particles f(1,2)
2299  /// return the input function with particles swapped g(1,2) = f(2,1)
2300  template <typename T, std::size_t NDIM>
2301  typename std::enable_if_t<NDIM%2==0, Function<T,NDIM>>
2303  // this could be done more efficiently for SVD, but it works decently
2304  std::vector<long> map(NDIM);
2305  constexpr std::size_t LDIM=NDIM/2;
2306  static_assert(LDIM*2==NDIM);
2307  for (auto d=0; d<LDIM; ++d) {
2308  map[d]=d+LDIM;
2309  map[d+LDIM]=d;
2310  }
2311 // map[0]=3;
2312 // map[1]=4;
2313 // map[2]=5; // 2 -> 1
2314 // map[3]=0;
2315 // map[4]=1;
2316 // map[5]=2; // 1 -> 2
2317  return mapdim(f,map);
2318  }
2319 
2320  /// symmetrize a function
2321 
2322  /// @param[in] symmetry possibilities are:
2323  /// (anti-) symmetric particle permutation ("sy_particle", "antisy_particle")
2324  /// symmetric mirror plane ("xy", "xz", "yz")
2325  /// @return a new function symmetrized according to the input parameter
2326  template <typename T, std::size_t NDIM>
2327  Function<T,NDIM>
2328  symmetrize(const Function<T,NDIM>& f, const std::string symmetry, bool fence=true) {
2329  Function<T,NDIM> result;
2330 
2331  MADNESS_ASSERT(NDIM==6); // works only for pair functions
2332  std::vector<long> map(NDIM);
2333 
2334  // symmetric particle permutation
2335  if (symmetry=="sy_particle") {
2336  map[0]=3; map[1]=4; map[2]=5;
2337  map[3]=0; map[4]=1; map[5]=2;
2338  } else if (symmetry=="cx") {
2339  map[0]=0; map[1]=2; map[2]=1;
2340  map[3]=3; map[4]=5; map[5]=4;
2341 
2342  } else if (symmetry=="cy") {
2343  map[0]=2; map[1]=1; map[2]=0;
2344  map[3]=5; map[4]=4; map[5]=3;
2345 
2346  } else if (symmetry=="cz") {
2347  map[0]=1; map[1]=0; map[2]=2;
2348  map[3]=4; map[4]=3; map[5]=5;
2349 
2350  } else {
2351  if (f.world().rank()==0) {
2352  print("unknown parameter in symmetrize:",symmetry);
2353  }
2354  MADNESS_EXCEPTION("unknown parameter in symmetrize",1);
2355  }
2356 
2357  result.mapdim(f,map,true); // need to fence here
2358  result.get_impl()->average(*f.get_impl());
2359 
2360  return result;
2361  }
2362 
2363 
2364 
2365  /// multiply a high-dimensional function with a low-dimensional function
2366 
2367  /// @param[in] f NDIM function of 2 particles: f=f(1,2)
2368  /// @param[in] g LDIM function of 1 particle: g=g(1) or g=g(2)
2369  /// @param[in] particle if g=g(1) or g=g(2)
2370  /// @return h(1,2) = f(1,2) * g(p)
2371  template<typename T, std::size_t NDIM, std::size_t LDIM>
2372  Function<T,NDIM> multiply(const Function<T,NDIM> f, const Function<T,LDIM> g, const int particle, const bool fence=true) {
2373 
2374  static_assert(LDIM+LDIM==NDIM);
2375  MADNESS_ASSERT(particle==1 or particle==2);
2376 
2377  Function<T,NDIM> result;
2378  result.set_impl(f, false);
2379 
2380 // Function<T,NDIM>& ff = const_cast< Function<T,NDIM>& >(f);
2381 // Function<T,LDIM>& gg = const_cast< Function<T,LDIM>& >(g);
2382 
2383  f.change_tree_state(redundant,false);
2384  g.change_tree_state(redundant);
2385  FunctionImpl<T,NDIM>* fimpl=f.get_impl().get();
2386  FunctionImpl<T,LDIM>* gimpl=g.get_impl().get();
2387 
2388  result.get_impl()->multiply(fimpl,gimpl,particle);
2389  result.world().gop.fence();
2390 
2391  f.change_tree_state(reconstructed,false);
2392  g.change_tree_state(reconstructed);
2393  return result;
2394  }
2395 
2396 
2397  template <typename T, std::size_t NDIM>
2398  Function<T,NDIM>
2402  bool fence=true)
2403  {
2404  PROFILE_FUNC;
2406  other.reconstruct();
2407  result.get_impl()->project(*other.get_impl(),fence);
2408  return result;
2409  }
2410 
2411 
2412  /// Computes the scalar/inner product between two functions
2413 
2414  /// In Maple this would be \c int(conjugate(f(x))*g(x),x=-infinity..infinity)
2415  template <typename T, typename R, std::size_t NDIM>
2417  PROFILE_FUNC;
2418  return f.inner(g);
2419  }
2420 
2421 
2422  /// Computes the partial scalar/inner product between two functions, returns a low-dim function
2423 
2424  /// syntax similar to the inner product in tensor.h
2425  /// e.g result=inner<3>(f,g),{0},{1}) : r(x,y) = int f(x1,x) g(y,x1) dx1
2426  /// @param[in] task 0: everything, 1; prepare only (fence), 2: work only (no fence), 3: finalize only (fence)
2427  template<std::size_t NDIM, typename T, std::size_t LDIM, typename R, std::size_t KDIM,
2428  std::size_t CDIM = (KDIM + LDIM - NDIM) / 2>
2429  std::vector<Function<TENSOR_RESULT_TYPE(T, R), NDIM>>
2430  innerXX(const Function<T, LDIM>& f, const std::vector<Function<R, KDIM>>& vg, const std::array<int, CDIM> v1,
2431  const std::array<int, CDIM> v2, int task=0) {
2432  bool prepare = ((task==0) or (task==1));
2433  bool work = ((task==0) or (task==2));
2434  bool finish = ((task==0) or (task==3));
2435 
2436  static_assert((KDIM + LDIM - NDIM) % 2 == 0, "faulty dimensions in inner (partial version)");
2437  static_assert(KDIM + LDIM - 2 * CDIM == NDIM, "faulty dimensions in inner (partial version)");
2438 
2439  // contraction indices must be contiguous and either in the beginning or at the end
2440  for (int i=0; i<CDIM-1; ++i) MADNESS_CHECK((v1[i]+1)==v1[i+1]);
2441  MADNESS_CHECK((v1[0]==0) or (v1[CDIM-1]==LDIM-1));
2442 
2443  for (int i=0; i<CDIM-1; ++i) MADNESS_CHECK((v2[i]+1)==v2[i+1]);
2444  MADNESS_CHECK((v2[0]==0) or (v2[CDIM-1]==KDIM-1));
2445 
2446  MADNESS_CHECK(f.is_initialized());
2447  MADNESS_CHECK(vg[0].is_initialized());
2448  MADNESS_CHECK(f.world().id() == vg[0].world().id());
2449  // this needs to be run in a single world, so that all coefficients are local.
2450  // Use macrotasks if run on multiple processes.
2451  World& world=f.world();
2452  MADNESS_CHECK(world.size() == 1);
2453 
2454  if (prepare) {
2455  f.change_tree_state(nonstandard);
2457  world.gop.fence();
2458  f.get_impl()->compute_snorm_and_dnorm(false);
2459  for (auto& g : vg) g.get_impl()->compute_snorm_and_dnorm(false);
2460  world.gop.fence();
2461  }
2462 
2463  typedef TENSOR_RESULT_TYPE(T, R) resultT;
2464  std::vector<Function<resultT,NDIM>> result(vg.size());
2465  if (work) {
2466  world.gop.set_forbid_fence(true);
2467  for (int i=0; i<vg.size(); ++i) {
2468  result[i]=FunctionFactory<resultT,NDIM>(world)
2469  .k(f.k()).thresh(f.thresh()).empty().nofence();
2470  result[i].get_impl()->partial_inner(*f.get_impl(),*(vg[i]).get_impl(),v1,v2);
2471  result[i].get_impl()->set_tree_state(nonstandard_after_apply);
2472  }
2473  world.gop.set_forbid_fence(false);
2474  }
2475 
2476  if (finish) {
2477 
2478  world.gop.fence();
2479 // result.get_impl()->reconstruct(true);
2480 
2482 // result.reconstruct();
2483  // restore initial state of g and h
2484  auto erase_list = [] (const auto& funcimpl) {
2485  typedef typename std::decay_t<decltype(funcimpl)>::keyT keyTT;
2486  std::list<keyTT> to_be_erased;
2487  for (auto it=funcimpl.get_coeffs().begin(); it!=funcimpl.get_coeffs().end(); ++it) {
2488  const auto& key=it->first;
2489  const auto& node=it->second;
2490  if (not node.has_children()) to_be_erased.push_back(key);
2491  }
2492  return to_be_erased;
2493  };
2494 
2495  FunctionImpl<T,LDIM>& f_nc=const_cast<FunctionImpl<T,LDIM>&>(*f.get_impl());
2496  for (auto& key : erase_list(f_nc)) f_nc.get_coeffs().erase(key);
2497  for (auto& g : vg) {
2498  FunctionImpl<R,KDIM>& g_nc=const_cast<FunctionImpl<R,KDIM>&>(*g.get_impl());
2499  for (auto& key : erase_list(g_nc)) g_nc.get_coeffs().erase(key);
2500  }
2501  world.gop.fence();
2503  f_nc.reconstruct(false);
2504  world.gop.fence();
2505 
2506  }
2507 
2508  return result;
2509  }
2510 
2511 
2512  /// Computes the partial scalar/inner product between two functions, returns a low-dim function
2513 
2514  /// syntax similar to the inner product in tensor.h
2515  /// e.g result=inner<3>(f,g),{0},{1}) : r(x,y) = int f(x1,x) g(y,x1) dx1
2516  /// @param[in] task 0: everything, 1; prepare only (fence), 2: work only (no fence), 3: finalize only (fence)
2517  template<std::size_t NDIM, typename T, std::size_t LDIM, typename R, std::size_t KDIM,
2518  std::size_t CDIM = (KDIM + LDIM - NDIM) / 2>
2519  Function<TENSOR_RESULT_TYPE(T, R), NDIM>
2520  innerXX(const Function<T, LDIM>& f, const Function<R, KDIM>& g, const std::array<int, CDIM> v1,
2521  const std::array<int, CDIM> v2, int task=0) {
2522  return innerXX<NDIM,T,LDIM,R,KDIM>(f,std::vector<Function<R,KDIM>>({g}),v1,v2,task)[0];
2523  }
2524 
2525  /// Computes the partial scalar/inner product between two functions, returns a low-dim function
2526 
2527  /// syntax similar to the inner product in tensor.h
2528  /// e.g result=inner<3>(f,g),{0},{1}) : r(x,y) = int f(x1,x) g(y,x1) dx1
2529  template <typename T, std::size_t LDIM, typename R, std::size_t KDIM>
2530  Function<TENSOR_RESULT_TYPE(T,R),KDIM+LDIM-2>
2531  inner(const Function<T,LDIM>& f, const Function<R,KDIM>& g, const std::tuple<int> v1, const std::tuple<int> v2) {
2532  return innerXX<KDIM+LDIM-2>(f,g,
2533  std::array<int,1>({std::get<0>(v1)}),
2534  std::array<int,1>({std::get<0>(v2)}));
2535  }
2536 
2537  /// Computes the partial scalar/inner product between two functions, returns a low-dim function
2538 
2539  /// syntax similar to the inner product in tensor.h
2540  /// e.g result=inner<3>(f,g),{0,1},{1,2}) : r(y) = int f(x1,x2) g(y,x1,x2) dx1 dx2
2541  template <typename T, std::size_t LDIM, typename R, std::size_t KDIM>
2542  Function<TENSOR_RESULT_TYPE(T,R),KDIM+LDIM-4>
2543  inner(const Function<T,LDIM>& f, const Function<R,KDIM>& g, const std::tuple<int,int> v1, const std::tuple<int,int> v2) {
2544  return innerXX<KDIM+LDIM-4>(f,g,
2545  std::array<int,2>({std::get<0>(v1),std::get<1>(v1)}),
2546  std::array<int,2>({std::get<0>(v2),std::get<1>(v2)}));
2547  }
2548 
2549  /// Computes the partial scalar/inner product between two functions, returns a low-dim function
2550 
2551  /// syntax similar to the inner product in tensor.h
2552  /// e.g result=inner<3>(f,g),{1},{2}) : r(x,y,z) = int f(x,x1) g(y,z,x1) dx1
2553  template <typename T, std::size_t LDIM, typename R, std::size_t KDIM>
2554  Function<TENSOR_RESULT_TYPE(T,R),KDIM+LDIM-6>
2555  inner(const Function<T,LDIM>& f, const Function<R,KDIM>& g, const std::tuple<int,int,int> v1, const std::tuple<int,int,int> v2) {
2556  return innerXX<KDIM+LDIM-6>(f,g,
2557  std::array<int,3>({std::get<0>(v1),std::get<1>(v1),std::get<2>(v1)}),
2558  std::array<int,3>({std::get<0>(v2),std::get<1>(v2),std::get<2>(v2)}));
2559  }
2560 
2561 
2562 
2563  /// Computes the scalar/inner product between an MRA function and an external functor
2564 
2565  /// Currently this defaults to inner_adaptive, which might be more expensive
2566  /// than inner_ext since it loops over all leaf nodes. If you feel inner_ext
2567  /// is more efficient you need to call it directly
2568  /// @param[in] f MRA function
2569  /// @param[in] g functor
2570  /// @result inner(f,g)
2571  template <typename T, typename opT, std::size_t NDIM>
2572  TENSOR_RESULT_TYPE(T,typename opT::value_type) inner(const Function<T,NDIM>& f, const opT& g) {
2573  PROFILE_FUNC;
2574  std::shared_ptr< FunctionFunctorInterface<double,3> > func(new opT(g));
2575  return f.inner_adaptive(func);
2576  }
2577 
2578  /// Computes the scalar/inner product between an MRA function and an external functor
2579 
2580  /// Currently this defaults to inner_adaptive, which might be more expensive
2581  /// than inner_ext since it loops over all leaf nodes. If you feel inner_ext
2582  /// is more efficient you need to call it directly
2583  /// @param[in] g functor
2584  /// @param[in] f MRA function
2585  /// @result inner(f,g)
2586  template <typename T, typename opT, std::size_t NDIM>
2587  TENSOR_RESULT_TYPE(T,typename opT::value_type) inner(const opT& g, const Function<T,NDIM>& f) {
2588  return inner(f,g);
2589  }
2590 
2591  template <typename T, typename R, std::size_t NDIM>
2592  typename IsSupported<TensorTypeData<R>, Function<TENSOR_RESULT_TYPE(T,R),NDIM> >::type
2594  return (f*R(1.0)).add_scalar(r);
2595  }
2596 
2597  template <typename T, typename R, std::size_t NDIM>
2598  typename IsSupported<TensorTypeData<R>, Function<TENSOR_RESULT_TYPE(T,R),NDIM> >::type
2600  return (f*R(1.0)).add_scalar(r);
2601  }
2602 
2603  template <typename T, typename R, std::size_t NDIM>
2604  typename IsSupported<TensorTypeData<R>, Function<TENSOR_RESULT_TYPE(T,R),NDIM> >::type
2606  return (f*R(1.0)).add_scalar(-r);
2607  }
2608 
2609  template <typename T, typename R, std::size_t NDIM>
2610  typename IsSupported<TensorTypeData<R>, Function<TENSOR_RESULT_TYPE(T,R),NDIM> >::type
2612  return (f*R(-1.0)).add_scalar(r);
2613  }
2614 
2615  namespace detail {
2616  template <std::size_t NDIM>
2617  struct realop {
2618  typedef double resultT;
2620  return real(t);
2621  }
2622 
2623  template <typename Archive> void serialize (Archive& ar) {}
2624  };
2625 
2626  template <std::size_t NDIM>
2627  struct imagop {
2628  typedef double resultT;
2630  return imag(t);
2631  }
2632 
2633  template <typename Archive> void serialize (Archive& ar) {}
2634  };
2635 
2636  template <std::size_t NDIM>
2637  struct abssqop {
2638  typedef double resultT;
2640  Tensor<double> r = abs(t);
2641  return r.emul(r);
2642  }
2643 
2644  template <typename Archive> void serialize (Archive& ar) {}
2645  };
2646 
2647  template <std::size_t NDIM>
2648  struct absop {
2649  typedef double resultT;
2651  Tensor<double> r = abs(t);
2652  return r;
2653  }
2654 
2655  template <typename Archive> void serialize (Archive& ar) {}
2656  };
2657 
2658  }
2659 
2660  /// Returns a new function that is the real part of the input
2661  template <std::size_t NDIM>
2663  return unary_op_coeffs(z, detail::realop<NDIM>(), fence);
2664  }
2665 
2666  /// Returns a new function that is the real part of the input
2667  template <std::size_t NDIM>
2668  Function<double,NDIM> real(const Function<double,NDIM>& z, bool fence=true) {
2669  return copy(z);
2670  }
2671 
2672  /// Returns a new function that is the imaginary part of the input
2673  template <std::size_t NDIM>
2675  return unary_op_coeffs(z, detail::imagop<NDIM>(), fence);
2676  }
2677 
2678 
2679  /// Create a new function that is the square of f - global comm only if not reconstructed
2680  template <typename T, std::size_t NDIM>
2681  Function<T,NDIM> square(const Function<T,NDIM>& f, bool fence=true) {
2682  PROFILE_FUNC;
2683  Function<T,NDIM> result = copy(f,true); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2684  return result.square(true); //fence); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2685  }
2686 
2687  /// Create a new function that is the abs of f - global comm only if not reconstructed
2688  template <typename T, std::size_t NDIM>
2689  Function<T,NDIM> abs(const Function<T,NDIM>& f, bool fence=true) {
2690  PROFILE_FUNC;
2691  Function<T,NDIM> result = copy(f,true); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2692  return result.abs(true); //fence); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2693  }
2694 
2695  /// Create a new function that is the abs_square of f - global comm only if not reconstructed
2696  template <typename T, std::size_t NDIM>
2697  typename std::enable_if<!TensorTypeData<T>::iscomplex, Function<T,NDIM> >::type
2698  abs_square(const Function<T,NDIM>& f, bool fence=true) {
2699  PROFILE_FUNC;
2700  Function<T,NDIM> result = copy(f,true); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2701  return result.abs_square(true); //fence); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2702  }
2703 
2704  /// Create a new function that is the abs_square of f - global comm only if not reconstructed
2705  template <typename T, std::size_t NDIM>
2706  typename std::enable_if<TensorTypeData<T>::iscomplex, Function<typename Tensor<T>::scalar_type,NDIM> >::type
2707  abs_square(const Function<T,NDIM>& f, bool fence=true) {
2708  return unary_op(f, detail::abssqop<NDIM>(), fence);
2709  }
2710 
2711  /// Returns a new function that is the square of the absolute value of the input
2712  template <std::size_t NDIM>
2714  return unary_op(z, detail::abssqop<NDIM>(), fence);
2715  }
2716 
2717  /// Returns a new function that is the absolute value of the input
2718  template <std::size_t NDIM>
2720  return unary_op(z, detail::absop<NDIM>(), fence);
2721  }
2722 
2723 }
2724 
2725 #include <madness/mra/funcplot.h>
2726 
2727 namespace madness {
2728  namespace archive {
2729  template <class archiveT, class T, std::size_t NDIM>
2731  static inline void load(const ParallelInputArchive<archiveT>& ar, Function<T,NDIM>& f) {
2732  f.load(*ar.get_world(), ar);
2733  }
2734  };
2735 
2736  template <class archiveT, class T, std::size_t NDIM>
2738  static inline void store(const ParallelOutputArchive<archiveT>& ar, const Function<T,NDIM>& f) {
2739  f.store(ar);
2740  }
2741  };
2742  }
2743 
2744  template <class T, std::size_t NDIM>
2745  void save(const Function<T,NDIM>& f, const std::string name) {
2747  ar2 & f;
2748  }
2749 
2750  template <class T, std::size_t NDIM>
2751  void load(Function<T,NDIM>& f, const std::string name) {
2753  ar2 & f;
2754  }
2755 
2756 }
2757 
2758 namespace madness {
2759  // type traits to check if a template parameter is a Function
2760  template<typename>
2761  struct is_madness_function : std::false_type {};
2762 
2763  template<typename T, std::size_t NDIM>
2764  struct is_madness_function<madness::Function<T, NDIM>> : std::true_type {};
2765 }
2766 
2767 
2768 /* @} */
2769 
2770 #include <madness/mra/derivative.h>
2771 #include <madness/mra/operator.h>
2772 #include <madness/mra/functypedefs.h>
2773 #include <madness/mra/vmra.h>
2774 // #include <madness/mra/mraimpl.h> !!!!!!!!!!!!! NOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO !!!!!!!!!!!!!!!!!!
2775 
2776 #endif // MADNESS_MRA_MRA_H__INCLUDED
double q(double t)
Definition: DKops.h:18
This header should include pretty much everything needed for the parallel runtime.
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
This class is used to specify boundary conditions for all operators.
Definition: funcdefaults.h:101
CompositeFunctorInterface implements a wrapper of holding several functions and functors.
Definition: function_interface.h:165
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
FunctionFactory & empty()
Definition: function_factory.h:239
FunctionFactory & fence(bool fence=true)
Definition: function_factory.h:269
FunctionFactory & nofence()
Definition: function_factory.h:275
virtual FunctionFactory & thresh(double thresh)
Definition: function_factory.h:191
virtual FunctionFactory & k(int k)
Definition: function_factory.h:186
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition: funcimpl.h:941
const dcT & get_coeffs() const
Definition: mraimpl.h:322
void reconstruct(bool fence)
reconstruct this tree – respects fence
Definition: mraimpl.h:1490
bool is_on_demand() const
Definition: mraimpl.h:262
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition: funcimpl.h:124
A multiresolution adaptive numerical function.
Definition: mra.h:122
void print_tree_json(std::ostream &os=std::cout) const
Definition: mra.h:864
TENSOR_RESULT_TYPE(T, R) inner(const Function< R
Returns the inner product.
Function< T, NDIM > & gaxpy_oop(T alpha, const Function< L, NDIM > &left, T beta, const Function< R, NDIM > &right, bool fence)
This is replaced with alpha*left + beta*right ... private.
Definition: mra.h:1594
Function< T, NDIM > & fill_cuspy_tree(const bool fence=true)
Special refinement on 6D boxes where the electrons come close (meet)
Definition: mra.h:1179
void gaxpy_ext(const Function< L, NDIM > &left, T(*f)(const coordT &), T alpha, T beta, double tol, bool fence=true) const
Definition: mra.h:1358
const Function< T, NDIM > & refine(bool fence=true) const
Inplace autorefines the function using same test as for squaring.
Definition: mra.h:830
T inner_adaptive(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > f, const bool leaf_refine=true) const
Definition: mra.h:1343
void unaryop_coeff(const opT &op, bool fence=true)
Unary operation applied inplace to the coefficients.
Definition: mra.h:914
bool is_compressed() const
Returns true if compressed, false otherwise. No communication.
Definition: mra.h:462
Function< T, NDIM > & fill_tree(const opT &op, bool fence=true)
With this being an on-demand function, fill the MRA tree according to different criteria.
Definition: mra.h:1143
Function< T, NDIM > & fill_cuspy_tree(const opT &op, const bool fence=true)
Definition: mra.h:1166
void set_impl(const Function< R, NDIM > &f, bool zero=true)
Replace current FunctionImpl with a new one using the same parameters & map as f.
Definition: mra.h:642
bool autorefine() const
Retunrs.
Definition: mra.h:548
void print_size(const std::string name) const
print some info about this
Definition: mra.h:498
void broaden(const BoundaryConditions< NDIM > &bc=FunctionDefaults< NDIM >::get_bc(), bool fence=true) const
Inplace broadens support in scaling function basis.
Definition: mra.h:836
void print_info() const
Print a summary of the load balancing info.
Definition: mra.h:880
Future< Level > evaldepthpt(const coordT &xuser) const
Definition: mra.h:252
void norm_tree(bool fence=true) const
Initializes information about the function norm at all length scales.
Definition: mra.h:701
void load(World &world, Archive &ar)
Replaces this function with one loaded from an archive using the default processor map.
Definition: mra.h:1439
double norm2sq_local() const
Returns the square of the norm of the local function ... no communication.
Definition: mra.h:678
void sum_down(bool fence=true) const
Sums scaling coeffs down tree restoring state with coeffs only at leaves. Optional fence....
Definition: mra.h:798
Level depthpt(const coordT &xuser) const
Definition: mra.h:406
Function< T, NDIM > & abs(bool fence=true)
Returns *this for chaining.
Definition: mra.h:1068
const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > & get_pmap() const
Returns a shared pointer to the process map.
Definition: mra.h:656
void set_autorefine(bool value, bool fence=true)
Sets the value of the autorefine flag. Optional global fence.
Definition: mra.h:558
Function< T, NDIM > & fill_tree(const Function< R, NDIM > &g, bool fence=true)
With this being an on-demand function, fill the MRA tree according to different criteria.
Definition: mra.h:1125
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition: mra.h:1099
void print_tree_graphviz(std::ostream &os=std::cout) const
Process 0 prints a graphviz-formatted output of all nodes in the tree (collective)
Definition: mra.h:870
double norm2() const
Returns the 2-norm of the function ... global sum ... works in either basis.
Definition: mra.h:688
void change_tensor_type(const TensorArgs &targs, bool fence=true)
change the tensor type of the coefficients in the FunctionNode
Definition: mra.h:1472
World & world() const
Returns the world.
Definition: mra.h:648
TreeState gstate
Definition: mra.h:1275
TENSOR_RESULT_TYPE(T, R) inner_on_demand(const Function< R
Returns the inner product for one on-demand function.
T operator()(const coordT &xuser) const
Evaluates the function at a point in user coordinates. Collective operation.
Definition: mra.h:371
Function< T, NDIM > & map_and_mirror(const Function< T, NDIM > &f, const std::vector< long > &map, const std::vector< long > &mirror, bool fence)
This is replaced with mirror(map(f)) ... private.
Definition: mra.h:1642
Function< T, NDIM > & operator=(const Function< T, NDIM > &f)
Assignment is shallow. No communication, works in either basis.
Definition: mra.h:171
Function< T, NDIM > & multiop_values(const opT &op, const std::vector< Function< T, NDIM > > &vf)
This is replaced with op(vector of functions) ... private.
Definition: mra.h:1502
MADNESS_ASSERT(g.is_compressed())
static std::vector< std::shared_ptr< FunctionImpl< Q, D > > > vimpl(const std::vector< Function< Q, D > > &v)
Returns vector of FunctionImpl pointers corresponding to vector of functions.
Definition: mra.h:1493
const Function< T, NDIM > & change_tree_state(const TreeState finalstate, bool fence=true) const
changes tree state to given state
Definition: mra.h:785
Function< T, NDIM > & chop_at_level(const int n, const bool fence=true)
remove all nodes with level higher than n
Definition: mra.h:1680
void verify_tree() const
Verifies the tree data structure ... global sync implied.
Definition: mra.h:453
IsSupported< TensorTypeData< Q >, Function< T, NDIM > >::type & operator*=(const Q q)
Inplace scaling by a constant.
Definition: mra.h:1049
std::pair< bool, T > eval_local_only(const Vector< double, NDIM > &xuser, Level maxlevel) const
Evaluate function only if point is local returning (true,value); otherwise return (false,...
Definition: mra.h:220
Function< TENSOR_RESULT_TYPE(T, R), NDIM-LDIM > project_out(const Function< R, LDIM > &g, const int dim) const
project this on the low-dim function g: h(x) = <f(x,y) | g(y)>
Definition: mra.h:1398
impl world gop fence()
int k() const
Returns the number of multiwavelets (k). No communication.
Definition: mra.h:586
double thresh() const
Returns value of truncation threshold. No communication.
Definition: mra.h:567
Function< T, NDIM > & add_scalar(T t, bool fence=true)
Inplace add scalar. No communication except for optional fence.
Definition: mra.h:963
void set_thresh(double value, bool fence=true)
Sets the value of the truncation threshold. Optional global fence.
Definition: mra.h:577
void distribute(std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > newmap) const
distribute this function according to newmap
Definition: mra.h:669
Function< T, NDIM > & mirror(const Function< T, NDIM > &f, const std::vector< long > &mirrormap, bool fence)
This is replaced with mirror(f) ... private.
Definition: mra.h:1623
Future< long > evalR(const coordT &xuser) const
Evaluates the function rank at a point in user coordinates. Possible non-blocking comm.
Definition: mra.h:289
void unaryop(const opT &op, bool fence=true)
Inplace unary operation on function values.
Definition: mra.h:904
void standard(bool fence=true)
Converts the function standard compressed form. Possible non-blocking comm.
Definition: mra.h:748
void unaryop_node(const opT &op, bool fence=true)
Unary operation applied inplace to the nodes.
Definition: mra.h:924
std::size_t size() const
Returns the number of coefficients in the function ... collective global sum.
Definition: mra.h:538
bool is_on_demand() const
Definition: mra.h:636
std::size_t max_nodes() const
Returns the max number of nodes on a processor.
Definition: mra.h:523
std::shared_ptr< FunctionImpl< T, NDIM > > impl
Definition: mra.h:129
Function< T, NDIM > conj(bool fence=true)
Inplace complex conjugate. No communication except for optional fence.
Definition: mra.h:942
Function< T, NDIM > & scale(const Q q, bool fence=true)
Inplace, scale the function by a constant. No communication except for optional fence.
Definition: mra.h:953
void replicate(bool fence=true) const
replicate this function, generating a unique pmap
Definition: mra.h:663
Function< typename opT::resultT, NDIM > & unary_op_coeffs(const Function< Q, NDIM > &func, const opT &op, bool fence)
This is replaced with left*right ... private.
Definition: mra.h:1480
T trace_local() const
Returns local contribution to int(f(x),x) ... no communication.
Definition: mra.h:1090
void mul_on_demand(const Function< L, NDIM > &f, const Function< R, NDIM > &g, bool fence=true)
Same as operator* but with optional fence and no automatic reconstruction.
Definition: mra.h:1565
Vector< double, NDIM > coordT
Type of vector holding coordinates.
Definition: mra.h:139
void store(Archive &ar) const
Stores the function to an archive.
Definition: mra.h:1460
Function< T, NDIM > & square(bool fence=true)
Inplace squaring of function ... global comm only if not reconstructed.
Definition: mra.h:1059
std::size_t max_local_depth() const
Returns the maximum local depth of the function tree ... no communications.
Definition: mra.h:515
Function< T, NDIM > & operator-=(const Function< Q, NDIM > &other)
Inplace subtraction of functions in the wavelet basis.
Definition: mra.h:1028
void vtransform(const std::vector< Function< R, NDIM > > &v, const Tensor< Q > &c, std::vector< Function< T, NDIM > > &vresult, double tol, bool fence=true)
sparse transformation of a vector of functions ... private
Definition: mra.h:1583
Function< T, NDIM > & fill_tree(bool fence=true)
With this being an on-demand function, fill the MRA tree according to different criteria.
Definition: mra.h:1154
Function< T, NDIM > & truncate(double tol=0.0, bool fence=true)
Truncate the function with optional fence. Compresses with fence if not compressed.
Definition: mra.h:602
impl world gop sum(local)
TENSOR_RESULT_TYPE(T, R) local
std::size_t max_depth() const
Returns the maximum depth of the function tree ... collective global sum.
Definition: mra.h:507
Function()
Default constructor makes uninitialized function. No communication.
Definition: mra.h:154
std::size_t tree_size() const
Returns the number of nodes in the function tree ... collective global sum.
Definition: mra.h:491
FunctionImpl< T, NDIM > implT
Definition: mra.h:136
void clear(bool fence=true)
Clears the function as if constructed uninitialized. Optional fence.
Definition: mra.h:847
void refine_general(const opT &op, bool fence=true) const
Inplace autorefines the function. Optional fence. Possible non-blocking comm.
Definition: mra.h:811
static void doconj(const Key< NDIM >, Tensor< T > &t)
Definition: mra.h:934
const std::shared_ptr< FunctionImpl< T, NDIM > > & get_impl() const
Returns a shared-pointer to the implementation.
Definition: mra.h:614
TreeState state
Definition: mra.h:1274
TENSOR_RESULT_TYPE(T, R) inner_local(const Function< R
Returns local part of inner product ... throws if both not compressed.
void set_functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > functor)
Replace the current functor with the provided new one.
Definition: mra.h:631
bool impl_initialized() const
Definition: mra.h:132
Future< T > eval(const coordT &xuser) const
Evaluates the function at a point in user coordinates. Possible non-blocking comm.
Definition: mra.h:186
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition: mra.h:775
return local
Definition: mra.h:1299
Function< T, NDIM > & fill_nuclear_cuspy_tree(const opT &op, const size_t particle, const bool fence=true)
Definition: mra.h:1194
auto func
Definition: mra.h:1377
~Function()
Destruction of any underlying implementation is deferred to next global fence.
Definition: mra.h:178
void do_hartree_product(const std::vector< std::shared_ptr< FunctionImpl< T, LDIM >>> left, const std::vector< std::shared_ptr< FunctionImpl< T, KDIM >>> right, const opT *op)
perform the hartree product of f*g, invoked by result
Definition: mra.h:1221
Function(const Function< T, NDIM > &f)
Copy constructor is shallow. No communication, works in either basis.
Definition: mra.h:165
void set_impl(const std::shared_ptr< FunctionImpl< T, NDIM > > &impl)
Replace current FunctionImpl with provided new one.
Definition: mra.h:621
T operator()(double x, double y=0, double z=0, double xx=0, double yy=0, double zz=0) const
Evaluates the function at a point in user coordinates. Collective operation.
Definition: mra.h:385
T inner_ext_local(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > f, const bool leaf_refine=true, const bool keep_redundant=false) const
Definition: mra.h:1310
std::size_t min_nodes() const
Returns the min number of nodes on a processor.
Definition: mra.h:530
constexpr std::size_t LDIM
Definition: mra.h:1376
return impl inner_local * g())
void do_hartree_product(const std::vector< std::shared_ptr< FunctionImpl< T, LDIM >>> left, const std::vector< std::shared_ptr< FunctionImpl< T, KDIM >>> right)
perform the hartree product of f*g, invoked by result
Definition: mra.h:1235
Function< T, NDIM > & operator+=(const Function< Q, NDIM > &other)
Inplace addition of functions in the wavelet basis.
Definition: mra.h:1008
bool is_nonstandard() const
Returns true if nonstandard-compressed, false otherwise. No communication.
Definition: mra.h:484
void verify() const
Asserts that the function is initialized.
Definition: mra.h:142
Function< T, NDIM/2 > dirac_convolution(const bool fence=true) const
Definition: mra.h:1422
double err(const funcT &func) const
Returns an estimate of the difference ||this-func|| ... global sum performed.
Definition: mra.h:440
Function< T, NDIM > & reduce_rank(const double thresh=0.0, const bool fence=true)
reduce the rank of the coefficient tensors
Definition: mra.h:1672
T inner_ext(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > f, const bool leaf_refine=true, const bool keep_redundant=false) const
Definition: mra.h:1326
void make_redundant(bool fence=true)
Converts the function to redundant form, i.e. sum coefficients on all levels.
Definition: mra.h:760
double check_symmetry() const
check symmetry of a function by computing the 2nd derivative
Definition: mra.h:1658
void multi_to_multi_op_values(const opT &op, const std::vector< Function< T, NDIM > > &vin, std::vector< Function< T, NDIM > > &vout, const bool fence=true)
apply op on the input vector yielding an output vector of functions
Definition: mra.h:1521
Function< T, NDIM > & fill_nuclear_cuspy_tree(const size_t particle, const bool fence=true)
Special refinement on 6D boxes for the nuclear potentials (regularized with cusp, non-regularized wit...
Definition: mra.h:1207
FunctionFactory< T, NDIM > factoryT
Definition: mra.h:138
bool is_initialized() const
Returns true if the function is initialized.
Definition: mra.h:147
change_tree_state(state, false)
bool is_reconstructed() const
Returns true if reconstructed, false otherwise. No communication.
Definition: mra.h:473
MADNESS_ASSERT(is_compressed())
void vmulXX(const Function< L, NDIM > &left, const std::vector< Function< R, NDIM > > &right, std::vector< Function< T, NDIM > > &result, double tol, bool fence)
Multiplication of function * vector of functions using recursive algorithm of mulxx.
Definition: mra.h:1542
Function< T, NDIM > & mapdim(const Function< T, NDIM > &f, const std::vector< long > &map, bool fence)
This is replaced with mapdim(f) ... private.
Definition: mra.h:1608
double errsq_local(const funcT &func) const
Returns an estimate of the difference ||this-func||^2 from local data.
Definition: mra.h:425
Function< T, NDIM > & abs_square(bool fence=true)
Returns *this for chaining.
Definition: mra.h:1077
void make_nonstandard(bool keepleaves, bool fence=true) const
Compresses the function retaining scaling function coeffs. Possible non-blocking comm.
Definition: mra.h:734
void print_tree(std::ostream &os=std::cout) const
Process 0 prints a summary of all nodes in the tree (collective)
Definition: mra.h:857
FunctionNode< T, NDIM > nodeT
Definition: mra.h:137
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition: mra.h:721
Tensor< T > eval_cube(const Tensor< double > &cell, const std::vector< long > &npt, bool eval_refine=false) const
Evaluates a cube/slice of points (probably for plotting) ... collective but no fence necessary.
Definition: mra.h:330
Function(const factoryT &factory)
Constructor from FunctionFactory provides named parameter idiom. Possible non-blocking communication.
Definition: mra.h:158
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
void unaryop(T(*f)(T))
Inplace unary operation on function values.
Definition: mra.h:895
A future is a possibly yet unevaluated value.
Definition: future.h:373
remote_refT remote_ref(World &world) const
Returns a structure used to pass references to another process.
Definition: future.h:675
Key is the index for a node of the 2^NDIM-tree.
Definition: key.h:66
Definition: leafop.h:406
Definition: leafop.h:275
Traits class to specify support of numeric types.
Definition: type_data.h:56
A tensor is a multidimension array.
Definition: tensor.h:317
Tensor< T > & conj()
Inplace complex conjugate.
Definition: tensor.h:716
Tensor< T > & emul(const Tensor< T > &t)
Inplace multiply by corresponding elements of argument Tensor.
Definition: tensor.h:1798
void erase(const keyT &key)
Erases entry from container (non-blocking comm if remote)
Definition: worlddc.h:1105
Interface to be provided by any process map.
Definition: worlddc.h:82
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
bool set_forbid_fence(bool value)
Set forbid_fence flag to new value and return old value.
Definition: worldgop.h:674
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition: world.h:328
unsigned long id() const
Definition: world.h:313
WorldGopInterface & gop
Global operations.
Definition: world.h:205
World * get_world() const
Returns a pointer to the world.
Definition: parallel_archive.h:130
An archive for storing local or parallel data, wrapping a BinaryFstreamInputArchive.
Definition: parallel_archive.h:366
An archive for storing local or parallel data wrapping a BinaryFstreamOutputArchive.
Definition: parallel_archive.h:321
Objects that implement their own parallel archive interface should derive from this class.
Definition: parallel_archive.h:58
static const double R
Definition: csqrt.cc:46
Declaration and initialization of tree traversal functions and generic derivative.
double(* f1)(const coord_3d &)
Definition: derivatives.cc:55
double(* f2)(const coord_3d &)
Definition: derivatives.cc:56
const double delta
Definition: dielectric_external_field.cc:119
Provides FunctionDefaults and utilities for coordinate transformation.
Provides FunctionCommonData, FunctionImpl and FunctionFactory.
Defines/implements plotting interface for functions.
Provides typedefs to hide use of templates and to increase interoperability.
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
const double beta
Definition: gygi_soltion.cc:62
static const double v
Definition: hatom_sf_dirac.cc:20
Provides IndexIterator.
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
Multidimension Key for MRA tree and associated iterators.
Implements (2nd generation) static load/data balancing for functions.
#define max(a, b)
Definition: lda.h:51
#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
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#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
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:210
Header to declare stuff which has not yet found a home.
static const bool VERIFY_TREE
Definition: mra.h:57
double norm(const T &t)
Definition: adquad.h:42
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
Function< TENSOR_RESULT_TYPE(L, R), NDIM > sub(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator- but with optional fence and no automatic compression.
Definition: mra.h:1971
Function< T, KDIM+LDIM > hartree_product(const std::vector< Function< T, KDIM >> &left, const std::vector< Function< T, LDIM >> &right)
Performs a Hartree/outer product on the two given low-dimensional function vectors.
Definition: mra.h:1832
std::vector< CCPairFunction< T, NDIM > > operator+(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1047
double abs(double x)
Definition: complexfun.h:48
Function< TENSOR_RESULT_TYPE(L, R), NDIM > add(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator+ but with optional fence and no automatic compression.
Definition: mra.h:1927
std::enable_if_t< NDIM%2==0, Function< T, NDIM > > swap_particles(const Function< T, NDIM > &f)
swap particles 1 and 2
Definition: mra.h:2302
Function< T, NDIM > mapdim(const Function< T, NDIM > &f, const std::vector< long > &map, bool fence=true)
Generate a new function by reordering dimensions ... optional fence.
Definition: mra.h:2259
Function< TENSOR_RESULT_TYPE(Q, T), NDIM > mul(const Q alpha, const Function< T, NDIM > &f, bool fence=true)
Returns new function equal to alpha*f(x) with optional fence.
Definition: mra.h:1701
std::vector< CCPairFunction< T, NDIM > > operator-(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1055
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
Function< double, NDIM > abssq(const Function< double_complex, NDIM > &z, bool fence=true)
Returns a new function that is the square of the absolute value of the input.
Definition: mra.h:2713
Function< typename opT::resultT, NDIM > unary_op_coeffs(const Function< Q, NDIM > &func, const opT &op, bool fence=true)
Out of place application of unary operation to scaling function coefficients with optional fence.
Definition: mra.h:1793
Function< typename TensorTypeData< Q >::scalar_type, NDIM > abs_square(const Function< Q, NDIM > &func)
Definition: complexfun.h:121
Function< T, NDIM > multiply(const Function< T, NDIM > f, const Function< T, LDIM > g, const int particle, const bool fence=true)
multiply a high-dimensional function with a low-dimensional function
Definition: mra.h:2372
Function< T, NDIM > multiop_values(const opT &op, const std::vector< Function< T, NDIM > > &vf)
Definition: mra.h:1691
Function< T, NDIM > project(const Function< T, NDIM > &other, int k=FunctionDefaults< NDIM >::get_k(), double thresh=FunctionDefaults< NDIM >::get_thresh(), bool fence=true)
Definition: mra.h:2399
const std::vector< Function< T, NDIM > > & change_tree_state(const std::vector< Function< T, NDIM >> &v, const TreeState finalstate, const bool fence=true)
change tree state of the functions
Definition: vmra.h:277
TreeState
Definition: funcdefaults.h:58
@ nonstandard_after_apply
s and d coeffs, state after operator application
Definition: funcdefaults.h:63
@ redundant_after_merge
s coeffs everywhere, must be summed up to yield the result
Definition: funcdefaults.h:65
@ reconstructed
s coeffs at the leaves only
Definition: funcdefaults.h:59
@ nonstandard
s and d coeffs in internal nodes
Definition: funcdefaults.h:61
@ unknown
Definition: funcdefaults.h:67
@ compressed
d coeffs in internal nodes, s and d coeffs at the root
Definition: funcdefaults.h:60
@ redundant
s coeffs everywhere
Definition: funcdefaults.h:64
@ nonstandard_with_leaves
like nonstandard, with s coeffs at the leaves
Definition: funcdefaults.h:62
static void user_to_sim(const Vector< double, NDIM > &xuser, Vector< double, NDIM > &xsim)
Convert user coords (cell[][]) to simulation coords ([0,1]^ndim)
Definition: funcdefaults.h:524
void standard(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Generates standard form of a vector of functions.
Definition: vmra.h:260
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
Function< Q, NDIM > convert(const Function< T, NDIM > &f, bool fence=true)
Type conversion implies a deep copy. No communication except for optional fence.
Definition: mra.h:2032
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition: mra.h:2681
std::vector< std::shared_ptr< FunctionImpl< T, NDIM > > > get_impl(const std::vector< Function< T, NDIM >> &v)
Definition: vmra.h:663
Function< T, NDIM > map_and_mirror(const Function< T, NDIM > &f, const std::vector< long > &map, const std::vector< long > &mirror, bool fence=true)
This is replaced with mirror(map(f)), optional fence.
Definition: mra.h:2288
int Level
Definition: key.h:55
std::vector< Function< TENSOR_RESULT_TYPE(L, R), D > > vmulXX(const Function< L, D > &left, const std::vector< Function< R, D > > &vright, double tol, bool fence=true)
Use the vmra/mul(...) interface instead.
Definition: mra.h:1807
std::string get_mra_data_dir()
Definition: startup.cc:208
Function< TENSOR_RESULT_TYPE(L, R), NDIM > gaxpy_oop(TENSOR_RESULT_TYPE(L, R) alpha, const Function< L, NDIM > &left, TENSOR_RESULT_TYPE(L, R) beta, const Function< R, NDIM > &right, bool fence=true)
Returns new function alpha*left + beta*right optional fence and no automatic compression.
Definition: mra.h:1917
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
Function< TENSOR_RESULT_TYPE(L, R), NDIM > binary_op(const Function< L, NDIM > &left, const Function< R, NDIM > &right, const opT &op, bool fence=true)
Generate new function = op(left,right) where op acts on the function values.
Definition: mra.h:1766
TENSOR_RESULT_TYPE(T, R) inner(const Function< T
Computes the scalar/inner product between two functions.
Definition: vmra.h:1059
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > innerXX(const Function< T, LDIM > &f, const std::vector< Function< R, KDIM >> &vg, const std::array< int, CDIM > v1, const std::array< int, CDIM > v2, int task=0)
Computes the partial scalar/inner product between two functions, returns a low-dim function.
Definition: mra.h:2430
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
@ TT_2D
Definition: gentensor.h:120
NDIM & f
Definition: mra.h:2416
Function< T, NDIM > symmetrize(const Function< T, NDIM > &f, const std::string symmetry, bool fence=true)
symmetrize a function
Definition: mra.h:2328
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition: mra.h:1742
NDIM const Function< R, NDIM > & g
Definition: mra.h:2416
std::vector< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition: ccpairfunction.h:1084
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
double imag(double x)
Definition: complexfun.h:56
void startup(World &world, int argc, char **argv, bool doprint=false, bool make_stdcout_nice_to_reals=true)
Definition: startup.cc:64
std::string type(const PairType &n)
Definition: PNOParameters.h:18
static bool print_timings
Definition: SCF.cc:104
Function< typename opT::resultT, NDIM > unary_op(const Function< Q, NDIM > &func, const opT &op, bool fence=true)
Out of place application of unary operation to function values with optional fence.
Definition: mra.h:1780
Function< T, NDIM > gaxpy_oop_reconstructed(const double alpha, const Function< T, NDIM > &left, const double beta, const Function< T, NDIM > &right, const bool fence=true)
Returns new function alpha*left + beta*right optional fence, having both addends reconstructed.
Definition: mra.h:1935
Function< T, NDIM > mirror(const Function< T, NDIM > &f, const std::vector< long > &mirrormap, bool fence=true)
Generate a new function by mirroring within the dimensions .. optional fence.
Definition: mra.h:2271
void load(Function< T, NDIM > &f, const std::string name)
Definition: mra.h:2751
double real(double x)
Definition: complexfun.h:52
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_only(const opT &op, const Function< R, NDIM > &f, bool fence=true)
Apply operator ONLY in non-standard form - required other steps missing !!
Definition: mra.h:2120
std::string name(const FuncType &type, const int ex=-1)
Definition: ccpairfunction.h:28
void save(const Function< T, NDIM > &f, const std::string name)
Definition: mra.h:2745
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_1d_realspace_push(const opT &op, const Function< R, NDIM > &f, int axis, bool fence=true)
Definition: mra.h:2231
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
Implements most functionality of separated operators.
Implements ParallelInputArchive and ParallelOutputArchive for parallel serialization of data.
double Q(double a)
Definition: relops.cc:20
static const double c
Definition: relops.cc:10
static const double L
Definition: rk.cc:46
static const double thresh
Definition: rk.cc:45
static const long k
Definition: rk.cc:44
Definition: test_ar.cc:204
Definition: leafop.h:142
void serialize(Archive &ar)
Definition: mra.h:891
T(* f)(T)
Definition: mra.h:886
void operator()(const Key< NDIM > &key, Tensor< T > &t) const
Definition: mra.h:888
SimpleUnaryOpWrapper(T(*f)(T))
Definition: mra.h:887
void serialize(Archive &ar)
Definition: mra.h:824
bool operator()(implT *impl, const Key< NDIM > &key, const nodeT &t) const
Definition: mra.h:820
Definition: type_data.h:146
Definition: leafop.h:188
Definition: leafop.h:62
TensorArgs holds the arguments for creating a LowRankTensor.
Definition: gentensor.h:134
static void load(const ParallelInputArchive< archiveT > &ar, Function< T, NDIM > &f)
Definition: mra.h:2731
Default load of an object via serialize(ar, t).
Definition: archive.h:666
static void store(const ParallelOutputArchive< archiveT > &ar, const Function< T, NDIM > &f)
Definition: mra.h:2738
Default store of an object via serialize(ar, t).
Definition: archive.h:611
Definition: mra.h:2648
double resultT
Definition: mra.h:2649
void serialize(Archive &ar)
Definition: mra.h:2655
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition: mra.h:2650
Definition: mra.h:2637
void serialize(Archive &ar)
Definition: mra.h:2644
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition: mra.h:2639
double resultT
Definition: mra.h:2638
Definition: mra.h:2627
void serialize(Archive &ar)
Definition: mra.h:2633
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition: mra.h:2629
double resultT
Definition: mra.h:2628
Definition: mra.h:2617
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition: mra.h:2619
double resultT
Definition: mra.h:2618
void serialize(Archive &ar)
Definition: mra.h:2623
Definition: mra.h:110
Definition: funcimpl.h:606
returns true if the result of a hartree_product is a leaf node (compute norm & error)
Definition: funcimpl.h:496
Definition: mra.h:2761
Definition: mra.h:95
Definition: mra.h:98
Definition: funcimpl.h:560
Definition: lowrankfunction.h:332
Defines and implements most of Tensor.
#define UNARY_OPTIMIZED_ITERATOR(X, x, exp)
Definition: tensor_macros.h:658
const double alpha
Definition: test_coulomb.cc:54
int task(int i)
Definition: test_runtime.cpp:4
void d()
Definition: test_sig.cc:79
void e()
Definition: test_sig.cc:75
static const std::size_t NDIM
Definition: testpdiff.cc:42
std::size_t axis
Definition: testpdiff.cc:59
Defines operations on vectors of Functions.
Implements WorldContainer.
#define PROFILE_FUNC
Definition: worldprofile.h:209
#define PROFILE_MEMBER_FUNC(classname)
Definition: worldprofile.h:210