MADNESS  0.10.1
vmra.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  $Id$
32 */
33 #ifndef MADNESS_MRA_VMRA_H__INCLUDED
34 #define MADNESS_MRA_VMRA_H__INCLUDED
35 
36 /*!
37  \file vmra.h
38  \brief Defines operations on vectors of Functions
39  \ingroup mra
40 
41  This file defines a number of operations on vectors of functions.
42  Assume v is a vector of NDIM-D functions of a certain type.
43 
44 
45  Operations on array of functions
46 
47  *) copying: deep copying of vectors of functions to vector of functions
48  \code
49  vector2 = copy(world, vector1,fence);
50  \endcode
51 
52  *) compress: convert multiwavelet representation to legendre representation
53  \code
54  compress(world, vector, fence);
55  \endcode
56 
57  *) reconstruct: convert representation to multiwavelets
58  \code
59  reconstruct(world, vector, fence);
60  \endcode
61 
62  *) make_nonstandard: convert to non-standard form
63  \code
64  make_nonstandard(world, v, fence);
65  \endcode
66 
67  *) standard: convert to standard form
68  \code
69  standard(world, v, fence);
70  \endcode
71 
72  *) truncate: truncating vectors of functions to desired precision
73  \code
74  truncate(world, v, tolerance, fence);
75  \endcode
76 
77 
78  *) zero function: create a vector of zero functions of length n
79  \code
80  v=zero(world, n);
81  \endcode
82 
83  *) transform: transform a representation from one basis to another
84  \code
85  transform(world, vector, tensor, tolerance, fence )
86  \endcode
87 
88  Setting thresh-hold for precision
89 
90  *) set_thresh: setting a finite thresh-hold for a vector of functions
91  \code
92  void set_thresh(World& world, std::vector< Function<T,NDIM> >& v, double thresh, bool fence=true);
93  \endcode
94 
95  Arithmetic Operations on arrays of functions
96 
97  *) conjugation: conjugate a vector of complex functions
98 
99  *) add
100  *) sub
101  *) mul
102  - mul_sparse
103  *) square
104  *) gaxpy
105  *) apply
106 
107  Norms, inner-products, blas-1 like operations on vectors of functions
108 
109  *) inner
110  *) matrix_inner
111  *) norm_tree
112  *) normalize
113  *) norm2
114  - norm2s
115  *) scale(world, v, alpha);
116 
117 
118 
119 
120 */
121 
122 #include <madness/mra/mra.h>
123 #include <madness/mra/derivative.h>
125 #include <cstdio>
126 
127 namespace madness {
128 
129 
130 
131  /// Compress a vector of functions
132  template <typename T, std::size_t NDIM>
133  void compress(World& world,
134  const std::vector< Function<T,NDIM> >& v,
135  bool fence=true) {
136 
137  PROFILE_BLOCK(Vcompress);
139 // bool must_fence = false;
140 // for (unsigned int i=0; i<v.size(); ++i) {
141 // if (!v[i].is_compressed()) {
142 // v[i].compress(false);
143 // must_fence = true;
144 // }
145 // }
146 //
147 // if (fence && must_fence) world.gop.fence();
148  }
149 
150 
151  /// reconstruct a vector of functions
152 
153  /// implies fence
154  /// return v for chaining
155  template <typename T, std::size_t NDIM>
156  const std::vector< Function<T,NDIM> >& reconstruct(const std::vector< Function<T,NDIM> >& v) {
158  }
159 
160  /// compress a vector of functions
161 
162  /// implies fence
163  /// return v for chaining
164  template <typename T, std::size_t NDIM>
165  const std::vector< Function<T,NDIM> >& compress(const std::vector< Function<T,NDIM> >& v) {
167  }
168 
169  /// Reconstruct a vector of functions
170  template <typename T, std::size_t NDIM>
171  void reconstruct(World& world,
172  const std::vector< Function<T,NDIM> >& v,
173  bool fence=true) {
174  PROFILE_BLOCK(Vreconstruct);
175 // bool must_fence = false;
177 // for (unsigned int i=0; i<v.size(); ++i) {
178 // if (v[i].is_compressed() or v[i].is_nonstandard()) {
179 // v[i].reconstruct(false);
180 // must_fence = true;
181 // }
182 // }
183 //
184 // if (fence && must_fence) world.gop.fence();
185  }
186 
187  /// change tree_state of a vector of functions to redundant
188  template <typename T, std::size_t NDIM>
189  void make_redundant(World& world,
190  const std::vector< Function<T,NDIM> >& v,
191  bool fence=true) {
192 
193  PROFILE_BLOCK(Vcompress);
195 // bool must_fence = false;
196 // for (unsigned int i=0; i<v.size(); ++i) {
197 // if (!v[i].get_impl()->is_redundant()) {
198 // v[i].get_impl()->make_redundant(false);
199 // must_fence = true;
200 // }
201 // }
202 //
203 // if (fence && must_fence) world.gop.fence();
204  }
205 
206  /// refine the functions according to the autorefine criteria
207  template <typename T, std::size_t NDIM>
208  void refine(World& world, const std::vector<Function<T,NDIM> >& vf,
209  bool fence=true) {
210  for (const auto& f : vf) f.refine(false);
211  if (fence) world.gop.fence();
212  }
213 
214  /// refine all functions to a common (finest) level
215 
216  /// if functions are not initialized (impl==NULL) they are ignored
217  template <typename T, std::size_t NDIM>
218  void refine_to_common_level(World& world, std::vector<Function<T,NDIM> >& vf,
219  bool fence=true) {
220 
221  reconstruct(world,vf);
223  std::vector<FunctionImpl<T,NDIM>*> v_ptr;
224 
225  // push initialized function pointers into the vector v_ptr
226  for (unsigned int i=0; i<vf.size(); ++i) {
227  if (vf[i].is_initialized()) v_ptr.push_back(vf[i].get_impl().get());
228  }
229 
230  // sort and remove duplicates to not confuse the refining function
231  std::sort(v_ptr.begin(),v_ptr.end());
232  typename std::vector<FunctionImpl<T, NDIM>*>::iterator it;
233  it = std::unique(v_ptr.begin(), v_ptr.end());
234  v_ptr.resize( std::distance(v_ptr.begin(),it) );
235 
236  std::vector< Tensor<T> > c(v_ptr.size());
237  v_ptr[0]->refine_to_common_level(v_ptr, c, key0);
238  if (fence) v_ptr[0]->world.gop.fence();
239  if (VERIFY_TREE)
240  for (unsigned int i=0; i<vf.size(); i++) vf[i].verify_tree();
241  }
242 
243  /// Generates non-standard form of a vector of functions
244  template <typename T, std::size_t NDIM>
245  void make_nonstandard(World& world,
246  std::vector< Function<T,NDIM> >& v,
247  bool fence= true) {
248  PROFILE_BLOCK(Vnonstandard);
250 // reconstruct(world, v);
251 // for (unsigned int i=0; i<v.size(); ++i) {
252 // v[i].make_nonstandard(false, false);
253 // }
254 // if (fence) world.gop.fence();
255  }
256 
257 
258  /// Generates standard form of a vector of functions
259  template <typename T, std::size_t NDIM>
260  void standard(World& world,
261  std::vector< Function<T,NDIM> >& v,
262  bool fence=true) {
263  PROFILE_BLOCK(Vstandard);
265 // for (unsigned int i=0; i<v.size(); ++i) {
266 // v[i].standard(false);
267 // }
268 // if (fence) world.gop.fence();
269  }
270 
271 
272  /// change tree state of the functions
273 
274  /// will respect fence
275  /// @return v for chaining
276  template <typename T, std::size_t NDIM>
277  const std::vector<Function<T,NDIM>>& change_tree_state(const std::vector<Function<T,NDIM>>& v,
278  const TreeState finalstate,
279  const bool fence=true) {
280  if (v.size()==0) return v;
281  // find initialized function with world
282  Function<T,NDIM> dummy;
283  for (const auto& f : v)
284  if (f.is_initialized()) {
285  dummy=f;
286  break;
287  }
288  if (not dummy.is_initialized()) return v;
289  World& world=dummy.world();
290 
291 
292  // if a tree state cannot directly be changed to finalstate, we need to go via intermediate
293  auto change_initial_to_intermediate =[](const std::vector<Function<T,NDIM>>& v,
294  const TreeState initialstate,
295  const TreeState intermediatestate) {
296  int must_fence=0;
297  for (auto& f : v) {
298  if (f.is_initialized() and f.get_impl()->get_tree_state()==initialstate) {
299  f.change_tree_state(intermediatestate,false);
300  must_fence=1;
301  }
302  }
303  return must_fence;
304  };
305 
306  int do_fence=0;
307  if (finalstate==compressed) {
308  do_fence+=change_initial_to_intermediate(v,redundant,TreeState::reconstructed);
309  }
310  if (finalstate==nonstandard) {
311  do_fence+=change_initial_to_intermediate(v,compressed,TreeState::reconstructed);
312  do_fence+=change_initial_to_intermediate(v,redundant,TreeState::reconstructed);
313  }
314  if (finalstate==nonstandard_with_leaves) {
315  do_fence+=change_initial_to_intermediate(v,compressed,TreeState::reconstructed);
316  do_fence+=change_initial_to_intermediate(v,nonstandard,TreeState::reconstructed);
317  do_fence+=change_initial_to_intermediate(v,redundant,TreeState::reconstructed);
318  }
319  if (finalstate==redundant) {
320  do_fence+=change_initial_to_intermediate(v,compressed,TreeState::reconstructed);
321  do_fence+=change_initial_to_intermediate(v,nonstandard,TreeState::reconstructed);
322  do_fence+=change_initial_to_intermediate(v,nonstandard_with_leaves,TreeState::reconstructed);
323  }
324  if (do_fence>0) world.gop.fence();
325 
326  for (unsigned int i=0; i<v.size(); ++i) v[i].change_tree_state(finalstate,fence);
327  if (fence) world.gop.fence();
328 
329  return v;
330  }
331 
332 
333  /// Truncates a vector of functions
334  template <typename T, std::size_t NDIM>
335  void truncate(World& world,
336  std::vector< Function<T,NDIM> >& v,
337  double tol=0.0,
338  bool fence=true) {
339  PROFILE_BLOCK(Vtruncate);
340 
341  // truncate in compressed form only for low-dimensional functions
342  // compression is very expensive if low-rank tensor approximations are used
343  if (NDIM<4) compress(world, v);
344 
345  for (unsigned int i=0; i<v.size(); ++i) {
346  v[i].truncate(tol, false);
347  }
348 
349  if (fence) world.gop.fence();
350  }
351 
352  /// Truncates a vector of functions
353 
354  /// @return the truncated vector for chaining
355  template <typename T, std::size_t NDIM>
356  std::vector< Function<T,NDIM> > truncate(std::vector< Function<T,NDIM> > v,
357  double tol=0.0, bool fence=true) {
358  if (v.size()>0) truncate(v[0].world(),v,tol,fence);
359  return v;
360  }
361 
362  /// reduces the tensor rank of the coefficient tensor (if applicable)
363 
364  /// @return the vector for chaining
365  template <typename T, std::size_t NDIM>
366  std::vector< Function<T,NDIM> > reduce_rank(std::vector< Function<T,NDIM> > v,
367  double thresh=0.0, bool fence=true) {
368  if (v.size()==0) return v;
369  for (auto& vv : v) vv.reduce_rank(thresh,false);
370  if (fence) v[0].world().gop.fence();
371  return v;
372  }
373 
374 
375  /// Applies a derivative operator to a vector of functions
376  template <typename T, std::size_t NDIM>
377  std::vector< Function<T,NDIM> >
378  apply(World& world,
379  const Derivative<T,NDIM>& D,
380  const std::vector< Function<T,NDIM> >& v,
381  bool fence=true)
382  {
383  reconstruct(world, v);
384  std::vector< Function<T,NDIM> > df(v.size());
385  for (unsigned int i=0; i<v.size(); ++i) {
386  df[i] = D(v[i],false);
387  }
388  if (fence) world.gop.fence();
389  return df;
390  }
391 
392  /// Generates a vector of zero functions (reconstructed)
393  template <typename T, std::size_t NDIM>
394  std::vector< Function<T,NDIM> >
395  zero_functions(World& world, int n, bool fence=true) {
396  PROFILE_BLOCK(Vzero_functions);
397  std::vector< Function<T,NDIM> > r(n);
398  for (int i=0; i<n; ++i)
399  r[i] = Function<T,NDIM>(FunctionFactory<T,NDIM>(world).fence(false));
400 
401  if (n && fence) world.gop.fence();
402 
403  return r;
404  }
405 
406  /// Generates a vector of zero functions (compressed)
407  template <typename T, std::size_t NDIM>
408  std::vector< Function<T,NDIM> >
409  zero_functions_compressed(World& world, int n, bool fence=true) {
410  PROFILE_BLOCK(Vzero_functions);
411  std::vector< Function<T,NDIM> > r(n);
412  for (int i=0; i<n; ++i)
413  r[i] = Function<T,NDIM>(FunctionFactory<T,NDIM>(world).fence(false).compressed(true).initial_level(1));
414  if (n && fence) world.gop.fence();
415  return r;
416  }
417 
418 
419  /// orthonormalize the vectors
420  template<typename T, std::size_t NDIM>
421  std::vector<Function<T,NDIM>> orthonormalize(const std::vector<Function<T,NDIM> >& vf_in) {
422  if (vf_in.size()==0) return std::vector<Function<T,NDIM>>();
423  World& world=vf_in.front().world();
424  auto vf=copy(world,vf_in);
425  normalize(world,vf);
426  if (vf.size()==1) return copy(world,vf_in);
427  double maxq;
428  double trantol=0.0;
429  auto Q2=[](const Tensor<T>& s) {
430  Tensor<T> Q = -0.5*s;
431  for (int i=0; i<s.dim(0); ++i) Q(i,i) += 1.5;
432  return Q;
433  };
434 
435  do {
436  Tensor<T> Q = Q2(matrix_inner(world, vf, vf));
437  maxq=0.0;
438  for (int i=0; i<Q.dim(0); ++i)
439  for (int j=0; j<i; ++j)
440  maxq = std::max(maxq,std::abs(Q(i,j)));
441 
442  vf = transform(world, vf, Q, trantol, true);
443  truncate(world, vf);
444 
445  } while (maxq>0.01);
446  normalize(world,vf);
447  return vf;
448  }
449 
450 
451  /// symmetric orthonormalization (see e.g. Szabo/Ostlund)
452 
453  /// @param[in] the vector to orthonormalize
454  /// @param[in] overlap matrix
455  template <typename T, std::size_t NDIM>
456  std::vector<Function<T,NDIM> > orthonormalize_symmetric(
457  const std::vector<Function<T,NDIM> >& v,
458  const Tensor<T>& ovlp) {
459  if(v.empty()) return v;
460 
461  Tensor<T> U;
463  syev(ovlp,U,s);
464 
465  // transform s to s^{-1}
466  for(size_t i=0;i<v.size();++i) s(i)=1.0/(sqrt(s(i)));
467 
468  // save Ut before U gets modified with s^{-1}
469  const Tensor<T> Ut=transpose(U);
470  for(size_t i=0;i<v.size();++i){
471  for(size_t j=0;j<v.size();++j){
472  U(i,j)=U(i,j)*s(j);
473  }
474  }
475 
476  Tensor<T> X=inner(U,Ut,1,0);
477 
478  World& world=v.front().world();
479  return transform(world,v,X);
480 
481  }
482  /// convenience routine for symmetric orthonormalization (see e.g. Szabo/Ostlund)
483  /// overlap matrix is calculated
484  /// @param[in] the vector to orthonormalize
485  template <typename T, std::size_t NDIM>
486  std::vector<Function<T,NDIM> > orthonormalize_symmetric(const std::vector<Function<T,NDIM> >& v){
487  if(v.empty()) return v;
488 
489 
490  World& world=v.front().world();
491  Tensor<T> ovlp = matrix_inner(world, v, v);
492 
493  return orthonormalize_symmetric(v,ovlp);
494  }
495 
496  /// canonical orthonormalization (see e.g. Szabo/Ostlund)
497  /// @param[in] the vector to orthonormalize
498  /// @param[in] overlap matrix
499  /// @param[in] lindep linear dependency threshold relative to largest eigenvalue
500  template <typename T, std::size_t NDIM>
501  std::vector<Function<T,NDIM> > orthonormalize_canonical(
502  const std::vector<Function<T,NDIM> >& v,
503  const Tensor<T>& ovlp,
504  double lindep) {
505 
506  if(v.empty()) return v;
507 
508  Tensor<T> U;
510  syev(ovlp,U,s);
511  lindep*=s(s.size()-1); // eigenvalues are in ascending order
512 
513  // transform s to s^{-1}
514  int rank=0,lo=0;
516  for(size_t i=0;i<v.size();++i) {
517  if (s(i)>lindep) {
518  sqrts(i)=1.0/(sqrt(s(i)));
519  rank++;
520  } else {
521  sqrts(i)=0.0;
522  lo++;
523  }
524  }
525  MADNESS_ASSERT(size_t(lo+rank)==v.size());
526 
527  for(size_t i=0;i<v.size();++i){
528  for(size_t j=0;j<v.size();++j){
529  U(i,j)=U(i,j)*(sqrts(j));
530  }
531  }
532  Tensor<T> X=U(_,Slice(lo,-1));
533 
534  World& world=v.front().world();
535  return transform(world,v,X);
536 
537  }
538 
539  /// convenience routine for canonical routine for symmetric orthonormalization (see e.g. Szabo/Ostlund)
540  /// overlap matrix is calculated
541  /// @param[in] the vector to orthonormalize
542  template <typename T, std::size_t NDIM>
543  std::vector<Function<T,NDIM> > orthonormalize_canonical(const std::vector<Function<T,NDIM> >& v,
544  const double lindep){
545  if(v.empty()) return v;
546 
547  World& world=v.front().world();
548  Tensor<T> ovlp = matrix_inner(world, v, v);
549 
550  return orthonormalize_canonical(v,ovlp,lindep);
551  }
552 
553  /// cholesky orthonormalization without pivoting
554  /// @param[in] the vector to orthonormalize
555  /// @param[in] overlap matrix, destroyed on return!
556  template <typename T, std::size_t NDIM>
557  std::vector<Function<T,NDIM> > orthonormalize_cd(
558  const std::vector<Function<T,NDIM> >& v,
559  Tensor<T>& ovlp) {
560 
561  if (v.empty()) return v;
562 
563  cholesky(ovlp); // destroys ovlp and gives back Upper ∆ Matrix from CD
564 
565  Tensor<T> L = transpose(ovlp);
566  Tensor<T> Linv = inverse(L);
567  Tensor<T> U = transpose(Linv);
568 
569  World& world=v.front().world();
570  return transform(world, v, U);
571 
572  }
573 
574  /// convenience routine for cholesky orthonormalization without pivoting
575  /// @param[in] the vector to orthonormalize
576  /// @param[in] overlap matrix
577  template <typename T, std::size_t NDIM>
578  std::vector<Function<T,NDIM> > orthonormalize_cd(const std::vector<Function<T,NDIM> >& v){
579  if(v.empty()) return v;
580 
581  World& world=v.front().world();
582  Tensor<T> ovlp = matrix_inner(world, v, v);
583 
584  return orthonormalize_cd(v,ovlp);
585  }
586 
587  /// @param[in] the vector to orthonormalize
588  /// @param[in] overlap matrix, will be destroyed on return!
589  /// @param[in] tolerance for numerical rank reduction
590  /// @param[out] pivoting vector, no allocation on input needed
591  /// @param[out] rank
592  /// @return orthonrormalized vector (may or may not be truncated)
593  template <typename T, std::size_t NDIM>
594  std::vector<Function<T,NDIM> > orthonormalize_rrcd(
595  const std::vector<Function<T,NDIM> >& v,
596  Tensor<T>& ovlp,
597  const double tol,
598  Tensor<integer>& piv,
599  int& rank) {
600 
601  if (v.empty()) {
602  return v;
603  }
604 
605  rr_cholesky(ovlp,tol,piv,rank); // destroys ovlp and gives back Upper ∆ Matrix from CCD
606 
607  // rearrange and truncate the functions according to the pivoting of the rr_cholesky
608  std::vector<Function<T,NDIM> > pv(rank);
609  for(integer i=0;i<rank;++i){
610  pv[i]=v[piv[i]];
611  }
612  ovlp=ovlp(Slice(0,rank-1),Slice(0,rank-1));
613 
614  Tensor<T> L = transpose(ovlp);
615  Tensor<T> Linv = inverse(L);
616  Tensor<T> U = transpose(Linv);
617 
618  World& world=v.front().world();
619  return transform(world, pv, U);
620  }
621 
622  /// convenience routine for orthonromalize_cholesky: orthonromalize_cholesky without information on pivoting and rank
623  /// @param[in] the vector to orthonormalize
624  /// @param[in] overlap matrix
625  /// @param[in] tolerance for numerical rank reduction
626  template <typename T, std::size_t NDIM>
627  std::vector<Function<T,NDIM> > orthonormalize_rrcd(const std::vector<Function<T,NDIM> >& v, Tensor<T> ovlp , const double tol) {
628  Tensor<integer> piv;
629  int rank;
630  return orthonormalize_rrcd(v,ovlp,tol,piv,rank);
631  }
632 
633  /// convenience routine for orthonromalize_cholesky: computes the overlap matrix and then calls orthonromalize_cholesky
634  /// @param[in] the vector to orthonormalize
635  /// @param[in] tolerance for numerical rank reduction
636  template <typename T, std::size_t NDIM>
637  std::vector<Function<T,NDIM> > orthonormalize_rrcd(const std::vector<Function<T,NDIM> >& v, const double tol) {
638  if (v.empty()) {
639  return v;
640  }
641  // compute overlap
642  World& world=v.front().world();
643  Tensor<T> ovlp = matrix_inner(world, v, v);
644  return orthonormalize_rrcd(v,ovlp,tol);
645  }
646 
647  /// combine two vectors
648  template <typename T, std::size_t NDIM>
649  std::vector<Function<T,NDIM> > append(const std::vector<Function<T,NDIM> > & lhs, const std::vector<Function<T,NDIM> > & rhs){
650  std::vector<Function<T,NDIM> > v=lhs;
651  for (std::size_t i = 0; i < rhs.size(); ++i) v.push_back(rhs[i]);
652  return v;
653  }
654 
655  template <typename T, std::size_t NDIM>
656  std::vector<Function<T,NDIM> > flatten(const std::vector< std::vector<Function<T,NDIM> > >& vv){
657  std::vector<Function<T,NDIM> >result;
658  for(const auto& x:vv) result=append(result,x);
659  return result;
660  }
661 
662  template<typename T, std::size_t NDIM>
663  std::vector<std::shared_ptr<FunctionImpl<T,NDIM>>> get_impl(const std::vector<Function<T,NDIM>>& v) {
664  std::vector<std::shared_ptr<FunctionImpl<T,NDIM>>> result;
665  for (auto& f : v) result.push_back(f.get_impl());
666  return result;
667  }
668 
669  template<typename T, std::size_t NDIM>
670  void set_impl(std::vector<Function<T,NDIM>>& v, const std::vector<std::shared_ptr<FunctionImpl<T,NDIM>>> vimpl) {
671  MADNESS_CHECK(vimpl.size()==v.size());
672  for (std::size_t i=0; i<vimpl.size(); ++i) v[i].set_impl(vimpl[i]);
673  }
674 
675  template<typename T, std::size_t NDIM>
676  std::vector<Function<T,NDIM>> impl2function(const std::vector<std::shared_ptr<FunctionImpl<T,NDIM>>> vimpl) {
677  std::vector<Function<T,NDIM>> v(vimpl.size());
678  for (std::size_t i=0; i<vimpl.size(); ++i) v[i].set_impl(vimpl[i]);
679  return v;
680  }
681 
682 
683  /// Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i]
684 
685  /// Uses sparsity in the transformation matrix --- set small elements to
686  /// zero to take advantage of this.
687  template <typename T, typename R, std::size_t NDIM>
688  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> >
689  transform(World& world,
690  const std::vector< Function<T,NDIM> >& v,
691  const Tensor<R>& c,
692  bool fence=true) {
693 
694  PROFILE_BLOCK(Vtransformsp);
695  typedef TENSOR_RESULT_TYPE(T,R) resultT;
696  int n = v.size(); // n is the old dimension
697  int m = c.dim(1); // m is the new dimension
698  MADNESS_ASSERT(n==c.dim(0));
699 
700  std::vector< Function<resultT,NDIM> > vc = zero_functions_compressed<resultT,NDIM>(world, m);
701  compress(world, v);
702 
703  for (int i=0; i<m; ++i) {
704  for (int j=0; j<n; ++j) {
705  if (c(j,i) != R(0.0)) vc[i].gaxpy(resultT(1.0),v[j],resultT(c(j,i)),false);
706  }
707  }
708 
709  if (fence) world.gop.fence();
710  return vc;
711  }
712 
713  /// Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i]
714 
715  /// all trees are in reconstructed state, final trees have to be summed down if no fence is present
716  template <typename T, typename R, std::size_t NDIM>
717  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> >
719  const std::vector< Function<T,NDIM> >& v,
720  const Tensor<R>& c,
721  bool fence=true) {
722 
723  PROFILE_BLOCK(Vtransformsp);
724  typedef TENSOR_RESULT_TYPE(T,R) resultT;
725  int n = v.size(); // n is the old dimension
726  int m = c.dim(1); // m is the new dimension
727  MADNESS_CHECK(n==c.dim(0));
728 
729  // if we fence set the right tree state here, otherwise it has to be correct from the start.
730  if (fence) change_tree_state(v,reconstructed);
731  for (const auto& vv : v) MADNESS_CHECK_THROW(
732  vv.get_impl()->get_tree_state()==reconstructed,"trees have to be reconstructed in transform_reconstructed");
733 
734  std::vector< Function<resultT,NDIM> > result = zero_functions<resultT,NDIM>(world, m);
735 
736  for (int i=0; i<m; ++i) {
737  result[i].get_impl()->set_tree_state(redundant_after_merge);
738  for (int j=0; j<n; ++j) {
739  if (c(j,i) != R(0.0)) v[j].get_impl()->accumulate_trees(*(result[i].get_impl()),resultT(c(j,i)),true);
740  }
741  }
742 
743  // if we fence we can as well finish the job here. Otherwise no harm done, as the tree state is well-defined.
744  if (fence) {
745  world.gop.fence();
746  // for (auto& r : vc) r.sum_down(false);
747  for (auto& r : result) r.get_impl()->finalize_sum();
748  world.gop.fence();
749  }
750  return result;
751  }
752 
753  /// this version of transform uses Function::vtransform and screens
754  /// using both elements of `c` and `v`
755  template <typename L, typename R, std::size_t NDIM>
756  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> >
757  transform(World& world, const std::vector< Function<L,NDIM> >& v,
758  const Tensor<R>& c, double tol, bool fence) {
759  PROFILE_BLOCK(Vtransform);
760  MADNESS_ASSERT(v.size() == (unsigned int)(c.dim(0)));
761 
762  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> > vresult
763  = zero_functions_compressed<TENSOR_RESULT_TYPE(L,R),NDIM>(world, c.dim(1));
764 
765  compress(world, v, true);
766  vresult[0].vtransform(v, c, vresult, tol, fence);
767  return vresult;
768  }
769 
770  template <typename T, typename R, std::size_t NDIM>
771  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> >
772  transform(World& world,
773  const std::vector< Function<T,NDIM> >& v,
774  const DistributedMatrix<R>& c,
775  bool fence=true) {
776  PROFILE_FUNC;
777 
778  typedef TENSOR_RESULT_TYPE(T,R) resultT;
779  long n = v.size(); // n is the old dimension
780  long m = c.rowdim(); // m is the new dimension
781  MADNESS_ASSERT(n==c.coldim());
782 
783  // new(i) = sum(j) old(j) c(j,i)
784 
785  Tensor<T> tmp(n,m);
786  c.copy_to_replicated(tmp); // for debugging
787  tmp = transpose(tmp);
788 
789  std::vector< Function<resultT,NDIM> > vc = zero_functions_compressed<resultT,NDIM>(world, m);
790  compress(world, v);
791 
792  for (int i=0; i<m; ++i) {
793  for (int j=0; j<n; ++j) {
794  if (tmp(j,i) != R(0.0)) vc[i].gaxpy(1.0,v[j],tmp(j,i),false);
795  }
796  }
797 
798  if (fence) world.gop.fence();
799  return vc;
800  }
801 
802 
803  /// Scales inplace a vector of functions by distinct values
804  template <typename T, typename Q, std::size_t NDIM>
805  void scale(World& world,
806  std::vector< Function<T,NDIM> >& v,
807  const std::vector<Q>& factors,
808  bool fence=true) {
809  PROFILE_BLOCK(Vscale);
810  for (unsigned int i=0; i<v.size(); ++i) v[i].scale(factors[i],false);
811  if (fence) world.gop.fence();
812  }
813 
814  /// Scales inplace a vector of functions by the same
815  template <typename T, typename Q, std::size_t NDIM>
816  void scale(World& world,
817  std::vector< Function<T,NDIM> >& v,
818  const Q factor,
819  bool fence=true) {
820  PROFILE_BLOCK(Vscale);
821  for (unsigned int i=0; i<v.size(); ++i) v[i].scale(factor,false);
822  if (fence) world.gop.fence();
823  }
824 
825  /// Computes the 2-norms of a vector of functions
826  template <typename T, std::size_t NDIM>
827  std::vector<double> norm2s(World& world,
828  const std::vector< Function<T,NDIM> >& v) {
829  PROFILE_BLOCK(Vnorm2);
830  std::vector<double> norms(v.size());
831  for (unsigned int i=0; i<v.size(); ++i) norms[i] = v[i].norm2sq_local();
832  world.gop.sum(&norms[0], norms.size());
833  for (unsigned int i=0; i<v.size(); ++i) norms[i] = sqrt(norms[i]);
834  world.gop.fence();
835  return norms;
836  }
837  /// Computes the 2-norms of a vector of functions
838  template <typename T, std::size_t NDIM>
839  Tensor<double> norm2s_T(World& world, const std::vector<Function<T, NDIM>>& v) {
840  PROFILE_BLOCK(Vnorm2);
841  Tensor<double> norms(v.size());
842  for (unsigned int i = 0; i < v.size(); ++i) norms[i] = v[i].norm2sq_local();
843  world.gop.sum(&norms[0], norms.size());
844  for (unsigned int i = 0; i < v.size(); ++i) norms[i] = sqrt(norms[i]);
845  world.gop.fence();
846  return norms;
847  }
848 
849  /// Computes the 2-norm of a vector of functions
850  template <typename T, std::size_t NDIM>
851  double norm2(World& world,const std::vector< Function<T,NDIM> >& v) {
852  PROFILE_BLOCK(Vnorm2);
853  if (v.size()==0) return 0.0;
854  std::vector<double> norms(v.size());
855  for (unsigned int i=0; i<v.size(); ++i) norms[i] = v[i].norm2sq_local();
856  world.gop.sum(&norms[0], norms.size());
857  for (unsigned int i=1; i<v.size(); ++i) norms[0] += norms[i];
858  world.gop.fence();
859  return sqrt(norms[0]);
860  }
861 
862  inline double conj(double x) {
863  return x;
864  }
865 
866  inline double conj(float x) {
867  return x;
868  }
869 
870 // !!! FIXME: this task is broken because FunctionImpl::inner_local forces a
871 // future on return from WorldTaskQueue::reduce, which will causes a deadlock if
872 // run inside a task. This behavior must be changed before this task can be used
873 // again.
874 //
875 // template <typename T, typename R, std::size_t NDIM>
876 // struct MatrixInnerTask : public TaskInterface {
877 // Tensor<TENSOR_RESULT_TYPE(T,R)> result; // Must be a copy
878 // const Function<T,NDIM>& f;
879 // const std::vector< Function<R,NDIM> >& g;
880 // long jtop;
881 //
882 // MatrixInnerTask(const Tensor<TENSOR_RESULT_TYPE(T,R)>& result,
883 // const Function<T,NDIM>& f,
884 // const std::vector< Function<R,NDIM> >& g,
885 // long jtop)
886 // : result(result), f(f), g(g), jtop(jtop) {}
887 //
888 // void run(World& world) {
889 // for (long j=0; j<jtop; ++j) {
890 // result(j) = f.inner_local(g[j]);
891 // }
892 // }
893 //
894 // private:
895 // /// Get the task id
896 //
897 // /// \param id The id to set for this task
898 // virtual void get_id(std::pair<void*,unsigned short>& id) const {
899 // PoolTaskInterface::make_id(id, *this);
900 // }
901 // }; // struct MatrixInnerTask
902 
903 
904 
905  template <typename T, std::size_t NDIM>
907  const std::vector< Function<T,NDIM> >& f,
908  const std::vector< Function<T,NDIM> >& g,
909  bool sym=false)
910  {
911  PROFILE_FUNC;
913  const int64_t n = A.coldim();
914  const int64_t m = A.rowdim();
915  MADNESS_ASSERT(int64_t(f.size()) == n && int64_t(g.size()) == m);
916 
917  // Assume we can always create an ichunk*jchunk matrix locally
918  const int ichunk = 1000;
919  const int jchunk = 1000; // 1000*1000*8 = 8 MBytes
920  for (int64_t ilo=0; ilo<n; ilo+=ichunk) {
921  int64_t ihi = std::min(ilo + ichunk, n);
922  std::vector< Function<T,NDIM> > ivec(f.begin()+ilo, f.begin()+ihi);
923  for (int64_t jlo=0; jlo<m; jlo+=jchunk) {
924  int64_t jhi = std::min(jlo + jchunk, m);
925  std::vector< Function<T,NDIM> > jvec(g.begin()+jlo, g.begin()+jhi);
926 
927  Tensor<T> P = matrix_inner(A.get_world(),ivec,jvec);
928  A.copy_from_replicated_patch(ilo, ihi-1, jlo, jhi-1, P);
929  }
930  }
931  return A;
932  }
933 
934  /// Computes the matrix inner product of two function vectors - q(i,j) = inner(f[i],g[j])
935 
936  /// For complex types symmetric is interpreted as Hermitian.
937  ///
938  /// The current parallel loop is non-optimal but functional.
939  template <typename T, typename R, std::size_t NDIM>
941  const std::vector< Function<T,NDIM> >& f,
942  const std::vector< Function<R,NDIM> >& g,
943  bool sym=false)
944  {
945  world.gop.fence();
946  compress(world, f);
947 // if ((void*)(&f) != (void*)(&g)) compress(world, g);
948  compress(world, g);
949 
950 
951  std::vector<const FunctionImpl<T,NDIM>*> left(f.size());
952  std::vector<const FunctionImpl<R,NDIM>*> right(g.size());
953  for (unsigned int i=0; i<f.size(); i++) left[i] = f[i].get_impl().get();
954  for (unsigned int i=0; i<g.size(); i++) right[i]= g[i].get_impl().get();
955 
957 
958  world.gop.fence();
959  world.gop.sum(r.ptr(),f.size()*g.size());
960 
961  return r;
962  }
963 
964  /// Computes the matrix inner product of two function vectors - q(i,j) = inner(f[i],g[j])
965 
966  /// For complex types symmetric is interpreted as Hermitian.
967  ///
968  /// The current parallel loop is non-optimal but functional.
969  template <typename T, typename R, std::size_t NDIM>
971  const std::vector< Function<T,NDIM> >& f,
972  const std::vector< Function<R,NDIM> >& g,
973  bool sym=false) {
974  PROFILE_BLOCK(Vmatrix_inner);
975  long n=f.size(), m=g.size();
976  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n,m);
977  if (sym) MADNESS_ASSERT(n==m);
978 
979  world.gop.fence();
980  compress(world, f);
981  if ((void*)(&f) != (void*)(&g)) compress(world, g);
982 
983  for (long i=0; i<n; ++i) {
984  long jtop = m;
985  if (sym) jtop = i+1;
986  for (long j=0; j<jtop; ++j) {
987  r(i,j) = f[i].inner_local(g[j]);
988  if (sym) r(j,i) = conj(r(i,j));
989  }
990  }
991 
992 // for (long i=n-1; i>=0; --i) {
993 // long jtop = m;
994 // if (sym) jtop = i+1;
995 // world.taskq.add(new MatrixInnerTask<T,R,NDIM>(r(i,_), f[i], g, jtop));
996 // }
997  world.gop.fence();
998  world.gop.sum(r.ptr(),n*m);
999 
1000 // if (sym) {
1001 // for (int i=0; i<n; ++i) {
1002 // for (int j=0; j<i; ++j) {
1003 // r(j,i) = conj(r(i,j));
1004 // }
1005 // }
1006 // }
1007  return r;
1008  }
1009 
1010  /// Computes the element-wise inner product of two function vectors - q(i) = inner(f[i],g[i])
1011  template <typename T, typename R, std::size_t NDIM>
1013  const std::vector< Function<T,NDIM> >& f,
1014  const std::vector< Function<R,NDIM> >& g) {
1015  PROFILE_BLOCK(Vinnervv);
1016  long n=f.size(), m=g.size();
1017  MADNESS_CHECK(n==m);
1018  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
1019 
1020  compress(world, f);
1021  compress(world, g);
1022 
1023  for (long i=0; i<n; ++i) {
1024  r(i) = f[i].inner_local(g[i]);
1025  }
1026 
1027  world.taskq.fence();
1028  world.gop.sum(r.ptr(),n);
1029  world.gop.fence();
1030  return r;
1031  }
1032 
1033 
1034  /// Computes the inner product of a function with a function vector - q(i) = inner(f,g[i])
1035  template <typename T, typename R, std::size_t NDIM>
1037  const Function<T,NDIM>& f,
1038  const std::vector< Function<R,NDIM> >& g) {
1039  PROFILE_BLOCK(Vinner);
1040  long n=g.size();
1041  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
1042 
1043  f.compress();
1044  compress(world, g);
1045 
1046  for (long i=0; i<n; ++i) {
1047  r(i) = f.inner_local(g[i]);
1048  }
1049 
1050  world.taskq.fence();
1051  world.gop.sum(r.ptr(),n);
1052  world.gop.fence();
1053  return r;
1054  }
1055 
1056  /// inner function with right signature for the nonlinear solver
1057  /// this is needed for the KAIN solvers and other functions
1058  template <typename T, typename R, std::size_t NDIM>
1059  TENSOR_RESULT_TYPE(T,R) inner( const std::vector< Function<T,NDIM> >& f,
1060  const std::vector< Function<R,NDIM> >& g){
1061  MADNESS_ASSERT(f.size()==g.size());
1062  if(f.empty()) return 0.0;
1063  else return inner(f[0].world(),f,g).sum();
1064  }
1065 
1066 
1067  /// Multiplies a function against a vector of functions --- q[i] = a * v[i]
1068  template <typename T, typename R, std::size_t NDIM>
1069  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1070  mul(World& world,
1071  const Function<T,NDIM>& a,
1072  const std::vector< Function<R,NDIM> >& v,
1073  bool fence=true) {
1074  PROFILE_BLOCK(Vmul);
1075  a.reconstruct(false);
1076  reconstruct(world, v, false);
1077  world.gop.fence();
1078  return vmulXX(a, v, 0.0, fence);
1079  }
1080 
1081  /// Multiplies a function against a vector of functions using sparsity of a and v[i] --- q[i] = a * v[i]
1082  template <typename T, typename R, std::size_t NDIM>
1083  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1085  const Function<T,NDIM>& a,
1086  const std::vector< Function<R,NDIM> >& v,
1087  double tol,
1088  bool fence=true) {
1089  PROFILE_BLOCK(Vmulsp);
1090  a.reconstruct(false);
1091  reconstruct(world, v, false);
1092  world.gop.fence();
1093  for (unsigned int i=0; i<v.size(); ++i) {
1094  v[i].norm_tree(false);
1095  }
1096  a.norm_tree();
1097  return vmulXX(a, v, tol, fence);
1098  }
1099 
1100 
1101  /// Outer product of a vector of functions with a vector of functions using sparsity
1102 
1103  /// \tparam T type parameter for first factor
1104  /// \tparam R type parameter for second factor
1105  /// \tparam NDIM dimension of first and second factors
1106  /// \param world the world
1107  /// \param f first vector of functions
1108  /// \param g second vector of functions
1109  /// \param tol threshold for multiplication
1110  /// \param fence force fence (will always fence if necessary)
1111  /// \param symm if true, only compute f(i) * g(j) for j<=i
1112  /// \return fg(i,j) = f(i) * g(j), as a vector of vectors
1113  template <typename T, typename R, std::size_t NDIM>
1114  std::vector<std::vector<Function<TENSOR_RESULT_TYPE(T, R), NDIM> > >
1116  const std::vector<Function<R, NDIM> > &f,
1117  const std::vector<Function<R, NDIM> > &g,
1118  double tol,
1119  bool fence = true,
1120  bool symm = false) {
1121  PROFILE_BLOCK(Vmulsp);
1122  bool same=(&f == &g);
1123  reconstruct(world, f, false);
1124  if (not same) reconstruct(world, g, false);
1125  world.gop.fence();
1126  for (auto& ff : f) ff.norm_tree(false);
1127  if (not same) for (auto& gg : g) gg.norm_tree(false);
1128  world.gop.fence();
1129 
1130  std::vector<std::vector<Function<R,NDIM> > >result(f.size());
1131  std::vector<Function<R,NDIM>> g_i;
1132  for (int64_t i=f.size()-1; i>=0; --i) {
1133  if (!symm)
1134  result[i]= vmulXX(f[i], g, tol, false);
1135  else {
1136  if (g_i.empty()) g_i = g;
1137  g_i.resize(i+1); // this shrinks g_i down to single function for i=0
1138  result[i]= vmulXX(f[i], g_i, tol, false);
1139  }
1140  }
1141  if (fence) world.gop.fence();
1142  return result;
1143  }
1144 
1145  /// Makes the norm tree for all functions in a vector
1146  template <typename T, std::size_t NDIM>
1147  void norm_tree(World& world,
1148  const std::vector< Function<T,NDIM> >& v,
1149  bool fence=true)
1150  {
1151  PROFILE_BLOCK(Vnorm_tree);
1152  for (unsigned int i=0; i<v.size(); ++i) {
1153  v[i].norm_tree(false);
1154  }
1155  if (fence) world.gop.fence();
1156  }
1157 
1158  /// Multiplies two vectors of functions q[i] = a[i] * b[i]
1159  template <typename T, typename R, std::size_t NDIM>
1160  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1161  mul(World& world,
1162  const std::vector< Function<T,NDIM> >& a,
1163  const std::vector< Function<R,NDIM> >& b,
1164  bool fence=true) {
1165  PROFILE_BLOCK(Vmulvv);
1166  reconstruct(world, a, true);
1167  reconstruct(world, b, true);
1168 // if (&a != &b) reconstruct(world, b, true); // fails if type(a) != type(b)
1169 
1170  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > q(a.size());
1171  for (unsigned int i=0; i<a.size(); ++i) {
1172  q[i] = mul(a[i], b[i], false);
1173  }
1174  if (fence) world.gop.fence();
1175  return q;
1176  }
1177 
1178 
1179  /// multiply a high-dimensional function with a low-dimensional function
1180 
1181  /// @param[in] f NDIM function of NDIM dimensions
1182  /// @param[in] g LDIM function of LDIM
1183  /// @param[in] v dimension indices of f to multiply
1184  /// @return h[i](0,1,2,3) = f(0,1,2,3) * g[i](1,2,3) for v={1,2,3}
1185  template<typename T, std::size_t NDIM, std::size_t LDIM>
1186  std::vector<Function<T,NDIM> > partial_mul(const Function<T,NDIM> f, const std::vector<Function<T,LDIM> > g,
1187  const int particle) {
1188 
1189  World& world=f.world();
1190  std::vector<Function<T,NDIM> > result(g.size());
1191  for (auto& r : result) r.set_impl(f, false);
1192 
1193  FunctionImpl<T,NDIM>* fimpl=f.get_impl().get();
1194 // fimpl->make_redundant(false);
1195  fimpl->change_tree_state(redundant,false);
1196  make_redundant(world,g,false);
1197  world.gop.fence();
1198 
1199  for (std::size_t i=0; i<result.size(); ++i) {
1200  FunctionImpl<T,LDIM>* gimpl=g[i].get_impl().get();
1201  result[i].get_impl()->multiply(fimpl,gimpl,particle); // stupid naming inconsistency
1202  }
1203  world.gop.fence();
1204 
1205  fimpl->undo_redundant(false);
1206  for (auto& ig : g) ig.get_impl()->undo_redundant(false);
1207  world.gop.fence();
1208  return result;
1209  }
1210 
1211  template<typename T, std::size_t NDIM, std::size_t LDIM>
1212  std::vector<Function<T,NDIM> > multiply(const Function<T,NDIM> f, const std::vector<Function<T,LDIM> > g,
1213  const std::tuple<int,int,int> v) {
1214  return partial_mul<T,NDIM,LDIM>(f,g,std::array<int,3>({std::get<0>(v),std::get<1>(v),std::get<2>(v)}));
1215  }
1216 
1217 
1218 /// Computes the square of a vector of functions --- q[i] = v[i]**2
1219  template <typename T, std::size_t NDIM>
1220  std::vector< Function<T,NDIM> >
1221  square(World& world,
1222  const std::vector< Function<T,NDIM> >& v,
1223  bool fence=true) {
1224  return mul<T,T,NDIM>(world, v, v, fence);
1225 // std::vector< Function<T,NDIM> > vsq(v.size());
1226 // for (unsigned int i=0; i<v.size(); ++i) {
1227 // vsq[i] = square(v[i], false);
1228 // }
1229 // if (fence) world.gop.fence();
1230 // return vsq;
1231  }
1232 
1233 
1234  /// Computes the square of a vector of functions --- q[i] = abs(v[i])**2
1235  template <typename T, std::size_t NDIM>
1236  std::vector< Function<typename Tensor<T>::scalar_type,NDIM> >
1237  abssq(World& world,
1238  const std::vector< Function<T,NDIM> >& v,
1239  bool fence=true) {
1240  typedef typename Tensor<T>::scalar_type scalartype;
1241  reconstruct(world,v);
1242  std::vector<Function<scalartype,NDIM> > result(v.size());
1243  for (size_t i=0; i<v.size(); ++i) result[i]=abs_square(v[i],false);
1244  if (fence) world.gop.fence();
1245  return result;
1246  }
1247 
1248 
1249  /// Sets the threshold in a vector of functions
1250  template <typename T, std::size_t NDIM>
1251  void set_thresh(World& world, std::vector< Function<T,NDIM> >& v, double thresh, bool fence=true) {
1252  for (unsigned int j=0; j<v.size(); ++j) {
1253  v[j].set_thresh(thresh,false);
1254  }
1255  if (fence) world.gop.fence();
1256  }
1257 
1258  /// Returns the complex conjugate of the vector of functions
1259  template <typename T, std::size_t NDIM>
1260  std::vector< Function<T,NDIM> >
1261  conj(World& world,
1262  const std::vector< Function<T,NDIM> >& v,
1263  bool fence=true) {
1264  PROFILE_BLOCK(Vconj);
1265  std::vector< Function<T,NDIM> > r = copy(world, v); // Currently don't have oop conj
1266  for (unsigned int i=0; i<v.size(); ++i) {
1267  r[i].conj(false);
1268  }
1269  if (fence) world.gop.fence();
1270  return r;
1271  }
1272 
1273  /// Returns a deep copy of a vector of functions
1274  template <typename T, typename R, std::size_t NDIM>
1275  std::vector< Function<R,NDIM> > convert(World& world,
1276  const std::vector< Function<T,NDIM> >& v, bool fence=true) {
1277  PROFILE_BLOCK(Vcopy);
1278  std::vector< Function<R,NDIM> > r(v.size());
1279  for (unsigned int i=0; i<v.size(); ++i) {
1280  r[i] = convert<T,R,NDIM>(v[i], false);
1281  }
1282  if (fence) world.gop.fence();
1283  return r;
1284  }
1285 
1286 
1287  /// Returns a deep copy of a vector of functions
1288  template <typename T, std::size_t NDIM>
1289  std::vector< Function<T,NDIM> >
1290  copy(World& world,
1291  const std::vector< Function<T,NDIM> >& v,
1292  bool fence=true) {
1293  PROFILE_BLOCK(Vcopy);
1294  std::vector< Function<T,NDIM> > r(v.size());
1295  for (unsigned int i=0; i<v.size(); ++i) {
1296  r[i] = copy(v[i], false);
1297  }
1298  if (fence) world.gop.fence();
1299  return r;
1300  }
1301 
1302 
1303  /// Returns a deep copy of a vector of functions
1304  template <typename T, std::size_t NDIM>
1305  std::vector< Function<T,NDIM> >
1306  copy(const std::vector< Function<T,NDIM> >& v, bool fence=true) {
1307  PROFILE_BLOCK(Vcopy);
1308  std::vector< Function<T,NDIM> > r(v.size());
1309  if (v.size()>0) r=copy(v.front().world(),v,fence);
1310  return r;
1311  }
1312 
1313  /// Returns a vector of `n` deep copies of a function
1314  template <typename T, std::size_t NDIM>
1315  std::vector< Function<T,NDIM> >
1316  copy(World& world,
1317  const Function<T,NDIM>& v,
1318  const unsigned int n,
1319  bool fence=true) {
1320  PROFILE_BLOCK(Vcopy1);
1321  std::vector< Function<T,NDIM> > r(n);
1322  for (unsigned int i=0; i<n; ++i) {
1323  r[i] = copy(v, false);
1324  }
1325  if (fence) world.gop.fence();
1326  return r;
1327  }
1328 
1329  /// Create a new copy of the function with different distribution and optional
1330  /// fence
1331 
1332  /// Works in either basis. Different distributions imply
1333  /// asynchronous communication and the optional fence is
1334  /// collective.
1335  //
1336  /// Returns a deep copy of a vector of functions
1337 
1338  template <typename T, std::size_t NDIM>
1339  std::vector<Function<T, NDIM>> copy(World& world,
1340  const std::vector<Function<T, NDIM>>& v,
1341  const std::shared_ptr<WorldDCPmapInterface<Key<NDIM>>>& pmap,
1342  bool fence = true) {
1343  PROFILE_BLOCK(Vcopy);
1344  std::vector<Function<T, NDIM>> r(v.size());
1345  for (unsigned int i = 0; i < v.size(); ++i) {
1346  r[i] = copy(v[i], pmap, false);
1347  }
1348  if (fence) world.gop.fence();
1349  return r;
1350  }
1351 
1352  /// Returns new vector of functions --- q[i] = a[i] + b[i]
1353  template <typename T, typename R, std::size_t NDIM>
1354  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1355  add(World& world,
1356  const std::vector< Function<T,NDIM> >& a,
1357  const std::vector< Function<R,NDIM> >& b,
1358  bool fence=true) {
1359  PROFILE_BLOCK(Vadd);
1360  MADNESS_ASSERT(a.size() == b.size());
1361  compress(world, a);
1362  compress(world, b);
1363 
1364  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
1365  for (unsigned int i=0; i<a.size(); ++i) {
1366  r[i] = add(a[i], b[i], false);
1367  }
1368  if (fence) world.gop.fence();
1369  return r;
1370  }
1371 
1372  /// Returns new vector of functions --- q[i] = a + b[i]
1373  template <typename T, typename R, std::size_t NDIM>
1374  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1375  add(World& world,
1376  const Function<T,NDIM> & a,
1377  const std::vector< Function<R,NDIM> >& b,
1378  bool fence=true) {
1379  PROFILE_BLOCK(Vadd1);
1380  a.compress();
1381  compress(world, b);
1382 
1383  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(b.size());
1384  for (unsigned int i=0; i<b.size(); ++i) {
1385  r[i] = add(a, b[i], false);
1386  }
1387  if (fence) world.gop.fence();
1388  return r;
1389  }
1390  template <typename T, typename R, std::size_t NDIM>
1391  inline std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1392  add(World& world,
1393  const std::vector< Function<R,NDIM> >& b,
1394  const Function<T,NDIM> & a,
1395  bool fence=true) {
1396  return add(world, a, b, fence);
1397  }
1398 
1399  /// Returns new vector of functions --- q[i] = a[i] - b[i]
1400  template <typename T, typename R, std::size_t NDIM>
1401  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1402  sub(World& world,
1403  const std::vector< Function<T,NDIM> >& a,
1404  const std::vector< Function<R,NDIM> >& b,
1405  bool fence=true) {
1406  PROFILE_BLOCK(Vsub);
1407  MADNESS_ASSERT(a.size() == b.size());
1408  compress(world, a);
1409  compress(world, b);
1410 
1411  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
1412  for (unsigned int i=0; i<a.size(); ++i) {
1413  r[i] = sub(a[i], b[i], false);
1414  }
1415  if (fence) world.gop.fence();
1416  return r;
1417  }
1418 
1419  /// Returns new function --- q = sum_i f[i]
1420  template <typename T, std::size_t NDIM>
1421  Function<T, NDIM> sum(World& world, const std::vector<Function<T,NDIM> >& f,
1422  bool fence=true) {
1423 
1424  compress(world, f);
1426 
1427  for (unsigned int i=0; i<f.size(); ++i) r.gaxpy(1.0,f[i],1.0,false);
1428  if (fence) world.gop.fence();
1429  return r;
1430  }
1431 
1432 
1433  /// Multiplies and sums two vectors of functions r = \sum_i a[i] * b[i]
1434  template <typename T, typename R, std::size_t NDIM>
1435  Function<TENSOR_RESULT_TYPE(T,R), NDIM>
1436  dot(World& world,
1437  const std::vector< Function<T,NDIM> >& a,
1438  const std::vector< Function<R,NDIM> >& b,
1439  bool fence=true) {
1440  MADNESS_CHECK(a.size()==b.size());
1441  return sum(world,mul(world,a,b,true),fence);
1442  }
1443 
1444 
1445  /// out-of-place gaxpy for two vectors: result[i] = alpha * a[i] + beta * b[i]
1446  template <typename T, typename Q, typename R, std::size_t NDIM>
1447  std::vector<Function<TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)),NDIM> >
1449  const std::vector< Function<T,NDIM> >& a,
1450  Q beta,
1451  const std::vector< Function<R,NDIM> >& b,
1452  bool fence=true) {
1453 
1454  MADNESS_ASSERT(a.size() == b.size());
1455  typedef TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)) resultT;
1456  if (a.size()==0) return std::vector<Function<resultT,NDIM> >();
1457 
1458  World& world=a[0].world();
1459  compress(world,a);
1460  compress(world,b);
1461  std::vector<Function<resultT,NDIM> > result(a.size());
1462  for (unsigned int i=0; i<a.size(); ++i) {
1463  result[i]=gaxpy_oop(alpha, a[i], beta, b[i], false);
1464  }
1465  if (fence) world.gop.fence();
1466  return result;
1467  }
1468 
1469 
1470  /// out-of-place gaxpy for a vectors and a function: result[i] = alpha * a[i] + beta * b
1471  template <typename T, typename Q, typename R, std::size_t NDIM>
1472  std::vector<Function<TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)),NDIM> >
1474  const std::vector< Function<T,NDIM> >& a,
1475  Q beta,
1476  const Function<R,NDIM>& b,
1477  bool fence=true) {
1478 
1479  typedef TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)) resultT;
1480  if (a.size()==0) return std::vector<Function<resultT,NDIM> >();
1481 
1482  World& world=a[0].world();
1483  compress(world,a);
1484  b.compress();
1485  std::vector<Function<resultT,NDIM> > result(a.size());
1486  for (unsigned int i=0; i<a.size(); ++i) {
1487  result[i]=gaxpy_oop(alpha, a[i], beta, b, false);
1488  }
1489  if (fence) world.gop.fence();
1490  return result;
1491  }
1492 
1493 
1494  /// Generalized A*X+Y for vectors of functions ---- a[i] = alpha*a[i] + beta*b[i]
1495  template <typename T, typename Q, typename R, std::size_t NDIM>
1496  void gaxpy(Q alpha, std::vector<Function<T,NDIM>>& a, Q beta, const std::vector<Function<R,NDIM>>& b, const bool fence) {
1497  if (a.size() == 0) return;
1498  World& world=a.front().world();
1499  gaxpy(world,alpha,a,beta,b,fence);
1500  }
1501 
1502  /// Generalized A*X+Y for vectors of functions ---- a[i] = alpha*a[i] + beta*b[i]
1503  template <typename T, typename Q, typename R, std::size_t NDIM>
1504  void gaxpy(World& world,
1505  Q alpha,
1506  std::vector< Function<T,NDIM> >& a,
1507  Q beta,
1508  const std::vector< Function<R,NDIM> >& b,
1509  bool fence=true) {
1510  PROFILE_BLOCK(Vgaxpy);
1511  MADNESS_ASSERT(a.size() == b.size());
1512  if (fence) {
1513  compress(a.front().world(), a);
1514  compress(b.front().world(), b);
1515  }
1516  for (const auto& aa : a) MADNESS_CHECK_THROW(aa.is_compressed(),"vector-gaxpy requires compressed functions");
1517  for (const auto& bb : b) MADNESS_CHECK_THROW(bb.is_compressed(),"vector-gaxpy requires compressed functions");
1518 
1519  for (unsigned int i=0; i<a.size(); ++i) {
1520  a[i].gaxpy(alpha, b[i], beta, false);
1521  }
1522  if (fence) world.gop.fence();
1523  }
1524 
1525 
1526  /// Applies a vector of operators to a vector of functions --- q[i] = apply(op[i],f[i])
1527  template <typename opT, typename R, std::size_t NDIM>
1528  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> >
1529  apply(World& world,
1530  const std::vector< std::shared_ptr<opT> >& op,
1531  const std::vector< Function<R,NDIM> > f) {
1532 
1533  PROFILE_BLOCK(Vapplyv);
1534  MADNESS_ASSERT(f.size()==op.size());
1535 
1536  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
1537 
1538 // reconstruct(world, f);
1539  make_nonstandard(world, ncf);
1540 
1541  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> > result(f.size());
1542  for (unsigned int i=0; i<f.size(); ++i) {
1543  result[i] = apply_only(*op[i], f[i], false);
1544  result[i].get_impl()->set_tree_state(nonstandard_after_apply);
1545  }
1546 
1547  world.gop.fence();
1548 
1549  standard(world, ncf, false); // restores promise of logical constness
1550  reconstruct(result);
1551  world.gop.fence();
1552 
1553  return result;
1554  }
1555 
1556 
1557  /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i])
1558  template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
1559  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1561  const std::vector< Function<R,NDIM> > f) {
1562  return apply(op.get_world(),op,f);
1563  }
1564 
1565 
1566  /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i])
1567  template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
1568  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1569  apply(World& world,
1571  const std::vector< Function<R,NDIM> > f) {
1572  PROFILE_BLOCK(Vapply);
1573 
1574  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
1575  bool print_timings=(NDIM==6) and (world.rank()==0) and op.print_timings;
1576 
1577  double wall0=wall_time();
1578 // reconstruct(world, f);
1579  make_nonstandard(world, ncf);
1580  double wall1=wall_time();
1581  if (print_timings) printf("timer: %20.20s %8.2fs\n", "make_nonstandard", wall1-wall0);
1582 
1583  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> > result(f.size());
1584  for (unsigned int i=0; i<f.size(); ++i) {
1585  result[i] = apply_only(op, f[i], false);
1586  }
1587 
1588  world.gop.fence();
1589 
1590  // restores promise of logical constness
1591  if (op.destructive()) {
1592  for (auto& ff : ncf) ff.clear(false);
1593  world.gop.fence();
1594  } else {
1595  reconstruct(world,f);
1596  }
1597 
1598  // svd-tensor requires some cleanup after apply
1599  if (result[0].get_impl()->get_tensor_type()==TT_2D) {
1600  for (auto& r : result) r.get_impl()->finalize_apply();
1601  }
1602 
1603  if (print_timings) {
1604  for (auto& r : result) r.get_impl()->print_timer();
1605  op.print_timer();
1606  }
1607  reconstruct(world, result);
1608 
1609  return result;
1610  }
1611 
1612  /// Normalizes a vector of functions --- v[i] = v[i].scale(1.0/v[i].norm2())
1613  template <typename T, std::size_t NDIM>
1614  void normalize(World& world, std::vector< Function<T,NDIM> >& v, bool fence=true) {
1615  PROFILE_BLOCK(Vnormalize);
1616  std::vector<double> nn = norm2s(world, v);
1617  for (unsigned int i=0; i<v.size(); ++i) v[i].scale(1.0/nn[i],false);
1618  if (fence) world.gop.fence();
1619  }
1620 
1621  template <typename T, std::size_t NDIM>
1622  void print_size(World &world, const std::vector<Function<T,NDIM> > &v, const std::string &msg = "vectorfunction" ){
1623  if(v.empty()){
1624  if(world.rank()==0) std::cout << "print_size: " << msg << " is empty" << std::endl;
1625  }else if(v.size()==1){
1626  v.front().print_size(msg);
1627  }else{
1628  for(auto x:v){
1629  x.print_size(msg);
1630  }
1631  }
1632  }
1633 
1634  // gives back the size in GB
1635  template <typename T, std::size_t NDIM>
1636  double get_size(World& world, const std::vector< Function<T,NDIM> >& v){
1637 
1638  if (v.empty()) return 0.0;
1639 
1640  const double d=sizeof(T);
1641  const double fac=1024*1024*1024;
1642 
1643  double size=0.0;
1644  for(unsigned int i=0;i<v.size();i++){
1645  if (v[i].is_initialized()) size+=v[i].size();
1646  }
1647 
1648  return size/fac*d;
1649 
1650  }
1651 
1652  // gives back the size in GB
1653  template <typename T, std::size_t NDIM>
1654  double get_size(const Function<T,NDIM> & f){
1655  const double d=sizeof(T);
1656  const double fac=1024*1024*1024;
1657  double size=f.size();
1658  return size/fac*d;
1659  }
1660 
1661  /// apply op on the input vector yielding an output vector of functions
1662 
1663  /// @param[in] op the operator working on vin
1664  /// @param[in] vin vector of input Functions; needs to be refined to common level!
1665  /// @return vector of output Functions vout = op(vin)
1666  template <typename T, typename opT, std::size_t NDIM>
1667  std::vector<Function<T,NDIM> > multi_to_multi_op_values(const opT& op,
1668  const std::vector< Function<T,NDIM> >& vin,
1669  const bool fence=true) {
1670  MADNESS_ASSERT(vin.size()>0);
1671  MADNESS_ASSERT(vin[0].is_initialized()); // might be changed
1672  World& world=vin[0].world();
1673  Function<T,NDIM> dummy;
1674  dummy.set_impl(vin[0], false);
1675  std::vector<Function<T,NDIM> > vout=zero_functions<T,NDIM>(world, op.get_result_size());
1676  for (auto& out : vout) out.set_impl(vin[0],false);
1677  dummy.multi_to_multi_op_values(op, vin, vout, fence);
1678  return vout;
1679  }
1680 
1681 
1682 
1683 
1684  // convenience operators
1685 
1686  /// result[i] = a[i] + b[i]
1687  template <typename T, std::size_t NDIM>
1688  std::vector<Function<T,NDIM> > operator+(const std::vector<Function<T,NDIM> >& lhs,
1689  const std::vector<Function<T,NDIM>>& rhs) {
1690  MADNESS_CHECK(lhs.size() == rhs.size());
1691  return gaxpy_oop(1.0,lhs,1.0,rhs);
1692  }
1693 
1694  /// result[i] = a[i] - b[i]
1695  template <typename T, std::size_t NDIM>
1696  std::vector<Function<T,NDIM> > operator-(const std::vector<Function<T,NDIM> >& lhs,
1697  const std::vector<Function<T,NDIM> >& rhs) {
1698  MADNESS_CHECK(lhs.size() == rhs.size());
1699  return gaxpy_oop(1.0,lhs,-1.0,rhs);
1700  }
1701 
1702  /// result[i] = a[i] + b
1703  template <typename T, std::size_t NDIM>
1704  std::vector<Function<T,NDIM> > operator+(const std::vector<Function<T,NDIM> >& lhs,
1705  const Function<T,NDIM>& rhs) {
1706  // MADNESS_CHECK(lhs.size() == rhs.size()); // no!!
1707  return gaxpy_oop(1.0,lhs,1.0,rhs);
1708  }
1709 
1710  /// result[i] = a[i] - b
1711  template <typename T, std::size_t NDIM>
1712  std::vector<Function<T,NDIM> > operator-(const std::vector<Function<T,NDIM> >& lhs,
1713  const Function<T,NDIM>& rhs) {
1714  // MADNESS_CHECK(lhs.size() == rhs.size()); // no
1715  return gaxpy_oop(1.0,lhs,-1.0,rhs);
1716  }
1717 
1718  /// result[i] = a + b[i]
1719  template <typename T, std::size_t NDIM>
1720  std::vector<Function<T,NDIM> > operator+(const Function<T,NDIM>& lhs,
1721  const std::vector<Function<T,NDIM> >& rhs) {
1722  // MADNESS_CHECK(lhs.size() == rhs.size()); // no
1723  return gaxpy_oop(1.0,rhs,1.0,lhs);
1724  }
1725 
1726  /// result[i] = a - b[i]
1727  template <typename T, std::size_t NDIM>
1728  std::vector<Function<T,NDIM> > operator-(const Function<T,NDIM>& lhs,
1729  const std::vector<Function<T,NDIM> >& rhs) {
1730 // MADNESS_CHECK(lhs.size() == rhs.size()); // no
1731  return gaxpy_oop(-1.0,rhs,1.0,lhs);
1732  }
1733 
1734 
1735  template <typename T, typename R, std::size_t NDIM>
1736  std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > operator*(const R fac,
1737  const std::vector<Function<T,NDIM> >& rhs) {
1738  if (rhs.size()>0) {
1739  std::vector<Function<T,NDIM> > tmp=copy(rhs[0].world(),rhs);
1740  scale(tmp[0].world(),tmp,TENSOR_RESULT_TYPE(T,R)(fac));
1741  return tmp;
1742  }
1743  return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1744  }
1745 
1746  template <typename T, typename R, std::size_t NDIM>
1747  std::vector<Function<T,NDIM> > operator*(const std::vector<Function<T,NDIM> >& rhs,
1748  const R fac) {
1749  if (rhs.size()>0) {
1750  std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > tmp=copy(rhs[0].world(),rhs);
1751  scale(tmp[0].world(),tmp,TENSOR_RESULT_TYPE(T,R)(fac));
1752  return tmp;
1753  }
1754  return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1755  }
1756 
1757  /// multiply a vector of functions with a function: r[i] = v[i] * a
1758  template <typename T, typename R, std::size_t NDIM>
1760  const std::vector<Function<R,NDIM> >& v) {
1761  if (v.size()>0) return mul(v[0].world(),a,v,true);
1762  return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1763  }
1764 
1765 
1766  /// multiply a vector of functions with a function: r[i] = a * v[i]
1767  template <typename T, typename R, std::size_t NDIM>
1768  std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > operator*(const std::vector<Function<T,NDIM> >& v,
1769  const Function<R,NDIM>& a) {
1770  if (v.size()>0) return mul(v[0].world(),a,v,true);
1771  return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1772  }
1773 
1774 
1775  template <typename T, std::size_t NDIM>
1776  std::vector<Function<T,NDIM> > operator+=(std::vector<Function<T,NDIM> >& lhs, const std::vector<Function<T,NDIM> >& rhs) {
1777  MADNESS_CHECK(lhs.size() == rhs.size());
1778  if (lhs.size() > 0) gaxpy(lhs.front().world(), 1.0, lhs, 1.0, rhs);
1779  return lhs;
1780  }
1781 
1782  template <typename T, std::size_t NDIM>
1783  std::vector<Function<T,NDIM> > operator-=(std::vector<Function<T,NDIM> >& lhs,
1784  const std::vector<Function<T,NDIM> >& rhs) {
1785  MADNESS_CHECK(lhs.size() == rhs.size());
1786  if (lhs.size() > 0) gaxpy(lhs.front().world(), 1.0, lhs, -1.0, rhs);
1787  return lhs;
1788  }
1789 
1790  /// return the real parts of the vector's function (if complex)
1791  template <typename T, std::size_t NDIM>
1792  std::vector<Function<typename Tensor<T>::scalar_type,NDIM> >
1793  real(const std::vector<Function<T,NDIM> >& v, bool fence=true) {
1794  std::vector<Function<typename Tensor<T>::scalar_type,NDIM> > result(v.size());
1795  for (std::size_t i=0; i<v.size(); ++i) result[i]=real(v[i],false);
1796  if (fence and result.size()>0) result[0].world().gop.fence();
1797  return result;
1798  }
1799 
1800  /// return the imaginary parts of the vector's function (if complex)
1801  template <typename T, std::size_t NDIM>
1802  std::vector<Function<typename Tensor<T>::scalar_type,NDIM> >
1803  imag(const std::vector<Function<T,NDIM> >& v, bool fence=true) {
1804  std::vector<Function<typename Tensor<T>::scalar_type,NDIM> > result(v.size());
1805  for (std::size_t i=0; i<v.size(); ++i) result[i]=imag(v[i],false);
1806  if (fence and result.size()>0) result[0].world().gop.fence();
1807  return result;
1808  }
1809 
1810  /// shorthand gradient operator
1811 
1812  /// returns the differentiated function f in all NDIM directions
1813  /// @param[in] f the function on which the grad operator works on
1814  /// @param[in] refine refinement before diff'ing makes the result more accurate
1815  /// @param[in] fence fence after completion; if reconstruction is needed always fence
1816  /// @return the vector \frac{\partial}{\partial x_i} f
1817  template <typename T, std::size_t NDIM>
1818  std::vector<Function<T,NDIM> > grad(const Function<T,NDIM>& f,
1819  bool refine=false, bool fence=true) {
1820 
1821  World& world=f.world();
1822  f.reconstruct();
1823  if (refine) f.refine(); // refine to make result more precise
1824 
1825  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1826  gradient_operator<T,NDIM>(world);
1827 
1828  std::vector<Function<T,NDIM> > result(NDIM);
1829  for (size_t i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1830  if (fence) world.gop.fence();
1831  return result;
1832  }
1833 
1834  // BLM first derivative
1835  template <typename T, std::size_t NDIM>
1836  std::vector<Function<T,NDIM> > grad_ble_one(const Function<T,NDIM>& f,
1837  bool refine=false, bool fence=true) {
1838 
1839  World& world=f.world();
1840  f.reconstruct();
1841  if (refine) f.refine(); // refine to make result more precise
1842 
1843  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1844  gradient_operator<T,NDIM>(world);
1845 
1846  // Read in new coeff for each operator
1847  for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_ble1();
1848 
1849  std::vector<Function<T,NDIM> > result(NDIM);
1850  for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1851  if (fence) world.gop.fence();
1852  return result;
1853  }
1854 
1855  // BLM second derivative
1856  template <typename T, std::size_t NDIM>
1857  std::vector<Function<T,NDIM> > grad_ble_two(const Function<T,NDIM>& f,
1858  bool refine=false, bool fence=true) {
1859 
1860  World& world=f.world();
1861  f.reconstruct();
1862  if (refine) f.refine(); // refine to make result more precise
1863 
1864  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1865  gradient_operator<T,NDIM>(world);
1866 
1867  // Read in new coeff for each operator
1868  for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_ble2();
1869 
1870  std::vector<Function<T,NDIM> > result(NDIM);
1871  for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1872  if (fence) world.gop.fence();
1873  return result;
1874  }
1875 
1876  // Bspline first derivative
1877  template <typename T, std::size_t NDIM>
1878  std::vector<Function<T,NDIM> > grad_bspline_one(const Function<T,NDIM>& f,
1879  bool refine=false, bool fence=true) {
1880 
1881  World& world=f.world();
1882  f.reconstruct();
1883  if (refine) f.refine(); // refine to make result more precise
1884 
1885  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1886  gradient_operator<T,NDIM>(world);
1887 
1888  // Read in new coeff for each operator
1889  for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline1();
1890 
1891  std::vector<Function<T,NDIM> > result(NDIM);
1892  for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1893  if (fence) world.gop.fence();
1894  return result;
1895  }
1896 
1897  // Bpsline second derivative
1898  template <typename T, std::size_t NDIM>
1899  std::vector<Function<T,NDIM> > grad_bpsline_two(const Function<T,NDIM>& f,
1900  bool refine=false, bool fence=true) {
1901 
1902  World& world=f.world();
1903  f.reconstruct();
1904  if (refine) f.refine(); // refine to make result more precise
1905 
1906  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1907  gradient_operator<T,NDIM>(world);
1908 
1909  // Read in new coeff for each operator
1910  for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline2();
1911 
1912  std::vector<Function<T,NDIM> > result(NDIM);
1913  for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1914  if (fence) world.gop.fence();
1915  return result;
1916  }
1917 
1918  // Bspline third derivative
1919  template <typename T, std::size_t NDIM>
1920  std::vector<Function<T,NDIM> > grad_bspline_three(const Function<T,NDIM>& f,
1921  bool refine=false, bool fence=true) {
1922 
1923  World& world=f.world();
1924  f.reconstruct();
1925  if (refine) f.refine(); // refine to make result more precise
1926 
1927  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1928  gradient_operator<T,NDIM>(world);
1929 
1930  // Read in new coeff for each operator
1931  for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline3();
1932 
1933  std::vector<Function<T,NDIM> > result(NDIM);
1934  for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1935  if (fence) world.gop.fence();
1936  return result;
1937  }
1938 
1939 
1940 
1941  /// shorthand div operator
1942 
1943  /// returns the dot product of nabla with a vector f
1944  /// @param[in] f the vector of functions on which the div operator works on
1945  /// @param[in] refine refinement before diff'ing makes the result more accurate
1946  /// @param[in] fence fence after completion; currently always fences
1947  /// @return the vector \frac{\partial}{\partial x_i} f
1948  /// TODO: add this to operator fusion
1949  template <typename T, std::size_t NDIM>
1950  Function<T,NDIM> div(const std::vector<Function<T,NDIM> >& v,
1951  bool do_refine=false, bool fence=true) {
1952 
1953  MADNESS_ASSERT(v.size()>0);
1954  World& world=v[0].world();
1955  reconstruct(world,v);
1956  if (do_refine) refine(world,v); // refine to make result more precise
1957 
1958  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1959  gradient_operator<T,NDIM>(world);
1960 
1961  std::vector<Function<T,NDIM> > result(NDIM);
1962  for (size_t i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),v[i],false);
1963  world.gop.fence();
1964  return sum(world,result,fence);
1965  }
1966 
1967  /// shorthand rot operator
1968 
1969  /// returns the cross product of nabla with a vector f
1970  /// @param[in] f the vector of functions on which the rot operator works on
1971  /// @param[in] refine refinement before diff'ing makes the result more accurate
1972  /// @param[in] fence fence after completion; currently always fences
1973  /// @return the vector \frac{\partial}{\partial x_i} f
1974  /// TODO: add this to operator fusion
1975  template <typename T, std::size_t NDIM>
1976  std::vector<Function<T,NDIM> > rot(const std::vector<Function<T,NDIM> >& v,
1977  bool do_refine=false, bool fence=true) {
1978 
1979  MADNESS_ASSERT(v.size()==3);
1980  World& world=v[0].world();
1981  reconstruct(world,v);
1982  if (do_refine) refine(world,v); // refine to make result more precise
1983 
1984  std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1985  gradient_operator<T,NDIM>(world);
1986 
1987  std::vector<Function<T,NDIM> > d(NDIM),dd(NDIM);
1988  d[0]=apply(*(grad[1]),v[2],false); // Dy z
1989  d[1]=apply(*(grad[2]),v[0],false); // Dz x
1990  d[2]=apply(*(grad[0]),v[1],false); // Dx y
1991  dd[0]=apply(*(grad[2]),v[1],false); // Dz y
1992  dd[1]=apply(*(grad[0]),v[2],false); // Dx z
1993  dd[2]=apply(*(grad[1]),v[0],false); // Dy x
1994  world.gop.fence();
1995 
1996  d[0].gaxpy(1.0,dd[0],-1.0,false);
1997  d[1].gaxpy(1.0,dd[1],-1.0,false);
1998  d[2].gaxpy(1.0,dd[2],-1.0,false);
1999 
2000  world.gop.fence();
2001  return d;
2002  }
2003 
2004  /// shorthand cross operator
2005 
2006  /// returns the cross product of vectors f and g
2007  /// @param[in] f the vector of functions on which the rot operator works on
2008  /// @param[in] g the vector of functions on which the rot operator works on
2009  /// @param[in] fence fence after completion; currently always fences
2010  /// @return the vector \frac{\partial}{\partial x_i} f
2011  /// TODO: add this to operator fusion
2012  template <typename T, typename R, std::size_t NDIM>
2013  std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > cross(const std::vector<Function<T,NDIM> >& f,
2014  const std::vector<Function<R,NDIM> >& g,
2015  bool do_refine=false, bool fence=true) {
2016 
2017  MADNESS_ASSERT(f.size()==3);
2018  MADNESS_ASSERT(g.size()==3);
2019  World& world=f[0].world();
2020  reconstruct(world,f,false);
2021  reconstruct(world,g);
2022 
2023  std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > d(f.size()),dd(f.size());
2024 
2025  d[0]=mul(f[1],g[2],false);
2026  d[1]=mul(f[2],g[0],false);
2027  d[2]=mul(f[0],g[1],false);
2028 
2029  dd[0]=mul(f[2],g[1],false);
2030  dd[1]=mul(f[0],g[2],false);
2031  dd[2]=mul(f[1],g[0],false);
2032  world.gop.fence();
2033 
2034  compress(world,d,false);
2035  compress(world,dd);
2036 
2037  d[0].gaxpy(1.0,dd[0],-1.0,false);
2038  d[1].gaxpy(1.0,dd[1],-1.0,false);
2039  d[2].gaxpy(1.0,dd[2],-1.0,false);
2040 
2041  world.gop.fence();
2042  return d;
2043  }
2044 
2045 
2046  /// load a vector of functions
2047  template<typename T, size_t NDIM>
2048  void load_function(World& world, std::vector<Function<T,NDIM> >& f,
2049  const std::string name) {
2050  if (world.rank()==0) print("loading vector of functions",name);
2052  std::size_t fsize=0;
2053  ar & fsize;
2054  f.resize(fsize);
2055  for (std::size_t i=0; i<fsize; ++i) ar & f[i];
2056  }
2057 
2058  /// save a vector of functions
2059  template<typename T, size_t NDIM>
2060  void save_function(const std::vector<Function<T,NDIM> >& f, const std::string name) {
2061  if (f.size()>0) {
2062  World& world=f.front().world();
2063  if (world.rank()==0) print("saving vector of functions",name);
2065  std::size_t fsize=f.size();
2066  ar & fsize;
2067  for (std::size_t i=0; i<fsize; ++i) ar & f[i];
2068  }
2069  }
2070 
2071 
2072 }
2073 #endif // MADNESS_MRA_VMRA_H__INCLUDED
double q(double t)
Definition: DKops.h:18
real_convolution_3d A(World &world)
Definition: DKops.h:230
Definition: test_ar.cc:118
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
Definition: distributed_matrix.h:68
Manages data associated with a row/column/block distributed array.
Definition: distributed_matrix.h:388
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
FunctionFactory & compressed(bool value=true)
Definition: function_factory.h:168
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition: funcimpl.h:941
void undo_redundant(const bool fence)
convert this from redundant to standard reconstructed form
Definition: mraimpl.h:1560
static Tensor< TENSOR_RESULT_TYPE(T, R) > inner_local(const std::vector< const FunctionImpl< T, NDIM > * > &left, const std::vector< const FunctionImpl< R, NDIM > * > &right, bool sym)
Definition: funcimpl.h:5722
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
A multiresolution adaptive numerical function.
Definition: mra.h:122
World & world() const
Returns the world.
Definition: mra.h:648
void set_impl(const std::shared_ptr< FunctionImpl< T, NDIM > > &impl)
Replace current FunctionImpl with provided new one.
Definition: mra.h:621
void multi_to_multi_op_values(const opT &op, const std::vector< Function< T, NDIM > > &vin, std::vector< Function< T, NDIM > > &vout, const bool fence=true)
apply op on the input vector yielding an output vector of functions
Definition: mra.h:1521
bool is_initialized() const
Returns true if the function is initialized.
Definition: mra.h:147
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
Key is the index for a node of the 2^NDIM-tree.
Definition: key.h:66
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
A tensor is a multidimension array.
Definition: tensor.h:317
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition: tensor.h:409
Interface to be provided by any process map.
Definition: worlddc.h:82
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:870
void fence()
Returns after all local tasks have completed.
Definition: world_task_queue.h:1384
A parallel world class.
Definition: world.h:132
WorldTaskQueue & taskq
Task queue.
Definition: world.h:204
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
An archive for storing local or parallel data, wrapping a BinaryFstreamInputArchive.
Definition: parallel_archive.h:366
An archive for storing local or parallel data wrapping a BinaryFstreamOutputArchive.
Definition: parallel_archive.h:321
int integer
Definition: crayio.c:25
static const double R
Definition: csqrt.cc:46
Declaration and initialization of tree traversal functions and generic derivative.
static double lo
Definition: dirac-hatom.cc:23
const double m
Definition: gfit.cc:199
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
const double beta
Definition: gygi_soltion.cc:62
static const double v
Definition: hatom_sf_dirac.cc:20
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
#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_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
Main include file for MADNESS and defines Function interface.
static const bool VERIFY_TREE
Definition: mra.h:57
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
void save_function(const std::vector< Function< T, NDIM > > &f, const std::string name)
save a vector of functions
Definition: vmra.h:2060
std::vector< Function< T, NDIM > > grad_bpsline_two(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1899
Function< TENSOR_RESULT_TYPE(L, R), NDIM > sub(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator- but with optional fence and no automatic compression.
Definition: mra.h:1971
std::vector< CCPairFunction< T, NDIM > > operator+(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1047
void rr_cholesky(Tensor< T > &A, typename Tensor< T >::scalar_type tol, Tensor< integer > &piv, int &rank)
Compute the rank-revealing Cholesky factorization.
Definition: lapack.cc:1203
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
Function< TENSOR_RESULT_TYPE(T, R), NDIM > dot(World &world, const std::vector< Function< T, NDIM > > &a, const std::vector< Function< R, NDIM > > &b, bool fence=true)
Multiplies and sums two vectors of functions r = \sum_i a[i] * b[i].
Definition: vmra.h:1436
std::vector< Function< T, NDIM > > orthonormalize_canonical(const std::vector< Function< T, NDIM > > &v, const Tensor< T > &ovlp, double lindep)
Definition: vmra.h:501
std::vector< Function< T, NDIM > > append(const std::vector< Function< T, NDIM > > &lhs, const std::vector< Function< T, NDIM > > &rhs)
combine two vectors
Definition: vmra.h:649
std::vector< Function< T, NDIM > > orthonormalize_rrcd(const std::vector< Function< T, NDIM > > &v, Tensor< T > &ovlp, const double tol, Tensor< integer > &piv, int &rank)
Definition: vmra.h:594
Function< TENSOR_RESULT_TYPE(L, R), NDIM > add(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator+ but with optional fence and no automatic compression.
Definition: mra.h:1927
Function< TENSOR_RESULT_TYPE(Q, T), NDIM > mul(const Q alpha, const Function< T, NDIM > &f, bool fence=true)
Returns new function equal to alpha*f(x) with optional fence.
Definition: mra.h:1701
Function< T, NDIM > div(const std::vector< Function< T, NDIM > > &v, bool do_refine=false, bool fence=true)
shorthand div operator
Definition: vmra.h:1950
void set_impl(std::vector< Function< T, NDIM >> &v, const std::vector< std::shared_ptr< FunctionImpl< T, NDIM >>> vimpl)
Definition: vmra.h:670
response_space scale(response_space a, double b)
std::vector< CCPairFunction< T, NDIM > > operator-(const std::vector< CCPairFunction< T, NDIM >> c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition: ccpairfunction.h:1055
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
Function< double, NDIM > abssq(const Function< double_complex, NDIM > &z, bool fence=true)
Returns a new function that is the square of the absolute value of the input.
Definition: mra.h:2713
Function< typename TensorTypeData< Q >::scalar_type, NDIM > abs_square(const Function< Q, NDIM > &func)
Definition: complexfun.h:121
std::vector< Function< T, NDIM > > grad_bspline_three(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1920
std::vector< Function< T, NDIM > > orthonormalize_symmetric(const std::vector< Function< T, NDIM > > &v, const Tensor< T > &ovlp)
symmetric orthonormalization (see e.g. Szabo/Ostlund)
Definition: vmra.h:456
std::vector< Function< T, NDIM > > orthonormalize(const std::vector< Function< T, NDIM > > &vf_in)
orthonormalize the vectors
Definition: vmra.h:421
std::vector< Function< T, NDIM > > grad_ble_two(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1857
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
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
tensorT Q2(const tensorT &s)
Given overlap matrix, return rotation with 2nd order error to orthonormalize the vectors.
Definition: SCF.cc:135
const std::vector< Function< T, NDIM > > & change_tree_state(const std::vector< Function< T, NDIM >> &v, const TreeState finalstate, const bool fence=true)
change tree state of the functions
Definition: vmra.h:277
TreeState
Definition: funcdefaults.h:58
@ nonstandard_after_apply
s and d coeffs, state after operator application
Definition: funcdefaults.h:63
@ redundant_after_merge
s coeffs everywhere, must be summed up to yield the result
Definition: funcdefaults.h:65
@ reconstructed
s coeffs at the leaves only
Definition: funcdefaults.h:59
@ nonstandard
s and d coeffs in internal nodes
Definition: funcdefaults.h:61
@ 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
void cholesky(Tensor< T > &A)
Compute the Cholesky factorization.
Definition: lapack.cc:1174
std::vector< Function< T, NDIM > > zero_functions(World &world, int n, bool fence=true)
Generates a vector of zero functions (reconstructed)
Definition: vmra.h:395
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
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:133
Function< Q, NDIM > convert(const Function< T, NDIM > &f, bool fence=true)
Type conversion implies a deep copy. No communication except for optional fence.
Definition: mra.h:2032
response_space transpose(response_space &f)
Definition: basic_operators.cc:10
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition: mra.h:2681
std::vector< std::shared_ptr< FunctionImpl< T, NDIM > > > get_impl(const std::vector< Function< T, NDIM >> &v)
Definition: vmra.h:663
void set_thresh(World &world, std::vector< Function< T, NDIM > > &v, double thresh, bool fence=true)
Sets the threshold in a vector of functions.
Definition: vmra.h:1251
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
static void verify_tree(World &world, const std::vector< Function< T, NDIM > > &v)
Definition: SCF.cc:72
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform_reconstructed(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:718
std::vector< Function< T, NDIM > > zero_functions_compressed(World &world, int n, bool fence=true)
Generates a vector of zero functions (compressed)
Definition: vmra.h:409
static const Slice _(0,-1, 1)
std::vector< Function< T, NDIM > > flatten(const std::vector< std::vector< Function< T, NDIM > > > &vv)
Definition: vmra.h:656
std::vector< Function< T, NDIM > > grad_ble_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1836
std::vector< Function< TENSOR_RESULT_TYPE(L, R), D > > vmulXX(const Function< L, D > &left, const std::vector< Function< R, D > > &vright, double tol, bool fence=true)
Use the vmra/mul(...) interface instead.
Definition: mra.h:1807
Function< TENSOR_RESULT_TYPE(L, R), NDIM > gaxpy_oop(TENSOR_RESULT_TYPE(L, R) alpha, const Function< L, NDIM > &left, TENSOR_RESULT_TYPE(L, R) beta, const Function< R, NDIM > &right, bool fence=true)
Returns new function alpha*left + beta*right optional fence and no automatic compression.
Definition: mra.h:1917
std::vector< Function< T, NDIM > > rot(const std::vector< Function< T, NDIM > > &v, bool do_refine=false, bool fence=true)
shorthand rot operator
Definition: vmra.h:1976
TENSOR_RESULT_TYPE(T, R) inner(const Function< T
Computes the scalar/inner product between two functions.
Definition: vmra.h:1059
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
std::vector< Function< T, NDIM > > orthonormalize_cd(const std::vector< Function< T, NDIM > > &v, Tensor< T > &ovlp)
Definition: vmra.h:557
@ TT_2D
Definition: gentensor.h:120
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
void print_size(World &world, const std::vector< Function< T, NDIM > > &v, const std::string &msg="vectorfunction")
Definition: vmra.h:1622
std::vector< CCPairFunction< T, NDIM > > & operator+=(std::vector< CCPairFunction< T, NDIM > > &lhs, const CCPairFunction< T, NDIM > &rhs)
Definition: ccpairfunction.h:1063
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
std::vector< CCPairFunction< T, NDIM > > & operator-=(std::vector< CCPairFunction< T, NDIM > > &rhs, const std::vector< CCPairFunction< T, NDIM > > &lhs)
Definition: ccpairfunction.h:1077
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition: mra.h:1742
std::vector< Function< T, NDIM > > partial_mul(const Function< T, NDIM > f, const std::vector< Function< T, LDIM > > g, const int particle)
multiply a high-dimensional function with a low-dimensional function
Definition: vmra.h:1186
std::vector< std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > > matrix_mul_sparse(World &world, const std::vector< Function< R, NDIM > > &f, const std::vector< Function< R, NDIM > > &g, double tol, bool fence=true, bool symm=false)
Outer product of a vector of functions with a vector of functions using sparsity.
Definition: vmra.h:1115
NDIM const Function< R, NDIM > & g
Definition: mra.h:2416
std::vector< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition: ccpairfunction.h:1084
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition: timers.cc:48
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition: vmra.h:1818
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition: lapack.cc:832
std::vector< Function< T, NDIM > > grad_bspline_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition: vmra.h:1878
std::vector< Function< T, NDIM > > reduce_rank(std::vector< Function< T, NDIM > > v, double thresh=0.0, bool fence=true)
reduces the tensor rank of the coefficient tensor (if applicable)
Definition: vmra.h:366
Tensor< TENSOR_RESULT_TYPE(T, R) > matrix_inner_old(World &world, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< R, NDIM > > &g, bool sym=false)
Computes the matrix inner product of two function vectors - q(i,j) = inner(f[i],g[j])
Definition: vmra.h:970
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
double imag(double x)
Definition: complexfun.h:56
std::vector< Function< T, NDIM > > impl2function(const std::vector< std::shared_ptr< FunctionImpl< T, NDIM >>> vimpl)
Definition: vmra.h:676
void load_function(World &world, std::vector< Function< T, NDIM > > &f, const std::string name)
load a vector of functions
Definition: vmra.h:2048
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
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition: vmra.h:1614
Tensor< double > norm2s_T(World &world, const std::vector< Function< T, NDIM >> &v)
Computes the 2-norms of a vector of functions.
Definition: vmra.h:839
double real(double x)
Definition: complexfun.h:52
std::vector< Function< T, NDIM > > multi_to_multi_op_values(const opT &op, const std::vector< Function< T, NDIM > > &vin, const bool fence=true)
apply op on the input vector yielding an output vector of functions
Definition: vmra.h:1667
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_only(const opT &op, const Function< R, NDIM > &f, bool fence=true)
Apply operator ONLY in non-standard form - required other steps missing !!
Definition: mra.h:2120
std::string name(const FuncType &type, const int ex=-1)
Definition: ccpairfunction.h:28
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
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition: distpm.cc:46
double get_size(World &world, const std::vector< Function< T, NDIM > > &v)
Definition: vmra.h:1636
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition: vmra.h:827
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:969
void make_nonstandard(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Generates non-standard form of a vector of functions.
Definition: vmra.h:245
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > cross(const std::vector< Function< T, NDIM > > &f, const std::vector< Function< R, NDIM > > &g, bool do_refine=false, bool fence=true)
shorthand cross operator
Definition: vmra.h:2013
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
int distance(const madness::Hash_private::HashIterator< hashT > &it, const madness::Hash_private::HashIterator< hashT > &jt)
Definition: worldhashmap.h:616
static long abs(long a)
Definition: tensor.h:218
static const double b
Definition: nonlinschro.cc:119
static const double a
Definition: nonlinschro.cc:118
double Q(double a)
Definition: relops.cc:20
static const double c
Definition: relops.cc:10
static const double L
Definition: rk.cc:46
static const double thresh
Definition: rk.cc:45
Definition: test_ar.cc:204
Definition: lowrankfunction.h:332
Definition: dirac-hatom.cc:108
int P
Definition: test_binsorter.cc:9
const double alpha
Definition: test_coulomb.cc:54
void d()
Definition: test_sig.cc:79
double aa
Definition: testbsh.cc:68
static const std::size_t NDIM
Definition: testpdiff.cc:42
#define PROFILE_FUNC
Definition: worldprofile.h:209
#define PROFILE_BLOCK(name)
Definition: worldprofile.h:208