MADNESS 0.10.1
mraimpl.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_MRAIMPL_H__INCLUDED
33#define MADNESS_MRA_MRAIMPL_H__INCLUDED
34
35#ifndef MPRAIMPLX
36#error "mraimpl.h should ONLY be included in one of the mraX.cc files (x=1..6)"
37#endif
38
39#include <memory>
40#include <math.h>
41#include <cmath>
46
49
50namespace std {
51 template <typename T>
52 bool isnan(const std::complex<T>& v) {
53 return ::std::isnan(v.real()) || ::std::isnan(v.imag());
54 }
55}
56
57/// \file mra/mraimpl.h
58/// \brief Declaration and initialization of static data, some implementation, some instantiation
59
60namespace madness {
61
62 // Definition and initialization of FunctionDefaults static members
63 // It cannot be an instance of FunctionFactory since we want to
64 // set the defaults independent of the data type.
65
66 template <typename T, std::size_t NDIM>
68 if (! two_scale_hg(k, &hg)) throw "failed to get twoscale coefficients";
69 hgT = copy(transpose(hg));
70
71 Slice sk(0,k-1), sk2(k,-1);
72 hgsonly = copy(hg(Slice(0,k-1),_));
73
74 h0 = copy(hg(sk,sk));
75 h1 = copy(hg(sk,sk2));
76 g0 = copy(hg(sk2,sk));
77 g1 = copy(hg(sk2,sk2));
78
79 h0T = copy(transpose(hg(sk,sk)));
80 h1T = copy(transpose(hg(sk,sk2)));
81 g0T = copy(transpose(hg(sk2,sk)));
82 g1T = copy(transpose(hg(sk2,sk2)));
83
84 }
85
86 template <typename T, std::size_t NDIM>
88 (int k, int npt, Tensor<double>& quad_x, Tensor<double>& quad_w,
89 Tensor<double>& quad_phi, Tensor<double>& quad_phiw, Tensor<double>& quad_phit) {
90 quad_x = Tensor<double>(npt); // point
91 quad_w = Tensor<double>(npt); // wheight
92 quad_phi = Tensor<double>(npt,k);
93 quad_phiw = Tensor<double>(npt,k);
94
95 gauss_legendre(npt,0.0,1.0,quad_x.ptr(),quad_w.ptr());
96 for (int mu=0; mu<npt; ++mu) {
97 double phi[200];
98 legendre_scaling_functions(quad_x(mu),k,phi);
99 for (int j=0; j<k; ++j) {
100 quad_phi(mu,j) = phi[j];
101 quad_phiw(mu,j) = quad_w(mu)*phi[j];
102 }
103 }
104 quad_phit = transpose(quad_phi);
105 }
106
107
108 template <typename T, std::size_t NDIM>
111 world.gop.fence(); // Make sure nothing is going on
112
113 // Verify consistency of compression status, existence and size of coefficients,
114 // and has_children() flag.
115 for (typename dcT::const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
116 const keyT& key = it->first;
117 const nodeT& node = it->second;
118 bool bad;
119
120 if (is_compressed()) {
121 if (node.has_children()) {
122 bad = (node.coeff().has_data()) and (node.coeff().dim(0) != 2*cdata.k);
123 }
124 else {
125 // bad = node.coeff().size() != 0;
126 bad = node.coeff().has_data();
128 }
129 else {
130 if (node.has_children()) {
131 // bad = node.coeff().size() != 0;
132 bad = node.coeff().has_data();
133 }
134 else {
135 bad = (node.coeff().has_data()) and ( node.coeff().dim(0) != cdata.k);
136 }
137 }
138
139 if (bad) {
140 print(world.rank(), "FunctionImpl: verify: INCONSISTENT TREE NODE, key =", key, ", node =", node,
141 ", dim[0] =",node.coeff().dim(0),", compressed =",is_compressed());
142 std::cout.flush();
143 MADNESS_EXCEPTION("FunctionImpl: verify: INCONSISTENT TREE NODE", 0);
144 }
145 }
146
147 // Ensure that parents and children exist appropriately
148 for (typename dcT::const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
149 const keyT& key = it->first;
150 const nodeT& node = it->second;
151
152 if (key.level() > 0) {
153 const keyT parent = key.parent();
154 typename dcT::const_iterator pit = coeffs.find(parent).get();
155 if (pit == coeffs.end()) {
156 print(world.rank(), "FunctionImpl: verify: MISSING PARENT",key,parent);
157 std::cout.flush();
158 MADNESS_EXCEPTION("FunctionImpl: verify: MISSING PARENT", 0);
159 }
160 const nodeT& pnode = pit->second;
161 if (!pnode.has_children()) {
162 print(world.rank(), "FunctionImpl: verify: PARENT THINKS IT HAS NO CHILDREN",key,parent);
163 std::cout.flush();
164 MADNESS_EXCEPTION("FunctionImpl: verify: PARENT THINKS IT HAS NO CHILDREN", 0);
165 }
166 }
167
168 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
169 typename dcT::const_iterator cit = coeffs.find(kit.key()).get();
170 if (cit == coeffs.end()) {
171 if (node.has_children()) {
172 print(world.rank(), "FunctionImpl: verify: MISSING CHILD",key,kit.key());
173 std::cout.flush();
174 MADNESS_EXCEPTION("FunctionImpl: verify: MISSING CHILD", 0);
175 }
176 }
177 else {
178 if (! node.has_children()) {
179 print(world.rank(), "FunctionImpl: verify: UNEXPECTED CHILD",key,kit.key());
180 std::cout.flush();
181 MADNESS_EXCEPTION("FunctionImpl: verify: UNEXPECTED CHILD", 0);
182 }
183 }
184 }
185 }
186
187 world.gop.fence();
188 }
189
190
191 template <typename T, std::size_t NDIM>
192 const std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >& FunctionImpl<T,NDIM>::get_pmap() const {
193 return coeffs.get_pmap();
194 }
195
196
197 /// perform: this= alpha*f + beta*g, invoked by result
198
199 /// f and g are reconstructed, so we can save on the compress operation,
200 /// walk down the joint tree, and add leaf coefficients; effectively refines
201 /// to common finest level.
202 /// @param[in] alpha prefactor for f
203 /// @param[in] f first addend
204 /// @param[in] beta prefactor for g
205 /// @param[in] g second addend
206 /// @return nothing, but leaves this's tree reconstructed and as sum of f and g
207 template <typename T, std::size_t NDIM>
209 const double beta, const implT& g, const bool fence) {
210
211 MADNESS_ASSERT(not f.is_compressed());
212 MADNESS_ASSERT(not g.is_compressed());
213
214 ProcessID owner = coeffs.owner(cdata.key0);
215 if (world.rank() == owner) {
216
219
220 typedef add_op coeff_opT;
221 coeff_opT coeff_op(ff,gg,alpha,beta);
222 typedef insert_op<T,NDIM> apply_opT;
223 apply_opT apply_op(this);
224
225 woT::task(world.rank(), &implT:: template forward_traverse<coeff_opT,apply_opT>,
226 coeff_op, apply_op, cdata.key0);
227
228 }
229 set_tree_state(reconstructed);
230 if (fence) world.gop.fence();
231 }
232
233 /// Returns true if the function is compressed.
234 template <typename T, std::size_t NDIM>
236 return (tree_state==compressed);
237 }
238
239 /// Returns true if the function is reconstructed.
240 template <typename T, std::size_t NDIM>
242 return (tree_state==reconstructed);
243 }
244
245 /// Returns true if the function is redundant.
246 template <typename T, std::size_t NDIM>
248 return (tree_state==redundant);
249 }
250
251 template <typename T, std::size_t NDIM>
253 return (tree_state==nonstandard);
254 }
255
256 template <typename T, std::size_t NDIM>
258 return (tree_state==nonstandard_with_leaves);
259 }
260
261 template <typename T, std::size_t NDIM>
263 return tree_state==on_demand;
264 }
265
266 template <typename T, std::size_t NDIM>
268 return (tree_state==nonstandard_with_leaves);
269 }
270
271 template <typename T, std::size_t NDIM>
273 set_tree_state(on_demand);
274// this->on_demand=true;
275 functor=functor1;
276 }
277
278 template <typename T, std::size_t NDIM>
279 std::shared_ptr<FunctionFunctorInterface<T,NDIM> > FunctionImpl<T,NDIM>::get_functor() {
280 MADNESS_ASSERT(this->functor);
281 return functor;
282 }
283
284 template <typename T, std::size_t NDIM>
285 std::shared_ptr<FunctionFunctorInterface<T,NDIM> > FunctionImpl<T,NDIM>::get_functor() const {
286 MADNESS_ASSERT(this->functor);
287 return functor;
288 }
289
290 template <typename T, std::size_t NDIM>
292// this->on_demand=false;
293 set_tree_state(unknown);
294 functor.reset();
295 }
296
297 template <typename T, std::size_t NDIM>
299
300 template <typename T, std::size_t NDIM>
302
303 template <typename T, std::size_t NDIM>
305
306 template <typename T, std::size_t NDIM>
308
309 template <typename T, std::size_t NDIM>
310 void FunctionImpl<T,NDIM>::set_thresh(double value) {thresh = value;}
311
312 template <typename T, std::size_t NDIM>
313 bool FunctionImpl<T,NDIM>::get_autorefine() const {return autorefine;}
314
315 template <typename T, std::size_t NDIM>
316 void FunctionImpl<T,NDIM>::set_autorefine(bool value) {autorefine = value;}
317
318 template <typename T, std::size_t NDIM>
319 int FunctionImpl<T,NDIM>::get_k() const {return k;}
320
321 template <typename T, std::size_t NDIM>
322 const typename FunctionImpl<T,NDIM>::dcT& FunctionImpl<T,NDIM>::get_coeffs() const {return coeffs;}
323
324 template <typename T, std::size_t NDIM>
326
327 template <typename T, std::size_t NDIM>
329
330 template <typename T, std::size_t NDIM>
331 void FunctionImpl<T,NDIM>::accumulate_timer(const double time) const {
332 timer_accumulate.accumulate(time);
333 }
334
335 template <typename T, std::size_t NDIM>
337 if (world.rank()==0) {
338 timer_accumulate.print("accumulate");
339 timer_target_driven.print("target_driven");
340 timer_lr_result.print("result2low_rank");
341 }
342 }
343
344 template <typename T, std::size_t NDIM>
346 if (world.rank()==0) {
347 timer_accumulate.reset();
348 timer_target_driven.reset();
349 timer_lr_result.reset();
350 }
351 }
352
353 /// Truncate according to the threshold with optional global fence
354
355 /// If thresh<=0 the default value of this->thresh is used
356 template <typename T, std::size_t NDIM>
357 void FunctionImpl<T,NDIM>::truncate(double tol, bool fence) {
358 // Cannot put tol into object since it would make a race condition
359 if (tol <= 0.0)
360 tol = thresh;
361 if (world.rank() == coeffs.owner(cdata.key0)) {
362 if (is_compressed()) {
363 truncate_spawn(cdata.key0,tol);
364 } else {
365 truncate_reconstructed_spawn(cdata.key0,tol);
366 }
367 }
368 if (fence)
369 world.gop.fence();
370 }
371
372 template <typename T, std::size_t NDIM>
374 return cdata.key0;
375 }
376
377 /// Print a plane ("xy", "xz", or "yz") containing the point x to file
378
379 /// works for all dimensions; we walk through the tree, and if a leaf node
380 /// inside the sub-cell touches the plane we print it in pstricks format
381 template <typename T, std::size_t NDIM>
382 void FunctionImpl<T,NDIM>::print_plane(const std::string filename, const int xaxis, const int yaxis, const coordT& el2) {
383
384 // get the local information
385 Tensor<double> localinfo=print_plane_local(xaxis,yaxis,el2);
386
387 // lump all the local information together, and gather on node0
388 std::vector<Tensor<double> > localinfo_vec(1,localinfo);
389 std::vector<Tensor<double> > printinfo=world.gop.concat0(localinfo_vec);
390 world.gop.fence();
391
392 // do the actual print
393 if (world.rank()==0) do_print_plane(filename,printinfo,xaxis,yaxis,el2);
394 }
395
396 /// collect the data for a plot of the MRA structure locally on each node
397
398 /// @param[in] xaxis the x-axis in the plot (can be any axis of the MRA box)
399 /// @param[in] yaxis the y-axis in the plot (can be any axis of the MRA box)
400 /// @param[in] el2
401 template <typename T, std::size_t NDIM>
402 Tensor<double> FunctionImpl<T,NDIM>::print_plane_local(const int xaxis, const int yaxis, const coordT& el2) {
403 coordT x_sim;
404 user_to_sim<NDIM>(el2,x_sim);
405 x_sim[0]+=1.e-10;
406
407 // dimensions are: (# boxes)(hue, x lo left, y lo left, x hi right, y hi right)
408 Tensor<double> plotinfo(coeffs.size(),5);
409 long counter=0;
410
411 // loop over local boxes, if the fit, add the info to the output tensor
412 typename dcT::const_iterator end = coeffs.end();
413 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
414 const keyT& key = it->first;
415 const nodeT& node = it->second;
416
417 // thisKeyContains ignores dim0 and dim1
418 if (key.thisKeyContains(x_sim,xaxis,yaxis) and node.is_leaf() and (node.has_coeff())) {
419
420 Level n=key.level();
422 // get the diametral edges of the node in the plotting plane
423 double scale=std::pow(0.5,double(n));
424 double xloleft = scale*l[xaxis];
425 double yloleft = scale*l[yaxis];
426 double xhiright = scale*(l[xaxis]+1);
427 double yhiright = scale*(l[yaxis]+1);
428
429 // convert back to user coordinates
430 Vector<double,4> user;
435
436
437 // if ((xloleft<-5.0) or (yloleft<-5.0) or (xhiright>5.0) or (yhiright>5.0)) continue;
438 if ((user[0]<-5.0) or (user[1]<-5.0) or (user[2]>5.0) or (user[3]>5.0)) continue;
439
440 // do rank or do error
441 double color=0.0;
442 if (1) {
443
444 const double maxrank=40;
445 do_convert_to_color hue(maxrank,false);
446 color=hue(node.coeff().rank());
447 } else {
448
449 // Make quadrature rule of higher order
450 const int npt = cdata.npt + 1;
451 Tensor<double> qx, qw, quad_phi, quad_phiw, quad_phit;
452 FunctionCommonData<T,NDIM>::_init_quadrature(k+1, npt, qx, qw, quad_phi, quad_phiw, quad_phit);
453 do_err_box< FunctionFunctorInterface<T,NDIM> > op(this, this->get_functor().get(), npt, qx, quad_phit, quad_phiw);
454
455 do_convert_to_color hue(1000.0,true);
456 double error=op(it);
457 error=sqrt(error);//*pow(2,key.level()*6);
458 color=hue(error);
459 }
460
461 plotinfo(counter,0)=color;
462 plotinfo(counter,1)=user[0];
463 plotinfo(counter,2)=user[1];
464 plotinfo(counter,3)=user[2];
465 plotinfo(counter,4)=user[3];
466 ++counter;
467 }
468 }
469
470 // shrink the info
471 if (counter==0) plotinfo=Tensor<double>();
472 else plotinfo=plotinfo(Slice(0,counter-1),Slice(_));
473 return plotinfo;
474 }
475
476 /// print the MRA structure
477 template <typename T, std::size_t NDIM>
478 void FunctionImpl<T,NDIM>::do_print_plane(const std::string filename, std::vector<Tensor<double> > plotinfo,
479 const int xaxis, const int yaxis, const coordT el2) {
480
481 // invoke only on master node
482 MADNESS_ASSERT(world.rank()==0);
483
484 // prepare file
485 FILE * pFile;
486 pFile = fopen(filename.c_str(), "w");
488
489
490 fprintf(pFile,"\\psset{unit=1cm}\n");
491 fprintf(pFile,"\\begin{pspicture}(%4.2f,%4.2f)(%4.2f,%4.2f)\n",
492 // cell(xaxis,0),cell(xaxis,1),cell(yaxis,0),cell(yaxis,1));
493 -5.0,-5.0,5.0,5.0);
494 fprintf(pFile,"\\pslinewidth=0.1pt\n");
495
496 for (std::vector<Tensor<double> >::const_iterator it=plotinfo.begin(); it!=plotinfo.end(); ++it) {
497
498 Tensor<double> localinfo=*it;
499 if (localinfo.has_data()) {
500
501 for (long i=0; i<localinfo.dim(0); ++i) {
502
503 fprintf(pFile,"\\newhsbcolor{mycolor}{%8.4f 1.0 0.7}\n",localinfo(i,0));
504 fprintf(pFile,"\\psframe["//linewidth=0.5pt,"
505 "fillstyle=solid,"
506 "fillcolor=mycolor]"
507 "(%12.8f,%12.8f)(%12.8f,%12.8f)\n",
508 localinfo(i,1),localinfo(i,2),localinfo(i,3),localinfo(i,4));
509 }
510 }
511 }
512
513
514 fprintf(pFile,"\\end{pspicture}\n");
515 fclose(pFile);
516 }
517
518 /// print the grid (the roots of the quadrature of each leaf box)
519 /// of this function in user xyz coordinates
520 template <typename T, std::size_t NDIM>
521 void FunctionImpl<T,NDIM>::print_grid(const std::string filename) const {
522
523 // get the local information
524 std::vector<keyT> local_keys=local_leaf_keys();
525
526 // lump all the local information together, and gather on node0
527 std::vector<keyT> all_keys=world.gop.concat0(local_keys);
528 world.gop.fence();
529
530 // do the actual print
531 if (world.rank()==0) do_print_grid(filename,all_keys);
532
533 }
534
535 /// return the keys of the local leaf boxes
536 template <typename T, std::size_t NDIM>
537 std::vector<typename FunctionImpl<T,NDIM>::keyT> FunctionImpl<T,NDIM>::local_leaf_keys() const {
538
539 // coeffs.size is maximum number of keys (includes internal keys)
540 std::vector<keyT> keys(coeffs.size());
541
542 // loop over local boxes, if they are leaf boxes add their quadrature roots
543 // to the output tensor
544 int i=0;
545 typename dcT::const_iterator end = coeffs.end();
546 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
547 const keyT& key = it->first;
548 const nodeT& node = it->second;
549 if (node.is_leaf()) keys[i++]=key;
550 }
551
552 // shrink the vector to number of leaf keys
553 keys.resize(i);
554 return keys;
555 }
556
557 /// print the grid in xyz format
558
559 /// the quadrature points and the key information will be written to file,
560 /// @param[in] filename where the quadrature points will be written to
561 /// @param[in] keys all leaf keys
562 template <typename T, std::size_t NDIM>
563 void FunctionImpl<T,NDIM>::do_print_grid(const std::string filename, const std::vector<keyT>& keys) const {
564 // invoke only on master node
565 MADNESS_ASSERT(world.rank()==0);
566
567 // the quadrature points in simulation coordinates of the root node
568 const Tensor<double> qx=cdata.quad_x;
569 const size_t npt = qx.dim(0);
570
571 // the number of coordinates (grid point tuples) per box ({x1},{x2},{x3},..,{xNDIM})
572 long npoints=power<NDIM>(npt);
573 // the number of boxes
574 long nboxes=keys.size();
575
576 // prepare file
577 FILE * pFile;
578 pFile = fopen(filename.c_str(), "w");
579
580 fprintf(pFile,"%ld\n",npoints*nboxes);
581 fprintf(pFile,"%ld points per box and %ld boxes \n",npoints,nboxes);
582
583 // loop over all leaf boxes
584 typename std::vector<keyT>::const_iterator key_it=keys.begin();
585 for (key_it=keys.begin(); key_it!=keys.end(); ++key_it) {
586
587 const keyT& key=*key_it;
588 fprintf(pFile,"# key: %8d",key.level());
589 for (size_t d=0; d<NDIM; d++) fprintf(pFile,"%8d",int(key.translation()[d]));
590 fprintf(pFile,"\n");
591
592 // this is borrowed from fcube
593 const Vector<Translation,NDIM>& l = key.translation();
594 const Level n = key.level();
595 const double h = std::pow(0.5,double(n));
596 coordT c; // will hold the point in user coordinates
597
600
601 if (NDIM == 3) {
602 for (size_t i=0; i<npt; ++i) {
603 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
604 for (size_t j=0; j<npt; ++j) {
605 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
606 for (size_t k=0; k<npt; ++k) {
607 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
608 // grid weights
609 // double scale = pow(0.5,0.5*NDIM*key.level())*
610 // sqrt(FunctionDefaults<NDIM>::get_cell_volume());
611 // double w=cdata.quad_phiw[i]*cdata.quad_phiw[j]*cdata.quad_phiw[k];
612
613 fprintf(pFile,"%18.12f %18.12f %18.12f\n",c[0],c[1],c[2]);
614 // fprintf(pFile,"%18.12e %18.12e %18.12e %18.12e\n",c[0],c[1],c[2],w*scale);
615 }
616 }
617 }
618 } else {
619 MADNESS_EXCEPTION("only NDIM=3 in print_grid",0);
620 }
621 }
622 fclose(pFile);
623 }
624
625
626 /// Returns the truncation threshold according to truncate_method
627 template <typename T, std::size_t NDIM>
628 double FunctionImpl<T,NDIM>::truncate_tol(double tol, const keyT& key) const {
629
630 // RJH ... introduced max level here to avoid runaway
631 // refinement due to truncation threshold going down to
632 // intrinsic numerical error
633 const int MAXLEVEL1 = 20; // 0.5**20 ~= 1e-6
634 const int MAXLEVEL2 = 10; // 0.25**10 ~= 1e-6
635
636 if (truncate_mode == 0) {
637 return tol;
638 }
639 else if (truncate_mode == 1) {
641 return tol*std::min(1.0,pow(0.5,double(std::min(key.level(),MAXLEVEL1)))*L);
642 }
643 else if (truncate_mode == 2) {
645 return tol*std::min(1.0,pow(0.25,double(std::min(key.level(),MAXLEVEL2)))*L*L);
646 }
647 else if (truncate_mode == 3) {
648 // similar to truncate mode 1, but with an additional factor to
649 // account for an increased number of boxes in higher dimensions
650
651 // here is our handwaving argument: this threshold will give each
652 // FunctionNode an error of less than tol. The total error can
653 // then be as high as sqrt(#nodes) * tol. Therefore in order to
654 // account for higher dimensions: divide tol by about the root of
655 // number of siblings (2^NDIM) that have a large error when we
656 // refine along a deep branch of the tree. FAB
657 //
658 // Nope ... it can easily be as high as #nodes * tol. The real
659 // fix for this is an end-to-end error analysis of the larger
660 // application and if desired to include this factor into the
661 // threshold selected by the application. RJH
662 const static double fac=1.0/std::pow(2,NDIM*0.5);
663 tol*=fac;
664
666 return tol*std::min(1.0,pow(0.5,double(std::min(key.level(),MAXLEVEL1)))*L);
667
668 } else {
669 MADNESS_EXCEPTION("truncate_mode invalid",truncate_mode);
670 }
671 }
672
673 /// Returns patch referring to coeffs of child in parent box
674 template <typename T, std::size_t NDIM>
675 std::vector<Slice> FunctionImpl<T,NDIM>::child_patch(const keyT& child) const {
676 std::vector<Slice> s(NDIM);
677 const Vector<Translation,NDIM>& l = child.translation();
678 for (std::size_t i=0; i<NDIM; ++i)
679 s[i] = cdata.s[l[i]&1]; // Lowest bit of translation
680 return s;
681 }
682
683 /// Directly project parent NS coeffs to child NS coeffs
684
685 template <typename T, std::size_t NDIM>
687 const keyT& child, const keyT& parent, const coeffT& coeff) const {
688
689 const implT* f=this;
690 // MADNESS_ASSERT(coeff.tensor_type()==TT_FULL);
691 coeffT result(f->cdata.v2k,coeff.tensor_type());
692
693 // if the node for child is existent in f, and it is an internal node, we
694 // automatically have the NS form; if it is a leaf node, we only have the
695 // sum coeffs, so we take zero difference coeffs
696 if (child==parent) {
697 if (coeff.dim(0)==2*f->get_k()) result=coeff; // internal node
698 else if (coeff.dim(0)==f->get_k()) { // leaf node
699 result(f->cdata.s0)+=coeff;
700 } else {
701 MADNESS_EXCEPTION("confused k in parent_to_child_NS",1);
702 }
703 } else if (child.level()>parent.level()) {
704
705 // parent and coeff should refer to a leaf node with sum coeffs only
706 // b/c tree should be compressed with leaves kept.
707 MADNESS_ASSERT(coeff.dim(0)==f->get_k());
708 const coeffT scoeff=f->parent_to_child(coeff,parent,child);
709 result(f->cdata.s0)+=scoeff;
710 } else {
711 MADNESS_EXCEPTION("confused keys in parent_to_child_NS",1);
712 }
713 return result;
714 }
715
716 /// truncate tree at a certain level
717 template <typename T, std::size_t NDIM>
718 void FunctionImpl<T,NDIM>::erase(const Level& max_level) {
719 this->make_redundant(true);
720
721 typename dcT::iterator end = coeffs.end();
722 for (typename dcT::iterator it= coeffs.begin(); it!=end; ++it) {
723 keyT key=it->first;
724 nodeT& node=it->second;
725 if (key.level()>max_level) coeffs.erase(key);
726 if (key.level()==max_level) node.set_has_children(false);
727 }
728 this->undo_redundant(true);
729 }
730
731
732 /// Returns some asymmetry measure ... no comms
733 template <typename T, std::size_t NDIM>
737 return world.taskq.reduce<double,rangeT,do_check_symmetry_local>(rangeT(coeffs.begin(),coeffs.end()),
739 }
740
741
742 /// Refine multiple functions down to the same finest level
743
744 /// @param[v] is the vector of functions we are refining.
745 /// @param[key] is the current node.
746 /// @param[c] is the vector of coefficients passed from above.
747 template <typename T, std::size_t NDIM>
749 const std::vector<tensorT>& c,
750 const keyT key) {
751 if (key == cdata.key0 && coeffs.owner(key)!=world.rank()) return;
752
753 // First insert coefficients from above ... also get write accessors here
754 std::unique_ptr<typename dcT::accessor[]> acc(new typename dcT::accessor[v.size()]);
755 for (unsigned int i=0; i<c.size(); i++) {
756 MADNESS_ASSERT(v[i]->coeffs.get_pmap() == coeffs.get_pmap());
757 MADNESS_ASSERT(v[i]->coeffs.owner(key) == world.rank());
758 bool exists = ! v[i]->coeffs.insert(acc[i],key);
759 if (c[i].size()) {
760 MADNESS_CHECK(!exists);
761 acc[i]->second = nodeT(coeffT(c[i],targs),false);
762 }
763 else {
764 MADNESS_ASSERT(exists);
765 }
766 }
767
768 // If everyone has coefficients we are done
769 bool done = true;
770 for (unsigned int i=0; i<v.size(); i++) {
771 done &= acc[i]->second.has_coeff();
772 }
773
774 if (!done) {
775 // Those functions with coefficients need to be refined down
776 std::vector<tensorT> d(v.size());
777 for (unsigned int i=0; i<v.size(); i++) {
778 if (acc[i]->second.has_coeff()) {
779 tensorT s(cdata.v2k);
780 // s(cdata.s0) = acc[i]->second.coeff()(___);
781 s(cdata.s0) = acc[i]->second.coeff().full_tensor_copy();
782 acc[i]->second.clear_coeff();
783 d[i] = unfilter(s);
784 acc[i]->second.set_has_children(true);
785 }
786 }
787
788 // Loop thru children and pass down
789 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
790 const keyT& child = kit.key();
791 std::vector<Slice> cp = child_patch(child);
792 std::vector<tensorT> childc(v.size());
793 for (unsigned int i=0; i<v.size(); i++) {
794 if (d[i].size()) childc[i] = copy(d[i](cp));
795 }
796 woT::task(coeffs.owner(child), &implT::refine_to_common_level, v, childc, child);
797 }
798 }
799 }
800
801 // horrifically non-scalable
802 template <typename T, std::size_t NDIM>
803 void FunctionImpl<T,NDIM>::put_in_box(ProcessID from, long nl, long ni) const {
804 if (world.size()> 1000)
805 throw "NO!";
806 box_leaf[from] = nl;
807 box_interior[from] = ni;
808 }
809
810 /// Prints summary of data distribution
811 template <typename T, std::size_t NDIM>
813 if (world.size() >= 1000)
814 return;
815 for (int i=0; i<world.size(); ++i)
816 box_leaf[i] = box_interior[i] == 0;
817 world.gop.fence();
818 long nleaf=0, ninterior=0;
819 typename dcT::const_iterator end = coeffs.end();
820 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
821 const nodeT& node = it->second;
822 if (node.is_leaf())
823 ++nleaf;
824 else
825 ++ninterior;
826 }
827 this->send(0, &implT::put_in_box, world.rank(), nleaf, ninterior);
828 world.gop.fence();
829 if (world.rank() == 0) {
830 for (int i=0; i<world.size(); ++i) {
831 printf("load: %5d %8ld %8ld\n", i, box_leaf[i], box_interior[i]);
832 }
833 }
834 world.gop.fence();
835 }
836
837 template <typename T, std::size_t NDIM>
838 bool FunctionImpl<T,NDIM>::noautorefine(const keyT& key, const tensorT& t) const {
839 return false;
840 }
841
842 /// Returns true if this block of coeffs needs autorefining
843 template <typename T, std::size_t NDIM>
844 bool FunctionImpl<T,NDIM>::autorefine_square_test(const keyT& key, const nodeT& t) const {
845 double lo, hi;
846 tnorm(t.coeff().full_tensor_copy(), &lo, &hi);
847 double test = 2*lo*hi + hi*hi;
848 //print("autoreftest",key,thresh,truncate_tol(thresh, key),lo,hi,test);
849 return test> truncate_tol(thresh, key);
850 }
851
852
853 /// is this the same as trickle_down() ?
854 template <typename T, std::size_t NDIM>
856 typename dcT::accessor acc;
857 coeffs.insert(acc,key);
858 nodeT& node = acc->second;
859 coeffT& c = node.coeff();
860
861 //print(key,"received",s.normf(),c.normf(),node.has_children());
862
863 if (s.size() > 0) {
864 if (c.size() > 0)
865 c.gaxpy(1.0,s,1.0);
866 else
867 c = s;
868 }
869
870 if (node.has_children()) {
871 coeffT d;
872 if (c.has_data()) {
873 d = coeffT(cdata.v2k,targs);
874 d(cdata.s0) += c;
875 d = unfilter(d);
876 node.clear_coeff();
877 }
878 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
879 coeffT ss;
880 const keyT& child = kit.key();
881 if (d.size() > 0) ss = copy(d(child_patch(child)));
882 //print(key,"sending",ss.normf(),"to",child);
883 woT::task(coeffs.owner(child), &implT::sum_down_spawn, child, ss);
884 }
885 }
886 else {
887 // Missing coeffs assumed to be zero
888 if (c.size() <= 0) c = coeffT(cdata.vk,targs);
889 }
890 }
891
892 /// After 1d push operator must sum coeffs down the tree to restore correct scaling function coefficients
893 template <typename T, std::size_t NDIM>
895 if (world.rank() == coeffs.owner(cdata.key0)) sum_down_spawn(cdata.key0, coeffT());
896
897 if (fence) world.gop.fence();
898 }
899
900
901 template <typename T, std::size_t NDIM>
903 const implT* f,
904 const keyT& key,
905 const std::pair<keyT,coeffT>& left,
906 const std::pair<keyT,coeffT>& center,
907 const std::pair<keyT,coeffT>& right) {
908 D->forward_do_diff1(f,this,key,left,center,right);
909 }
910
911
912 template <typename T, std::size_t NDIM>
914 const implT* f,
915 const keyT& key,
916 const std::pair<keyT,coeffT>& left,
917 const std::pair<keyT,coeffT>& center,
918 const std::pair<keyT,coeffT>& right) {
919 D->do_diff1(f,this,key,left,center,right);
920 }
921
922
923 // Called by result function to differentiate f
924 template <typename T, std::size_t NDIM>
925 void FunctionImpl<T,NDIM>::diff(const DerivativeBase<T,NDIM>* D, const implT* f, bool fence) {
926 typedef std::pair<keyT,coeffT> argT;
927 for (const auto& [key, node]: f->coeffs) {
928 if (node.has_coeff()) {
929 Future<argT> left = D->find_neighbor(f, key,-1);
930 argT center(key,node.coeff());
931 Future<argT> right = D->find_neighbor(f, key, 1);
932 world.taskq.add(*this, &implT::do_diff1, D, f, key, left, center, right, TaskAttributes::hipri());
933 }
934 else {
935 coeffs.replace(key,nodeT(coeffT(),true)); // Empty internal node
936 }
937 }
938 if (fence) world.gop.fence();
939 }
940
941
942 /// return the a std::pair<key, node>, which MUST exist
943 template <typename T, std::size_t NDIM>
945 MADNESS_ASSERT(coeffs.probe(key));
946 ShallowNode<T,NDIM> snode(coeffs.find(key).get()->second);
947 return std::pair<Key<NDIM>,ShallowNode<T,NDIM> >(key,snode);
948 }
949
950 /// multiply the ket with a one-electron potential rr(1,2)= f(1,2)*g(1)
951
952 /// @param[in] val_ket function values of f(1,2)
953 /// @param[in] val_pot function values of g(1)
954 /// @param[in] particle if 0 then g(1), if 1 then g(2)
955 /// @return the resulting function values
956 template <typename T, std::size_t NDIM>
958 const coeffT& val_pot, int particle) const {
959
962 MADNESS_ASSERT(val_ket.is_svd_tensor());
963
964 std::vector<long> vkhalf=std::vector<long>(NDIM/2,cdata.vk[0]);
965 tensorT ones=tensorT(vkhalf);
966 ones=1.0;
967
968 TensorArgs targs(-1.0,val_ket.tensor_type());
969 coeffT pot12;
970 if (particle==0) pot12=outer(val_pot.full_tensor(),ones,targs);
971 else if (particle==1) pot12=outer(ones,val_pot.full_tensor(),targs);
972
973 coeffT result=copy(val_ket);
974 result.emul(pot12);
975
976 return result;
977 }
978
979
980 /// given several coefficient tensors, assemble a result tensor
981
982 /// the result looks like: (v(1,2) + v(1) + v(2)) |ket(1,2)>
983 /// or (v(1,2) + v(1) + v(2)) |p(1) p(2)>
984 /// i.e. coefficients for the ket and coefficients for the two particles are
985 /// mutually exclusive. All potential terms are optional, just pass in empty coeffs.
986 /// @param[in] key the key of the FunctionNode to which these coeffs belong
987 /// @param[in] cket coefficients of the ket
988 /// @param[in] vpotential1 function values of the potential for particle 1
989 /// @param[in] vpotential2 function values of the potential for particle 2
990 /// @param[in] veri function values for the 2-particle potential
991 template <typename T, std::size_t NDIM>
993 const keyT& key, const coeffT& coeff_ket, const coeffT& vpotential1,
994 const coeffT& vpotential2, const tensorT& veri) const {
995
996 // take a shortcut if we are already done
997 bool ket_only=(not (vpotential1.has_data() or vpotential2.has_data() or veri.has_data()));
998 if (ket_only) return coeff_ket;
999
1000 // switch to values instead of coefficients
1001 coeffT val_ket=coeffs2values(key,coeff_ket);
1002
1003 // the result tensor
1004 coeffT val_result=coeffT(val_ket.ndim(),val_ket.dims(),this->get_tensor_args().tt);
1005 coeffT coeff_result;
1006
1007 // potential for particles 1 and 2, must be done in TT_2D
1008 if (vpotential1.has_data() or vpotential2.has_data()) {
1009 val_ket=val_ket.convert(TensorArgs(-1.0,TT_2D));
1010 }
1011 if (vpotential1.has_data()) val_result+=multiply(val_ket,vpotential1,0);
1012 if (vpotential2.has_data()) val_result+=multiply(val_ket,vpotential2,1);
1013
1014 // values for eri: this must be done in full rank...
1015 if (veri.has_data()) {
1016 tensorT val_ket2=val_ket.full_tensor_copy().emul(veri);
1017 if (val_result.has_data()) val_ket2+=val_result.full_tensor_copy();
1018 // values2coeffs expensive (30%), coeffT() (relatively) cheap (8%)
1019 coeff_result=coeffT(values2coeffs(key,val_ket2),this->get_tensor_args());
1020
1021 } else {
1022
1023 // convert back to original tensor type
1024 val_ket=val_ket.convert(get_tensor_args());
1025 MADNESS_ASSERT(val_result.has_data());
1026 coeff_result=values2coeffs(key,val_result);
1027 coeff_result.reduce_rank(this->get_tensor_args().thresh);
1028 }
1029
1030 return coeff_result;
1031
1032 }
1033
1034 /// Permute the dimensions of f according to map, result on this
1035 template <typename T, std::size_t NDIM>
1036 void FunctionImpl<T,NDIM>::mapdim(const implT& f, const std::vector<long>& map, bool fence) {
1037
1039 const_cast<implT*>(&f)->flo_unary_op_node_inplace(do_mapdim(map,*this),fence);
1040
1041 }
1042
1043 /// mirror the dimensions of f according to mirror, result on this
1044 template <typename T, std::size_t NDIM>
1045 void FunctionImpl<T,NDIM>::mirror(const implT& f, const std::vector<long>& mirrormap, bool fence) {
1047 const_cast<implT*>(&f)->flo_unary_op_node_inplace(do_mirror(mirrormap,*this),fence);
1048 }
1049
1050 /// map and mirror the translation index and the coefficients, result on this
1051
1052 /// first map the dimensions, the mirror!
1053 /// this = mirror(map(f))
1054 template <typename T, std::size_t NDIM>
1055 void FunctionImpl<T,NDIM>::map_and_mirror(const implT& f, const std::vector<long>& map,
1056 const std::vector<long>& mirror, bool fence) {
1058 const_cast<implT*>(&f)->flo_unary_op_node_inplace(do_map_and_mirror(map,mirror,*this),fence);
1059 }
1060
1061
1062
1063 /// take the average of two functions, similar to: this=0.5*(this+rhs)
1064
1065 /// works in either basis and also in nonstandard form
1066 template <typename T, std::size_t NDIM>
1068
1069 rhs.flo_unary_op_node_inplace(do_average(*this),true);
1070 this->scale_inplace(0.5,true);
1071 flo_unary_op_node_inplace(do_reduce_rank(targs),true);
1072 }
1073
1074 /// change the tensor type of the coefficients in the FunctionNode
1075
1076 /// @param[in] targs target tensor arguments (threshold and full/low rank)
1077 template <typename T, std::size_t NDIM>
1079 flo_unary_op_node_inplace(do_change_tensor_type(targs,*this),fence);
1080 }
1081
1082 /// reduce the rank of the coefficients tensors
1083
1084 /// @param[in] targs target tensor arguments (threshold and full/low rank)
1085 template <typename T, std::size_t NDIM>
1086 void FunctionImpl<T,NDIM>::reduce_rank(const double thresh, bool fence) {
1087 flo_unary_op_node_inplace(do_reduce_rank(thresh),fence);
1088 }
1089
1090 /// reduce the rank of the coefficients tensors
1091
1092 /// @param[in] targs target tensor arguments (threshold and full/low rank)
1093 template <typename T, std::size_t NDIM>
1094 void FunctionImpl<T,NDIM>::chop_at_level(const int n, bool fence) {
1095 std::list<keyT> to_be_erased;
1096 for (auto it=coeffs.begin(); it!=coeffs.end(); ++it) {
1097 const keyT& key=it->first;
1098 nodeT& node=it->second;
1099 if (key.level()==n) node.set_is_leaf(true);
1100 if (key.level()>n) to_be_erased.push_back(key);
1101 }
1102 for (auto& key : to_be_erased) coeffs.erase(key);
1103 }
1104
1105
1106/// compute norm of s and d coefficients for all nodes
1107
1108 /// @param[in] targs target tensor arguments (threshold and full/low rank)
1109 template <typename T, std::size_t NDIM>
1111 //const auto& data=FunctionCommonData<T,NDIM>::get(get_k());
1112 flo_unary_op_node_inplace(
1113 do_compute_snorm_and_dnorm(cdata),fence);
1114 }
1115
1116
1117/// Transform sum coefficients at level n to sums+differences at level n-1
1118
1119/// Given scaling function coefficients s[n][l][i] and s[n][l+1][i]
1120/// return the scaling function and wavelet coefficients at the
1121/// coarser level. I.e., decompose Vn using Vn = Vn-1 + Wn-1.
1122 /// \code
1123 /// s_i = sum(j) h0_ij*s0_j + h1_ij*s1_j
1124 /// d_i = sum(j) g0_ij*s0_j + g1_ij*s1_j
1125 // \endcode
1126 /// Returns a new tensor and has no side effects. Works for any
1127 /// number of dimensions.
1128 ///
1129 /// No communication involved.
1130 template <typename T, std::size_t NDIM>
1132 tensorT r(cdata.v2k,false);
1133 tensorT w(cdata.v2k,false);
1134 return fast_transform(s,cdata.hgT,r,w);
1135 //return transform(s,cdata.hgT);
1136 }
1137
1138 template <typename T, std::size_t NDIM>
1140 coeffT result=transform(s,cdata.hgT);
1141 return result;
1142 }
1143
1144 /// Transform sums+differences at level n to sum coefficients at level n+1
1145
1146 /// Given scaling function and wavelet coefficients (s and d)
1147 /// returns the scaling function coefficients at the next finer
1148 /// level. I.e., reconstruct Vn using Vn = Vn-1 + Wn-1.
1149 /// \code
1150 /// s0 = sum(j) h0_ji*s_j + g0_ji*d_j
1151 /// s1 = sum(j) h1_ji*s_j + g1_ji*d_j
1152 /// \endcode
1153 /// Returns a new tensor and has no side effects
1154 ///
1155 /// If (sonly) ... then ss is only the scaling function coeff (and
1156 /// assume the d are zero). Works for any number of dimensions.
1157 ///
1158 /// No communication involved.
1159 template <typename T, std::size_t NDIM>
1161 tensorT r(cdata.v2k,false);
1162 tensorT w(cdata.v2k,false);
1163 return fast_transform(s,cdata.hg,r,w);
1164 //return transform(s, cdata.hg);
1165 }
1166
1167 template <typename T, std::size_t NDIM>
1169 return transform(s,cdata.hg);
1170 }
1171
1172 /// downsample the sum coefficients of level n+1 to sum coeffs on level n
1173
1174 /// specialization of the filter method, will yield only the sum coefficients
1175 /// @param[in] key key of level n
1176 /// @param[in] v vector of sum coefficients of level n+1
1177 /// @param[in] args TensorArguments for possible low rank approximations
1178 /// @return sum coefficients on level n in full tensor format
1179 template <typename T, std::size_t NDIM>
1181
1182 tensorT result(cdata.vk);
1183
1184 // the twoscale coefficients: for downsampling use h0/h1; see Alpert Eq (3.34a)
1185 const tensorT h[2] = {cdata.h0T, cdata.h1T};
1186 tensorT matrices[NDIM];
1187
1188 // loop over all child nodes, transform and accumulate
1189 long i=0;
1190 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1191
1192 // get the appropriate twoscale coefficients for each dimension
1193 for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=h[kit.key().translation()[ii]%2];
1194
1195 // transform and accumulate on the result
1196 result+=general_transform(v[i].get(),matrices).full_tensor_copy();
1197
1198 }
1199 return result;
1200 }
1201
1202 /// upsample the sum coefficients of level 1 to sum coeffs on level n+1
1203
1204 /// specialization of the unfilter method, will transform only the sum coefficients
1205 /// @param[in] key key of level n+1
1206 /// @param[in] coeff sum coefficients of level n (does NOT belong to key!!)
1207 /// @param[in] args TensorArguments for possible low rank approximations
1208 /// @return sum coefficients on level n+1
1209 template <typename T, std::size_t NDIM>
1211
1212 // the twoscale coefficients: for upsampling use h0/h1; see Alpert Eq (3.35a/b)
1213 // note there are no difference coefficients; if you want that use unfilter
1214 const tensorT h[2] = {cdata.h0, cdata.h1};
1215 tensorT matrices[NDIM];
1216
1217 // get the appropriate twoscale coefficients for each dimension
1218 for (size_t ii=0; ii<NDIM; ++ii) matrices[ii]=h[key.translation()[ii]%2];
1219
1220 // transform and accumulate on the result
1221 const coeffT result=general_transform(coeff,matrices);
1222 return result;
1223 }
1224
1225
1226 /// Projects old function into new basis (only in reconstructed form)
1227 template <typename T, std::size_t NDIM>
1228 void FunctionImpl<T,NDIM>::project(const implT& old, bool fence) {
1229 long kmin = std::min(cdata.k,old.cdata.k);
1230 std::vector<Slice> s(NDIM,Slice(0,kmin-1));
1231 typename dcT::const_iterator end = old.coeffs.end();
1232 for (typename dcT::const_iterator it=old.coeffs.begin(); it!=end; ++it) {
1233 const keyT& key = it->first;
1234 const nodeT& node = it->second;
1235 if (node.has_coeff()) {
1236 coeffT c(cdata.vk,targs);
1237 c(s) += node.coeff()(s);
1238 coeffs.replace(key,nodeT(c,false));
1239 }
1240 else {
1241 coeffs.replace(key,nodeT(coeffT(),true));
1242 }
1243 }
1244 if (fence)
1245 world.gop.fence();
1246 }
1247
1248 template <typename T, std::size_t NDIM>
1250 return coeffs.probe(key) && coeffs.find(key).get()->second.has_children();
1251 }
1253 template <typename T, std::size_t NDIM>
1255 return coeffs.probe(key) && (not coeffs.find(key).get()->second.has_children());
1256 }
1258
1259 template <typename T, std::size_t NDIM>
1260 void FunctionImpl<T,NDIM>::broaden_op(const keyT& key, const std::vector< Future <bool> >& v) {
1261 for (unsigned int i=0; i<v.size(); ++i) {
1262 if (v[i]) {
1263 refine_op(true_refine_test(), key);
1264 break;
1265 }
1266 }
1267 }
1268
1269 // For each local node sets value of norm tree, snorm and dnorm to 0.0
1270 template <typename T, std::size_t NDIM>
1272 typename dcT::iterator end = coeffs.end();
1273 for (typename dcT::iterator it=coeffs.begin(); it!=end; ++it) {
1274 it->second.set_norm_tree(0.0);
1275 it->second.set_snorm(0.0);
1276 it->second.set_dnorm(0.0);
1278 }
1279
1280 // Broaden tree
1281 template <typename T, std::size_t NDIM>
1282 void FunctionImpl<T,NDIM>::broaden(std::vector<bool> is_periodic, bool fence) {
1283 typename dcT::iterator end = coeffs.end();
1284 for (typename dcT::iterator it=coeffs.begin(); it!=end; ++it) {
1285 const keyT& key = it->first;
1286 typename dcT::accessor acc;
1287 const auto found = coeffs.find(acc,key);
1288 MADNESS_CHECK(found);
1289 nodeT& node = acc->second;
1290 if (node.has_coeff() &&
1291 node.get_norm_tree() != -1.0 &&
1292 node.coeff().normf() >= truncate_tol(thresh,key)) {
1294 node.set_norm_tree(-1.0); // Indicates already broadened or result of broadening/refining
1296 //int ndir = std::pow(3,NDIM);
1297 int ndir = static_cast<int>(std::pow(static_cast<double>(3), static_cast<int>(NDIM)));
1298 std::vector< Future <bool> > v = future_vector_factory<bool>(ndir);
1299 keyT neigh;
1300 int i=0;
1301 for (HighDimIndexIterator it(NDIM,3); it; ++it) {
1303 for (std::size_t d=0; d<NDIM; ++d) {
1304 const int odd = key.translation()[d] & 0x1L; // 1 if odd, 0 if even
1305 l[d] -= 1; // (0,1,2) --> (-1,0,1)
1306 if (l[d] == -1)
1307 l[d] = -1-odd;
1308 else if (l[d] == 1)
1309 l[d] = 2 - odd;
1310 }
1311 keyT neigh = neighbor(key, keyT(key.level(),l), is_periodic);
1312
1313 if (neigh.is_valid()) {
1314 v[i++] = this->task(coeffs.owner(neigh), &implT::exists_and_has_children, neigh);
1315 }
1316 else {
1317 v[i++].set(false);
1319 }
1320 woT::task(world.rank(), &implT::broaden_op, key, v);
1321 }
1322 }
1323 // Reset value of norm tree so that can repeat broadening
1324 if (fence) {
1325 world.gop.fence();
1326 zero_norm_tree();
1327 world.gop.fence();
1328 }
1329 }
1330
1331 /// sum all the contributions from all scales after applying an operator in mod-NS form
1332 template <typename T, std::size_t NDIM>
1334 set_tree_state(reconstructed);
1335 if (world.rank() == coeffs.owner(cdata.key0))
1336 woT::task(world.rank(), &implT::trickle_down_op, cdata.key0,coeffT());
1337 if (fence) world.gop.fence();
1338 }
1339
1340 /// sum all the contributions from all scales after applying an operator in mod-NS form
1341
1342 /// cf reconstruct_op
1343 template <typename T, std::size_t NDIM>
1345 // Note that after application of an integral operator not all
1346 // siblings may be present so it is necessary to check existence
1347 // and if absent insert an empty leaf node.
1348 //
1349 // If summing the result of an integral operator (i.e., from
1350 // non-standard form) there will be significant scaling function
1351 // coefficients at all levels and possibly difference coefficients
1352 // in leaves, hence the tree may refine as a result.
1353 typename dcT::iterator it = coeffs.find(key).get();
1354 if (it == coeffs.end()) {
1355 coeffs.replace(key,nodeT(coeffT(),false));
1356 it = coeffs.find(key).get();
1357 }
1358 nodeT& node = it->second;
1359
1360 // The integral operator will correctly connect interior nodes
1361 // to children but may leave interior nodes without coefficients
1362 // ... but they still need to sum down so just give them zeros
1363 if (node.coeff().has_no_data()) node.coeff()=coeffT(cdata.vk,targs);
1365 // if (node.has_children() || node.has_coeff()) { // Must allow for inconsistent state from transform, etc.
1366 if (node.has_children()) { // Must allow for inconsistent state from transform, etc.
1367 coeffT d = node.coeff();
1368 if (key.level() > 0) d += s; // -- note accumulate for NS summation
1369 node.clear_coeff();
1370 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
1371 const keyT& child = kit.key();
1372 coeffT ss= upsample(child,d);
1373 ss.reduce_rank(thresh);
1374 PROFILE_BLOCK(recon_send);
1375 woT::task(coeffs.owner(child), &implT::trickle_down_op, child, ss);
1376 }
1378 else {
1379 node.coeff()+=s;
1380 node.coeff().reduce_rank(thresh);
1381 }
1382 }
1383
1384 /// change the tree state of this function, might or might not respect fence!
1385 template <typename T, std::size_t NDIM>
1386 void FunctionImpl<T,NDIM>::change_tree_state(const TreeState finalstate, bool fence) {
1387
1388 TreeState current_state=get_tree_state();
1389 if (current_state==finalstate) return;
1390
1391 // very special case
1392 if (get_tree_state()==nonstandard_after_apply) {
1393 MADNESS_CHECK(finalstate==reconstructed);
1394 reconstruct(fence);
1395 return;
1396 }
1397 MADNESS_CHECK_THROW(current_state!=TreeState::nonstandard_after_apply,"unknown tree state");
1398 bool must_fence=false;
1399
1400 if (finalstate==reconstructed) {
1401 if (current_state==reconstructed) return;
1402 if (current_state==compressed) reconstruct(fence);
1403 if (current_state==nonstandard) reconstruct(fence);
1404 if (current_state==nonstandard_with_leaves) remove_internal_coefficients(fence);
1405 if (current_state==redundant) remove_internal_coefficients(fence);
1406 set_tree_state(reconstructed);
1407 } else if (finalstate==compressed) {
1408 if (current_state==reconstructed) compress(compressed,fence);
1409 if (current_state==compressed) return;
1410 if (current_state==nonstandard) standard(fence);
1411 if (current_state==nonstandard_with_leaves) standard(fence);
1412 if (current_state==redundant) {
1413 remove_internal_coefficients(true);
1414 must_fence=true;
1415 set_tree_state(reconstructed);
1416 compress(compressed,fence);
1417 }
1418 set_tree_state(compressed);
1419 } else if (finalstate==nonstandard) {
1420 if (current_state==reconstructed) compress(nonstandard,fence);
1421 if (current_state==compressed) {
1422 reconstruct(true);
1423 must_fence=true;
1424 compress(nonstandard,fence);
1425 }
1426 if (current_state==nonstandard) return;
1427 if (current_state==nonstandard_with_leaves) remove_leaf_coefficients(fence);
1428 if (current_state==redundant) {
1429 remove_internal_coefficients(true);
1430 must_fence=true;
1431 set_tree_state(reconstructed);
1432 compress(nonstandard,fence);
1434 set_tree_state(nonstandard);
1435 } else if (finalstate==nonstandard_with_leaves) {
1437 if (current_state==compressed) {
1438 reconstruct(true);
1439 must_fence=true;
1441 }
1442 if (current_state==nonstandard) {
1444 must_fence=true;
1445 reconstruct(true);
1447 }
1448 if (current_state==nonstandard_with_leaves) return;
1449 if (current_state==redundant) {
1450 remove_internal_coefficients(true);
1451 must_fence=true;
1452 set_tree_state(reconstructed);
1454 }
1455 set_tree_state(nonstandard_with_leaves);
1456 } else if (finalstate==redundant) {
1457 if (current_state==reconstructed) make_redundant(fence);
1458 if (current_state==compressed) {
1459 reconstruct(true);
1460 must_fence=true;
1461 make_redundant(fence);
1462 }
1463 if (current_state==nonstandard) {
1464 standard(true);
1465 must_fence=true;
1466 reconstruct(true);
1467 make_redundant(fence);
1468 }
1469 if (current_state==nonstandard_with_leaves) {
1470 remove_internal_coefficients(true);
1471 must_fence=true;
1472 set_tree_state(reconstructed);
1473 make_redundant(fence);
1474 }
1475 if (current_state==redundant) return;
1476 set_tree_state(redundant);
1477 } else {
1478 MADNESS_EXCEPTION("unknown/unsupported final tree state",1);
1479 }
1480 if (must_fence and world.rank()==0) {
1481 print("could not respect fence in change_tree_state");
1482 }
1483 if (fence && VERIFY_TREE) verify_tree(); // Must be after in case nonstandard
1484 return;
1485 }
1486
1487
1488
1489 template <typename T, std::size_t NDIM>
1491
1492 if (is_reconstructed()) return;
1493
1494 if (is_redundant() or is_nonstandard_with_leaves()) {
1495 set_tree_state(reconstructed);
1496 this->remove_internal_coefficients(fence);
1497 } else if (is_compressed() or tree_state==nonstandard_after_apply) {
1498 // Must set true here so that successive calls without fence do the right thing
1499 set_tree_state(reconstructed);
1500 if (world.rank() == coeffs.owner(cdata.key0))
1501 woT::task(world.rank(), &implT::reconstruct_op, cdata.key0,coeffT(), true);
1502 } else if (is_nonstandard()) {
1503 // Must set true here so that successive calls without fence do the right thing
1504 set_tree_state(reconstructed);
1505 if (world.rank() == coeffs.owner(cdata.key0))
1506 woT::task(world.rank(), &implT::reconstruct_op, cdata.key0,coeffT(), false);
1507 } else {
1508 MADNESS_EXCEPTION("cannot reconstruct this tree",1);
1509 }
1510 if (fence) world.gop.fence();
1511
1512 }
1513
1514 /// compress the wave function
1515
1516 /// after application there will be sum coefficients at the root level,
1517 /// and difference coefficients at all other levels; furthermore:
1518 /// @param[in] nonstandard keep sum coeffs at all other levels, except leaves
1519 /// @param[in] keepleaves keep sum coeffs (but no diff coeffs) at leaves
1520 /// @param[in] redundant keep only sum coeffs at all levels, discard difference coeffs
1521 template <typename T, std::size_t NDIM>
1522 void FunctionImpl<T,NDIM>::compress(const TreeState newstate, bool fence) {
1523 MADNESS_CHECK_THROW(is_reconstructed(),"impl::compress wants a reconstructe tree");
1524 // Must set true here so that successive calls without fence do the right thing
1525 set_tree_state(newstate);
1526 bool keepleaves1=(tree_state==nonstandard_with_leaves) or (tree_state==redundant);
1527 bool nonstandard1=(tree_state==nonstandard) or (tree_state==nonstandard_with_leaves);
1528 bool redundant1=(tree_state==redundant);
1529
1530 if (world.rank() == coeffs.owner(cdata.key0)) {
1531
1532 compress_spawn(cdata.key0, nonstandard1, keepleaves1, redundant1);
1533 }
1534 if (fence)
1535 world.gop.fence();
1536 }
1537
1538 template <typename T, std::size_t NDIM>
1540 flo_unary_op_node_inplace(remove_internal_coeffs(),fence);
1541 }
1542
1543 template <typename T, std::size_t NDIM>
1545 flo_unary_op_node_inplace(remove_leaf_coeffs(),fence);
1546 }
1547
1548 /// convert this to redundant, i.e. have sum coefficients on all levels
1549 template <typename T, std::size_t NDIM>
1551
1552 // fast return if possible
1553 if (is_redundant()) return;
1554 MADNESS_CHECK_THROW(is_reconstructed(),"impl::make_redundant() wants a reconstructed tree");
1555 compress(redundant,fence);
1556 }
1557
1558 /// convert this from redundant to standard reconstructed form
1559 template <typename T, std::size_t NDIM>
1561 MADNESS_CHECK_THROW(is_redundant(),"impl::undo_redundant() wants a redundant tree");
1562 set_tree_state(reconstructed);
1563 flo_unary_op_node_inplace(remove_internal_coeffs(),fence);
1564 }
1565
1566
1567 /// compute for each FunctionNode the norm of the function inside that node
1568 template <typename T, std::size_t NDIM>
1570 if (world.rank() == coeffs.owner(cdata.key0))
1571 norm_tree_spawn(cdata.key0);
1572 if (fence)
1573 world.gop.fence();
1574 }
1575
1576 template <typename T, std::size_t NDIM>
1577 double FunctionImpl<T,NDIM>::norm_tree_op(const keyT& key, const std::vector< Future<double> >& v) {
1578 //PROFILE_MEMBER_FUNC(FunctionImpl);
1579 double sum = 0.0;
1580 int i=0;
1581 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1582 double value = v[i].get();
1583 sum += value*value;
1584 }
1585 sum = sqrt(sum);
1586 coeffs.task(key, &nodeT::set_norm_tree, sum); // why a task? because send is deprecated to keep comm thread free
1587 //if (key.level() == 0) std::cout << "NORM_TREE_TOP " << sum << "\n";
1588 return sum;
1589 }
1590
1591 template <typename T, std::size_t NDIM>
1593 nodeT& node = coeffs.find(key).get()->second;
1594 if (node.has_children()) {
1595 std::vector< Future<double> > v = future_vector_factory<double>(1<<NDIM);
1596 int i=0;
1597 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1598 v[i] = woT::task(coeffs.owner(kit.key()), &implT::norm_tree_spawn, kit.key());
1599 }
1600 return woT::task(world.rank(),&implT::norm_tree_op, key, v);
1601 }
1602 else {
1603 // return Future<double>(node.coeff().normf());
1604 const double norm=node.coeff().normf();
1605 // invoked locally anyways
1606 node.set_norm_tree(norm);
1607 return Future<double>(norm);
1608 }
1609 }
1610
1611 /// truncate using a tree in reconstructed form
1612
1613 /// must be invoked where key is local
1614 template <typename T, std::size_t NDIM>
1616 MADNESS_ASSERT(coeffs.probe(key));
1617 nodeT& node = coeffs.find(key).get()->second;
1618
1619 // if this is a leaf node just return the sum coefficients
1620 if (not node.has_children()) return Future<coeffT>(node.coeff());
1621
1622 // if this is an internal node, wait for all the children's sum coefficients
1623 // and use them to determine if the children can be removed
1624 std::vector<Future<coeffT> > v = future_vector_factory<coeffT>(1<<NDIM);
1625 int i=0;
1626 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1627 v[i] = woT::task(coeffs.owner(kit.key()), &implT::truncate_reconstructed_spawn, kit.key(),tol,TaskAttributes::hipri());
1628 }
1629
1630 // will return (possibly empty) sum coefficients
1631 return woT::task(world.rank(),&implT::truncate_reconstructed_op,key,v,tol,TaskAttributes::hipri());
1632
1633 }
1634
1635 /// given the sum coefficients of all children, truncate or not
1636
1637 /// @return new sum coefficients (empty if internal, not empty, if new leaf); might delete its children
1638 template <typename T, std::size_t NDIM>
1639 typename FunctionImpl<T,NDIM>::coeffT FunctionImpl<T,NDIM>::truncate_reconstructed_op(const keyT& key, const std::vector< Future<coeffT > >& v, const double tol) {
1640
1641 MADNESS_ASSERT(coeffs.probe(key));
1642
1643 // the sum coefficients might be empty, which means they come from an internal node
1644 // and we must not truncate; so just return empty coeffs again
1645 for (size_t i=0; i<v.size(); ++i) if (v[i].get().has_no_data()) return coeffT();
1646
1647 // do not truncate below level 1
1648 if (key.level()<2) return coeffT();
1649
1650 // compute the wavelet coefficients from the child nodes
1651 typename dcT::accessor acc;
1652 const auto found = coeffs.find(acc, key);
1654 int i=0;
1655 tensorT d(cdata.v2k);
1656 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1657 // d(child_patch(kit.key())) += v[i].get();
1658 d(child_patch(kit.key())) += v[i].get().full_tensor_copy();
1659 }
1660
1661 d = filter(d);
1662 tensorT s=copy(d(cdata.s0));
1663 d(cdata.s0) = 0.0;
1664 const double error=d.normf();
1665
1666 nodeT& node = coeffs.find(key).get()->second;
1667
1668 if (error < truncate_tol(tol,key)) {
1669 node.set_has_children(false);
1670 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
1671 coeffs.erase(kit.key());
1672 }
1673 // "replace" children with new sum coefficients
1674 coeffT ss=coeffT(s,targs);
1675 acc->second.set_coeff(ss);
1676 return ss;
1677 } else {
1678 return coeffT();
1679 }
1680 }
1681
1682 /// calculate the wavelet coefficients using the sum coefficients of all child nodes
1683
1684 /// @param[in] key this's key
1685 /// @param[in] v sum coefficients of the child nodes
1686 /// @param[in] nonstandard keep the sum coefficients with the wavelet coefficients
1687 /// @param[in] redundant keep only the sum coefficients, discard the wavelet coefficients
1688 /// @return the sum coefficients
1689 template <typename T, std::size_t NDIM>
1690 std::pair<typename FunctionImpl<T,NDIM>::coeffT,double> FunctionImpl<T,NDIM>::compress_op(const keyT& key,
1691 const std::vector< Future<std::pair<coeffT,double> > >& v, bool nonstandard1) {
1692 //PROFILE_MEMBER_FUNC(FunctionImpl);
1693
1694 double cpu0=cpu_time();
1695 // Copy child scaling coeffs into contiguous block
1696 tensorT d(cdata.v2k);
1697 // coeffT d(cdata.v2k,targs);
1698 int i=0;
1699 double norm_tree2=0.0;
1700 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1701 // d(child_patch(kit.key())) += v[i].get();
1702 d(child_patch(kit.key())) += v[i].get().first.full_tensor_copy();
1703 norm_tree2+=v[i].get().second*v[i].get().second;
1704 }
1705
1706 d = filter(d);
1707 double cpu1=cpu_time();
1708 timer_filter.accumulate(cpu1-cpu0);
1709 cpu0=cpu1;
1711 typename dcT::accessor acc;
1712 const auto found = coeffs.find(acc, key);
1713 MADNESS_CHECK(found);
1714 MADNESS_CHECK_THROW(!acc->second.has_coeff(),"compress_op: existing coeffs where there should be none");
1715
1716 // tighter thresh for internal nodes
1717 TensorArgs targs2=targs;
1718 targs2.thresh*=0.1;
1719
1720 // need the deep copy for contiguity
1721 coeffT ss=coeffT(copy(d(cdata.s0)));
1722 double snorm=ss.normf();
1723
1724 if (key.level()> 0 && !nonstandard1) d(cdata.s0) = 0.0;
1725
1726 coeffT dd=coeffT(d,targs2);
1727 double dnorm=dd.normf();
1728 double norm_tree=sqrt(norm_tree2);
1729
1730 acc->second.set_snorm(snorm);
1731 acc->second.set_dnorm(dnorm);
1732 acc->second.set_norm_tree(norm_tree);
1733
1734 acc->second.set_coeff(dd);
1735 cpu1=cpu_time();
1736 timer_compress_svd.accumulate(cpu1-cpu0);
1737
1738 // return sum coefficients
1739 return std::make_pair(ss,snorm);
1740 }
1741
1742 /// similar to compress_op, but insert only the sum coefficients in the tree
1743
1744 /// also sets snorm, dnorm and norm_tree for all nodes
1745 /// @param[in] key this's key
1746 /// @param[in] v sum coefficients of the child nodes
1747 /// @return the sum coefficients
1748 template <typename T, std::size_t NDIM>
1749 std::pair<typename FunctionImpl<T,NDIM>::coeffT,double>
1750 FunctionImpl<T,NDIM>::make_redundant_op(const keyT& key, const std::vector< Future<std::pair<coeffT,double> > >& v) {
1751
1752 tensorT d(cdata.v2k);
1753 int i=0;
1754 double norm_tree2=0.0;
1755 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
1756 d(child_patch(kit.key())) += v[i].get().first.full_tensor_copy();
1757 norm_tree2+=v[i].get().second*v[i].get().second;
1758 }
1759 d = filter(d);
1760 double norm_tree=sqrt(norm_tree2);
1761
1762 // tighter thresh for internal nodes
1763 TensorArgs targs2=targs;
1764 targs2.thresh*=0.1;
1765
1766 // need the deep copy for contiguity
1767 coeffT s=coeffT(copy(d(cdata.s0)),targs2);
1768 d(cdata.s0)=0.0;
1769 double dnorm=d.normf();
1770 double snorm=s.normf();
1771
1772 typename dcT::accessor acc;
1773 const auto found = coeffs.find(acc, key);
1774 MADNESS_CHECK(found);
1775
1776 acc->second.set_coeff(s);
1777 acc->second.set_dnorm(dnorm);
1778 acc->second.set_snorm(snorm);
1779 acc->second.set_norm_tree(norm_tree);
1780
1781 // return sum coefficients
1782 return std::make_pair(s,norm_tree);
1783 }
1784
1785 /// Changes non-standard compressed form to standard compressed form
1786 template <typename T, std::size_t NDIM>
1788
1789 if (is_compressed()) return;
1790 set_tree_state(compressed);
1791 flo_unary_op_node_inplace(do_standard(this),fence);
1792// make_nonstandard = false;
1793 }
1794
1795
1796 /// after apply we need to do some cleanup;
1797
1798 /// forces fence
1799 template <typename T, std::size_t NDIM>
1801 bool print_timings=false;
1802 bool printme=(world.rank()==0 and print_timings);
1803 TensorArgs tight_args(targs);
1804 tight_args.thresh*=0.01;
1805 double begin=wall_time();
1806 double begin1=wall_time();
1807 flo_unary_op_node_inplace(do_consolidate_buffer(tight_args),true);
1808 double end1=wall_time();
1809 if (printme) printf("time in consolidate_buffer %8.4f\n",end1-begin1);
1810
1811
1812 // reduce the rank of the final nodes, leave full tensors unchanged
1813 // flo_unary_op_node_inplace(do_reduce_rank(tight_args.thresh),true);
1814 begin1=wall_time();
1815 flo_unary_op_node_inplace(do_reduce_rank(targs),true);
1816 end1=wall_time();
1817 if (printme) printf("time in do_reduce_rank %8.4f\n",end1-begin1);
1818
1819 // change TT_FULL to low rank
1820 begin1=wall_time();
1821 flo_unary_op_node_inplace(do_change_tensor_type(targs,*this),true);
1822 end1=wall_time();
1823 if (printme) printf("time in do_change_tensor_type %8.4f\n",end1-begin1);
1824
1825 // truncate leaf nodes to avoid excessive tree refinement
1826 begin1=wall_time();
1827 flo_unary_op_node_inplace(do_truncate_NS_leafs(this),true);
1828 end1=wall_time();
1829 if (printme) printf("time in do_truncate_NS_leafs %8.4f\n",end1-begin1);
1830
1831 double end=wall_time();
1832 double elapsed=end-begin;
1833 set_tree_state(nonstandard_after_apply);
1834 world.gop.fence();
1835 return elapsed;
1836 }
1837
1838
1839 /// after summing up we need to do some cleanup;
1840
1841 /// forces fence
1842 template <typename T, std::size_t NDIM>
1844 world.gop.fence();
1845 flo_unary_op_node_inplace(do_consolidate_buffer(get_tensor_args()), true);
1846 sum_down(true);
1847 set_tree_state(reconstructed);
1848 }
1849
1850 /// Returns the square of the local norm ... no comms
1851 template <typename T, std::size_t NDIM>
1855 return world.taskq.reduce<double,rangeT,do_norm2sq_local>(rangeT(coeffs.begin(),coeffs.end()),
1857 }
1858
1859
1860
1861
1862 /// Returns the maximum local depth of the tree ... no communications.
1863 template <typename T, std::size_t NDIM>
1865 std::size_t maxdepth = 0;
1866 typename dcT::const_iterator end = coeffs.end();
1867 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
1868 std::size_t N = (std::size_t) it->first.level();
1869 if (N> maxdepth)
1870 maxdepth = N;
1871 }
1872 return maxdepth;
1873 }
1874
1875
1876 /// Returns the maximum depth of the tree ... collective ... global sum/broadcast
1877 template <typename T, std::size_t NDIM>
1879 std::size_t maxdepth = max_local_depth();
1880 world.gop.max(maxdepth);
1881 return maxdepth;
1882 }
1883
1884 /// Returns the max number of nodes on a processor
1885 template <typename T, std::size_t NDIM>
1887 std::size_t maxsize = 0;
1888 maxsize = coeffs.size();
1889 world.gop.max(maxsize);
1890 return maxsize;
1891 }
1892
1893 /// Returns the min number of nodes on a processor
1894 template <typename T, std::size_t NDIM>
1896 std::size_t minsize = 0;
1897 minsize = coeffs.size();
1898 world.gop.min(minsize);
1899 return minsize;
1900 }
1901
1902 /// Returns the size of the tree structure of the function ... collective global sum
1903 template <typename T, std::size_t NDIM>
1905 std::size_t sum = 0;
1906 sum = coeffs.size();
1907 world.gop.sum(sum);
1908 return sum;
1909 }
1910
1911 /// Returns the number of coefficients in the function ... collective global sum
1912 template <typename T, std::size_t NDIM>
1913 std::size_t FunctionImpl<T,NDIM>::size() const {
1914 std::size_t sum = 0;
1915#if 1
1916 typename dcT::const_iterator end = coeffs.end();
1917 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
1918 const nodeT& node = it->second;
1919 if (node.has_coeff())
1920 sum+=node.size();
1921 }
1922 // print("proc",world.rank(),sum);
1923#else
1924 typename dcT::const_iterator end = coeffs.end();
1925 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
1926 const nodeT& node = it->second;
1927 if (node.has_coeff())
1928 ++sum;
1929 }
1930 if (is_compressed())
1931 for (std::size_t i=0; i<NDIM; ++i)
1932 sum *= 2*cdata.k;
1933 else
1934 for (std::size_t i=0; i<NDIM; ++i)
1935 sum *= cdata.k;
1936#endif
1937 world.gop.sum(sum);
1938
1939 return sum;
1940 }
1941
1942 /// Returns the number of coefficients in the function ... collective global sum
1943 template <typename T, std::size_t NDIM>
1945 std::size_t sum = coeffs.size() * (sizeof(keyT) + sizeof(nodeT));
1946 typename dcT::const_iterator end = coeffs.end();
1947 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
1948 const nodeT& node = it->second;
1949 if (node.has_coeff()) sum+=node.coeff().real_size();
1950 }
1951 world.gop.sum(sum);
1952 return sum;
1953 }
1954
1955 /// Returns the number of coefficients in the function ... collective global sum
1956 template <typename T, std::size_t NDIM>
1957 std::size_t FunctionImpl<T,NDIM>::nCoeff() const {
1958 std::size_t sum = coeffs.size() * (sizeof(keyT) + sizeof(nodeT));
1959 typename dcT::const_iterator end = coeffs.end();
1960 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
1961 const nodeT& node = it->second;
1962 if (node.has_coeff()) sum+=node.coeff().nCoeff();
1963 }
1964 world.gop.sum(sum);
1965 return sum;
1966 }
1967
1968
1969 /// print tree size and size
1970 template <typename T, std::size_t NDIM>
1971 void FunctionImpl<T,NDIM>::print_size(const std::string name) const {
1972 const size_t tsize=this->tree_size();
1973// const size_t size=this->size();
1974 const size_t ncoeff=this->nCoeff();
1975 const double wall=wall_time();
1976 const double d=sizeof(T);
1977 const double fac=1024*1024*1024;
1978
1979 double norm=0.0;
1980 {
1981 double local = norm2sq_local();
1982 this->world.gop.sum(local);
1983 this->world.gop.fence();
1984 norm=sqrt(local);
1985 }
1986
1987 if (this->world.rank()==0) {
1988
1989 constexpr std::size_t bufsize=128;
1990 char buf[bufsize];
1991 snprintf(buf, bufsize, "%40s at time %.1fs: norm/tree/#coeff/size: %7.5f %zu, %6.3f m, %6.3f GByte",
1992 (name.c_str()), wall, norm, tsize,double(ncoeff)*1.e-6,double(ncoeff)/fac*d);
1993 print(std::string(buf));
1994 }
1995 }
1996
1997 /// print the number of configurations per node
1998 template <typename T, std::size_t NDIM>
2000 if (this->targs.tt==TT_FULL) return;
2001 int dim=NDIM/2;
2002 int k0=k;
2003 if (is_compressed()) k0=2*k;
2004 Tensor<long> n(int(std::pow(double(k0),double(dim))+1));
2005 long n_full=0;
2006 long n_large=0;
2007
2008 if (world.rank()==0) print("n.size(),k0,dim",n.size(),k0,dim);
2009 typename dcT::const_iterator end = coeffs.end();
2010 for (typename dcT::const_iterator it=coeffs.begin(); it!=end; ++it) {
2011 const nodeT& node = it->second;
2012 if (node.has_coeff()) {
2013 if (node.coeff().rank()>long(n.size())) {
2014 ++n_large;
2015 } else if (node.coeff().rank()==-1) {
2016 ++n_full;
2017 } else if (node.coeff().rank()<0) {
2018 print("small rank",node.coeff().rank());
2019 } else {
2020 n[node.coeff().rank()]++;
2021 }
2022 }
2023 }
2024
2025 world.gop.sum(n.ptr(), n.size());
2026
2027 if (world.rank()==0) {
2028 print("configurations number of nodes");
2029 print(" full rank ",n_full);
2030 for (unsigned int i=0; i<n.size(); i++) {
2031 print(" ",i," ",n[i]);
2032 }
2033 print(" large rank ",n_large);
2034
2035 // repeat for logarithmic scale: <3, <10, <30, <100, ..
2036 Tensor<long> nlog(6);
2037 nlog=0;
2038 for (unsigned int i=0; i<std::min(3l,n.size()); i++) nlog[0]+=n[i];
2039 for (unsigned int i=3; i<std::min(10l,n.size()); i++) nlog[1]+=n[i];
2040 for (unsigned int i=10; i<std::min(30l,n.size()); i++) nlog[2]+=n[i];
2041 for (unsigned int i=30; i<std::min(100l,n.size()); i++) nlog[3]+=n[i];
2042 for (unsigned int i=100; i<std::min(300l,n.size()); i++) nlog[4]+=n[i];
2043 for (unsigned int i=300; i<std::min(1000l,n.size()); i++) nlog[5]+=n[i];
2044
2045 std::vector<std::string> slog={"3","10","30","100","300","1000"};
2046 for (unsigned int i=0; i<nlog.size(); i++) {
2047 print(" < ",slog[i]," ",nlog[i]);
2048 }
2049 print(" large rank ",n_large);
2050
2051 }
2052 }
2053
2054 template <typename T, std::size_t NDIM>
2057 const int k = cdata.k;
2058 double px[NDIM][k];
2059 T sum = T(0.0);
2060
2061 for (std::size_t i=0; i<NDIM; ++i) legendre_scaling_functions(x[i],k,px[i]);
2062
2063 if (NDIM == 1) {
2064 for (int p=0; p<k; ++p)
2065 sum += c(p)*px[0][p];
2066 }
2067 else if (NDIM == 2) {
2068 for (int p=0; p<k; ++p)
2069 for (int q=0; q<k; ++q)
2070 sum += c(p,q)*px[0][p]*px[1][q];
2071 }
2072 else if (NDIM == 3) {
2073 for (int p=0; p<k; ++p)
2074 for (int q=0; q<k; ++q)
2075 for (int r=0; r<k; ++r)
2076 sum += c(p,q,r)*px[0][p]*px[1][q]*px[2][r];
2077 }
2078 else if (NDIM == 4) {
2079 for (int p=0; p<k; ++p)
2080 for (int q=0; q<k; ++q)
2081 for (int r=0; r<k; ++r)
2082 for (int s=0; s<k; ++s)
2083 sum += c(p,q,r,s)*px[0][p]*px[1][q]*px[2][r]*px[3][s];
2084 }
2085 else if (NDIM == 5) {
2086 for (int p=0; p<k; ++p)
2087 for (int q=0; q<k; ++q)
2088 for (int r=0; r<k; ++r)
2089 for (int s=0; s<k; ++s)
2090 for (int t=0; t<k; ++t)
2091 sum += c(p,q,r,s,t)*px[0][p]*px[1][q]*px[2][r]*px[3][s]*px[4][t];
2092 }
2093 else if (NDIM == 6) {
2094 for (int p=0; p<k; ++p)
2095 for (int q=0; q<k; ++q)
2096 for (int r=0; r<k; ++r)
2097 for (int s=0; s<k; ++s)
2098 for (int t=0; t<k; ++t)
2099 for (int u=0; u<k; ++u)
2100 sum += c(p,q,r,s,t,u)*px[0][p]*px[1][q]*px[2][r]*px[3][s]*px[4][t]*px[5][u];
2101 }
2102 else {
2103 MADNESS_EXCEPTION("FunctionImpl:eval_cube:NDIM?",NDIM);
2104 }
2105 return sum*pow(2.0,0.5*NDIM*n)/sqrt(FunctionDefaults<NDIM>::get_cell_volume());
2106 }
2107
2108 template <typename T, std::size_t NDIM>
2109 void FunctionImpl<T,NDIM>::reconstruct_op(const keyT& key, const coeffT& s, const bool accumulate_NS) {
2110 //PROFILE_MEMBER_FUNC(FunctionImpl);
2111 // Note that after application of an integral operator not all
2112 // siblings may be present so it is necessary to check existence
2113 // and if absent insert an empty leaf node.
2114 //
2115 // If summing the result of an integral operator (i.e., from
2116 // non-standard form) there will be significant scaling function
2117 // coefficients at all levels and possibly difference coefficients
2118 // in leaves, hence the tree may refine as a result.
2119 typename dcT::iterator it = coeffs.find(key).get();
2120 if (it == coeffs.end()) {
2121 coeffs.replace(key,nodeT(coeffT(),false));
2122 it = coeffs.find(key).get();
2123 }
2124 nodeT& node = it->second;
2125
2126 // The integral operator will correctly connect interior nodes
2127 // to children but may leave interior nodes without coefficients
2128 // ... but they still need to sum down so just give them zeros
2129 if (node.has_children() && !node.has_coeff()) {
2130 node.set_coeff(coeffT(cdata.v2k,targs));
2131 }
2132
2133 if (node.has_children() || node.has_coeff()) { // Must allow for inconsistent state from transform, etc.
2134 coeffT d = node.coeff();
2135 if (!d.has_data()) d = coeffT(cdata.v2k,targs);
2136 if (accumulate_NS and (key.level() > 0)) d(cdata.s0) += s; // -- note accumulate for NS summation
2137 if (d.dim(0)==2*get_k()) { // d might be pre-truncated if it's a leaf
2138 d = unfilter(d);
2139 node.clear_coeff();
2140 node.set_has_children(true);
2141 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2142 const keyT& child = kit.key();
2143 coeffT ss = copy(d(child_patch(child)));
2144 ss.reduce_rank(thresh);
2145 //PROFILE_BLOCK(recon_send); // Too fine grain for routine profiling
2146 woT::task(coeffs.owner(child), &implT::reconstruct_op, child, ss, accumulate_NS);
2147 }
2148 } else {
2149 MADNESS_ASSERT(node.is_leaf());
2150 // node.coeff()+=s;
2151 node.coeff().reduce_rank(targs.thresh);
2152 }
2153 }
2154 else {
2155 coeffT ss=s;
2156 if (s.has_no_data()) ss=coeffT(cdata.vk,targs);
2157 if (key.level()) node.set_coeff(copy(ss));
2158 else node.set_coeff(ss);
2159 }
2161
2162 template <typename T, std::size_t NDIM>
2163 Tensor<T> fcube(const Key<NDIM>& key, T (*f)(const Vector<double,NDIM>&), const Tensor<double>& qx) {
2164 // fcube(key,typename FunctionFactory<T,NDIM>::FunctorInterfaceWrapper(f) , qx, fval);
2165 std::vector<long> npt(NDIM,qx.dim(0));
2166 Tensor<T> fval(npt);
2167 fcube(key,ElementaryInterface<T,NDIM>(f) , qx, fval);
2168 return fval;
2169 }
2170
2171 template <typename T, std::size_t NDIM>
2173 // fcube(key,typename FunctionFactory<T,NDIM>::FunctorInterfaceWrapper(f) , qx, fval);
2174 std::vector<long> npt(NDIM,qx.dim(0));
2175 Tensor<T> fval(npt);
2176 fcube(key, f, qx, fval);
2177 return fval;
2178 }
2179
2180 template <typename T, std::size_t NDIM>
2181 // void FunctionImpl<T,NDIM>::fcube(const keyT& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, tensorT& fval) const {
2182 void fcube(const Key<NDIM>& key, const FunctionFunctorInterface<T,NDIM>& f, const Tensor<double>& qx, Tensor<T>& fval) {
2183 //~ template <typename T, std::size_t NDIM> template< typename FF>
2184 //~ void FunctionImpl<T,NDIM>::fcube(const keyT& key, const FF& f, const Tensor<double>& qx, tensorT& fval) const {
2186 //PROFILE_MEMBER_FUNC(FunctionImpl);
2187 const Vector<Translation,NDIM>& l = key.translation();
2188 const Level n = key.level();
2189 const double h = std::pow(0.5,double(n));
2190 coordT c; // will hold the point in user coordinates
2191 const int npt = qx.dim(0);
2192
2195
2196 // Do pre-screening of the FunctionFunctorInterface, f, before calculating f(r) at quadrature points
2197 coordT c1, c2;
2198 for (std::size_t i = 0; i < NDIM; i++) {
2199 c1[i] = cell(i,0) + h*cell_width[i]*(l[i] + qx((long)0));
2200 c2[i] = cell(i,0) + h*cell_width[i]*(l[i] + qx(npt-1));
2201 }
2202 if (f.screened(c1, c2)) {
2203 fval(___) = 0.0;
2204 return;
2205 }
2206
2207 Tensor<double> vqx;
2208 bool vectorized = f.supports_vectorized();
2209 if (vectorized) {
2210 T* fvptr = fval.ptr();
2211 if (NDIM == 1) {
2212 double* x1 = new double[npt];
2213 int idx = 0;
2214 for (int i=0; i<npt; ++i, ++idx) {
2215 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2216 x1[idx] = c[0];
2217 }
2218 Vector<double*,1> xvals {x1};
2219 f(xvals, fvptr, npt);
2220 delete [] x1;
2221 }
2222 else if (NDIM == 2) {
2223 double* x1 = new double[npt*npt];
2224 double* x2 = new double[npt*npt];
2225 int idx = 0;
2226 for (int i=0; i<npt; ++i) {
2227 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2228 for (int j=0; j<npt; ++j, ++idx) {
2229 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2230 x1[idx] = c[0];
2231 x2[idx] = c[1];
2232 }
2233 }
2234 Vector<double*,2> xvals {x1, x2};
2235 f(xvals, fvptr, npt*npt);
2236 delete [] x1;
2237 delete [] x2;
2238 }
2239 else if (NDIM == 3) {
2240 double* x1 = new double[npt*npt*npt];
2241 double* x2 = new double[npt*npt*npt];
2242 double* x3 = new double[npt*npt*npt];
2243 int idx = 0;
2244 for (int i=0; i<npt; ++i) {
2245 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2246 for (int j=0; j<npt; ++j) {
2247 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2248 for (int k=0; k<npt; ++k, ++idx) {
2249 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2250 x1[idx] = c[0];
2251 x2[idx] = c[1];
2252 x3[idx] = c[2];
2253 }
2254 }
2255 }
2256 Vector<double*,3> xvals {x1, x2, x3};
2257 f(xvals, fvptr, npt*npt*npt);
2258 delete [] x1;
2259 delete [] x2;
2260 delete [] x3;
2261 }
2262 else if (NDIM == 4) {
2263 double* x1 = new double[npt*npt*npt*npt];
2264 double* x2 = new double[npt*npt*npt*npt];
2265 double* x3 = new double[npt*npt*npt*npt];
2266 double* x4 = new double[npt*npt*npt*npt];
2267 int idx = 0;
2268 for (int i=0; i<npt; ++i) {
2269 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2270 for (int j=0; j<npt; ++j) {
2271 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2272 for (int k=0; k<npt; ++k) {
2273 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2274 for (int m=0; m<npt; ++m, ++idx) {
2275 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2276 x1[idx] = c[0];
2277 x2[idx] = c[1];
2278 x3[idx] = c[2];
2279 x4[idx] = c[3];
2280 }
2281 }
2282 }
2283 }
2284 Vector<double*,4> xvals {x1, x2, x3, x4};
2285 f(xvals, fvptr, npt*npt*npt*npt);
2286 delete [] x1;
2287 delete [] x2;
2288 delete [] x3;
2289 delete [] x4;
2290 }
2291 else if (NDIM == 5) {
2292 double* x1 = new double[npt*npt*npt*npt*npt];
2293 double* x2 = new double[npt*npt*npt*npt*npt];
2294 double* x3 = new double[npt*npt*npt*npt*npt];
2295 double* x4 = new double[npt*npt*npt*npt*npt];
2296 double* x5 = new double[npt*npt*npt*npt*npt];
2297 int idx = 0;
2298 for (int i=0; i<npt; ++i) {
2299 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2300 for (int j=0; j<npt; ++j) {
2301 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2302 for (int k=0; k<npt; ++k) {
2303 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2304 for (int m=0; m<npt; ++m) {
2305 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2306 for (int n=0; n<npt; ++n, ++idx) {
2307 c[4] = cell(4,0) + h*cell_width[4]*(l[4] + qx(n)); // yy
2308 x1[idx] = c[0];
2309 x2[idx] = c[1];
2310 x3[idx] = c[2];
2311 x4[idx] = c[3];
2312 x5[idx] = c[4];
2313 }
2314 }
2315 }
2316 }
2317 }
2318 Vector<double*,5> xvals {x1, x2, x3, x4, x5};
2319 f(xvals, fvptr, npt*npt*npt*npt*npt);
2320 delete [] x1;
2321 delete [] x2;
2322 delete [] x3;
2323 delete [] x4;
2324 delete [] x5;
2325 }
2326 else if (NDIM == 6) {
2327 double* x1 = new double[npt*npt*npt*npt*npt*npt];
2328 double* x2 = new double[npt*npt*npt*npt*npt*npt];
2329 double* x3 = new double[npt*npt*npt*npt*npt*npt];
2330 double* x4 = new double[npt*npt*npt*npt*npt*npt];
2331 double* x5 = new double[npt*npt*npt*npt*npt*npt];
2332 double* x6 = new double[npt*npt*npt*npt*npt*npt];
2333 int idx = 0;
2334 for (int i=0; i<npt; ++i) {
2335 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2336 for (int j=0; j<npt; ++j) {
2337 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2338 for (int k=0; k<npt; ++k) {
2339 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2340 for (int m=0; m<npt; ++m) {
2341 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2342 for (int n=0; n<npt; ++n) {
2343 c[4] = cell(4,0) + h*cell_width[4]*(l[4] + qx(n)); // yy
2344 for (int p=0; p<npt; ++p, ++idx) {
2345 c[5] = cell(5,0) + h*cell_width[5]*(l[5] + qx(p)); // zz
2346 x1[idx] = c[0];
2347 x2[idx] = c[1];
2348 x3[idx] = c[2];
2349 x4[idx] = c[3];
2350 x5[idx] = c[4];
2351 x6[idx] = c[5];
2352 }
2353 }
2354 }
2355 }
2356 }
2357 }
2358 Vector<double*,6> xvals {x1, x2, x3, x4, x5, x6};
2359 f(xvals, fvptr, npt*npt*npt*npt*npt*npt);
2360 delete [] x1;
2361 delete [] x2;
2362 delete [] x3;
2363 delete [] x4;
2364 delete [] x5;
2365 delete [] x6;
2366 }
2367 else {
2368 MADNESS_EXCEPTION("FunctionImpl: fcube: confused about NDIM?",NDIM);
2369 }
2370 }
2371 else {
2372 if (NDIM == 1) {
2373 for (int i=0; i<npt; ++i) {
2374 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2375 fval(i) = f(c);
2376 MADNESS_ASSERT(!std::isnan(fval(i)));
2377 }
2378 }
2379 else if (NDIM == 2) {
2380 for (int i=0; i<npt; ++i) {
2381 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2382 for (int j=0; j<npt; ++j) {
2383 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2384 fval(i,j) = f(c);
2385 MADNESS_ASSERT(!std::isnan(fval(i,j)));
2386 }
2387 }
2388 }
2389 else if (NDIM == 3) {
2390 for (int i=0; i<npt; ++i) {
2391 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2392 for (int j=0; j<npt; ++j) {
2393 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2394 for (int k=0; k<npt; ++k) {
2395 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2396 fval(i,j,k) = f(c);
2397 MADNESS_ASSERT(!std::isnan(fval(i,j,k)));
2398 }
2399 }
2400 }
2401 }
2402 else if (NDIM == 4) {
2403 for (int i=0; i<npt; ++i) {
2404 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2405 for (int j=0; j<npt; ++j) {
2406 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2407 for (int k=0; k<npt; ++k) {
2408 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2409 for (int m=0; m<npt; ++m) {
2410 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2411 fval(i,j,k,m) = f(c);
2412 MADNESS_ASSERT(!std::isnan(fval(i,j,k,m)));
2413 }
2414 }
2415 }
2416 }
2417 }
2418 else if (NDIM == 5) {
2419 for (int i=0; i<npt; ++i) {
2420 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2421 for (int j=0; j<npt; ++j) {
2422 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2423 for (int k=0; k<npt; ++k) {
2424 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2425 for (int m=0; m<npt; ++m) {
2426 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2427 for (int n=0; n<npt; ++n) {
2428 c[4] = cell(4,0) + h*cell_width[4]*(l[4] + qx(n)); // yy
2429 fval(i,j,k,m,n) = f(c);
2430 MADNESS_ASSERT(!std::isnan(fval(i,j,k,m,n)));
2431 }
2432 }
2433 }
2434 }
2435 }
2436 }
2437 else if (NDIM == 6) {
2438 for (int i=0; i<npt; ++i) {
2439 c[0] = cell(0,0) + h*cell_width[0]*(l[0] + qx(i)); // x
2440 for (int j=0; j<npt; ++j) {
2441 c[1] = cell(1,0) + h*cell_width[1]*(l[1] + qx(j)); // y
2442 for (int k=0; k<npt; ++k) {
2443 c[2] = cell(2,0) + h*cell_width[2]*(l[2] + qx(k)); // z
2444 for (int m=0; m<npt; ++m) {
2445 c[3] = cell(3,0) + h*cell_width[3]*(l[3] + qx(m)); // xx
2446 for (int n=0; n<npt; ++n) {
2447 c[4] = cell(4,0) + h*cell_width[4]*(l[4] + qx(n)); // yy
2448 for (int p=0; p<npt; ++p) {
2449 c[5] = cell(5,0) + h*cell_width[5]*(l[5] + qx(p)); // zz
2450 fval(i,j,k,m,n,p) = f(c);
2451 MADNESS_ASSERT(!std::isnan(fval(i,j,k,m,n,p)));
2452 }
2453 }
2454 }
2455 }
2456 }
2457 }
2458 }
2459 else {
2460 MADNESS_EXCEPTION("FunctionImpl: fcube: confused about NDIM?",NDIM);
2461 }
2462 }
2463 }
2464
2465 template <typename T, std::size_t NDIM>
2466 void FunctionImpl<T,NDIM>::fcube(const keyT& key, T (*f)(const coordT&), const Tensor<double>& qx, tensorT& fval) const {
2467 // fcube(key,typename FunctionFactory<T,NDIM>::FunctorInterfaceWrapper(f) , qx, fval);
2469 }
2470
2471 template <typename T, std::size_t NDIM>
2473 madness::fcube(key,f,qx,fval);
2474 }
2475
2476
2477 /// project the functor into this functionimpl, and "return" a tree in reconstructed,
2478 /// rank-reduced form.
2479
2480 /// @param[in] key current FunctionNode
2481 /// @param[in] do_refine
2482 /// @param[in] specialpts in case these are very spiky functions -- don't undersample
2483 template <typename T, std::size_t NDIM>
2485 bool do_refine,
2486 const std::vector<Vector<double,NDIM> >& specialpts) {
2487 //PROFILE_MEMBER_FUNC(FunctionImpl);
2488 if (do_refine && key.level() < max_refine_level) {
2489
2490 // Restrict special points to this box
2491 std::vector<Vector<double,NDIM> > newspecialpts;
2492 if (key.level() < functor->special_level() && specialpts.size() > 0) {
2494 std::vector<bool> bperiodic = bc.is_periodic();
2495 for (unsigned int i = 0; i < specialpts.size(); ++i) {
2496 coordT simpt;
2497 user_to_sim(specialpts[i], simpt);
2498 Key<NDIM> specialkey = simpt2key(simpt, key.level());
2499 if (specialkey.is_neighbor_of(key,bperiodic)) {
2500 newspecialpts.push_back(specialpts[i]);
2501 }
2502 }
2503 }
2504
2505 // If refining compute scaling function coefficients and
2506 // norm of difference coefficients
2507 tensorT r, s0;
2508 double dnorm = 0.0;
2509 //////////////////////////if (newspecialpts.size() == 0)
2510 {
2511 // Make in r child scaling function coeffs at level n+1
2512 r = tensorT(cdata.v2k);
2513 for (KeyChildIterator<NDIM> it(key); it; ++it) {
2514 const keyT& child = it.key();
2515 r(child_patch(child)) = project(child);
2516 }
2517 // Filter then test difference coeffs at level n
2518 tensorT d = filter(r);
2519 if (truncate_on_project) s0 = copy(d(cdata.s0));
2520 d(cdata.s0) = T(0);
2521 dnorm = d.normf();
2522 }
2523
2524 // If have special points always refine. If don't have special points
2525 // refine if difference norm is big
2526 if (newspecialpts.size() > 0 || dnorm >=truncate_tol(thresh,key.level())) {
2527 coeffs.replace(key,nodeT(coeffT(),true)); // Insert empty node for parent
2528 for (KeyChildIterator<NDIM> it(key); it; ++it) {
2529 const keyT& child = it.key();
2530 ProcessID p;
2532 p = world.random_proc();
2533 }
2534 else {
2535 p = coeffs.owner(child);
2536 }
2537 //PROFILE_BLOCK(proj_refine_send); // Too fine grain for routine profiling
2538 woT::task(p, &implT::project_refine_op, child, do_refine, newspecialpts);
2539 }
2540 }
2541 else {
2542 if (truncate_on_project) {
2544 coeffs.replace(key,nodeT(s,false));
2545 }
2546 else {
2547 coeffs.replace(key,nodeT(coeffT(),true)); // Insert empty node for parent
2548 for (KeyChildIterator<NDIM> it(key); it; ++it) {
2549 const keyT& child = it.key();
2550 coeffT s(r(child_patch(child)),thresh,FunctionDefaults<NDIM>::get_tensor_type());
2551 coeffs.replace(child,nodeT(s,false));
2552 }
2553 }
2554 }
2555 }
2556 else {
2557 coeffs.replace(key,nodeT(coeffT(project(key),targs),false));
2558 }
2559 }
2560
2561 template <typename T, std::size_t NDIM>
2563 std::vector<long> v0(NDIM,0L);
2564 std::vector<long> v1(NDIM,1L);
2565 std::vector<Slice> s(NDIM,Slice(0,0));
2566 const TensorArgs full_args(-1.0,TT_FULL);
2567 if (is_compressed()) {
2568 if (world.rank() == coeffs.owner(cdata.key0)) {
2569 typename dcT::iterator it = coeffs.find(cdata.key0).get();
2570 MADNESS_ASSERT(it != coeffs.end());
2571 nodeT& node = it->second;
2572 MADNESS_ASSERT(node.has_coeff());
2573 // node.node_to_full_rank();
2574 // node.full_tensor_reference()(v0) += t*sqrt(FunctionDefaults<NDIM>::get_cell_volume());
2575 // node.node_to_low_rank();
2576 change_tensor_type(node.coeff(),full_args);
2578 change_tensor_type(node.coeff(),targs);
2579 }
2580 }
2581 else {
2582 for (typename dcT::iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
2583 Level n = it->first.level();
2584 nodeT& node = it->second;
2585 if (node.has_coeff()) {
2586 // this looks funny, but is necessary for GenTensor, since you can't access a
2587 // single matrix element. Therefore make a (1^NDIM) tensor, convert to GenTensor, then
2588 // add to the original one by adding a slice.
2589 tensorT ttt(v1);
2590 ttt=t*sqrt(FunctionDefaults<NDIM>::get_cell_volume()*pow(0.5,double(NDIM*n)));
2591 coeffT tt(ttt,get_tensor_args());
2592 node.coeff()(s) += tt;
2593 // this was the original line:
2594 // node.coeff().full_tensor()(v0) += t*sqrt(FunctionDefaults<NDIM>::get_cell_volume()*pow(0.5,double(NDIM*n)));
2595
2596 }
2597 }
2598 }
2599 if (fence) world.gop.fence();
2600 }
2601
2602 template <typename T, std::size_t NDIM>
2605 if (is_compressed()) initial_level = std::max(initial_level,1); // Otherwise zero function is confused
2606 if (coeffs.is_local(key)) {
2607 if (is_compressed()) {
2608 if (key.level() == initial_level) {
2609 coeffs.replace(key, nodeT(coeffT(), false));
2610 }
2611 else {
2612 coeffs.replace(key, nodeT(coeffT(cdata.v2k,targs), true));
2613 }
2614 }
2615 else {
2616 if (key.level()<initial_level) {
2617 coeffs.replace(key, nodeT(coeffT(), true));
2618 }
2619 else {
2620 coeffs.replace(key, nodeT(coeffT(cdata.vk,targs), false));
2621 }
2622 }
2623 }
2624 if (key.level() < initial_level) {
2625 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2626 insert_zero_down_to_initial_level(kit.key());
2627 }
2628 }
2629
2630 }
2631
2632
2633 template <typename T, std::size_t NDIM>
2635 //PROFILE_MEMBER_FUNC(FunctionImpl);
2636 typename dcT::iterator it = coeffs.find(key).get();
2637 if (it == coeffs.end()) {
2638 // In a standard tree all children would exist but some ops (transform)
2639 // can leave the tree in a messy state. Just make the missing node as an
2640 // empty leaf.
2641 coeffs.replace(key,nodeT());
2642 it = coeffs.find(key).get();
2643 }
2644 nodeT& node = it->second;
2645 if (node.has_children()) {
2646 std::vector< Future<bool> > v = future_vector_factory<bool>(1<<NDIM);
2647 int i=0;
2648 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
2649 v[i] = woT::task(coeffs.owner(kit.key()), &implT::truncate_spawn, kit.key(), tol, TaskAttributes::generator());
2650 }
2651 return woT::task(world.rank(),&implT::truncate_op, key, tol, v);
2652 }
2653 else {
2654 // In compressed form leaves should not have coeffs ... however the
2655 // transform op could leave the tree with leaves that do have coeffs
2656 // in which case we want something sensible to happen
2657 //MADNESS_ASSERT(!node.has_coeff());
2658 if (node.has_coeff() && key.level()>1) {
2659 double dnorm = node.coeff().normf();
2660 if (dnorm < truncate_tol(tol,key)) {
2661 node.clear_coeff();
2662 }
2663 }
2664 return Future<bool>(node.has_coeff());
2665 }
2666 }
2667
2668
2669 template <typename T, std::size_t NDIM>
2670 bool FunctionImpl<T,NDIM>::truncate_op(const keyT& key, double tol, const std::vector< Future<bool> >& v) {
2671 //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling
2672 // If any child has coefficients, a parent cannot truncate
2673 for (int i=0; i<(1<<NDIM); ++i) if (v[i].get()) return true;
2674 nodeT& node = coeffs.find(key).get()->second;
2675
2676 // Interior nodes should always have coeffs but transform might
2677 // leave empty interior nodes ... hence just force no coeffs to
2678 // be zero coeff unless it is a leaf.
2679 if (node.has_children() && !node.has_coeff()) node.set_coeff(coeffT(cdata.v2k,targs));
2680
2681 if (key.level() > 1) { // >1 rather >0 otherwise reconstruct might get confused
2682 double dnorm = node.coeff().normf();
2683 if (dnorm < truncate_tol(tol,key)) {
2684 node.clear_coeff();
2685 if (node.has_children()) {
2686 node.set_has_children(false);
2687 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2688 coeffs.erase(kit.key());
2689 }
2690 }
2691 }
2692 }
2693 return node.has_coeff();
2694 }
2695
2696
2697 template <typename T, std::size_t NDIM>
2698 void FunctionImpl<T,NDIM>::print_tree(std::ostream& os, Level maxlevel) const {
2699 if (world.rank() == 0) do_print_tree(cdata.key0, os, maxlevel);
2700 world.gop.fence();
2701 if (world.rank() == 0) os.flush();
2702 world.gop.fence();
2703 }
2704
2705
2706 template <typename T, std::size_t NDIM>
2707 void FunctionImpl<T,NDIM>::do_print_tree(const keyT& key, std::ostream& os, Level maxlevel) const {
2708 typename dcT::const_iterator it = coeffs.find(key).get();
2709 if (it == coeffs.end()) {
2710 //MADNESS_EXCEPTION("FunctionImpl: do_print_tree: null node pointer",0);
2711 for (int i=0; i<key.level(); ++i) os << " ";
2712 os << key << " missing --> " << coeffs.owner(key) << "\n";
2713 }
2714 else {
2715 const nodeT& node = it->second;
2716 for (int i=0; i<key.level(); ++i) os << " ";
2717 os << key << " " << node << " --> " << coeffs.owner(key) << "\n";
2718 if (key.level() < maxlevel && node.has_children()) {
2719 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2720 do_print_tree(kit.key(),os,maxlevel);
2721 }
2723 }
2724 }
2725
2726 template <typename T, std::size_t NDIM>
2727 void FunctionImpl<T,NDIM>::print_tree_json(std::ostream& os, Level maxlevel) const {
2728 std::multimap<Level, std::tuple<tranT, std::string>> data;
2729 if (world.rank() == 0) do_print_tree_json(cdata.key0, data, maxlevel);
2730 world.gop.fence();
2731 if (world.rank() == 0) {
2732 for (Level level = 0; level != maxlevel; ++level) {
2733 if (data.count(level) == 0)
2734 break;
2735 else {
2736 if (level > 0)
2737 os << ",";
2738 os << "\"" << level << "\":{";
2739 os << "\"level\": " << level << ",";
2740 os << "\"nodes\":{";
2741 auto range = data.equal_range(level);
2742 for (auto it = range.first; it != range.second; ++it) {
2743 os << "\"" << std::get<0>(it->second) << "\":"
2744 << std::get<1>(it->second);
2745 if (std::next(it) != range.second)
2746 os << ",";
2747 }
2748 os << "}}";
2749 }
2750 }
2751 os.flush();
2752 }
2753 world.gop.fence();
2754 }
2755
2756
2757 template <typename T, std::size_t NDIM>
2758 void FunctionImpl<T,NDIM>::do_print_tree_json(const keyT& key, std::multimap<Level, std::tuple<tranT, std::string>>& data, Level maxlevel) const {
2759 typename dcT::const_iterator it = coeffs.find(key).get();
2760 if (it == coeffs.end()) {
2761 MADNESS_EXCEPTION("FunctionImpl: do_print_tree_json: null node pointer",0);
2762 }
2763 else {
2764 const nodeT& node = it->second;
2765 std::ostringstream oss;
2766 oss << "{";
2767 node.print_json(oss);
2768 oss << ",\"owner\": " << coeffs.owner(key) << "}";
2769 auto node_json_str = oss.str();
2770 data.insert(std::make_pair(key.level(), std::make_tuple(key.translation(), node_json_str)));
2771 if (key.level() < maxlevel && node.has_children()) {
2772 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2773 do_print_tree_json(kit.key(),data, maxlevel);
2774 }
2775 }
2776 }
2777 }
2778
2779 template <typename T, std::size_t NDIM>
2780 void FunctionImpl<T,NDIM>::print_tree_graphviz(std::ostream& os, Level maxlevel) const {
2781 // aggregate data by level, thus collect data first, then dump
2782 if (world.rank() == 0) do_print_tree_graphviz(cdata.key0, os, maxlevel);
2783 world.gop.fence();
2784 if (world.rank() == 0) os.flush();
2785 world.gop.fence();
2786 }
2787
2788 template <typename T, std::size_t NDIM>
2789 void FunctionImpl<T,NDIM>::do_print_tree_graphviz(const keyT& key, std::ostream& os, Level maxlevel) const {
2790
2791 struct uniqhash {
2792 static int64_t value(const keyT& key) {
2793 int64_t result = 0;
2794 for (int64_t j = 0; j <= key.level()-1; ++j) {
2795 result += (1 << j*NDIM);
2796 }
2797 result += key.translation()[0];
2798 return result;
2799 }
2800 };
2801
2802 typename dcT::const_iterator it = coeffs.find(key).get();
2803 if (it != coeffs.end()) {
2804 const nodeT& node = it->second;
2805 if (key.level() < maxlevel && node.has_children()) {
2806 for (KeyChildIterator<NDIM> kit(key); kit; ++kit) {
2807 os << uniqhash::value(key) << " -> " << uniqhash::value(kit.key()) << "\n";
2808 do_print_tree_graphviz(kit.key(),os,maxlevel);
2809 }
2810 }
2811 }
2812 }
2813
2814 template <typename T, std::size_t NDIM>
2816 //PROFILE_MEMBER_FUNC(FunctionImpl);
2817
2818 if (not functor) MADNESS_EXCEPTION("FunctionImpl: project: confusion about function?",0);
2819
2820 // if functor provides coeffs directly, awesome; otherwise use compute by yourself
2821 if (functor->provides_coeff()) return functor->coeff(key).full_tensor_copy();
2822
2823 MADNESS_ASSERT(cdata.npt == cdata.k); // only necessary due to use of fast transform
2824 tensorT fval(cdata.vq,false); // this will be the returned result
2825 tensorT work(cdata.vk,false); // initially evaluate the function in here
2826 tensorT workq(cdata.vq,false); // initially evaluate the function in here
2827
2828 // compute the values of the functor at the quadrature points and scale appropriately
2829 madness::fcube(key,*functor,cdata.quad_x,work);
2830 work.scale(sqrt(FunctionDefaults<NDIM>::get_cell_volume()*pow(0.5,double(NDIM*key.level()))));
2831 //return transform(work,cdata.quad_phiw);
2832 return fast_transform(work,cdata.quad_phiw,fval,workq);
2833 }
2834
2835 template <typename T, std::size_t NDIM>
2837 if (coeffs.probe(key)) {
2838 return Future<double>(coeffs.find(key).get()->second.get_norm_tree());
2839 }
2840 MADNESS_ASSERT(key.level());
2841 keyT parent = key.parent();
2842 return woT::task(coeffs.owner(parent), &implT::get_norm_tree_recursive, parent, TaskAttributes::hipri());
2843 }
2844
2845
2846 template <typename T, std::size_t NDIM>
2848 const RemoteReference< FutureImpl< std::pair<keyT,coeffT> > >& ref) const {
2849 //PROFILE_MEMBER_FUNC(FunctionImpl);
2850 if (coeffs.probe(key)) {
2851 const nodeT& node = coeffs.find(key).get()->second;
2852 Future< std::pair<keyT,coeffT> > result(ref);
2853 if (node.has_coeff()) {
2854 //madness::print("sock found it with coeff",key);
2855 result.set(std::pair<keyT,coeffT>(key,node.coeff()));
2856 }
2857 else {
2858 //madness::print("sock found it without coeff",key);
2859 result.set(std::pair<keyT,coeffT>(key,coeffT()));
2860 }
2861 }
2862 else {
2863 keyT parent = key.parent();
2864 //madness::print("sock forwarding to parent",key,parent);
2865 //PROFILE_BLOCK(sitome_send); // Too fine grain for routine profiling
2866 if (coeffs.is_local(parent))
2867 woT::send(coeffs.owner(parent), &FunctionImpl<T,NDIM>::sock_it_to_me, parent, ref);
2868 else
2869 woT::task(coeffs.owner(parent), &FunctionImpl<T,NDIM>::sock_it_to_me, parent, ref, TaskAttributes::hipri());
2870 }
2871 }
2872
2873 // like sock_it_to_me, but it replaces empty node with averaged coeffs from further down the tree
2874 template <typename T, std::size_t NDIM>
2876 const RemoteReference< FutureImpl< std::pair<keyT,coeffT> > >& ref) const {
2878 if (coeffs.probe(key)) {
2879 const nodeT& node = coeffs.find(key).get()->second;
2880 Future< std::pair<keyT,coeffT> > result(ref);
2881 if (node.has_coeff()) {
2882 result.set(std::pair<keyT,coeffT>(key,node.coeff()));
2883 }
2884 else {
2885 result.set(std::pair<keyT,coeffT>(key,nodeT(coeffT(project(key),targs),false).coeff()));
2886 }
2887 }
2888 else {
2889 keyT parent = key.parent();
2890 //PROFILE_BLOCK(sitome2_send); // Too fine grain for routine profiling
2891 woT::task(coeffs.owner(parent), &FunctionImpl<T,NDIM>::sock_it_to_me_too, parent, ref, TaskAttributes::hipri());
2892 }
2893 }
2894
2895
2896 template <typename T, std::size_t NDIM>
2898 const keyT& keyin,
2899 const typename Future<T>::remote_refT& ref) {
2900
2902 // This is ugly. We must figure out a clean way to use
2903 // owner computes rule from the container.
2904 Vector<double,NDIM> x = xin;
2905 keyT key = keyin;
2907 ProcessID me = world.rank();
2908 while (1) {
2909 ProcessID owner = coeffs.owner(key);
2910 if (owner != me) {
2911 //PROFILE_BLOCK(eval_send); // Too fine grain for routine profiling
2912 woT::task(owner, &implT::eval, x, key, ref, TaskAttributes::hipri());
2913 return;
2914 }
2915 else {
2916 typename dcT::futureT fut = coeffs.find(key);
2917 typename dcT::iterator it = fut.get();
2918 nodeT& node = it->second;
2919 if (node.has_coeff()) {
2920 Future<T>(ref).set(eval_cube(key.level(), x, node.coeff().full_tensor_copy()));
2921 return;
2922 }
2923 else {
2924 for (std::size_t i=0; i<NDIM; ++i) {
2925 double xi = x[i]*2.0;
2926 int li = int(xi);
2927 if (li == 2) li = 1;
2928 x[i] = xi - li;
2929 l[i] = 2*l[i] + li;
2930 }
2931 key = keyT(key.level()+1,l);
2932 }
2933 }
2934 }
2935 //MADNESS_EXCEPTION("should not be here",0);
2936 }
2937
2938
2939 template <typename T, std::size_t NDIM>
2940 std::pair<bool,T>
2942 Vector<double,NDIM> x = xin;
2943 keyT key(0);
2945 const ProcessID me = world.rank();
2946 while (key.level() <= maxlevel) {
2947 if (coeffs.owner(key) == me) {
2948 typename dcT::futureT fut = coeffs.find(key);
2949 typename dcT::iterator it = fut.get();
2950 if (it != coeffs.end()) {
2951 nodeT& node = it->second;
2952 if (node.has_coeff()) {
2953 return std::pair<bool,T>(true,eval_cube(key.level(), x, node.coeff().full_tensor_copy()));
2954 }
2955 }
2956 }
2957 for (std::size_t i=0; i<NDIM; ++i) {
2958 double xi = x[i]*2.0;
2959 int li = int(xi);
2960 if (li == 2) li = 1;
2961 x[i] = xi - li;
2962 l[i] = 2*l[i] + li;
2963 }
2964 key = keyT(key.level()+1,l);
2965 }
2966 return std::pair<bool,T>(false,0.0);
2967 }
2968
2969 template <typename T, std::size_t NDIM>
2971 const keyT& keyin,
2972 const typename Future<Level>::remote_refT& ref) {
2973
2975 // This is ugly. We must figure out a clean way to use
2976 // owner computes rule from the container.
2977 Vector<double,NDIM> x = xin;
2978 keyT key = keyin;
2980 ProcessID me = world.rank();
2981 while (1) {
2982 ProcessID owner = coeffs.owner(key);
2983 if (owner != me) {
2984 //PROFILE_BLOCK(eval_send); // Too fine grain for routine profiling
2985 woT::task(owner, &implT::evaldepthpt, x, key, ref, TaskAttributes::hipri());
2986 return;
2987 }
2988 else {
2989 typename dcT::futureT fut = coeffs.find(key);
2990 typename dcT::iterator it = fut.get();
2991 nodeT& node = it->second;
2992 if (node.has_coeff()) {
2993 Future<Level>(ref).set(key.level());
2994 return;
2995 }
2996 else {
2997 for (std::size_t i=0; i<NDIM; ++i) {
2998 double xi = x[i]*2.0;
2999 int li = int(xi);
3000 if (li == 2) li = 1;
3001 x[i] = xi - li;
3002 l[i] = 2*l[i] + li;
3003 }
3004 key = keyT(key.level()+1,l);
3005 }
3006 }
3007 }
3008 //MADNESS_EXCEPTION("should not be here",0);
3009 }
3010
3011 template <typename T, std::size_t NDIM>
3013 const keyT& keyin,
3014 const typename Future<long>::remote_refT& ref) {
3015
3017 // This is ugly. We must figure out a clean way to use
3018 // owner computes rule from the container.
3019 Vector<double,NDIM> x = xin;
3020 keyT key = keyin;
3022 ProcessID me = world.rank();
3023 while (1) {
3024 ProcessID owner = coeffs.owner(key);
3025 if (owner != me) {
3026 //PROFILE_BLOCK(eval_send); // Too fine grain for routine profiling
3027 woT::task(owner, &implT::evalR, x, key, ref, TaskAttributes::hipri());
3028 return;
3029 }
3030 else {
3031 typename dcT::futureT fut = coeffs.find(key);
3032 typename dcT::iterator it = fut.get();
3033 nodeT& node = it->second;
3034 if (node.has_coeff()) {
3035 Future<long>(ref).set(node.coeff().rank());
3036 return;
3037 }
3038 else {
3039 for (std::size_t i=0; i<NDIM; ++i) {
3040 double xi = x[i]*2.0;
3041 int li = int(xi);
3042 if (li == 2) li = 1;
3043 x[i] = xi - li;
3044 l[i] = 2*l[i] + li;
3045 }
3046 key = keyT(key.level()+1,l);
3047 }
3048 }
3049 }
3050 //MADNESS_EXCEPTION("should not be here",0);
3051 }
3052
3053
3054 template <typename T, std::size_t NDIM>
3055 void FunctionImpl<T,NDIM>::tnorm(const tensorT& t, double* lo, double* hi) {
3056 //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling
3057 auto& cdata=FunctionCommonData<T,NDIM>::get(t.dim(0));
3058 tensorT work = copy(t);
3059 tensorT tlo = work(cdata.sh);
3060 *lo = tlo.normf();
3061 tlo.fill(0.0);
3062 *hi = work.normf();
3063 }
3064
3065 template <typename T, std::size_t NDIM>
3066 void FunctionImpl<T,NDIM>::tnorm(const GenTensor<T>& t, double* lo, double* hi) {
3067 auto& cdata=FunctionCommonData<T,NDIM>::get(t.dim(0));
3068 coeffT shalf=t(cdata.sh);
3069 *lo=shalf.normf();
3070 coeffT sfull=copy(t);
3071 sfull(cdata.sh)-=shalf;
3072 *hi=sfull.normf();
3073 }
3074
3075 template <typename T, std::size_t NDIM>
3076 void FunctionImpl<T,NDIM>::tnorm(const SVDTensor<T>& t, double* lo, double* hi,
3077 const int particle) {
3078 *lo=0.0;
3079 *hi=0.0;
3080 auto& cdata=FunctionCommonData<T,NDIM>::get(t.dim(0));
3081 if (t.rank()==0) return;
3082 const tensorT vec=t.flat_vector(particle-1);
3083 for (long i=0; i<t.rank(); ++i) {
3084 double lo1,hi1;
3085 tensorT c=vec(Slice(i,i),_).reshape(cdata.vk);
3086 tnorm(c, &lo1, &hi1); // note we use g instead of h, since g is 3D
3087 *lo+=lo1*t.weights(i);
3088 *hi+=hi1*t.weights(i);
3089 }
3090 }
3091
3092
3093 namespace detail {
3094 template <typename A, typename B>
3095 struct noop {
3096 void operator()(const A& a, const B& b) const {};
3097
3098 template <typename Archive> void serialize(Archive& ar) {}
3099 };
3100
3101 template <typename T, std::size_t NDIM>
3105 // G++ 4.1.2 ICEs on BGP ... scaleinplace(T q) : q(q) {}
3106 scaleinplace(T q) {this->q = q;}
3107 void operator()(const Key<NDIM>& key, Tensor<T>& t) const {
3108 t.scale(q);
3109 }
3110 void operator()(const Key<NDIM>& key, FunctionNode<T,NDIM>& node) const {
3111 node.coeff().scale(q);
3112 }
3113 template <typename Archive> void serialize(Archive& ar) {
3114 ar & q;
3115 }
3116 };
3117
3118 template <typename T, std::size_t NDIM>
3120 void operator()(const Key<NDIM>& key, Tensor<T>& t) const {
3121 t.emul(t);
3122 }
3123 template <typename Archive> void serialize(Archive& ar) {}
3124 };
3125
3126 template <typename T, std::size_t NDIM>
3127 struct absinplace {
3128 void operator()(const Key<NDIM>& key, Tensor<T>& t) const {t=abs(t);}
3129 template <typename Archive> void serialize(Archive& ar) {}
3130 };
3131
3132 template <typename T, std::size_t NDIM>
3134 void operator()(const Key<NDIM>& key, Tensor<T>& t) const {abs(t.emul(t));}
3135 template <typename Archive> void serialize(Archive& ar) {}
3136 };
3137
3138 }
3139
3140template <typename T, std::size_t NDIM>
3141 void FunctionImpl<T,NDIM>::scale_inplace(const T q, bool fence) {
3142 // unary_op_coeff_inplace(detail::scaleinplace<T,NDIM>(q), fence);
3143 unary_op_node_inplace(detail::scaleinplace<T,NDIM>(q), fence);
3144 }
3145
3146 template <typename T, std::size_t NDIM>
3148 //unary_op_value_inplace(&implT::autorefine_square_test, detail::squareinplace<T,NDIM>(), fence);
3149 unary_op_value_inplace(detail::squareinplace<T,NDIM>(), fence);
3150 }
3151
3152 template <typename T, std::size_t NDIM>
3154 unary_op_value_inplace(detail::absinplace<T,NDIM>(), fence);
3155 }
3156
3157 template <typename T, std::size_t NDIM>
3159 unary_op_value_inplace(detail::abssquareinplace<T,NDIM>(), fence);
3160 }
3161
3162 template <typename T, std::size_t NDIM>
3164 //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling
3165 double p[200];
3166 double scale = pow(2.0,double(np-nc));
3167 for (int mu=0; mu<cdata.npt; ++mu) {
3168 double xmu = scale*(cdata.quad_x(mu)+lc) - lp;
3169 MADNESS_ASSERT(xmu>-1e-15 && xmu<(1+1e-15));
3170 legendre_scaling_functions(xmu,cdata.k,p);
3171 for (int i=0; i<k; ++i) phi(i,mu) = p[i];
3172 }
3173 phi.scale(pow(2.0,0.5*np));
3174 }
3175
3176 template <typename T, std::size_t NDIM>
3177
3178 const GenTensor<T> FunctionImpl<T,NDIM>::parent_to_child(const coeffT& s, const keyT& parent, const keyT& child) const {
3179 //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling
3180 // An invalid parent/child means that they are out of the box
3181 // and it is the responsibility of the caller to worry about that
3182 // ... most likely the coefficients (s) are zero to reflect
3183 // zero B.C. so returning s makes handling this easy.
3184 if (parent == child || parent.is_invalid() || child.is_invalid()) return s;
3185
3186 coeffT result = fcube_for_mul<T>(child, parent, s);
3187 result.scale(sqrt(FunctionDefaults<NDIM>::get_cell_volume()*pow(0.5,double(NDIM*child.level()))));
3188 result = transform(result,cdata.quad_phiw);
3189
3190 return result;
3191 }
3192
3193
3194 template <typename T, std::size_t NDIM>
3197 std::vector<long> v0(NDIM,0);
3198 T sum = 0.0;
3199 if (is_compressed()) {
3200 if (world.rank() == coeffs.owner(cdata.key0)) {
3201 typename dcT::const_iterator it = coeffs.find(cdata.key0).get();
3202 if (it != coeffs.end()) {
3203 const nodeT& node = it->second;
3204 if (node.has_coeff()) sum = node.coeff().full_tensor_copy()(v0);
3205 }
3206 }
3207 }
3208 else {
3209 for (typename dcT::const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
3210 const keyT& key = it->first;
3211 const nodeT& node = it->second;
3212 if (node.has_coeff()) sum += node.coeff().full_tensor_copy()(v0)*pow(0.5,NDIM*key.level()*0.5);
3213 }
3214 }
3216 }
3217
3218
3219 static inline bool enforce_bc(bool is_periodic, Level n, Translation& l) {
3220 Translation two2n = 1ul << n;
3221 if (l < 0) {
3222 if (is_periodic)
3223 l += two2n; // Periodic BC
3224 else
3225 return false; // Zero BC
3226 }
3227 else if (l >= two2n) {
3228 if (is_periodic)
3229 l -= two2n; // Periodic BC
3230 else
3231 return false; // Zero BC
3232 }
3233 return true;
3234 }
3235
3236 template <typename T, std::size_t NDIM>
3237 Key<NDIM> FunctionImpl<T,NDIM>::neighbor(const keyT& key, const Key<NDIM>& disp, const std::vector<bool>& is_periodic) const {
3239
3240 for (std::size_t axis=0; axis<NDIM; ++axis) {
3241 l[axis] += disp.translation()[axis];
3243 //if (!enforce_bc(bc(axis,0), bc(axis,1), key.level(), l[axis])) {
3244 if (!enforce_bc(is_periodic[axis], key.level(), l[axis])) {
3245 return keyT::invalid();
3246 }
3247 }
3248 return keyT(key.level(),l);
3249 }
3250
3251
3252 template <typename T, std::size_t NDIM>
3253 Future< std::pair< Key<NDIM>, GenTensor<T> > >
3255 //PROFILE_MEMBER_FUNC(FunctionImpl); // Too fine grain for routine profiling
3256 typedef std::pair< Key<NDIM>,coeffT > argT;
3257 Future<argT> result;
3258 //PROFILE_BLOCK(find_me_send); // Too fine grain for routine profiling
3259 woT::task(coeffs.owner(key), &implT::sock_it_to_me_too, key, result.remote_ref(world), TaskAttributes::hipri());
3260 return result;
3261 }
3262
3263
3264 /// will insert
3265 /// @return s coefficient and norm_tree for key
3266 template <typename T, std::size_t NDIM>
3268 bool nonstandard1, bool keepleaves, bool redundant1) {
3269 if (!coeffs.probe(key)) print("missing node",key);
3270 MADNESS_ASSERT(coeffs.probe(key));
3271
3272 // get fetches remote data (here actually local)
3273 nodeT& node = coeffs.find(key).get()->second;
3275 // internal node -> continue recursion
3276 if (node.has_children()) {
3277 std::vector< Future<std::pair<coeffT,double> > > v = future_vector_factory<std::pair<coeffT,double> >(1<<NDIM);
3278 int i=0;
3279 for (KeyChildIterator<NDIM> kit(key); kit; ++kit,++i) {
3280 //PROFILE_BLOCK(compress_send); // Too fine grain for routine profiling
3281 // readily available
3282 v[i] = woT::task(coeffs.owner(kit.key()), &implT::compress_spawn, kit.key(),
3283 nonstandard1, keepleaves, redundant1, TaskAttributes::hipri());
3284 }
3285 if (redundant1) return woT::task(world.rank(),&implT::make_redundant_op, key, v);
3286 return woT::task(world.rank(),&implT::compress_op, key, v, nonstandard1);
3287 }
3288
3289 // leaf node -> remove coefficients here and pass them back to parent for filtering
3290 // insert snorm, dnorm=0.0, normtree (=snorm)
3291 else {
3292 // special case: tree has only root node: keep sum coeffs and make zero diff coeffs
3293 if (key.level()==0) {
3294 coeffT result(node.coeff());
3295 coeffT sdcoeff(cdata.v2k,this->get_tensor_type());
3296 sdcoeff(cdata.s0)+=node.coeff();
3297 node.coeff()=sdcoeff;
3298 double snorm=node.coeff().normf();
3299 node.set_dnorm(0.0);
3300 node.set_snorm(snorm);
3301 node.set_norm_tree(snorm);
3302 return Future< std::pair<GenTensor<T>,double> >(std::make_pair(result,node.coeff().normf()));
3303
3304 } else { // this is a leaf node
3305 Future<coeffT > result(node.coeff());
3306 if (not keepleaves) node.clear_coeff();
3307
3308 auto snorm=(keepleaves) ? node.coeff().normf() : 0.0;
3309 node.set_norm_tree(snorm);
3310 node.set_snorm(snorm);
3311 node.set_dnorm(0.0);
3313 return Future< std::pair<GenTensor<T>,double> >(std::make_pair(result,snorm));
3314 }
3315 }
3316 }
3317
3318 template <typename T, std::size_t NDIM>
3320 const keyT& key,
3321 const coordT& plotlo, const coordT& plothi, const std::vector<long>& npt,
3322 bool eval_refine) const {
3324 Tensor<T>& r = *ptr;
3325
3326 coordT h; // Increment between points in each dimension
3327 for (std::size_t i=0; i<NDIM; ++i) {
3328 if (npt[i] > 1) {
3329 h[i] = (plothi[i]-plotlo[i])/(npt[i]-1);
3330 }
3331 else {
3332 MADNESS_ASSERT(plotlo[i] == plothi[i]);
3333 h[i] = 0.0;
3335 }
3336
3337 const Level n = key.level();
3338 const Vector<Translation,NDIM>& l = key.translation();
3339 const double twon = pow(2.0,double(n));
3340 const tensorT& coeff = coeffs.find(key).get()->second.coeff().full_tensor_copy(); // Ugh!
3341 // const tensorT coeff = coeffs.find(key).get()->second.full_tensor_copy(); // Ugh!
3342 long ind[NDIM];
3343 coordT x;
3344
3345 coordT boxlo, boxhi;
3346 Vector<int,NDIM> boxnpt;
3347 double fac = pow(0.5,double(key.level()));
3348 int npttotal = 1;
3349 for (std::size_t d=0; d<NDIM; ++d) {
3350 // Coords of box
3351 boxlo[d] = fac*key.translation()[d];
3352 boxhi[d] = boxlo[d]+fac;
3354 if (boxlo[d] > plothi[d] || boxhi[d] < plotlo[d]) {
3355 // Discard boxes out of the plot range
3356 npttotal = boxnpt[d] = 0;
3357 //print("OO range?");
3358 break;
3359 }
3360 else if (npt[d] == 1) {
3361 // This dimension is only a single point
3362 boxlo[d] = boxhi[d] = plotlo[d];
3363 boxnpt[d] = 1;
3365 else {
3366 // Restrict to plot range
3367 boxlo[d] = std::max(boxlo[d],plotlo[d]);
3368 boxhi[d] = std::min(boxhi[d],plothi[d]);
3369
3370 // Round lo up to next plot point; round hi down
3371 double xlo = long((boxlo[d]-plotlo[d])/h[d])*h[d] + plotlo[d];
3372 if (xlo < boxlo[d]) xlo += h[d];
3373 boxlo[d] = xlo;
3374 double xhi = long((boxhi[d]-plotlo[d])/h[d])*h[d] + plotlo[d];
3375 if (xhi > boxhi[d]) xhi -= h[d];
3376 // MADNESS_ASSERT(xhi >= xlo); // nope
3377 boxhi[d] = xhi;
3378 boxnpt[d] = long(round((boxhi[d] - boxlo[d])/h[d])) + 1;
3379 }
3380 npttotal *= boxnpt[d];
3382 //print(" box", boxlo, boxhi, boxnpt, npttotal);
3383 if (npttotal > 0) {
3384 for (IndexIterator it(boxnpt); it; ++it) {
3385 for (std::size_t d=0; d<NDIM; ++d) {
3386 double xd = boxlo[d] + it[d]*h[d]; // Sim. coords of point
3387 x[d] = twon*xd - l[d]; // Offset within box
3388 MADNESS_ASSERT(x[d]>=0.0 && x[d] <=1.0); // sanity
3389 if (npt[d] > 1) {
3390 ind[d] = long(round((xd-plotlo[d])/h[d])); // Index of plot point
3391 }
3392 else {
3393 ind[d] = 0;
3394 }
3395 MADNESS_ASSERT(ind[d]>=0 && ind[d]<npt[d]); // sanity
3396 }
3397 if (eval_refine) {
3398 r(ind) = n;
3399 }
3400 else {
3401 T tmp = eval_cube(n, x, coeff);
3402 r(ind) = tmp;
3403 //print(" eval", ind, tmp, r(ind));
3404 }
3405 }
3406 }
3407 }
3408
3409 /// Set plot_refine=true to get a plot of the refinment levels of
3410 /// the given function (defaulted to false in prototype).
3411 template <typename T, std::size_t NDIM>
3413 const coordT& plothi,
3414 const std::vector<long>& npt,
3415 const bool eval_refine) const {
3417 Tensor<T> r(NDIM, &npt[0]);
3418 //r(___) = 99.0;
3419 MADNESS_ASSERT(is_reconstructed());
3420
3421 for (typename dcT::const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {
3422 const keyT& key = it->first;
3423 const nodeT& node = it->second;
3424 if (node.has_coeff()) {
3425 woT::task(world.rank(), &implT::plot_cube_kernel,
3426 archive::archive_ptr< Tensor<T> >(&r), key, plotlo, plothi, npt, eval_refine);
3427 }
3428 }
3429
3430 // ITERATOR(r, if (r(IND) == 99.0) {print("BAD", IND); error("bad",0);});
3431
3432 world.taskq.fence();
3433 world.gop.sum(r.ptr(), r.size());
3434 world.gop.fence();
3435
3436 return r;
3437 }
3438
3439 static inline void dxprintvalue(FILE* f, const double t) {
3440 fprintf(f,"%.6e\n",t);
3441 }
3442
3443 static inline void dxprintvalue(FILE* f, const double_complex& t) {
3444 fprintf(f,"%.6e %.6e\n", t.real(), t.imag());
3445 }
3446
3447 template <typename T, std::size_t NDIM>
3449 const char* filename,
3450 const Tensor<double>& cell,
3451 const std::vector<long>& npt,
3452 bool binary) {
3454 MADNESS_ASSERT(NDIM<=6);
3455 const char* element[6] = {"lines","quads","cubes","cubes4D","cubes5D","cubes6D"};
3456
3457 function.verify();
3458 World& world = const_cast< Function<T,NDIM>& >(function).world();
3459 FILE *f=0;
3460 if (world.rank() == 0) {
3461 f = fopen(filename, "w");
3462 if (!f) MADNESS_EXCEPTION("plotdx: failed to open the plot file", 0);
3463
3464 fprintf(f,"object 1 class gridpositions counts ");
3465 for (std::size_t d=0; d<NDIM; ++d) fprintf(f," %ld",npt[d]);
3466 fprintf(f,"\n");
3467
3468 fprintf(f,"origin ");
3469 for (std::size_t d=0; d<NDIM; ++d) fprintf(f, " %.6e", cell(d,0));
3470 fprintf(f,"\n");
3471
3472 for (std::size_t d=0; d<NDIM; ++d) {
3473 fprintf(f,"delta ");
3474 for (std::size_t c=0; c<d; ++c) fprintf(f, " 0");
3475 double h = 0.0;
3476 if (npt[d]>1) h = (cell(d,1)-cell(d,0))/(npt[d]-1);
3477 fprintf(f," %.6e", h);
3478 for (std::size_t c=d+1; c<NDIM; ++c) fprintf(f, " 0");
3479 fprintf(f,"\n");
3480 }
3481 fprintf(f,"\n");
3482
3483 fprintf(f,"object 2 class gridconnections counts ");
3484 for (std::size_t d=0; d<NDIM; ++d) fprintf(f," %ld",npt[d]);
3485 fprintf(f,"\n");
3486 fprintf(f, "attribute \"element type\" string \"%s\"\n", element[NDIM-1]);
3487 fprintf(f, "attribute \"ref\" string \"positions\"\n");
3488 fprintf(f,"\n");
3489
3490 int npoint = 1;
3491 for (std::size_t d=0; d<NDIM; ++d) npoint *= npt[d];
3492 const char* iscomplex = "";
3493 if (TensorTypeData<T>::iscomplex) iscomplex = "category complex";
3494 const char* isbinary = "";
3495 if (binary) isbinary = "binary";
3496 fprintf(f,"object 3 class array type double %s rank 0 items %d %s data follows\n",
3497 iscomplex, npoint, isbinary);
3498 }
3499
3500 world.gop.fence();
3501 Tensor<T> r = function.eval_cube(cell, npt);
3502
3503 if (world.rank() == 0) {
3504 if (binary) {
3505 // This assumes that the values are double precision
3506 fflush(f);
3507 fwrite((void *) r.ptr(), sizeof(T), r.size(), f);
3508 fflush(f);
3509 }
3510 else {
3511 for (IndexIterator it(npt); it; ++it) {
3512 //fprintf(f,"%.6e\n",r(*it));
3513 dxprintvalue(f,r(*it));
3514 }
3515 }
3516 fprintf(f,"\n");
3517
3518 fprintf(f,"object \"%s\" class field\n",filename);
3519 fprintf(f,"component \"positions\" value 1\n");
3520 fprintf(f,"component \"connections\" value 2\n");
3521 fprintf(f,"component \"data\" value 3\n");
3522 fprintf(f,"\nend\n");
3523 fclose(f);
3524 }
3525 world.gop.fence();
3526 }
3527
3528 template <std::size_t NDIM>
3530 k = 6;
3531 thresh = 1e-4;
3532 initial_level = 2;
3533 special_level = 3;
3534 max_refine_level = 30;
3535 truncate_mode = 0;
3536 refine = true;
3537 autorefine = true;
3538 debug = false;
3539 truncate_on_project = true;
3540 apply_randomize = false;
3541 project_randomize = false;
3543 tt = TT_FULL;
3544 cell = make_default_cell();
3545 recompute_cell_info();
3546 set_default_pmap(world);
3547 }
3548
3549 template <std::size_t NDIM>
3551 //pmap = std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >(new WorldDCDefaultPmap< Key<NDIM> >(world));
3552 pmap = std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >(new madness::LevelPmap< Key<NDIM> >(world));
3553 //pmap = std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > >(new SimplePmap< Key<NDIM> >(world));
3554 }
3555
3556
3557 template <std::size_t NDIM>
3559 std::cout << "Function Defaults:" << std::endl;
3560 std::cout << " Dimension " << ": " << NDIM << std::endl;
3561 std::cout << " k" << ": " << k << std::endl;
3562 std::cout << " thresh" << ": " << thresh << std::endl;
3563 std::cout << " initial_level" << ": " << initial_level << std::endl;
3564 std::cout << " special_level" << ": " << special_level << std::endl;
3565 std::cout << " max_refine_level" << ": " << max_refine_level << std::endl;
3566 std::cout << " truncate_mode" << ": " << truncate_mode << std::endl;
3567 std::cout << " refine" << ": " << refine << std::endl;
3568 std::cout << " autorefine" << ": " << autorefine << std::endl;
3569 std::cout << " debug" << ": " << debug << std::endl;
3570 std::cout << " truncate_on_project" << ": " << truncate_on_project << std::endl;
3571 std::cout << " apply_randomize" << ": " << apply_randomize << std::endl;
3572 std::cout << " project_randomize" << ": " << project_randomize << std::endl;
3573 std::cout << " bc" << ": " << bc << std::endl;
3574 std::cout << " tt" << ": " << tt << std::endl;
3575 std::cout << " cell" << ": " << cell << std::endl;
3576 }
3577
3578 template <typename T, std::size_t NDIM>
3579 const FunctionCommonData<T,NDIM>* FunctionCommonData<T,NDIM>::data[MAXK] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
3580
3581 // default values match those in FunctionDefaults::set_defaults(world)
3582 template <std::size_t NDIM> int FunctionDefaults<NDIM>::k = 6;
3583 template <std::size_t NDIM> double FunctionDefaults<NDIM>::thresh = 1e-4;
3584 template <std::size_t NDIM> int FunctionDefaults<NDIM>::initial_level = 2;
3585 template <std::size_t NDIM> int FunctionDefaults<NDIM>::special_level = 3;
3586 template <std::size_t NDIM> int FunctionDefaults<NDIM>::max_refine_level = 30;
3587 template <std::size_t NDIM> int FunctionDefaults<NDIM>::truncate_mode = 0;
3588 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::refine = true;
3589 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::autorefine = true;
3590 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::debug = false;
3591 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::truncate_on_project = true;
3592 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::apply_randomize = false;
3593 template <std::size_t NDIM> bool FunctionDefaults<NDIM>::project_randomize = false;
3594 template <std::size_t NDIM> BoundaryConditions<NDIM> FunctionDefaults<NDIM>::bc = BoundaryConditions<NDIM>(BC_FREE);
3595 template <std::size_t NDIM> TensorType FunctionDefaults<NDIM>::tt = TT_FULL;
3596 template <std::size_t NDIM> Tensor<double> FunctionDefaults<NDIM>::cell = FunctionDefaults<NDIM>::make_default_cell();
3597 template <std::size_t NDIM> Tensor<double> FunctionDefaults<NDIM>::cell_width = FunctionDefaults<NDIM>::make_default_cell_width();
3598 template <std::size_t NDIM> Tensor<double> FunctionDefaults<NDIM>::rcell_width = FunctionDefaults<NDIM>::make_default_cell_width();
3599 template <std::size_t NDIM> double FunctionDefaults<NDIM>::cell_volume = 1.;
3600 template <std::size_t NDIM> double FunctionDefaults<NDIM>::cell_min_width = 1.;
3601 template <std::size_t NDIM> std::shared_ptr< WorldDCPmapInterface< Key<NDIM> > > FunctionDefaults<NDIM>::pmap;
3602
3603 template <std::size_t NDIM> std::vector< Key<NDIM> > Displacements<NDIM>::disp;
3604 template <std::size_t NDIM> std::vector< Key<NDIM> > Displacements<NDIM>::disp_periodicsum[64];
3605
3606}
3607
3608#endif // MADNESS_MRA_MRAIMPL_H__INCLUDED
double w(double t, double eps)
Definition DKops.h:22
double q(double t)
Definition DKops.h:18
std::complex< double > double_complex
Definition cfft.h:14
Definition test_ar.cc:118
Definition test_ar.cc:141
Definition test_tree.cc:78
long dim(int i) const
Returns the size of dimension i.
Definition basetensor.h:147
long size() const
Returns the number of elements in the tensor.
Definition basetensor.h:138
This class is used to specify boundary conditions for all operators.
Definition funcdefaults.h:101
std::vector< bool > is_periodic() const
Convenience for application of integral operators.
Definition funcdefaults.h:166
a class to track where relevant (parent) coeffs are
Definition funcimpl.h:788
Tri-diagonal operator traversing tree primarily for derivative operator.
Definition derivative.h:73
static std::vector< Key< NDIM > > disp
Definition displacements.h:41
static std::vector< Key< NDIM > > disp_periodicsum[64]
Definition displacements.h:42
ElementaryInterface (formerly FunctorInterfaceWrapper) interfaces a c-function.
Definition function_interface.h:266
FunctionCommonData holds all Function data common for given k.
Definition function_common_data.h:52
static const FunctionCommonData< T, NDIM > & get(int k)
Definition function_common_data.h:111
static void _init_quadrature(int k, int npt, Tensor< double > &quad_x, Tensor< double > &quad_w, Tensor< double > &quad_phi, Tensor< double > &quad_phiw, Tensor< double > &quad_phit)
Initialize the quadrature information.
Definition mraimpl.h:88
void _init_twoscale()
Private. Initialize the twoscale coefficients.
Definition mraimpl.h:67
FunctionDefaults holds default paramaters as static class members.
Definition funcdefaults.h:204
static Tensor< double > make_default_cell_width()
Definition funcdefaults.h:239
static bool truncate_on_project
If true initial projection inserts at n-1 not n.
Definition funcdefaults.h:218
static bool apply_randomize
If true use randomization for load balancing in apply integral operator.
Definition funcdefaults.h:219
static double cell_volume
Volume of simulation cell.
Definition funcdefaults.h:225
static int k
Wavelet order.
Definition funcdefaults.h:209
static int truncate_mode
Truncation method.
Definition funcdefaults.h:214
static void set_default_pmap(World &world)
Sets the default process map.
Definition mraimpl.h:3550
static Tensor< double > cell_width
Width of simulation cell in each dimension.
Definition funcdefaults.h:223
static BoundaryConditions< NDIM > bc
Default boundary conditions.
Definition funcdefaults.h:221
static double thresh
Truncation threshold.
Definition funcdefaults.h:210
static bool debug
Controls output of debug info.
Definition funcdefaults.h:217
static int special_level
Minimum level for fine scale projection of special boxes.
Definition funcdefaults.h:212
static double cell_min_width
Size of smallest dimension.
Definition funcdefaults.h:226
static std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > pmap
Default mapping of keys to processes.
Definition funcdefaults.h:228
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition funcdefaults.h:468
static const BoundaryConditions< NDIM > & get_bc()
Returns the default boundary conditions.
Definition funcdefaults.h:413
static Tensor< double > cell
cell[NDIM][2] Simulation cell, cell(0,0)=xlo, cell(0,1)=xhi, ...
Definition funcdefaults.h:222
static int initial_level
Initial level for fine scale projection.
Definition funcdefaults.h:211
static bool autorefine
Whether to autorefine in multiplication, etc.
Definition funcdefaults.h:216
static Tensor< double > make_default_cell()
Definition funcdefaults.h:230
static void print()
Definition mraimpl.h:3558
static Tensor< double > rcell_width
Reciprocal of width.
Definition funcdefaults.h:224
static double get_cell_min_width()
Returns the minimum width of any user cell dimension.
Definition funcdefaults.h:478
static bool project_randomize
If true use randomization for load balancing in project/refine.
Definition funcdefaults.h:220
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition funcdefaults.h:446
static void set_defaults(World &world)
Used to set defaults to k=7, thresh=1-5, for a unit cube [0,1].
Definition mraimpl.h:3529
static int max_refine_level
Level at which to stop refinement.
Definition funcdefaults.h:213
static TensorType tt
structure of the tensor in FunctionNode
Definition funcdefaults.h:227
static bool refine
Whether to refine new functions.
Definition funcdefaults.h:215
Abstract base class interface required for functors used as input to Functions.
Definition function_interface.h:68
Definition funcimpl.h:5321
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition funcimpl.h:942
bool is_nonstandard() const
Definition mraimpl.h:252
T eval_cube(Level n, coordT &x, const tensorT &c) const
Definition mraimpl.h:2055
void do_print_tree_graphviz(const keyT &key, std::ostream &os, Level maxlevel) const
Functor for the do_print_tree method (using GraphViz)
Definition mraimpl.h:2789
void change_tensor_type1(const TensorArgs &targs, bool fence)
change the tensor type of the coefficients in the FunctionNode
Definition mraimpl.h:1078
void evaldepthpt(const Vector< double, NDIM > &xin, const keyT &keyin, const typename Future< Level >::remote_refT &ref)
Get the depth of the tree at a point in simulation coordinates.
Definition mraimpl.h:2970
void scale_inplace(const T q, bool fence)
In-place scale by a constant.
Definition mraimpl.h:3141
void gaxpy_oop_reconstructed(const double alpha, const implT &f, const double beta, const implT &g, const bool fence)
perform: this= alpha*f + beta*g, invoked by result
Definition mraimpl.h:208
bool is_redundant() const
Returns true if the function is redundant.
Definition mraimpl.h:247
FunctionNode< Q, NDIM > nodeT
Type of node.
Definition funcimpl.h:952
void print_size(const std::string name) const
print tree size and size
Definition mraimpl.h:1971
void print_info() const
Prints summary of data distribution.
Definition mraimpl.h:812
void abs_inplace(bool fence)
Definition mraimpl.h:3153
void print_timer() const
Definition mraimpl.h:336
void evalR(const Vector< double, NDIM > &xin, const keyT &keyin, const typename Future< long >::remote_refT &ref)
Get the rank of leaf box of the tree at a point in simulation coordinates.
Definition mraimpl.h:3012
const FunctionCommonData< T, NDIM > & cdata
Definition funcimpl.h:980
void do_print_grid(const std::string filename, const std::vector< keyT > &keys) const
print the grid in xyz format
Definition mraimpl.h:563
std::size_t nCoeff() const
Returns the number of coefficients in the function ... collective global sum.
Definition mraimpl.h:1957
void compress(const TreeState newstate, bool fence)
compress the wave function
Definition mraimpl.h:1522
std::pair< coeffT, double > compress_op(const keyT &key, const std::vector< Future< std::pair< coeffT, double > > > &v, bool nonstandard)
calculate the wavelet coefficients using the sum coefficients of all child nodes
Definition mraimpl.h:1690
Future< bool > truncate_spawn(const keyT &key, double tol)
Returns true if after truncation this node has coefficients.
Definition mraimpl.h:2634
Future< double > norm_tree_spawn(const keyT &key)
Definition mraimpl.h:1592
std::vector< keyT > local_leaf_keys() const
return the keys of the local leaf boxes
Definition mraimpl.h:537
void do_print_tree(const keyT &key, std::ostream &os, Level maxlevel) const
Functor for the do_print_tree method.
Definition mraimpl.h:2707
void unset_functor()
Definition mraimpl.h:291
void do_print_plane(const std::string filename, std::vector< Tensor< double > > plotinfo, const int xaxis, const int yaxis, const coordT el2)
print the MRA structure
Definition mraimpl.h:478
std::pair< Key< NDIM >, ShallowNode< T, NDIM > > find_datum(keyT key) const
return the a std::pair<key, node>, which MUST exist
Definition mraimpl.h:944
void set_functor(const std::shared_ptr< FunctionFunctorInterface< T, NDIM > > functor1)
Definition mraimpl.h:272
void broaden(std::vector< bool > is_periodic, bool fence)
Definition mraimpl.h:1282
const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > & get_pmap() const
Definition mraimpl.h:192
void sock_it_to_me(const keyT &key, const RemoteReference< FutureImpl< std::pair< keyT, coeffT > > > &ref) const
Walk up the tree returning pair(key,node) for first node with coefficients.
Definition mraimpl.h:2847
double get_thresh() const
Definition mraimpl.h:307
void trickle_down(bool fence)
sum all the contributions from all scales after applying an operator in mod-NS form
Definition mraimpl.h:1333
std::pair< coeffT, double > make_redundant_op(const keyT &key, const std::vector< Future< std::pair< coeffT, double > > > &v)
similar to compress_op, but insert only the sum coefficients in the tree
Definition mraimpl.h:1750
void set_autorefine(bool value)
Definition mraimpl.h:316
tensorT filter(const tensorT &s) const
Transform sum coefficients at level n to sums+differences at level n-1.
Definition mraimpl.h:1131
void chop_at_level(const int n, const bool fence=true)
remove all nodes with level higher than n
Definition mraimpl.h:1094
void print_tree_json(std::ostream &os=std::cout, Level maxlevel=10000) const
Definition mraimpl.h:2727
coeffT parent_to_child_NS(const keyT &child, const keyT &parent, const coeffT &coeff) const
Directly project parent NS coeffs to child NS coeffs.
Definition mraimpl.h:686
void mapdim(const implT &f, const std::vector< long > &map, bool fence)
Permute the dimensions of f according to map, result on this.
Definition mraimpl.h:1036
bool is_compressed() const
Returns true if the function is compressed.
Definition mraimpl.h:235
void mirror(const implT &f, const std::vector< long > &mirror, bool fence)
mirror the dimensions of f according to map, result on this
Definition mraimpl.h:1045
void print_stats() const
print the number of configurations per node
Definition mraimpl.h:1999
coeffT truncate_reconstructed_op(const keyT &key, const std::vector< Future< coeffT > > &v, const double tol)
given the sum coefficients of all children, truncate or not
Definition mraimpl.h:1639
void fcube(const keyT &key, const FunctionFunctorInterface< T, NDIM > &f, const Tensor< double > &qx, tensorT &fval) const
Evaluate function at quadrature points in the specified box.
Definition mraimpl.h:2472
void forward_do_diff1(const DerivativeBase< T, NDIM > *D, const implT *f, const keyT &key, const std::pair< keyT, coeffT > &left, const std::pair< keyT, coeffT > &center, const std::pair< keyT, coeffT > &right)
Definition mraimpl.h:902
std::vector< Slice > child_patch(const keyT &child) const
Returns patch referring to coeffs of child in parent box.
Definition mraimpl.h:675
void print_tree_graphviz(std::ostream &os=std::cout, Level maxlevel=10000) const
Definition mraimpl.h:2780
std::size_t min_nodes() const
Returns the min number of nodes on a processor.
Definition mraimpl.h:1895
void make_redundant(const bool fence)
convert this to redundant, i.e. have sum coefficients on all levels
Definition mraimpl.h:1550
std::size_t max_nodes() const
Returns the max number of nodes on a processor.
Definition mraimpl.h:1886
coeffT upsample(const keyT &key, const coeffT &coeff) const
upsample the sum coefficients of level 1 to sum coeffs on level n+1
Definition mraimpl.h:1210
void flo_unary_op_node_inplace(const opT &op, bool fence)
Definition funcimpl.h:2137
void plot_cube_kernel(archive::archive_ptr< Tensor< T > > ptr, const keyT &key, const coordT &plotlo, const coordT &plothi, const std::vector< long > &npt, bool eval_refine) const
Definition mraimpl.h:3319
T trace_local() const
Returns int(f(x),x) in local volume.
Definition mraimpl.h:3195
void print_grid(const std::string filename) const
Definition mraimpl.h:521
Future< std::pair< coeffT, double > > compress_spawn(const keyT &key, bool nonstandard, bool keepleaves, bool redundant1)
Invoked on node where key is local.
Definition mraimpl.h:3267
bool get_autorefine() const
Definition mraimpl.h:313
void phi_for_mul(Level np, Translation lp, Level nc, Translation lc, Tensor< double > &phi) const
Compute the Legendre scaling functions for multiplication.
Definition mraimpl.h:3163
Future< std::pair< keyT, coeffT > > find_me(const keyT &key) const
find_me. Called by diff_bdry to get coefficients of boundary function
Definition mraimpl.h:3254
TensorType get_tensor_type() const
Definition mraimpl.h:298
void remove_leaf_coefficients(const bool fence)
Definition mraimpl.h:1544
void insert_zero_down_to_initial_level(const keyT &key)
Initialize nodes to zero function at initial_level of refinement.
Definition mraimpl.h:2603
void do_diff1(const DerivativeBase< T, NDIM > *D, const implT *f, const keyT &key, const std::pair< keyT, coeffT > &left, const std::pair< keyT, coeffT > &center, const std::pair< keyT, coeffT > &right)
Definition mraimpl.h:913
void do_print_tree_json(const keyT &key, std::multimap< Level, std::tuple< tranT, std::string > > &data, Level maxlevel) const
Functor for the do_print_tree_json method.
Definition mraimpl.h:2758
void finalize_sum()
after summing up we need to do some cleanup;
Definition mraimpl.h:1843
dcT coeffs
The coefficients.
Definition funcimpl.h:985
bool exists_and_is_leaf(const keyT &key) const
Definition mraimpl.h:1254
void verify_tree() const
Verify tree is properly constructed ... global synchronization involved.
Definition mraimpl.h:109
void set_tensor_args(const TensorArgs &t)
Definition mraimpl.h:304
Range< typename dcT::const_iterator > rangeT
Definition funcimpl.h:5450
std::size_t real_size() const
Returns the number of coefficients in the function ... collective global sum.
Definition mraimpl.h:1944
bool exists_and_has_children(const keyT &key) const
Definition mraimpl.h:1249
void sum_down_spawn(const keyT &key, const coeffT &s)
is this the same as trickle_down() ?
Definition mraimpl.h:855
void norm_tree(bool fence)
compute for each FunctionNode the norm of the function inside that node
Definition mraimpl.h:1569
bool has_leaves() const
Definition mraimpl.h:267
void reconstruct_op(const keyT &key, const coeffT &s, const bool accumulate_NS=true)
Definition mraimpl.h:2109
const coeffT parent_to_child(const coeffT &s, const keyT &parent, const keyT &child) const
Directly project parent coeffs to child coeffs.
Definition mraimpl.h:3178
void undo_redundant(const bool fence)
convert this from redundant to standard reconstructed form
Definition mraimpl.h:1560
GenTensor< Q > coeffT
Type of tensor used to hold coeffs.
Definition funcimpl.h:953
const keyT & key0() const
Returns cdata.key0.
Definition mraimpl.h:373
double finalize_apply()
after apply we need to do some cleanup;
Definition mraimpl.h:1800
const dcT & get_coeffs() const
Definition mraimpl.h:322
double norm2sq_local() const
Returns the square of the local norm ... no comms.
Definition mraimpl.h:1852
const FunctionCommonData< T, NDIM > & get_cdata() const
Definition mraimpl.h:328
void sum_down(bool fence)
After 1d push operator must sum coeffs down the tree to restore correct scaling function coefficients...
Definition mraimpl.h:894
bool noautorefine(const keyT &key, const tensorT &t) const
Always returns false (for when autorefine is not wanted)
Definition mraimpl.h:838
double truncate_tol(double tol, const keyT &key) const
Returns the truncation threshold according to truncate_method.
Definition mraimpl.h:628
bool autorefine_square_test(const keyT &key, const nodeT &t) const
Returns true if this block of coeffs needs autorefining.
Definition mraimpl.h:844
void erase(const Level &max_level)
truncate tree at a certain level
Definition mraimpl.h:718
void reconstruct(bool fence)
reconstruct this tree – respects fence
Definition mraimpl.h:1490
void multiply(const implT *f, const FunctionImpl< T, LDIM > *g, const int particle)
multiply f (a pair function of NDIM) with an orbital g (LDIM=NDIM/2)
Definition funcimpl.h:3555
coeffT assemble_coefficients(const keyT &key, const coeffT &coeff_ket, const coeffT &vpotential1, const coeffT &vpotential2, const tensorT &veri) const
given several coefficient tensors, assemble a result tensor
Definition mraimpl.h:992
static void tnorm(const tensorT &t, double *lo, double *hi)
Computes norm of low/high-order polyn. coeffs for autorefinement test.
Definition mraimpl.h:3055
std::pair< bool, T > eval_local_only(const Vector< double, NDIM > &xin, Level maxlevel)
Evaluate function only if point is local returning (true,value); otherwise return (false,...
Definition mraimpl.h:2941
std::size_t max_depth() const
Returns the maximum depth of the tree ... collective ... global sum/broadcast.
Definition mraimpl.h:1878
std::size_t size() const
Returns the number of coefficients in the function ... collective global sum.
Definition mraimpl.h:1913
void reduce_rank(const double thresh, bool fence)
reduce the rank of the coefficients tensors
Definition mraimpl.h:1086
std::shared_ptr< FunctionFunctorInterface< T, NDIM > > get_functor()
Definition mraimpl.h:279
tensorT unfilter(const tensorT &s) const
Transform sums+differences at level n to sum coefficients at level n+1.
Definition mraimpl.h:1160
Tensor< T > eval_plot_cube(const coordT &plotlo, const coordT &plothi, const std::vector< long > &npt, const bool eval_refine=false) const
Definition mraimpl.h:3412
void change_tree_state(const TreeState finalstate, bool fence=true)
change the tree state of this function, might or might not respect fence!
Definition mraimpl.h:1386
Future< coeffT > truncate_reconstructed_spawn(const keyT &key, const double tol)
truncate using a tree in reconstructed form
Definition mraimpl.h:1615
void map_and_mirror(const implT &f, const std::vector< long > &map, const std::vector< long > &mirror, bool fence)
map and mirror the translation index and the coefficients, result on this
Definition mraimpl.h:1055
void truncate(double tol, bool fence)
Truncate according to the threshold with optional global fence.
Definition mraimpl.h:357
bool is_reconstructed() const
Returns true if the function is compressed.
Definition mraimpl.h:241
double norm_tree_op(const keyT &key, const std::vector< Future< double > > &v)
Definition mraimpl.h:1577
void reset_timer()
Definition mraimpl.h:345
void refine_to_common_level(const std::vector< FunctionImpl< T, NDIM > * > &v, const std::vector< tensorT > &c, const keyT key)
Refine multiple functions down to the same finest level.
Definition mraimpl.h:748
int get_k() const
Definition mraimpl.h:319
void eval(const Vector< double, NDIM > &xin, const keyT &keyin, const typename Future< T >::remote_refT &ref)
Evaluate the function at a point in simulation coordinates.
Definition mraimpl.h:2897
bool truncate_op(const keyT &key, double tol, const std::vector< Future< bool > > &v)
Definition mraimpl.h:2670
void zero_norm_tree()
Definition mraimpl.h:1271
std::size_t max_local_depth() const
Returns the maximum local depth of the tree ... no communications.
Definition mraimpl.h:1864
tensorT project(const keyT &key) const
Definition mraimpl.h:2815
double check_symmetry_local() const
Returns some asymmetry measure ... no comms.
Definition mraimpl.h:734
Future< double > get_norm_tree_recursive(const keyT &key) const
Definition mraimpl.h:2836
Key< NDIM > keyT
Type of key.
Definition funcimpl.h:951
bool is_on_demand() const
Definition mraimpl.h:262
void accumulate_timer(const double time) const
Definition mraimpl.h:331
void trickle_down_op(const keyT &key, const coeffT &s)
sum all the contributions from all scales after applying an operator in mod-NS form
Definition mraimpl.h:1344
void set_thresh(double value)
Definition mraimpl.h:310
Tensor< double > print_plane_local(const int xaxis, const int yaxis, const coordT &el2)
collect the data for a plot of the MRA structure locally on each node
Definition mraimpl.h:402
void sock_it_to_me_too(const keyT &key, const RemoteReference< FutureImpl< std::pair< keyT, coeffT > > > &ref) const
Definition mraimpl.h:2875
void broaden_op(const keyT &key, const std::vector< Future< bool > > &v)
Definition mraimpl.h:1260
void print_plane(const std::string filename, const int xaxis, const int yaxis, const coordT &el2)
Print a plane ("xy", "xz", or "yz") containing the point x to file.
Definition mraimpl.h:382
void print_tree(std::ostream &os=std::cout, Level maxlevel=10000) const
Definition mraimpl.h:2698
void project_refine_op(const keyT &key, bool do_refine, const std::vector< Vector< double, NDIM > > &specialpts)
Definition mraimpl.h:2484
std::size_t tree_size() const
Returns the size of the tree structure of the function ... collective global sum.
Definition mraimpl.h:1904
void add_scalar_inplace(T t, bool fence)
Adds a constant to the function. Local operation, optional fence.
Definition mraimpl.h:2562
tensorT downsample(const keyT &key, const std::vector< Future< coeffT > > &v) const
downsample the sum coefficients of level n+1 to sum coeffs on level n
Definition mraimpl.h:1180
void abs_square_inplace(bool fence)
Definition mraimpl.h:3158
void put_in_box(ProcessID from, long nl, long ni) const
Definition mraimpl.h:803
TensorArgs get_tensor_args() const
Definition mraimpl.h:301
void average(const implT &rhs)
take the average of two functions, similar to: this=0.5*(this+rhs)
Definition mraimpl.h:1067
void diff(const DerivativeBase< T, NDIM > *D, const implT *f, bool fence)
Definition mraimpl.h:925
keyT neighbor(const keyT &key, const keyT &disp, const std::vector< bool > &is_periodic) const
Returns key of general neighbor enforcing BC.
Definition mraimpl.h:3237
void square_inplace(bool fence)
Pointwise squaring of function with optional global fence.
Definition mraimpl.h:3147
void remove_internal_coefficients(const bool fence)
Definition mraimpl.h:1539
void compute_snorm_and_dnorm(bool fence=true)
compute norm of s and d coefficients for all nodes
Definition mraimpl.h:1110
void standard(bool fence)
Changes non-standard compressed form to standard compressed form.
Definition mraimpl.h:1787
bool is_nonstandard_with_leaves() const
Definition mraimpl.h:257
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition funcimpl.h:124
bool has_coeff() const
Returns true if there are coefficients in this node.
Definition funcimpl.h:197
size_t size() const
Returns the number of coefficients in this node.
Definition funcimpl.h:239
void clear_coeff()
Clears the coefficients (has_coeff() will subsequently return false)
Definition funcimpl.h:292
bool is_leaf() const
Returns true if this does not have children.
Definition funcimpl.h:210
void set_has_children(bool flag)
Sets has_children attribute to value of flag.
Definition funcimpl.h:251
void set_snorm(const double sn)
set the precomputed norm of the (virtual) s coefficients
Definition funcimpl.h:318
void set_is_leaf(bool flag)
Sets has_children attribute to value of !flag.
Definition funcimpl.h:277
void print_json(std::ostream &s) const
Definition funcimpl.h:463
bool has_children() const
Returns true if this node has children.
Definition funcimpl.h:204
void set_coeff(const coeffT &coeffs)
Takes a shallow copy of the coeff — same as this->coeff()=coeff.
Definition funcimpl.h:282
void set_dnorm(const double dn)
set the precomputed norm of the (virtual) d coefficients
Definition funcimpl.h:323
coeffT & coeff()
Returns a non-const reference to the tensor containing the coeffs.
Definition funcimpl.h:224
void set_norm_tree(double norm_tree)
Sets the value of norm_tree.
Definition funcimpl.h:303
A multiresolution adaptive numerical function.
Definition mra.h:122
Implements the functionality of futures.
Definition future.h:74
A future is a possibly yet unevaluated value.
Definition future.h:373
T & get(bool dowork=true) &
Gets the value, waiting if necessary.
Definition future.h:574
remote_refT remote_ref(World &world) const
Returns a structure used to pass references to another process.
Definition future.h:675
void set(const Future< T > &other)
A.set(B), where A and B are futures ensures A has/will have the same value as B.
Definition future.h:508
RemoteReference< FutureImpl< T > > remote_refT
Definition future.h:398
Definition lowranktensor.h:59
GenTensor convert(const TensorArgs &targs) const
Definition gentensor.h:198
GenTensor full_tensor() const
Definition gentensor.h:200
long dim(const int i) const
return the number of entries in dimension i
Definition lowranktensor.h:391
Tensor< T > full_tensor_copy() const
Definition gentensor.h:206
long ndim() const
Definition lowranktensor.h:386
constexpr bool is_full_tensor() const
Definition gentensor.h:224
size_t real_size() const
Definition gentensor.h:214
GenTensor< T > & emul(const GenTensor< T > &other)
Inplace multiply by corresponding elements of argument Tensor.
Definition lowranktensor.h:631
float_scalar_type normf() const
Definition lowranktensor.h:406
void reduce_rank(const double &eps)
Definition gentensor.h:217
long rank() const
Definition gentensor.h:212
size_t nCoeff() const
Definition gentensor.h:215
TensorType tensor_type() const
Definition gentensor.h:221
bool has_data() const
Definition gentensor.h:210
const long * dims() const
return the number of entries in dimension i
Definition lowranktensor.h:397
GenTensor & gaxpy(const T alpha, const GenTensor &other, const T beta)
Definition lowranktensor.h:580
IsSupported< TensorTypeData< Q >, GenTensor< T > & >::type scale(Q fac)
Inplace multiplication by scalar of supported type (legacy name)
Definition lowranktensor.h:426
constexpr bool is_svd_tensor() const
Definition gentensor.h:222
Definition indexit.h:141
Definition indexit.h:55
Iterates in lexical order thru all children of a key.
Definition key.h:374
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:66
Level level() const
Definition key.h:159
bool thisKeyContains(const Vector< double, NDIM > &x, const unsigned int &dim0, const unsigned int &dim1) const
check if this MultiIndex contains point x, disregarding these two dimensions
Definition key.h:245
bool is_invalid() const
Checks if a key is invalid.
Definition key.h:109
Key parent(int generation=1) const
Returns the key of the parent.
Definition key.h:187
const Vector< Translation, NDIM > & translation() const
Definition key.h:164
bool is_neighbor_of(const Key &key, const std::vector< bool > &bperiodic) const
Assuming keys are at the same level, returns true if displaced by no more than 1 in any direction.
Definition key.h:222
A pmap that locates children on odd levels with their even level parents.
Definition funcimpl.h:102
Range, vaguely a la Intel TBB, to encapsulate a random-access, STL-like start and end iterator with c...
Definition range.h:64
Simple structure used to manage references/pointers to remote instances.
Definition worldref.h:395
double weights(const unsigned int &i) const
return the weight
Definition srconf.h:671
const Tensor< T > flat_vector(const unsigned int &idim) const
return shallow copy of a slice of one of the vectors, flattened to (r,kVec)
Definition srconf.h:545
Definition SVDTensor.h:42
long rank() const
Definition SVDTensor.h:77
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
static TaskAttributes hipri()
Definition thread.h:450
static TaskAttributes generator()
Definition thread.h:442
Traits class to specify support of numeric types.
Definition type_data.h:56
A tensor is a multidimension array.
Definition tensor.h:317
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
Tensor< T > & fill(T x)
Inplace fill with a scalar (legacy name)
Definition tensor.h:562
T sum() const
Returns the sum of all elements of the tensor.
Definition tensor.h:1662
Tensor< T > reshape(int ndimnew, const long *d)
Returns new view/tensor reshaping size/number of dimensions to conforming tensor.
Definition tensor.h:1384
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1824
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition tensor.h:686
Tensor< T > & emul(const Tensor< T > &t)
Inplace multiply by corresponding elements of argument Tensor.
Definition tensor.h:1798
bool has_data() const
Definition tensor.h:1886
iterator begin()
Returns an iterator to the beginning of the local data (no communication)
Definition worlddc.h:1070
implT::const_iterator const_iterator
Definition worlddc.h:872
iterator end()
Returns an iterator past the end of the local data (no communication)
Definition worlddc.h:1084
Future< iterator > futureT
Definition worlddc.h:875
implT::iterator iterator
Definition worlddc.h:871
implT::accessor accessor
Definition worlddc.h:873
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition worldgop.cc:161
A parallel world class.
Definition world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:318
WorldGopInterface & gop
Global operations.
Definition world.h:205
Wrapper for an opaque pointer for serialization purposes.
Definition archive.h:850
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
const std::size_t bufsize
Definition derivatives.cc:16
static double lo
Definition dirac-hatom.cc:23
static bool debug
Definition dirac-hatom.cc:16
Provides FunctionCommonData, FunctionImpl and FunctionFactory.
static double function(const coord_3d &r)
Normalized gaussian.
Definition functionio.cc:100
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
Tensor< TENSOR_RESULT_TYPE(T, Q) > & fast_transform(const Tensor< T > &t, const Tensor< Q > &c, Tensor< TENSOR_RESULT_TYPE(T, Q) > &result, Tensor< TENSOR_RESULT_TYPE(T, Q) > &workspace)
Restricted but heavily optimized form of transform()
Definition tensor.h:2443
const double beta
Definition gygi_soltion.cc:62
static const double v
Definition hatom_sf_dirac.cc:20
static double u(double r, double c)
Definition he.cc:20
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
static double pow(const double *a, const double *b)
Definition lda.h:74
#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
Vector< double, 3 > coordT
Definition mcpfit.cc:48
static const bool VERIFY_TREE
Definition mra.h:57
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
bool two_scale_hg(int k, Tensor< double > *hg)
Definition twoscale.cc:151
@ BC_FREE
Definition funcdefaults.h:56
void make_redundant(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
change tree_state of a vector of functions to redundant
Definition vmra.h:189
static const char * filename
Definition legendre.cc:96
static const std::vector< Slice > ___
Entire dimension.
Definition slice.h:128
double abs(double x)
Definition complexfun.h:48
static double cpu_time()
Returns the cpu time in seconds relative to an arbitrary origin.
Definition timers.h:127
Vector< double, 3 > coordT
Definition corepotential.cc:54
GenTensor< TENSOR_RESULT_TYPE(R, Q)> general_transform(const GenTensor< R > &t, const Tensor< Q > c[])
Definition gentensor.h:274
response_space scale(response_space a, double b)
void legendre_scaling_functions(double x, long k, double *p)
Evaluate the first k Legendre scaling functions.
Definition legendre.cc:85
void norm_tree(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Makes the norm tree for all functions in a vector.
Definition vmra.h:1147
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition vmra.h:689
static Key< NDIM > simpt2key(const Vector< T, NDIM > &pt, Level n)
Definition funcdefaults.h:533
TreeState
Definition funcdefaults.h:58
@ nonstandard_after_apply
s and d coeffs, state after operator application
Definition funcdefaults.h:63
@ on_demand
no coeffs anywhere, but a functor providing if necessary
Definition funcdefaults.h:66
@ reconstructed
s coeffs at the leaves only
Definition funcdefaults.h:59
@ nonstandard
s and d coeffs in internal nodes
Definition funcdefaults.h:61
@ unknown
Definition funcdefaults.h:67
@ compressed
d coeffs in internal nodes, s and d coeffs at the root
Definition funcdefaults.h:60
@ redundant
s coeffs everywhere
Definition funcdefaults.h:64
@ nonstandard_with_leaves
like nonstandard, with s coeffs at the leaves
Definition funcdefaults.h:62
static void user_to_sim(const Vector< double, NDIM > &xuser, Vector< double, NDIM > &xsim)
Convert user coords (cell[][]) to simulation coords ([0,1]^ndim)
Definition funcdefaults.h:524
void standard(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Generates standard form of a vector of functions.
Definition vmra.h:260
Tensor< double > tensorT
Definition distpm.cc:21
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition vmra.h:133
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
response_space transpose(response_space &f)
Definition basic_operators.cc:10
int64_t Translation
Definition key.h:54
void plotdx(const Function< T, NDIM > &f, const char *filename, const Tensor< double > &cell=FunctionDefaults< NDIM >::get_cell(), const std::vector< long > &npt=std::vector< long >(NDIM, 201L), bool binary=true)
Writes an OpenDX format file with a cube/slice of points on a uniform grid.
Definition mraimpl.h:3448
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
static void verify_tree(World &world, const std::vector< Function< T, NDIM > > &v)
Definition SCF.cc:72
static const Slice _(0,-1, 1)
void change_tensor_type(GenTensor< T > &t, const TensorArgs &targs)
change representation to targ.tt
Definition gentensor.h:284
int Level
Definition key.h:55
std::enable_if< std::is_base_of< ProjectorBase, projT >::value, OuterProjector< projT, projQ > >::type outer(const projT &p0, const projQ &p1)
Definition projector.h:457
bool gauss_legendre(int n, double xlo, double xhi, double *x, double *w)
Definition legendre.cc:226
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
Tensor< T > fcube(const Key< NDIM > &, T(*f)(const Vector< double, NDIM > &), const Tensor< double > &)
Definition mraimpl.h:2163
TensorType
low rank representations of tensors (see gentensor.h)
Definition gentensor.h:120
@ TT_2D
Definition gentensor.h:120
@ TT_FULL
Definition gentensor.h:120
static void dxprintvalue(FILE *f, const double t)
Definition mraimpl.h:3439
void refine(World &world, const std::vector< Function< T, NDIM > > &vf, bool fence=true)
refine the functions according to the autorefine criteria
Definition vmra.h:208
NDIM & f
Definition mra.h:2416
void error(const char *msg)
Definition world.cc:139
NDIM const Function< R, NDIM > & g
Definition mra.h:2416
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition timers.cc:48
static bool print_timings
Definition SCF.cc:104
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
static std::vector< double > ttt
Definition SCF.cc:105
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
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
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
Vector< T, sizeof...(Ts)+1 > vec(T t, Ts... ts)
Factory function for creating a madness::Vector.
Definition vector.h:711
static const int MAXK
The maximum wavelet order presently supported.
Definition funcdefaults.h:51
static bool enforce_bc(bool is_periodic, Level n, Translation &l)
Definition mraimpl.h:3219
Definition mraimpl.h:50
bool isnan(const std::complex< T > &v)
Definition mraimpl.h:52
const double mu
Definition navstokes_cosines.cc:95
static const double b
Definition nonlinschro.cc:119
static const double d
Definition nonlinschro.cc:121
static const double a
Definition nonlinschro.cc:118
static const double c
Definition relops.cc:10
static const double m
Definition relops.cc:9
static const double L
Definition rk.cc:46
static const double thresh
Definition rk.cc:45
static const long k
Definition rk.cc:44
const double xi
Exponent for delta function approx.
Definition siam_example.cc:60
Definition test_ar.cc:204
Definition test_ccpairfunction.cc:22
add two functions f and g: result=alpha * f + beta * g
Definition funcimpl.h:3494
"put" this on g
Definition funcimpl.h:2562
change representation of nodes' coeffs to low rank, optional fence
Definition funcimpl.h:2595
check symmetry wrt particle exchange
Definition funcimpl.h:2268
compute the norm of the wavelet coefficients
Definition funcimpl.h:4386
mirror dimensions of this, write result on f
Definition funcimpl.h:2496
map this on f
Definition funcimpl.h:2416
mirror dimensions of this, write result on f
Definition funcimpl.h:2446
Definition funcimpl.h:5375
reduce the rank of the nodes, optional fence
Definition funcimpl.h:2242
Changes non-standard compressed form to standard compressed form.
Definition funcimpl.h:4607
given an NS tree resulting from a convolution, truncate leafs if appropriate
Definition funcimpl.h:2163
remove all coefficients of internal nodes
Definition funcimpl.h:2188
remove all coefficients of leaf nodes
Definition funcimpl.h:2205
shallow-copy, pared-down version of FunctionNode, for special purpose only
Definition funcimpl.h:746
TensorArgs holds the arguments for creating a LowRankTensor.
Definition gentensor.h:134
double thresh
Definition gentensor.h:135
Definition mraimpl.h:3127
void operator()(const Key< NDIM > &key, Tensor< T > &t) const
Definition mraimpl.h:3128
void serialize(Archive &ar)
Definition mraimpl.h:3129
Definition mraimpl.h:3133
void operator()(const Key< NDIM > &key, Tensor< T > &t) const
Definition mraimpl.h:3134
void serialize(Archive &ar)
Definition mraimpl.h:3135
Definition mraimpl.h:3095
void operator()(const A &a, const B &b) const
Definition mraimpl.h:3096
void serialize(Archive &ar)
Definition mraimpl.h:3098
Definition mraimpl.h:3102
void operator()(const Key< NDIM > &key, FunctionNode< T, NDIM > &node) const
Definition mraimpl.h:3110
void serialize(Archive &ar)
Definition mraimpl.h:3113
T q
Definition mraimpl.h:3103
scaleinplace()
Definition mraimpl.h:3104
scaleinplace(T q)
Definition mraimpl.h:3106
void operator()(const Key< NDIM > &key, Tensor< T > &t) const
Definition mraimpl.h:3107
Definition mraimpl.h:3119
void serialize(Archive &ar)
Definition mraimpl.h:3123
void operator()(const Key< NDIM > &key, Tensor< T > &t) const
Definition mraimpl.h:3120
insert/replaces the coefficients into the function
Definition funcimpl.h:689
Definition lowrankfunction.h:332
int np
Definition tdse1d.cc:165
static const double s0
Definition tdse4.cc:83
AtomicInt sum
Definition test_atomicint.cc:46
int me
Definition test_binsorter.cc:10
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
#define N
Definition testconv.cc:37
static const double alpha
Definition testcosine.cc:10
static const int truncate_mode
Definition testcosine.cc:14
double h(const coord_1d &r)
Definition testgconv.cc:68
static const std::size_t NDIM
Definition testpdiff.cc:42
std::size_t axis
Definition testpdiff.cc:59
double k0
Definition testperiodic.cc:66
Defines and implements WorldObject.
Implements WorldContainer.
Defines and implements a concurrent hashmap.
#define PROFILE_FUNC
Definition worldprofile.h:209
#define PROFILE_MEMBER_FUNC(classname)
Definition worldprofile.h:210
#define PROFILE_BLOCK(name)
Definition worldprofile.h:208
int ProcessID
Used to clearly identify process number/rank.
Definition worldtypes.h:43
void test()
Definition y.cc:696