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>
125#include <cstdio>
126
127namespace 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
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) {
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>
190 const std::vector< Function<T,NDIM> >& v,
191 bool fence=true) {
192
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>
246 std::vector< Function<T,NDIM> >& v,
247 bool fence= true) {
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) {
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
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,
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) {
309 }
310 if (finalstate==nonstandard) {
313 }
318 }
319 if (finalstate==redundant) {
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) {
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) {
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) {
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
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
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
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,
597 const double tol,
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
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) {
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> >
690 const std::vector< Function<T,NDIM> >& v,
691 const Tensor<R>& c,
692 bool fence=true) {
693
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
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.
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) {
760 MADNESS_ASSERT(v.size() == (unsigned int)(c.dim(0)));
761
762 std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> > vresult
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> >
773 const std::vector< Function<T,NDIM> >& v,
774 const DistributedMatrix<R>& c,
775 bool fence=true) {
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) {
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) {
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) {
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) {
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) {
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 {
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) {
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) {
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) {
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>
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) {
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) {
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) {
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 {
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) {
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> >
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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);
1425 Function<T,NDIM> r=FunctionFactory<T,NDIM>(world).compressed();
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 template <typename T, std::size_t NDIM>
1434 const std::vector<Function<T, NDIM>>& f,
1435 const std::vector<Function<T, NDIM>>& g,
1436 bool sym=false)
1437 {
1440 const int64_t n = A.coldim();
1441 const int64_t m = A.rowdim();
1442 MADNESS_ASSERT(int64_t(f.size()) == n && int64_t(g.size()) == m);
1443
1444 // Assume we can always create an ichunk*jchunk matrix locally
1445 const int ichunk = 1000;
1446 const int jchunk = 1000; // 1000*1000*8 = 8 MBytes
1447 for (int64_t ilo = 0; ilo < n; ilo += ichunk) {
1448 int64_t ihi = std::min(ilo + ichunk, n);
1449 std::vector<Function<T, NDIM>> ivec(f.begin() + ilo, f.begin() + ihi);
1450 for (int64_t jlo = 0; jlo < m; jlo += jchunk) {
1451 int64_t jhi = std::min(jlo + jchunk, m);
1452 std::vector<Function<T, NDIM>> jvec(g.begin() + jlo, g.begin() + jhi);
1453
1454 Tensor<T> P = matrix_dot(A.get_world(), ivec, jvec, sym);
1455 A.copy_from_replicated_patch(ilo, ihi - 1, jlo, jhi - 1, P);
1456 }
1457 }
1458 return A;
1459 }
1460
1461 /// Computes the matrix inner product of two function vectors - q(i,j) = inner(f[i],g[j])
1462
1463 /// For complex types symmetric is interpreted as Hermitian.
1464 ///
1465 /// The current parallel loop is non-optimal but functional.
1466 template <typename T, typename R, std::size_t NDIM>
1468 const std::vector<Function<T, NDIM>>& f,
1469 const std::vector<Function<R, NDIM>>& g,
1470 bool sym=false)
1471 {
1472 world.gop.fence();
1473 compress(world, f);
1474 // if ((void*)(&f) != (void*)(&g)) compress(world, g);
1475 compress(world, g);
1476
1477 std::vector<const FunctionImpl<T, NDIM>*> left(f.size());
1478 std::vector<const FunctionImpl<R, NDIM>*> right(g.size());
1479 for (unsigned int i = 0; i < f.size(); i++) left[i] = f[i].get_impl().get();
1480 for (unsigned int i = 0; i < g.size(); i++) right[i] = g[i].get_impl().get();
1481
1483
1484 world.gop.fence();
1485 world.gop.sum(r.ptr(), f.size() * g.size());
1486
1487 return r;
1488 }
1489
1490 /// Multiplies and sums two vectors of functions r = \sum_i a[i] * b[i]
1491 template <typename T, typename R, std::size_t NDIM>
1492 Function<TENSOR_RESULT_TYPE(T,R), NDIM>
1493 dot(World& world,
1494 const std::vector< Function<T,NDIM> >& a,
1495 const std::vector< Function<R,NDIM> >& b,
1496 bool fence=true) {
1497 MADNESS_CHECK(a.size()==b.size());
1498 return sum(world,mul(world,a,b,true),fence);
1499 }
1500
1501
1502 /// out-of-place gaxpy for two vectors: result[i] = alpha * a[i] + beta * b[i]
1503 template <typename T, typename Q, typename R, std::size_t NDIM>
1504 std::vector<Function<TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)),NDIM> >
1506 const std::vector< Function<T,NDIM> >& a,
1507 Q beta,
1508 const std::vector< Function<R,NDIM> >& b,
1509 bool fence=true) {
1510
1511 MADNESS_ASSERT(a.size() == b.size());
1512 typedef TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)) resultT;
1513 if (a.size()==0) return std::vector<Function<resultT,NDIM> >();
1514
1515 World& world=a[0].world();
1516 std::vector<Function<resultT,NDIM> > result(a.size());
1517 if (NDIM<=3) {
1518 compress(world,a);
1519 compress(world,b);
1520 for (unsigned int i=0; i<a.size(); ++i) {
1521 result[i]=gaxpy_oop(alpha, a[i], beta, b[i], false);
1522 }
1523 } else {
1524 reconstruct(world,a);
1525 reconstruct(world,b);
1526 for (unsigned int i=0; i<a.size(); ++i) {
1527 result[i]=gaxpy_oop_reconstructed(alpha, a[i], beta, b[i], false);
1528 }
1529 }
1530
1531 if (fence) world.gop.fence();
1532 return result;
1533 }
1534
1535
1536 /// out-of-place gaxpy for a vectors and a function: result[i] = alpha * a[i] + beta * b
1537 template <typename T, typename Q, typename R, std::size_t NDIM>
1538 std::vector<Function<TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)),NDIM> >
1540 const std::vector< Function<T,NDIM> >& a,
1541 Q beta,
1542 const Function<R,NDIM>& b,
1543 bool fence=true) {
1544
1545 typedef TENSOR_RESULT_TYPE(Q,TENSOR_RESULT_TYPE(T,R)) resultT;
1546 if (a.size()==0) return std::vector<Function<resultT,NDIM> >();
1547
1548 World& world=a[0].world();
1549 compress(world,a);
1550 b.compress();
1551 std::vector<Function<resultT,NDIM> > result(a.size());
1552 for (unsigned int i=0; i<a.size(); ++i) {
1553 result[i]=gaxpy_oop(alpha, a[i], beta, b, false);
1554 }
1555 if (fence) world.gop.fence();
1556 return result;
1557 }
1558
1559
1560 /// Generalized A*X+Y for vectors of functions ---- a[i] = alpha*a[i] + beta*b[i]
1561 template <typename T, typename Q, typename R, std::size_t NDIM>
1562 void gaxpy(Q alpha, std::vector<Function<T,NDIM>>& a, Q beta, const std::vector<Function<R,NDIM>>& b, const bool fence) {
1563 if (a.size() == 0) return;
1564 World& world=a.front().world();
1565 gaxpy(world,alpha,a,beta,b,fence);
1566 }
1567
1568 /// Generalized A*X+Y for vectors of functions ---- a[i] = alpha*a[i] + beta*b[i]
1569 template <typename T, typename Q, typename R, std::size_t NDIM>
1570 void gaxpy(World& world,
1571 Q alpha,
1572 std::vector< Function<T,NDIM> >& a,
1573 Q beta,
1574 const std::vector< Function<R,NDIM> >& b,
1575 bool fence=true) {
1577 MADNESS_ASSERT(a.size() == b.size());
1578 if (fence) {
1579 compress(a.front().world(), a);
1580 compress(b.front().world(), b);
1581 }
1582 for (const auto& aa : a) MADNESS_CHECK_THROW(aa.is_compressed(),"vector-gaxpy requires compressed functions");
1583 for (const auto& bb : b) MADNESS_CHECK_THROW(bb.is_compressed(),"vector-gaxpy requires compressed functions");
1584
1585 for (unsigned int i=0; i<a.size(); ++i) {
1586 a[i].gaxpy(alpha, b[i], beta, false);
1587 }
1588 if (fence) world.gop.fence();
1589 }
1590
1591
1592 /// Applies a vector of operators to a vector of functions --- q[i] = apply(op[i],f[i])
1593 template <typename opT, typename R, std::size_t NDIM>
1594 std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> >
1595 apply(World& world,
1596 const std::vector< std::shared_ptr<opT> >& op,
1597 const std::vector< Function<R,NDIM> > f) {
1598
1600 MADNESS_ASSERT(f.size()==op.size());
1601
1602 std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
1603
1604// reconstruct(world, f);
1605 make_nonstandard(world, ncf);
1606
1607 std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> > result(f.size());
1608 for (unsigned int i=0; i<f.size(); ++i) {
1609 result[i] = apply_only(*op[i], f[i], false);
1610 result[i].get_impl()->set_tree_state(nonstandard_after_apply);
1611 }
1612
1613 world.gop.fence();
1614
1615 standard(world, ncf, false); // restores promise of logical constness
1616 reconstruct(result);
1617 world.gop.fence();
1618
1619 return result;
1620 }
1621
1622
1623 /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i])
1624 template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
1625 std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1627 const std::vector< Function<R,NDIM> > f) {
1628 return apply(op.get_world(),op,f);
1629 }
1630
1631
1632 /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i])
1633 template <typename T, typename R, std::size_t NDIM, std::size_t KDIM>
1634 std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
1635 apply(World& world,
1637 const std::vector< Function<R,NDIM> > f) {
1639
1640 std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
1641 bool print_timings=(NDIM==6) and (world.rank()==0) and op.print_timings;
1642
1643 double wall0=wall_time();
1644// reconstruct(world, f);
1645 make_nonstandard(world, ncf);
1646 double wall1=wall_time();
1647 if (print_timings) printf("timer: %20.20s %8.2fs\n", "make_nonstandard", wall1-wall0);
1648
1649 std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> > result(f.size());
1650 for (unsigned int i=0; i<f.size(); ++i) {
1651 result[i] = apply_only(op, f[i], false);
1652 }
1653
1654 world.gop.fence();
1655
1656 // restores promise of logical constness
1657 if (op.destructive()) {
1658 for (auto& ff : ncf) ff.clear(false);
1659 world.gop.fence();
1660 } else {
1661 reconstruct(world,f);
1662 }
1663
1664 // svd-tensor requires some cleanup after apply
1665 if (result[0].get_impl()->get_tensor_type()==TT_2D) {
1666 for (auto& r : result) r.get_impl()->finalize_apply();
1667 }
1668
1669 if (print_timings) {
1670 for (auto& r : result) r.get_impl()->print_timer();
1671 op.print_timer();
1672 }
1673 reconstruct(world, result);
1674
1675 return result;
1676 }
1677
1678 /// Normalizes a vector of functions --- v[i] = v[i].scale(1.0/v[i].norm2())
1679 template <typename T, std::size_t NDIM>
1680 void normalize(World& world, std::vector< Function<T,NDIM> >& v, bool fence=true) {
1682 std::vector<double> nn = norm2s(world, v);
1683 for (unsigned int i=0; i<v.size(); ++i) v[i].scale(1.0/nn[i],false);
1684 if (fence) world.gop.fence();
1685 }
1686
1687 template <typename T, std::size_t NDIM>
1688 void print_size(World &world, const std::vector<Function<T,NDIM> > &v, const std::string &msg = "vectorfunction" ){
1689 if(v.empty()){
1690 if(world.rank()==0) std::cout << "print_size: " << msg << " is empty" << std::endl;
1691 }else if(v.size()==1){
1692 v.front().print_size(msg);
1693 }else{
1694 for(auto x:v){
1695 x.print_size(msg);
1696 }
1697 }
1698 }
1699
1700 /// return the size of a vector of functions for each rank
1701 template <typename T, std::size_t NDIM>
1702 double get_size_local(World& world, const std::vector< Function<T,NDIM> >& v){
1703 double size=0.0;
1704 for(auto x:v){
1705 if (x.is_initialized()) size+=x.size_local();
1706 }
1707 const double d=sizeof(T);
1708 const double fac=1024*1024*1024;
1709 return size/fac*d;
1710 }
1711
1712 /// return the size of a function for each rank
1713 template <typename T, std::size_t NDIM>
1715 return get_size_local(f.world(),std::vector<Function<T,NDIM> >(1,f));
1716 }
1717
1718
1719 // gives back the size in GB
1720 template <typename T, std::size_t NDIM>
1721 double get_size(World& world, const std::vector< Function<T,NDIM> >& v){
1722
1723 if (v.empty()) return 0.0;
1724
1725 const double d=sizeof(T);
1726 const double fac=1024*1024*1024;
1727
1728 double size=0.0;
1729 for(unsigned int i=0;i<v.size();i++){
1730 if (v[i].is_initialized()) size+=v[i].size();
1731 }
1732
1733 return size/fac*d;
1734
1735 }
1736
1737 // gives back the size in GB
1738 template <typename T, std::size_t NDIM>
1739 double get_size(const Function<T,NDIM> & f){
1740 const double d=sizeof(T);
1741 const double fac=1024*1024*1024;
1742 double size=f.size();
1743 return size/fac*d;
1744 }
1745
1746 /// apply op on the input vector yielding an output vector of functions
1747
1748 /// @param[in] op the operator working on vin
1749 /// @param[in] vin vector of input Functions; needs to be refined to common level!
1750 /// @return vector of output Functions vout = op(vin)
1751 template <typename T, typename opT, std::size_t NDIM>
1752 std::vector<Function<T,NDIM> > multi_to_multi_op_values(const opT& op,
1753 const std::vector< Function<T,NDIM> >& vin,
1754 const bool fence=true) {
1755 MADNESS_ASSERT(vin.size()>0);
1756 MADNESS_ASSERT(vin[0].is_initialized()); // might be changed
1757 World& world=vin[0].world();
1759 dummy.set_impl(vin[0], false);
1760 std::vector<Function<T,NDIM> > vout=zero_functions<T,NDIM>(world, op.get_result_size());
1761 for (auto& out : vout) out.set_impl(vin[0],false);
1762 dummy.multi_to_multi_op_values(op, vin, vout, fence);
1763 return vout;
1764 }
1765
1766
1767
1768
1769 // convenience operators
1770
1771 /// result[i] = a[i] + b[i]
1772 template <typename T, std::size_t NDIM>
1773 std::vector<Function<T,NDIM> > operator+(const std::vector<Function<T,NDIM> >& lhs,
1774 const std::vector<Function<T,NDIM>>& rhs) {
1775 MADNESS_CHECK(lhs.size() == rhs.size());
1776 return gaxpy_oop(1.0,lhs,1.0,rhs);
1777 }
1778
1779 /// result[i] = a[i] - b[i]
1780 template <typename T, std::size_t NDIM>
1781 std::vector<Function<T,NDIM> > operator-(const std::vector<Function<T,NDIM> >& lhs,
1782 const std::vector<Function<T,NDIM> >& rhs) {
1783 MADNESS_CHECK(lhs.size() == rhs.size());
1784 return gaxpy_oop(1.0,lhs,-1.0,rhs);
1785 }
1786
1787 /// result[i] = a[i] + b
1788 template <typename T, std::size_t NDIM>
1789 std::vector<Function<T,NDIM> > operator+(const std::vector<Function<T,NDIM> >& lhs,
1790 const Function<T,NDIM>& rhs) {
1791 // MADNESS_CHECK(lhs.size() == rhs.size()); // no!!
1792 return gaxpy_oop(1.0,lhs,1.0,rhs);
1793 }
1794
1795 /// result[i] = a[i] - b
1796 template <typename T, std::size_t NDIM>
1797 std::vector<Function<T,NDIM> > operator-(const std::vector<Function<T,NDIM> >& lhs,
1798 const Function<T,NDIM>& rhs) {
1799 // MADNESS_CHECK(lhs.size() == rhs.size()); // no
1800 return gaxpy_oop(1.0,lhs,-1.0,rhs);
1801 }
1802
1803 /// result[i] = a + b[i]
1804 template <typename T, std::size_t NDIM>
1805 std::vector<Function<T,NDIM> > operator+(const Function<T,NDIM>& lhs,
1806 const std::vector<Function<T,NDIM> >& rhs) {
1807 // MADNESS_CHECK(lhs.size() == rhs.size()); // no
1808 return gaxpy_oop(1.0,rhs,1.0,lhs);
1809 }
1810
1811 /// result[i] = a - b[i]
1812 template <typename T, std::size_t NDIM>
1813 std::vector<Function<T,NDIM> > operator-(const Function<T,NDIM>& lhs,
1814 const std::vector<Function<T,NDIM> >& rhs) {
1815// MADNESS_CHECK(lhs.size() == rhs.size()); // no
1816 return gaxpy_oop(-1.0,rhs,1.0,lhs);
1817 }
1818
1819
1820 template <typename T, typename R, std::size_t NDIM>
1821 std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > operator*(const R fac,
1822 const std::vector<Function<T,NDIM> >& rhs) {
1823 if (rhs.size()>0) {
1824 std::vector<Function<T,NDIM> > tmp=copy(rhs[0].world(),rhs);
1825 scale(tmp[0].world(),tmp,TENSOR_RESULT_TYPE(T,R)(fac));
1826 return tmp;
1827 }
1828 return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1829 }
1830
1831 template <typename T, typename R, std::size_t NDIM>
1832 std::vector<Function<T,NDIM> > operator*(const std::vector<Function<T,NDIM> >& rhs,
1833 const R fac) {
1834 if (rhs.size()>0) {
1835 std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > tmp=copy(rhs[0].world(),rhs);
1836 scale(tmp[0].world(),tmp,TENSOR_RESULT_TYPE(T,R)(fac));
1837 return tmp;
1838 }
1839 return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1840 }
1841
1842 /// multiply a vector of functions with a function: r[i] = v[i] * a
1843 template <typename T, typename R, std::size_t NDIM>
1845 const std::vector<Function<R,NDIM> >& v) {
1846 if (v.size()>0) return mul(v[0].world(),a,v,true);
1847 return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1848 }
1849
1850
1851 /// multiply a vector of functions with a function: r[i] = a * v[i]
1852 template <typename T, typename R, std::size_t NDIM>
1853 std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > operator*(const std::vector<Function<T,NDIM> >& v,
1854 const Function<R,NDIM>& a) {
1855 if (v.size()>0) return mul(v[0].world(),a,v,true);
1856 return std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> >();
1857 }
1858
1859
1860 template <typename T, std::size_t NDIM>
1861 std::vector<Function<T,NDIM> > operator+=(std::vector<Function<T,NDIM> >& lhs, const std::vector<Function<T,NDIM> >& rhs) {
1862 MADNESS_CHECK(lhs.size() == rhs.size());
1863 if (lhs.size() > 0) gaxpy(lhs.front().world(), 1.0, lhs, 1.0, rhs);
1864 return lhs;
1865 }
1866
1867 template <typename T, std::size_t NDIM>
1868 std::vector<Function<T,NDIM> > operator-=(std::vector<Function<T,NDIM> >& lhs,
1869 const std::vector<Function<T,NDIM> >& rhs) {
1870 MADNESS_CHECK(lhs.size() == rhs.size());
1871 if (lhs.size() > 0) gaxpy(lhs.front().world(), 1.0, lhs, -1.0, rhs);
1872 return lhs;
1873 }
1874
1875 /// return the real parts of the vector's function (if complex)
1876 template <typename T, std::size_t NDIM>
1877 std::vector<Function<typename Tensor<T>::scalar_type,NDIM> >
1878 real(const std::vector<Function<T,NDIM> >& v, bool fence=true) {
1879 std::vector<Function<typename Tensor<T>::scalar_type,NDIM> > result(v.size());
1880 for (std::size_t i=0; i<v.size(); ++i) result[i]=real(v[i],false);
1881 if (fence and result.size()>0) result[0].world().gop.fence();
1882 return result;
1883 }
1884
1885 /// return the imaginary parts of the vector's function (if complex)
1886 template <typename T, std::size_t NDIM>
1887 std::vector<Function<typename Tensor<T>::scalar_type,NDIM> >
1888 imag(const std::vector<Function<T,NDIM> >& v, bool fence=true) {
1889 std::vector<Function<typename Tensor<T>::scalar_type,NDIM> > result(v.size());
1890 for (std::size_t i=0; i<v.size(); ++i) result[i]=imag(v[i],false);
1891 if (fence and result.size()>0) result[0].world().gop.fence();
1892 return result;
1893 }
1894
1895 /// shorthand gradient operator
1896
1897 /// returns the differentiated function f in all NDIM directions
1898 /// @param[in] f the function on which the grad operator works on
1899 /// @param[in] refine refinement before diff'ing makes the result more accurate
1900 /// @param[in] fence fence after completion; if reconstruction is needed always fence
1901 /// @return the vector \frac{\partial}{\partial x_i} f
1902 template <typename T, std::size_t NDIM>
1903 std::vector<Function<T,NDIM> > grad(const Function<T,NDIM>& f,
1904 bool refine=false, bool fence=true) {
1905
1906 World& world=f.world();
1907 f.reconstruct();
1908 if (refine) f.refine(); // refine to make result more precise
1909
1910 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1912
1913 std::vector<Function<T,NDIM> > result(NDIM);
1914 for (size_t i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1915 if (fence) world.gop.fence();
1916 return result;
1917 }
1918
1919 // BLM first derivative
1920 template <typename T, std::size_t NDIM>
1921 std::vector<Function<T,NDIM> > grad_ble_one(const Function<T,NDIM>& f,
1922 bool refine=false, bool fence=true) {
1923
1924 World& world=f.world();
1925 f.reconstruct();
1926 if (refine) f.refine(); // refine to make result more precise
1927
1928 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1930
1931 // Read in new coeff for each operator
1932 for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_ble1();
1933
1934 std::vector<Function<T,NDIM> > result(NDIM);
1935 for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1936 if (fence) world.gop.fence();
1937 return result;
1938 }
1939
1940 // BLM second derivative
1941 template <typename T, std::size_t NDIM>
1942 std::vector<Function<T,NDIM> > grad_ble_two(const Function<T,NDIM>& f,
1943 bool refine=false, bool fence=true) {
1944
1945 World& world=f.world();
1946 f.reconstruct();
1947 if (refine) f.refine(); // refine to make result more precise
1948
1949 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1951
1952 // Read in new coeff for each operator
1953 for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_ble2();
1954
1955 std::vector<Function<T,NDIM> > result(NDIM);
1956 for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1957 if (fence) world.gop.fence();
1958 return result;
1959 }
1960
1961 // Bspline first derivative
1962 template <typename T, std::size_t NDIM>
1963 std::vector<Function<T,NDIM> > grad_bspline_one(const Function<T,NDIM>& f,
1964 bool refine=false, bool fence=true) {
1965
1966 World& world=f.world();
1967 f.reconstruct();
1968 if (refine) f.refine(); // refine to make result more precise
1969
1970 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1972
1973 // Read in new coeff for each operator
1974 for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline1();
1975
1976 std::vector<Function<T,NDIM> > result(NDIM);
1977 for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1978 if (fence) world.gop.fence();
1979 return result;
1980 }
1981
1982 // Bpsline second derivative
1983 template <typename T, std::size_t NDIM>
1984 std::vector<Function<T,NDIM> > grad_bpsline_two(const Function<T,NDIM>& f,
1985 bool refine=false, bool fence=true) {
1986
1987 World& world=f.world();
1988 f.reconstruct();
1989 if (refine) f.refine(); // refine to make result more precise
1990
1991 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
1993
1994 // Read in new coeff for each operator
1995 for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline2();
1996
1997 std::vector<Function<T,NDIM> > result(NDIM);
1998 for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
1999 if (fence) world.gop.fence();
2000 return result;
2001 }
2002
2003 // Bspline third derivative
2004 template <typename T, std::size_t NDIM>
2005 std::vector<Function<T,NDIM> > grad_bspline_three(const Function<T,NDIM>& f,
2006 bool refine=false, bool fence=true) {
2007
2008 World& world=f.world();
2009 f.reconstruct();
2010 if (refine) f.refine(); // refine to make result more precise
2011
2012 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
2014
2015 // Read in new coeff for each operator
2016 for (unsigned int i=0; i<NDIM; ++i) (*grad[i]).set_bspline3();
2017
2018 std::vector<Function<T,NDIM> > result(NDIM);
2019 for (unsigned int i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),f,false);
2020 if (fence) world.gop.fence();
2021 return result;
2022 }
2023
2024
2025
2026 /// shorthand div operator
2027
2028 /// returns the dot product of nabla with a vector f
2029 /// @param[in] f the vector of functions on which the div operator works on
2030 /// @param[in] refine refinement before diff'ing makes the result more accurate
2031 /// @param[in] fence fence after completion; currently always fences
2032 /// @return the vector \frac{\partial}{\partial x_i} f
2033 /// TODO: add this to operator fusion
2034 template <typename T, std::size_t NDIM>
2036 bool do_refine=false, bool fence=true) {
2037
2038 MADNESS_ASSERT(v.size()>0);
2039 World& world=v[0].world();
2040 reconstruct(world,v);
2041 if (do_refine) refine(world,v); // refine to make result more precise
2042
2043 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
2045
2046 std::vector<Function<T,NDIM> > result(NDIM);
2047 for (size_t i=0; i<NDIM; ++i) result[i]=apply(*(grad[i]),v[i],false);
2048 world.gop.fence();
2049 return sum(world,result,fence);
2050 }
2051
2052 /// shorthand rot operator
2053
2054 /// returns the cross product of nabla with a vector f
2055 /// @param[in] f the vector of functions on which the rot operator works on
2056 /// @param[in] refine refinement before diff'ing makes the result more accurate
2057 /// @param[in] fence fence after completion; currently always fences
2058 /// @return the vector \frac{\partial}{\partial x_i} f
2059 /// TODO: add this to operator fusion
2060 template <typename T, std::size_t NDIM>
2061 std::vector<Function<T,NDIM> > rot(const std::vector<Function<T,NDIM> >& v,
2062 bool do_refine=false, bool fence=true) {
2063
2064 MADNESS_ASSERT(v.size()==3);
2065 World& world=v[0].world();
2066 reconstruct(world,v);
2067 if (do_refine) refine(world,v); // refine to make result more precise
2068
2069 std::vector< std::shared_ptr< Derivative<T,NDIM> > > grad=
2071
2072 std::vector<Function<T,NDIM> > d(NDIM),dd(NDIM);
2073 d[0]=apply(*(grad[1]),v[2],false); // Dy z
2074 d[1]=apply(*(grad[2]),v[0],false); // Dz x
2075 d[2]=apply(*(grad[0]),v[1],false); // Dx y
2076 dd[0]=apply(*(grad[2]),v[1],false); // Dz y
2077 dd[1]=apply(*(grad[0]),v[2],false); // Dx z
2078 dd[2]=apply(*(grad[1]),v[0],false); // Dy x
2079 world.gop.fence();
2080
2081 d[0].gaxpy(1.0,dd[0],-1.0,false);
2082 d[1].gaxpy(1.0,dd[1],-1.0,false);
2083 d[2].gaxpy(1.0,dd[2],-1.0,false);
2084
2085 world.gop.fence();
2086 return d;
2087 }
2088
2089 /// shorthand cross operator
2090
2091 /// returns the cross product of vectors f and g
2092 /// @param[in] f the vector of functions on which the rot operator works on
2093 /// @param[in] g the vector of functions on which the rot operator works on
2094 /// @param[in] fence fence after completion; currently always fences
2095 /// @return the vector \frac{\partial}{\partial x_i} f
2096 /// TODO: add this to operator fusion
2097 template <typename T, typename R, std::size_t NDIM>
2098 std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > cross(const std::vector<Function<T,NDIM> >& f,
2099 const std::vector<Function<R,NDIM> >& g,
2100 bool do_refine=false, bool fence=true) {
2101
2102 MADNESS_ASSERT(f.size()==3);
2103 MADNESS_ASSERT(g.size()==3);
2104 World& world=f[0].world();
2105 reconstruct(world,f,false);
2106 reconstruct(world,g);
2107
2108 std::vector<Function<TENSOR_RESULT_TYPE(T,R),NDIM> > d(f.size()),dd(f.size());
2109
2110 d[0]=mul(f[1],g[2],false);
2111 d[1]=mul(f[2],g[0],false);
2112 d[2]=mul(f[0],g[1],false);
2113
2114 dd[0]=mul(f[2],g[1],false);
2115 dd[1]=mul(f[0],g[2],false);
2116 dd[2]=mul(f[1],g[0],false);
2117 world.gop.fence();
2118
2119 compress(world,d,false);
2120 compress(world,dd);
2121
2122 d[0].gaxpy(1.0,dd[0],-1.0,false);
2123 d[1].gaxpy(1.0,dd[1],-1.0,false);
2124 d[2].gaxpy(1.0,dd[2],-1.0,false);
2125
2126 world.gop.fence();
2127 return d;
2128 }
2129
2130 template<typename T, std::size_t NDIM>
2131 void load_balance(World& world, std::vector<Function<T,NDIM> >& vf) {
2132
2133 struct LBCost {
2134 LBCost() = default;
2135 double operator()(const Key<NDIM>& key, const FunctionNode<T,NDIM>& node) const {
2136 return node.coeff().size();
2137 }
2138 };
2139
2140 LoadBalanceDeux<6> lb(world);
2141 for (const auto& f : vf) lb.add_tree(f, LBCost());
2143
2144 }
2145
2146 /// load a vector of functions
2147 template<typename T, size_t NDIM>
2148 void load_function(World& world, std::vector<Function<T,NDIM> >& f,
2149 const std::string name) {
2150 if (world.rank()==0) print("loading vector of functions",name);
2152 std::size_t fsize=0;
2153 ar & fsize;
2154 f.resize(fsize);
2155 for (std::size_t i=0; i<fsize; ++i) ar & f[i];
2156 }
2157
2158 /// save a vector of functions
2159 template<typename T, size_t NDIM>
2160 void save_function(const std::vector<Function<T,NDIM> >& f, const std::string name) {
2161 if (f.size()>0) {
2162 World& world=f.front().world();
2163 if (world.rank()==0) print("saving vector of functions",name);
2165 std::size_t fsize=f.size();
2166 ar & fsize;
2167 for (std::size_t i=0; i<fsize; ++i) ar & f[i];
2168 }
2169 }
2170
2171
2172}
2173#endif // MADNESS_MRA_VMRA_H__INCLUDED
double q(double t)
Definition DKops.h:18
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
static void redistribute(World &world, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &newpmap)
Sets the default process map and redistributes all functions using the old map.
Definition funcdefaults.h:414
FunctionFactory implements the named-parameter idiom for Function.
Definition function_factory.h:86
FunctionImpl holds all Function state to facilitate shallow copy semantics.
Definition funcimpl.h:945
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:5924
void undo_redundant(const bool fence)
convert this from redundant to standard reconstructed form
Definition mraimpl.h:1560
void multiply(const implT *f, const FunctionImpl< T, LDIM > *g, const int particle)
multiply f (a pair function of NDIM) with an orbital g (LDIM=NDIM/2)
Definition funcimpl.h:3560
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
static Tensor< TENSOR_RESULT_TYPE(T, R)> dot_local(const std::vector< const FunctionImpl< T, NDIM > * > &left, const std::vector< const FunctionImpl< R, NDIM > * > &right, bool sym)
Definition funcimpl.h:5976
FunctionNode holds the coefficients, etc., at each node of the 2^NDIM-tree.
Definition funcimpl.h:127
A multiresolution adaptive numerical function.
Definition mra.h:139
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:1006
void set_impl(const std::shared_ptr< FunctionImpl< T, NDIM > > &impl)
Replace current FunctionImpl with provided new one.
Definition mra.h:646
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:68
Definition lbdeux.h:228
std::shared_ptr< WorldDCPmapInterface< keyT > > load_balance(double fac=1.0, bool printstuff=false)
Actually does the partitioning of the tree.
Definition lbdeux.h:320
void add_tree(const Function< T, NDIM > &f, const costT &costfn, bool fence=false)
Accumulates cost from a function.
Definition lbdeux.h:289
Convolutions in separated form (including Gaussian)
Definition operator.h:139
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
A tensor is a multidimensional array.
Definition tensor.h:317
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition tensor.h:409
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1824
A simple, fixed dimension vector.
Definition vector.h:64
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:206
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition world.h:320
WorldGopInterface & gop
Global operations.
Definition world.h:207
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
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 rot(x, k)
Definition lookup3.c:72
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
#define MADNESS_CHECK_THROW(condition, msg)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:207
Main include file for MADNESS and defines Function interface.
static const bool VERIFY_TREE
Definition mra.h:57
Namespace for all elements and tools of MADNESS.
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:2160
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
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< 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:2745
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:1949
Function< typename TensorTypeData< Q >::scalar_type, NDIM > abs_square(const Function< Q, NDIM > &func)
Definition complexfun.h:121
const std::vector< Function< T, NDIM > > & change_tree_state(const std::vector< Function< T, NDIM > > &v, const TreeState finalstate, const bool fence=true)
change tree state of the functions
Definition vmra.h:277
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition mra.h:2713
response_space scale(response_space a, double b)
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:2003
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< 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
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
std::vector< Function< T, NDIM > > grad_bspline_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:1963
void set_impl(std::vector< Function< T, NDIM > > &v, const std::vector< std::shared_ptr< FunctionImpl< T, NDIM > > > vimpl)
Definition vmra.h:670
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:2064
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:1733
void truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
std::vector< std::shared_ptr< FunctionImpl< T, NDIM > > > get_impl(const std::vector< Function< T, NDIM > > &v)
Definition vmra.h:663
Function< T, NDIM > div(const std::vector< Function< T, NDIM > > &v, bool do_refine=false, bool fence=true)
shorthand div operator
Definition vmra.h:2035
std::vector< Function< T, NDIM > > orthonormalize_cd(const std::vector< Function< T, NDIM > > &v, Tensor< T > &ovlp)
Definition vmra.h:557
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
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition vmra.h:689
std::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:2098
TreeState
Definition funcdefaults.h:59
@ nonstandard_after_apply
s and d coeffs, state after operator application
Definition funcdefaults.h:64
@ redundant_after_merge
s coeffs everywhere, must be summed up to yield the result
Definition funcdefaults.h:66
@ reconstructed
s coeffs at the leaves only
Definition funcdefaults.h:60
@ nonstandard
s and d coeffs in internal nodes
Definition funcdefaults.h:62
@ compressed
d coeffs in internal nodes, s and d coeffs at the root
Definition funcdefaults.h:61
@ redundant
s coeffs everywhere
Definition funcdefaults.h:65
@ nonstandard_with_leaves
like nonstandard, with s coeffs at the leaves
Definition funcdefaults.h:63
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:2078
void cholesky(Tensor< T > &A)
Compute the Cholesky factorization.
Definition lapack.cc:1174
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
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
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition vmra.h:133
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
response_space transpose(response_space &f)
Definition basic_operators.cc:10
std::vector< Function< T, NDIM > > impl2function(const std::vector< std::shared_ptr< FunctionImpl< T, NDIM > > > vimpl)
Definition vmra.h:676
std::vector< Function< T, NDIM > > grad_bpsline_two(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:1984
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
std::vector< Function< T, NDIM > > flatten(const std::vector< std::vector< Function< T, NDIM > > > &vv)
Definition vmra.h:656
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< CCPairFunction< T, NDIM > > operator*(const double fac, const std::vector< CCPairFunction< T, NDIM > > &arg)
Definition ccpairfunction.h:1089
static void verify_tree(World &world, const std::vector< Function< T, NDIM > > &v)
Definition SCF.cc:72
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:1752
static const Slice _(0,-1, 1)
Tensor< T > inverse(const Tensor< T > &a_in)
invert general square matrix A
Definition lapack.cc:832
void load_balance(const real_function_6d &f, const bool leaf)
do some load-balancing
Definition madness/chem/mp2.cc:70
std::vector< CCPairFunction< T, NDIM > > operator-(const std::vector< CCPairFunction< T, NDIM > > c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition ccpairfunction.h:1060
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:1774
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
Function< T, NDIM > gaxpy_oop_reconstructed(const double alpha, const Function< T, NDIM > &left, const double beta, const Function< T, NDIM > &right, const bool fence=true)
Returns new function alpha*left + beta*right optional fence, having both addends reconstructed.
Definition mra.h:1967
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
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition print.h:225
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
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
@ 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:1688
NDIM & f
Definition mra.h:2448
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:1959
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:1493
std::vector< CCPairFunction< T, NDIM > > & operator-=(std::vector< CCPairFunction< T, NDIM > > &rhs, const std::vector< CCPairFunction< T, NDIM > > &lhs)
Definition ccpairfunction.h:1082
std::vector< Function< T, NDIM > > orthonormalize(const std::vector< Function< T, NDIM > > &vf_in)
orthonormalize the vectors
Definition vmra.h:421
NDIM const Function< R, NDIM > & g
Definition mra.h:2448
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition timers.cc:48
std::vector< Function< T, NDIM > > grad_ble_one(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:1921
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:2152
std::vector< Function< T, NDIM > > grad_bspline_three(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:2005
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
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
double imag(double x)
Definition complexfun.h:56
void load_function(World &world, std::vector< Function< T, NDIM > > &f, const std::string name)
load a vector of functions
Definition vmra.h:2148
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:1839
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
std::vector< Function< T, NDIM > > grad_ble_two(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
Definition vmra.h:1942
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:1680
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
std::vector< CCPairFunction< T, NDIM > > operator+(const std::vector< CCPairFunction< T, NDIM > > c1, const std::vector< CCPairFunction< T, NDIM > > &c2)
Definition ccpairfunction.h:1052
std::vector< Function< T, NDIM > > grad(const Function< T, NDIM > &f, bool refine=false, bool fence=true)
shorthand gradient operator
Definition vmra.h:1903
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:2404
DistributedMatrix< T > matrix_dot(const DistributedMatrixDistribution &d, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition vmra.h:1433
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
double real(double x)
Definition complexfun.h:52
std::vector< CCPairFunction< T, NDIM > > & operator+=(std::vector< CCPairFunction< T, NDIM > > &lhs, const CCPairFunction< T, NDIM > &rhs)
Definition ccpairfunction.h:1068
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
std::string name(const FuncType &type, const int ex=-1)
Definition ccpairfunction.h:28
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:1721
double get_size_local(World &world, const std::vector< Function< T, NDIM > > &v)
return the size of a vector of functions for each rank
Definition vmra.h:1702
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:2034
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
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
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
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:246
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 d
Definition nonlinschro.cc:121
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 m
Definition relops.cc:9
static const double L
Definition rk.cc:46
static const double thresh
Definition rk.cc:45
Definition test_ar.cc:204
Definition mp2.h:63
double operator()(const Key< 6 > &key, const FunctionNode< double, 6 > &node) const
Definition mp2.h:70
Definition lowrankfunction.h:332
Definition dirac-hatom.cc:108
AtomicInt sum
Definition test_atomicint.cc:46
int P
Definition test_binsorter.cc:9
double aa
Definition testbsh.cc:68
static const double alpha
Definition testcosine.cc:10
constexpr std::size_t NDIM
Definition testgconv.cc:54
#define TENSOR_RESULT_TYPE(L, R)
This macro simplifies access to TensorResultType.
Definition type_data.h:205
#define PROFILE_FUNC
Definition worldprofile.h:209
#define PROFILE_BLOCK(name)
Definition worldprofile.h:208