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
45#include <madness/misc/misc.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
57static const bool VERIFY_TREE = false; //true
58
59
60namespace 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>
68#include <madness/mra/indexit.h>
73#include <madness/mra/lbdeux.h>
75
76// some forward declarations
77namespace 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>
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
115namespace 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();
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();
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();
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);
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();
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();
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();
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.
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();
957 impl->scale_inplace(q,fence);
958 return *this;
959 }
960
961
962 /// Inplace add scalar. No communication except for optional fence.
965 verify();
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.
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();
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)
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>
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) {
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) {
1703 f.verify();
1704 if (VERIFY_TREE) f.verify_tree();
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) {
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) {
1744 left.verify();
1745 right.verify();
1747 if (VERIFY_TREE) left.verify_tree();
1748 if (VERIFY_TREE) right.verify_tree();
1749
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) {
1768 if (!left.is_reconstructed()) left.reconstruct();
1769 if (!right.is_reconstructed()) right.reconstruct();
1770
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();
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>
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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);
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) {
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) {
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) {
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) {
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);
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 {
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>
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) {
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>
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) {
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) {
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) {
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
2727namespace 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
2758namespace 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>
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 & nofence()
Definition function_factory.h:275
virtual FunctionFactory & thresh(double thresh)
Definition function_factory.h:191
FunctionFactory & fence(bool fence=true)
Definition function_factory.h:269
virtual FunctionFactory & k(int k)
Definition function_factory.h:186
FunctionFactory & empty()
Definition function_factory.h:239
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:942
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< 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
TENSOR_RESULT_TYPE(T, R) inner(const Function< R
Returns the inner product.
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
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
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/2 > dirac_convolution(const bool fence=true) const
Definition mra.h:1422
return impl inner_local * g())
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
Function< T, NDIM > & add_scalar(T t, bool fence=true)
Inplace add scalar. No communication except for optional fence.
Definition mra.h:963
void print_size(const std::string name) const
print some info about this
Definition mra.h:498
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
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
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 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
Function< T, NDIM > & abs_square(bool fence=true)
Returns *this for chaining.
Definition mra.h:1077
double norm2sq_local() const
Returns the square of the norm of the local function ... no communication.
Definition mra.h:678
IsSupported< TensorTypeData< Q >, Function< T, NDIM > >::type & operator*=(const Q q)
Inplace scaling by a constant.
Definition mra.h:1049
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 > & operator+=(const Function< Q, NDIM > &other)
Inplace addition of functions in the wavelet basis.
Definition mra.h:1008
Function< T, NDIM > & operator=(const Function< T, NDIM > &f)
Assignment is shallow. No communication, works in either basis.
Definition mra.h:171
void set_autorefine(bool value, bool fence=true)
Sets the value of the autorefine flag. Optional global fence.
Definition mra.h:558
World & world() const
Returns the world.
Definition mra.h:648
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition mra.h:1099
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
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< 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
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
Function< T, NDIM > & fill_cuspy_tree(const opT &op, const bool fence=true)
Definition mra.h:1166
void change_tensor_type(const TensorArgs &targs, bool fence=true)
change the tensor type of the coefficients in the FunctionNode
Definition mra.h:1472
const Function< T, NDIM > & refine(bool fence=true) const
Inplace autorefines the function using same test as for squaring.
Definition mra.h:830
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
MADNESS_ASSERT(g.is_compressed())
Function< T, NDIM > & square(bool fence=true)
Inplace squaring of function ... global comm only if not reconstructed.
Definition mra.h:1059
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
void verify_tree() const
Verifies the tree data structure ... global sync implied.
Definition mra.h:453
Future< Level > evaldepthpt(const coordT &xuser) const
Definition mra.h:252
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 > & mapdim(const Function< T, NDIM > &f, const std::vector< long > &map, bool fence)
This is replaced with mapdim(f) ... private.
Definition mra.h:1608
impl world gop fence()
int k() const
Returns the number of multiwavelets (k). No communication.
Definition mra.h:586
const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > & get_pmap() const
Returns a shared pointer to the process map.
Definition mra.h:656
double thresh() const
Returns value of truncation threshold. No communication.
Definition mra.h:567
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
void set_thresh(double value, bool fence=true)
Sets the value of the truncation threshold. Optional global fence.
Definition mra.h:577
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
void distribute(std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > newmap) const
distribute this function according to newmap
Definition mra.h:669
void unaryop(const opT &op, bool fence=true)
Inplace unary operation on function values.
Definition mra.h:904
const std::shared_ptr< FunctionImpl< T, NDIM > > & get_impl() const
Returns a shared-pointer to the implementation.
Definition mra.h:614
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
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
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
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
Function< T, NDIM > & operator-=(const Function< Q, NDIM > &other)
Inplace subtraction of functions in the wavelet basis.
Definition mra.h:1028
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
Function< T, NDIM > & abs(bool fence=true)
Returns *this for chaining.
Definition mra.h:1068
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
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
std::size_t max_nodes() const
Returns the max number of nodes on a processor.
Definition mra.h:523
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
std::shared_ptr< FunctionImpl< T, NDIM > > impl
Definition mra.h:129
void replicate(bool fence=true) const
replicate this function, generating a unique pmap
Definition mra.h:663
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
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
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
std::size_t max_local_depth() const
Returns the maximum local depth of the function tree ... no communications.
Definition mra.h:515
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
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
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
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
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
TreeState state
Definition mra.h:1274
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
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
Function< T, NDIM > & fill_nuclear_cuspy_tree(const opT &op, const size_t particle, const bool fence=true)
Definition mra.h:1194
return local
Definition mra.h:1299
auto func
Definition mra.h:1377
~Function()
Destruction of any underlying implementation is deferred to next global fence.
Definition mra.h:178
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
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
double err(const funcT &func) const
Returns an estimate of the difference ||this-func|| ... global sum performed.
Definition mra.h:440
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
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
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
FunctionFactory< T, NDIM > factoryT
Definition mra.h:138
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition mra.h:721
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
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
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
double errsq_local(const funcT &func) const
Returns an estimate of the difference ||this-func||^2 from local data.
Definition mra.h:425
void make_nonstandard(bool keepleaves, bool fence=true) const
Compresses the function retaining scaling function coeffs. Possible non-blocking comm.
Definition mra.h:734
Future< T > eval(const coordT &xuser) const
Evaluates the function at a point in user coordinates. Possible non-blocking comm.
Definition mra.h:186
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
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
Function(const factoryT &factory)
Constructor from FunctionFactory provides named parameter idiom. Possible non-blocking communication.
Definition mra.h:158
void unaryop(T(*f)(T))
Inplace unary operation on function values.
Definition mra.h:895
Function< T, NDIM > conj(bool fence=true)
Inplace complex conjugate. No communication except for optional fence.
Definition mra.h:942
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 > & emul(const Tensor< T > &t)
Inplace multiply by corresponding elements of argument Tensor.
Definition tensor.h:1798
Tensor< T > & conj()
Inplace complex conjugate.
Definition tensor.h:716
A simple, fixed dimension vector.
Definition vector.h:64
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 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:182
#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:207
Header to declare stuff which has not yet found a home.
static const bool VERIFY_TREE
Definition mra.h:57
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
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
double abs(double x)
Definition complexfun.h:48
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< 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
Function< typename TensorTypeData< Q >::scalar_type, NDIM > abs_square(const Function< Q, NDIM > &func)
Definition complexfun.h:121
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
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
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
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
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< 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
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< 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< std::shared_ptr< FunctionImpl< T, NDIM > > > get_impl(const std::vector< Function< T, NDIM > > &v)
Definition vmra.h:663
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
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
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
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
Function< T, NDIM > multiop_values(const opT &op, const std::vector< Function< T, NDIM > > &vf)
Definition mra.h:1691
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
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
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
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
std::vector< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition ccpairfunction.h:1084
std::shared_ptr< FunctionFunctorInterface< double, 3 > > func(new opT(g))
int Level
Definition key.h:55
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
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< 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
std::string get_mra_data_dir()
Definition startup.cc:208
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
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
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
@ TT_2D
Definition gentensor.h:120
NDIM & f
Definition mra.h:2416
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
Function< T, NDIM > symmetrize(const Function< T, NDIM > &f, const std::string symmetry, bool fence=true)
symmetrize a function
Definition mra.h:2328
NDIM const Function< R, NDIM > & g
Definition mra.h:2416
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
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
double imag(double x)
Definition complexfun.h:56
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
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
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
static bool print_timings
Definition SCF.cc:104
std::vector< CCPairFunction< T, NDIM > > operator+(const std::vector< CCPairFunction< T, NDIM > > c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition ccpairfunction.h:1047
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
void load(Function< T, NDIM > &f, const std::string name)
Definition mra.h:2751
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
double real(double x)
Definition complexfun.h:52
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< 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
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
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
static const double d
Definition nonlinschro.cc:121
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
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition mra.h:2650
double resultT
Definition mra.h:2649
void serialize(Archive &ar)
Definition mra.h:2655
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
Tensor< double > operator()(const Key< NDIM > &key, const Tensor< double_complex > &t) const
Definition mra.h:2629
void serialize(Archive &ar)
Definition mra.h:2633
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:607
returns true if the result of a hartree_product is a leaf node (compute norm & error)
Definition funcimpl.h:497
Definition mra.h:2761
Definition mra.h:95
Definition mra.h:98
Definition funcimpl.h:561
Definition lowrankfunction.h:332
Defines and implements most of Tensor.
#define UNARY_OPTIMIZED_ITERATOR(X, x, exp)
Definition tensor_macros.h:658
AtomicInt sum
Definition test_atomicint.cc:46
double norm(const T i1)
Definition test_cloud.cc:72
int task(int i)
Definition test_runtime.cpp:4
void e()
Definition test_sig.cc:75
static const double alpha
Definition testcosine.cc:10
static const std::size_t NDIM
Definition testpdiff.cc:42
std::size_t axis
Definition testpdiff.cc:59
#define TENSOR_RESULT_TYPE(L, R)
This macro simplifies access to TensorResultType.
Definition type_data.h:205
Defines operations on vectors of Functions.
Implements WorldContainer.
#define PROFILE_FUNC
Definition worldprofile.h:209
#define PROFILE_MEMBER_FUNC(classname)
Definition worldprofile.h:210