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>
43 #include <madness/world/worlddc.h>
46 
47 #include <madness/mra/funcimpl.h>
49 
50 namespace 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 
60 namespace 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();
127  }
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 
217  CoeffTracker<T,NDIM> ff(&f);
218  CoeffTracker<T,NDIM> gg(&g);
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>
307  double FunctionImpl<T,NDIM>::get_thresh() const {return thresh;}
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;
431  user[0]=xloleft*FunctionDefaults<NDIM>::get_cell_width()[xaxis] + FunctionDefaults<NDIM>::get_cell()(xaxis,0);
432  user[2]=xhiright*FunctionDefaults<NDIM>::get_cell_width()[xaxis] + FunctionDefaults<NDIM>::get_cell()(xaxis,0);
433  user[1]=yloleft*FunctionDefaults<NDIM>::get_cell_width()[yaxis] + FunctionDefaults<NDIM>::get_cell()(yaxis,0);
434  user[3]=yhiright*FunctionDefaults<NDIM>::get_cell_width()[yaxis] + FunctionDefaults<NDIM>::get_cell()(yaxis,0);
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()),
738  do_check_symmetry_local(*this));
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>
855  void FunctionImpl<T,NDIM>::sum_down_spawn(const keyT& key, const coeffT& s) {
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>
944  std::pair<Key<NDIM>,ShallowNode<T,NDIM> > FunctionImpl<T,NDIM>::find_datum(keyT key) const {
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 
960  MADNESS_ASSERT(particle==0 or particle==1);
961  MADNESS_ASSERT(val_pot.is_full_tensor());
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>
1180  typename FunctionImpl<T,NDIM>::tensorT FunctionImpl<T,NDIM>::downsample(const keyT& key, const std::vector< Future<coeffT > >& v) const {
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>
1210  typename FunctionImpl<T,NDIM>::coeffT FunctionImpl<T,NDIM>::upsample(const keyT& key, const coeffT& coeff) const {
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  }
1252 
1253  template <typename T, std::size_t NDIM>
1255  return coeffs.probe(key) && (not coeffs.find(key).get()->second.has_children());
1256  }
1257 
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);
1277  }
1278  }
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);
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)) {
1293 
1294  node.set_norm_tree(-1.0); // Indicates already broadened or result of broadening/refining
1295 
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);
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);
1318  }
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  }
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);
1364 
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  }
1377  }
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);
1433  }
1434  set_tree_state(nonstandard);
1435  } else if (finalstate==nonstandard_with_leaves) {
1436  if (current_state==reconstructed) compress(nonstandard_with_leaves,fence);
1437  if (current_state==compressed) {
1438  reconstruct(true);
1439  must_fence=true;
1441  }
1442  if (current_state==nonstandard) {
1443  standard(true);
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>
1550  void FunctionImpl<T,NDIM>::make_redundant(const bool fence) {
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>
1560  void FunctionImpl<T,NDIM>::undo_redundant(const bool fence) {
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);
1653  MADNESS_CHECK(found);
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;
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;
1710 
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()),
1856  do_norm2sq_local());
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>
1878  std::size_t FunctionImpl<T,NDIM>::max_depth() const {
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>
1886  std::size_t FunctionImpl<T,NDIM>::max_nodes() const {
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>
1895  std::size_t FunctionImpl<T,NDIM>::min_nodes() const {
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>
1904  std::size_t FunctionImpl<T,NDIM>::tree_size() const {
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>
1944  std::size_t FunctionImpl<T,NDIM>::real_size() const {
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  }
2160  }
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 {
2185  typedef Vector<double,NDIM> coordT;
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);
2468  madness::fcube(key,ElementaryInterface<T,NDIM>(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);
2577  node.coeff().full_tensor()(v0) += t*sqrt(FunctionDefaults<NDIM>::get_cell_volume());
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  }
2722  }
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>
3102  struct scaleinplace {
3103  T q;
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>
3119  struct squareinplace {
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 
3140 template <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  }
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];
3242 
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>
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;
3274 
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()));
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);
3312 
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 {
3323 
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;
3334  }
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;
3353 
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;
3364  }
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];
3381  }
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>
3448  void plotdx(const Function<T,NDIM>& function,
3449  const char* filename,
3450  const Tensor<double>& cell,
3451  const std::vector<long>& npt,
3452  bool binary) {
3453  PROFILE_FUNC;
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:787
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 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
static const FunctionCommonData< T, NDIM > & get(int k)
Definition: function_common_data.h:111
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
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 Tensor< double > make_default_cell()
Definition: funcdefaults.h:230
static double thresh
Truncation threshold.
Definition: funcdefaults.h:210
static Tensor< double > make_default_cell_width()
Definition: funcdefaults.h:239
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 BoundaryConditions< NDIM > & get_bc()
Returns the default boundary conditions.
Definition: funcdefaults.h:413
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition: funcdefaults.h:446
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 const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
static bool autorefine
Whether to autorefine in multiplication, etc.
Definition: funcdefaults.h:216
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 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:5319
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition: funcimpl.h:941
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
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 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:979
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
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:2136
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 finalize_sum()
after summing up we need to do some cleanup;
Definition: mraimpl.h:1843
dcT coeffs
The coefficients.
Definition: funcimpl.h:984
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
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
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:3554
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::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
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
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:196
size_t size() const
Returns the number of coefficients in this node.
Definition: funcimpl.h:238
void clear_coeff()
Clears the coefficients (has_coeff() will subsequently return false)
Definition: funcimpl.h:291
bool is_leaf() const
Returns true if this does not have children.
Definition: funcimpl.h:209
void set_has_children(bool flag)
Sets has_children attribute to value of flag.
Definition: funcimpl.h:250
coeffT & coeff()
Returns a non-const reference to the tensor containing the coeffs.
Definition: funcimpl.h:223
double get_norm_tree() const
Gets the value of norm_tree.
Definition: funcimpl.h:307
void set_snorm(const double sn)
set the precomputed norm of the (virtual) s coefficients
Definition: funcimpl.h:317
void set_is_leaf(bool flag)
Sets has_children attribute to value of !flag.
Definition: funcimpl.h:276
void print_json(std::ostream &s) const
Definition: funcimpl.h:462
bool has_children() const
Returns true if this node has children.
Definition: funcimpl.h:203
void set_coeff(const coeffT &coeffs)
Takes a shallow copy of the coeff — same as this->coeff()=coeff.
Definition: funcimpl.h:281
void set_dnorm(const double dn)
set the precomputed norm of the (virtual) d coefficients
Definition: funcimpl.h:322
void set_norm_tree(double norm_tree)
Sets the value of norm_tree.
Definition: funcimpl.h:302
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
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
T & get(bool dowork=true) &
Gets the value, waiting if necessary.
Definition: future.h:574
Definition: lowranktensor.h:59
GenTensor< T > & emul(const GenTensor< T > &other)
Inplace multiply by corresponding elements of argument Tensor.
Definition: lowranktensor.h:631
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
long ndim() const
Definition: lowranktensor.h:386
constexpr bool is_full_tensor() const
Definition: gentensor.h:224
IsSupported< TensorTypeData< Q >, GenTensor< T > & >::type scale(Q fac)
Inplace multiplication by scalar of supported type (legacy name)
Definition: lowranktensor.h:426
size_t real_size() const
Definition: gentensor.h:214
Tensor< T > full_tensor_copy() const
Definition: gentensor.h:206
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
constexpr bool is_svd_tensor() const
Definition: gentensor.h:222
Definition: worldhashmap.h:330
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
const Vector< Translation, NDIM > & translation() const
Definition: key.h:164
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
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
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
double weights(const unsigned int &i) const
return the weight
Definition: srconf.h:671
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
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
T sum() const
Returns the sum of all elements of the tensor.
Definition: tensor.h:1662
Tensor< T > & fill(T x)
Inplace fill with a scalar (legacy name)
Definition: tensor.h:562
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
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
const TensorIterator< T > & end() const
End point for forward iteration.
Definition: tensor.h:1876
Iterator for distributed container wraps the local iterator.
Definition: worlddc.h:244
iterator begin()
Returns an iterator to the beginning of the local data (no communication)
Definition: worlddc.h:1070
iterator end()
Returns an iterator past the end of the local data (no communication)
Definition: worlddc.h:1084
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
void test(World &world, bool doloadbal=false)
Definition: dataloadbal.cc:224
const std::size_t bufsize
Definition: derivatives.cc:16
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition: derivatives.cc:72
static double lo
Definition: dirac-hatom.cc:23
static bool debug
Definition: dirac-hatom.cc:16
Provides FunctionCommonData, FunctionImpl and FunctionFactory.
const double m
Definition: gfit.cc:199
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
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 max(a, b)
Definition: lda.h:51
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:190
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition: madness_exception.h:210
Vector< double, 3 > coordT
Definition: mcpfit.cc:48
static const bool VERIFY_TREE
Definition: mra.h:57
double norm(const T &t)
Definition: adquad.h:42
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
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
static Key< NDIM > simpt2key(const Vector< T, NDIM > &pt, Level n)
Definition: funcdefaults.h:533
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
Function< T, NDIM > multiply(const Function< T, NDIM > f, const Function< T, LDIM > g, const int particle, const bool fence=true)
multiply a high-dimensional function with a low-dimensional function
Definition: mra.h:2372
void 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
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
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
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
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
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
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
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
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
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
void refine_to_common_level(World &world, std::vector< Function< T, NDIM > > &vf, bool fence=true)
refine all functions to a common (finest) level
Definition: vmra.h:218
static bool print_timings
Definition: SCF.cc:104
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 std::vector< double > ttt
Definition: SCF.cc:105
std::string name(const FuncType &type, const int ex=-1)
Definition: ccpairfunction.h:28
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
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
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
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
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 a
Definition: nonlinschro.cc:118
static const double c
Definition: relops.cc:10
static const double L
Definition: rk.cc:46
static const double thresh
Definition: rk.cc:45
static const long k
Definition: rk.cc:44
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:3493
"put" this on g
Definition: funcimpl.h:2561
change representation of nodes' coeffs to low rank, optional fence
Definition: funcimpl.h:2594
check symmetry wrt particle exchange
Definition: funcimpl.h:2267
compute the norm of the wavelet coefficients
Definition: funcimpl.h:4385
mirror dimensions of this, write result on f
Definition: funcimpl.h:2495
map this on f
Definition: funcimpl.h:2415
mirror dimensions of this, write result on f
Definition: funcimpl.h:2445
Definition: funcimpl.h:5373
reduce the rank of the nodes, optional fence
Definition: funcimpl.h:2241
Changes non-standard compressed form to standard compressed form.
Definition: funcimpl.h:4606
given an NS tree resulting from a convolution, truncate leafs if appropriate
Definition: funcimpl.h:2162
remove all coefficients of internal nodes
Definition: funcimpl.h:2187
remove all coefficients of leaf nodes
Definition: funcimpl.h:2204
shallow-copy, pared-down version of FunctionNode, for special purpose only
Definition: funcimpl.h:745
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:688
Definition: lowrankfunction.h:332
int np
Definition: tdse1d.cc:165
static const double s0
Definition: tdse4.cc:83
int me
Definition: test_binsorter.cc:10
const double alpha
Definition: test_coulomb.cc:54
int task(int i)
Definition: test_runtime.cpp:4
void d()
Definition: test_sig.cc:79
void e()
Definition: test_sig.cc:75
#define N
Definition: testconv.cc:37
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
double u(const double x, const double expnt)
Definition: testperiodic.cc:56
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