MADNESS  0.10.1
vmra1.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 
34 #ifndef MADNESS_MRA_VMRA_H__INCLUDED
35 #define MADNESS_MRA_VMRA_H__INCLUDED
36 
37 /*!
38  \file vmra.h
39  \brief Defines operations on vectors of Functions
40  \ingroup mra
41 
42  This file defines a number of operations on vectors of functions.
43  Assume v is a vector of NDIM-D functions of a certain type.
44 
45 
46  Operations on array of functions
47 
48  *) copying: deep copying of vectors of functions to vector of functions
49  \code
50  vector2 = copy(world, vector1,fence);
51  \endcode
52 
53  *) compress: convert multiwavelet representation to legendre representation
54  \code
55  compress(world, vector, fence);
56  \endcode
57 
58  *) reconstruct: convert representation to multiwavelets
59  \code
60  reconstruct(world, vector, fence);
61  \endcode
62 
63  *) make_nonstandard: convert to non-standard form
64  \code
65  make_nonstandard(world, v, fence);
66  \endcode
67 
68  *) standard: convert to standard form
69  \code
70  standard(world, v, fence);
71  \endcode
72 
73  *) truncate: truncating vectors of functions to desired precision
74  \code
75  truncate(world, v, tolerance, fence);
76  \endcode
77 
78 
79  *) zero function: create a vector of zero functions of length n
80  \code
81  v=zero(world, n);
82  \endcode
83 
84  *) transform: transform a representation from one basis to another
85  \code
86  transform(world, vector, tensor, tolerance, fence )
87  \endcode
88 
89  Setting thresh-hold for precision
90 
91  *) set_thresh: setting a finite thresh-hold for a vector of functions
92  \code
93  void set_thresh(World& world, std::vector< Function<T,NDIM> >& v, double thresh, bool fence=true);
94  \endcode
95 
96  Arithmetic Operations on arrays of functions
97 
98  *) conjugation: conjugate a vector of complex functions
99 
100  *) add
101  *) sub
102  *) mul
103  - mul_sparse
104  *) square
105  *) gaxpy
106  *) apply
107 
108  Norms, inner-products, blas-1 like operations on vectors of functions
109 
110  *) inner
111  *) matrix_inner
112  *) norm_tree
113  *) normalize
114  *) norm2
115  - norm2s
116  *) scale(world, v, alpha);
117 
118 
119 
120 
121 */
122 
123 #include <madness/mra/mra.h>
124 #include <madness/mra/derivative.h>
125 #include <cstdio>
126 
127 namespace madness {
128 
129  /// Compress a vector of functions
130  template <typename T, std::size_t NDIM>
131  void compress(World& world,
132  const std::vector< Function<T,NDIM> >& v,
133  unsigned int blk=1,
134  bool fence=true){
135 
136  PROFILE_BLOCK(Vcompress);
137  bool must_fence = false;
138  unsigned int vvsize = v.size();
139  for (unsigned int i=0; i<vvsize; i+= blk) {
140  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
141  if (!v[j].is_compressed()) {
142  v[j].compress(false);
143  must_fence = true;
144  }
145  }
146  if ( blk!=1 && must_fence && fence) world.gop.fence();
147  }
148 
149  if (fence && must_fence) world.gop.fence();
150  }
151 
152 
153  /// Reconstruct a vector of functions
154  template <typename T, std::size_t NDIM>
155  void reconstruct(World& world,
156  const std::vector< Function<T,NDIM> >& v,
157  unsigned int blk=1,
158  bool fence=true){ // reconstr
159  PROFILE_BLOCK(Vreconstruct);
160  bool must_fence = false;
161  unsigned int vvsize = v.size();
162  for (unsigned int i=0; i<vvsize; i+= blk) {
163  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
164  if (v[j].is_compressed()) {
165  v[j].reconstruct(false);
166  must_fence = true;
167  }
168  }
169  if ( blk!=1 && must_fence && fence) world.gop.fence();
170  }
171 
172  if (fence && must_fence) world.gop.fence();
173  } // reconstr
174 
175  /// Generates non-standard form of a vector of functions
176  template <typename T, std::size_t NDIM>
177  void nonstandard(World& world,
178  std::vector< Function<T,NDIM> >& v,
179  unsigned int blk=1,
180  bool fence=true) { // nonstand
181  PROFILE_BLOCK(Vnonstandard);
182  unsigned int vvsize = v.size();
183  reconstruct(world, v, blk);
184  for (unsigned int i=0; i<vvsize; i+= blk) {
185  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
186  v[j].make_nonstandard(false, false);
187  }
188  if ( blk!=1 && fence) world.gop.fence();
189  }
190  if (fence) world.gop.fence();
191  } //nonstand
192 
193 
194  /// Generates standard form of a vector of functions
195 
196  template <typename T, std::size_t NDIM>
197  void standard(World& world,
198  std::vector< Function<T,NDIM> >& v,
199  unsigned int blk=1,
200  bool fence=true){ // standard
201  PROFILE_BLOCK(Vstandard);
202  unsigned int vvsize = v.size();
203  for (unsigned int i=0; i<vvsize; i+= blk) {
204  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
205  v[j].standard(false);
206  }
207  if ( blk!=1 && fence) world.gop.fence();
208  }
209  if (fence) world.gop.fence();
210  } // standard
211 
212 
213  /// Truncates a vector of functions
214  template <typename T, std::size_t NDIM>
215  void truncate(World& world,
216  std::vector< Function<T,NDIM> >& v,
217  double tol=0.0,
218  unsigned int blk=1,
219  bool fence=true){ // truncate
220  PROFILE_BLOCK(Vtruncate);
221 
222  compress(world, v, blk);
223 
224  unsigned int vvsize = v.size();
225  for (unsigned int i=0; i<vvsize; i+= blk) {
226  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
227  v[j].truncate(tol, false);
228  }
229  if ( blk!=1 && fence) world.gop.fence();
230  }
231  if (fence) world.gop.fence();
232  } //truncate
233 
234  /// Applies a derivative operator to a vector of functions
235 
236  template <typename T, std::size_t NDIM>
237  std::vector< Function<T,NDIM> >
238  apply(World& world,
239  const Derivative<T,NDIM>& D,
240  const std::vector< Function<T,NDIM> >& v,
241  const unsigned int blk=1,
242  const bool fence=true)
243  {
244  reconstruct(world, v, blk);
245  std::vector< Function<T,NDIM> > df(v.size());
246 
247  unsigned int vvsize = v.size();
248 
249  for (unsigned int i=0; i<vvsize; i+= blk) {
250  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
251  df[j] = D(v[j],false);
252  }
253  if (blk!= 1 && fence) world.gop.fence();
254  }
255  if (fence) world.gop.fence();
256 
257  return df;
258  }
259 
260  /// Generates a vector of zero functions
261 
262  template <typename T, std::size_t NDIM>
263  std::vector< Function<T,NDIM> >
264  zero_functions(World& world, int n) {
265  PROFILE_BLOCK(Vzero_functions);
266  std::vector< Function<T,NDIM> > r(n);
267  for (int i=0; i<n; ++i)
268  r[i] = Function<T,NDIM>(FunctionFactory<T,NDIM>(world));
269 
270  return r;
271  }
272 
273 
274  /// Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i]
275 
276  /// Uses sparsity in the transformation matrix --- set small elements to
277  /// zero to take advantage of this.
278 
279  template <typename T, typename R, std::size_t NDIM>
280  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> >
281  transform(World& world,
282  const std::vector< Function<T,NDIM> >& v,
283  const Tensor<R>& c,
284  unsigned int blki=1,
285  unsigned int blkj=1,
286  const bool fence=true){
287 
288  PROFILE_BLOCK(Vtransformsp);
289 
290  typedef TENSOR_RESULT_TYPE(T,R) resultT;
291 
292  unsigned int blk = min(blki, blkj);
293  unsigned int n = v.size(); // n is the old dimension
294  unsigned int m = c.dim(1); // m is the new dimension
295  MADNESS_ASSERT(n==c.dim(0));
296 
297  std::vector< Function<resultT,NDIM> > vc = zero_functions_compressed<resultT,NDIM>(world, m);
298  compress(world, v, blk);
299 
300  for (unsigned int i=0; i<m; i+= blki) {
301  for (unsigned int ii=i; ii<std::min(m,(i+1)*blki); ii++) {
302  for (unsigned int j=0; j<n; j+= blkj) {
303  for (unsigned int jj=j; jj<std::min(n, (j+1)*blkj); jj++)
304  if (c(jj,ii) != R(0.0)) vc[ii].gaxpy(1.0,v[jj],c(jj,ii),false);
305  if (fence && (blkj!=1)) world.gop.fence();
306  }
307  }
308  if (fence && (blki!=1)) world.gop.fence(); // a bit conservative
309  }
310 
311 
312  // for (unsigned int i=0; i<m; ++i) {
313  // for (unsigned int j=0; j<n; ++j) {
314  // if (c(j,i) != R(0.0)) vc[i].gaxpy(1.0,v[j],c(j,i),false);
315  // }
316  // }
317 
318  if (fence) world.gop.fence();
319  return vc;
320  }
321 
322 
323  template <typename L, typename R, std::size_t NDIM>
324  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> >
325  transform(World& world,
326  const std::vector< Function<L,NDIM> >& v,
327  const Tensor<R>& c,
328  const double tol,
329  const unsigned int blki=1,
330  const bool fence)
331  {
332  PROFILE_BLOCK(Vtransform);
333  MADNESS_ASSERT(v.size() == (unsigned int)(c.dim(0)));
334 
335  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> > vresult(c.dim(1));
336  unsigned int m=c.dim(1);
337 
338  for (unsigned int i=0; i<m; i+= blki) {
339  for (unsigned int ii=i; ii<std::min(m,(i+1)*blki); ii++) {
340  vresult[ii] = Function<TENSOR_RESULT_TYPE(L,R),NDIM>(FunctionFactory<TENSOR_RESULT_TYPE(L,R),NDIM>(world));
341  }
342  if (fence && (blki!=1)) world.gop.fence(); // a bit conservative
343  }
344 
345 
346  // for (unsigned int i=0; i<c.dim(1); ++i) {
347  // vresult[i] = Function<TENSOR_RESULT_TYPE(L,R),NDIM>(FunctionFactory<TENSOR_RESULT_TYPE(L,R),NDIM>(world));
348  // }
349  compress(world, v, blki, false);
350  compress(world, vresult, blki, false);
351  world.gop.fence();
352  vresult[0].vtransform(v, c, vresult, tol, fence);
353  return vresult;
354  }
355 
356  /// Scales inplace a vector of functions by distinct values
357 
358  template <typename T, typename Q, std::size_t NDIM>
359  void scale(World& world,
360  std::vector< Function<T,NDIM> >& v,
361  const std::vector<Q>& factors,
362  const unsigned int blk=1,
363  const bool fence=true)
364  {
365  PROFILE_BLOCK(Vscale);
366 
367  unsigned int vvsize = v.size();
368  for (unsigned int i=0; i<vvsize; i+= blk) {
369  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
370  v[j].scale(factors[j],false);
371  }
372  if (fence && blk!=1 ) world.gop.fence();
373  }
374  if (fence) world.gop.fence();
375  }
376 
377  /// Scales inplace a vector of functions by the same
378 
379  template <typename T, typename Q, std::size_t NDIM>
380  void scale(World& world,
381  std::vector< Function<T,NDIM> >& v,
382  const Q factor,
383  const unsigned int blk=1,
384  const bool fence=true){
385 
386  PROFILE_BLOCK(Vscale); // shouldn't need blocking since it is local
387 
388  unsigned int vvsize = v.size();
389  for (unsigned int i=0; i<vvsize; i+= blk) {
390  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
391  v[j].scale(factor,false);
392  }
393  if (fence && blk!=1 ) world.gop.fence();
394  }
395  if (fence) world.gop.fence();
396  }
397 
398  /// Computes the 2-norms of a vector of functions
399  template <typename T, std::size_t NDIM>
400  std::vector<double> norm2s(World& world,
401  const std::vector< Function<T,NDIM> >& v,
402  const unsigned int blk=1,
403  const bool fence=true){
404 
405  PROFILE_BLOCK(Vnorm2);
406  unsigned int vvsize = v.size();
407  std::vector<double> norms(vvsize);
408 
409  for (unsigned int i=0; i<vvsize; i+= blk) {
410  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
411  norms[j] = v[j].norm2sq_local();
412  }
413  if (fence && (blk!=1)) world.gop.fence();
414  }
415  if (fence ) world.gop.fence();
416 
417  world.gop.sum(&norms[0], norms.size());
418 
419  for (unsigned int i=0; i<vvsize; i+= blk) {
420  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
421  norms[j] = sqrt(norms[j]);
422  if (fence && (blk!=1)) world.gop.fence();
423  }
424 
425  world.gop.fence();
426  return norms;
427  }
428 
429  /// Computes the 2-norm of a vector of functions
430  // should be local; norms[0] contains the result
431 
432  template <typename T, std::size_t NDIM>
433  double norm2(World& world,
434  const std::vector< Function<T,NDIM> >& v)
435  {
436 
437  PROFILE_BLOCK(Vnorm2);
438 
439  std::vector<double> norms(v.size());
440 
441  for (unsigned int i=0; i<v.size(); ++i)
442  norms[i] = v[i].norm2sq_local();
443 
444  world.gop.sum(&norms[0], norms.size());
445 
446  for (unsigned int i=1; i<v.size(); ++i)
447  norms[0] += norms[i];
448 
449  world.gop.fence();
450  return sqrt(norms[0]);
451  }
452 
453  inline double conj(double x) {
454  return x;
455  }
456 
457  inline double conj(float x) {
458  return x;
459  }
460 
461 
462  template <typename T, typename R, std::size_t NDIM>
463  struct MatrixInnerTask : public TaskInterface {
464  Tensor<TENSOR_RESULT_TYPE(T,R)> result; // Must be a copy
465  const Function<T,NDIM>& f;
466  const std::vector< Function<R,NDIM> >& g;
467  long jtop;
468 
469  MatrixInnerTask(const Tensor<TENSOR_RESULT_TYPE(T,R)>& result,
470  const Function<T,NDIM>& f,
471  const std::vector< Function<R,NDIM> >& g,
472  long jtop)
473  : result(result), f(f), g(g), jtop(jtop) {}
474 
475  void run(World& world) {
476  for (long j=0; j<jtop; ++j) {
477  result(j) = f.inner_local(g[j]);
478  }
479  }
480 
481  private:
482  /// Get the task id
483 
484  /// \param id The id to set for this task
485  virtual void get_id(std::pair<void*,unsigned short>& id) const {
486  PoolTaskInterface::make_id(id, *this);
487  }
488  }; // struct MatrixInnerTask
489 
490  /// Computes the matrix inner product of two function vectors - q(i,j) = inner(f[i],g[j])
491 
492  /// For complex types symmetric is interpreted as Hermitian.
493  ///
494  /// The current parallel loop is non-optimal but functional.
495 
496  template <typename T, typename R, std::size_t NDIM>
497  Tensor< TENSOR_RESULT_TYPE(T,R) > matrix_inner(World& world,
498  const std::vector< Function<T,NDIM> >& f,
499  const std::vector< Function<R,NDIM> >& g,
500  bool sym=false) {
501  PROFILE_BLOCK(Vmatrix_inner);
502  unsigned int n=f.size(), m=g.size();
503  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n,m);
504  if (sym) MADNESS_ASSERT(n==m);
505 
506  world.gop.fence();
507  compress(world, f);
508  if (&f != &g) compress(world, g);
509 
510  // for (long i=0; i<n; ++i) {
511  // long jtop = m;
512  // if (sym) jtop = i+1;
513  // for (long j=0; j<jtop; ++j) {
514  // r(i,j) = f[i].inner_local(g[j]);
515  // if (sym) r(j,i) = conj(r(i,j));
516  // }
517  // }
518 
519  for (unsigned int i=n-1; i>=0; --i) {
520  unsigned int jtop = m;
521  if (sym) jtop = i+1;
522  world.taskq.add(new MatrixInnerTask<T,R,NDIM>(r(i,_), f[i], g, jtop));
523  }
524  world.gop.fence();
525  world.gop.sum(r.ptr(),n*m);
526 
527  if (sym) {
528  for (unsigned int i=0; i<n; ++i) {
529  for (unsigned int j=0; j<i; ++j) {
530  r(j,i) = conj(r(i,j));
531  }
532  }
533  }
534  return r;
535  }
536 
537  /// Computes the element-wise inner product of two function vectors - q(i) = inner(f[i],g[i])
538 
539  template <typename T, typename R, std::size_t NDIM>
540  Tensor< TENSOR_RESULT_TYPE(T,R) > inner(World& world,
541  const std::vector< Function<T,NDIM> >& f,
542  const std::vector< Function<R,NDIM> >& g) {
543  PROFILE_BLOCK(Vinnervv);
544  long n=f.size(), m=g.size();
545  MADNESS_ASSERT(n==m);
546  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
547 
548  compress(world, f);
549  compress(world, g);
550 
551  for (long i=0; i<n; ++i) {
552  r(i) = f[i].inner_local(g[i]);
553  }
554 
555  world.taskq.fence();
556  world.gop.sum(r.ptr(),n);
557  world.gop.fence();
558  return r;
559  }
560 
561 
562  /// Computes the inner product of a function with a function vector - q(i) = inner(f,g[i])
563 
564  template <typename T, typename R, std::size_t NDIM>
565  Tensor< TENSOR_RESULT_TYPE(T,R) > inner(World& world,
566  const Function<T,NDIM>& f,
567  const std::vector< Function<R,NDIM> >& g) {
568  PROFILE_BLOCK(Vinner);
569  long n=g.size();
570  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
571 
572  f.compress();
573  compress(world, g);
574 
575  for (long i=0; i<n; ++i) {
576  r(i) = f.inner_local(g[i]);
577  }
578 
579  world.taskq.fence();
580  world.gop.sum(r.ptr(),n);
581  world.gop.fence();
582  return r;
583  }
584 
585 
586  /// Multiplies a function against a vector of functions --- q[i] = a * v[i]
587 
588  template <typename T, typename R, std::size_t NDIM>
589  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
590  mul(World& world,
591  const Function<T,NDIM>& a,
592  const std::vector< Function<R,NDIM> >& v,
593  const unsigned int blk=1,
594  const bool fence=true) {
595 
596  PROFILE_BLOCK(Vmul);
597  a.reconstruct(false);
598  reconstruct(world, v, blk, false);
599  world.gop.fence();
600  return vmulXX(a, v, 0.0, fence);
601  }
602 
603  /// Multiplies a function against a vector of functions using sparsity of a and v[i] --- q[i] = a * v[i]
604  template <typename T, typename R, std::size_t NDIM>
605  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
606  mul_sparse(World& world,
607  const Function<T,NDIM>& a,
608  const std::vector< Function<R,NDIM> >& v,
609  const double tol,
610  const bool fence=true,
611  const unsigned int blk=1)
612  {
613  PROFILE_BLOCK(Vmulsp);
614  a.reconstruct(false);
615  reconstruct(world, v, blk, false);
616  world.gop.fence();
617 
618  unsigned int vvsize = v.size();
619  for (unsigned int i=0; i<vvsize; i+= blk) {
620  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
621  v[j].norm_tree(false);
622  if ( fence && (blk == 1)) world.gop.fence();
623  }
624  a.norm_tree();
625  return vmulXX(a, v, tol, fence);
626  }
627 
628  /// Makes the norm tree for all functions in a vector
629  template <typename T, std::size_t NDIM>
630  void norm_tree(World& world,
631  const std::vector< Function<T,NDIM> >& v,
632  bool fence=true,
633  unsigned int blk=1){
634  PROFILE_BLOCK(Vnorm_tree);
635 
636  unsigned int vvsize = v.size();
637  for (unsigned int i=0; i<vvsize; i+= blk) {
638  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
639  v[j].norm_tree(false);
640  if (fence && blk!=1 ) world.gop.fence();
641  }
642  if (fence) world.gop.fence();
643  }
644 
645  /// Multiplies two vectors of functions q[i] = a[i] * b[i]
646 
647  template <typename T, typename R, std::size_t NDIM>
648  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
649  mul(World& world,
650  const std::vector< Function<T,NDIM> >& a,
651  const std::vector< Function<R,NDIM> >& b,
652  bool fence=true,
653  unsigned int blk=1){
654  PROFILE_BLOCK(Vmulvv);
655  reconstruct(world, a, blk, false);
656  if (&a != &b) reconstruct(world, b, blk, false);
657  world.gop.fence();
658 
659  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > q(a.size());
660 
661  unsigned int vvsize = a.size();
662  for (unsigned int i=0; i<vvsize; i+= blk) {
663  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
664  q[j] = mul(a[j], b[j], false);
665  if (fence && (blk !=1 )) world.gop.fence();
666  }
667  if (fence) world.gop.fence();
668  return q;
669  }
670 
671 
672  /// Computes the square of a vector of functions --- q[i] = v[i]**2
673 
674  template <typename T, std::size_t NDIM>
675  std::vector< Function<T,NDIM> >
676  square(World& world,
677  const std::vector< Function<T,NDIM> >& v,
678  bool fence=true) {
679  return mul<T,T,NDIM>(world, v, v, fence);
680  // std::vector< Function<T,NDIM> > vsq(v.size());
681  // for (unsigned int i=0; i<v.size(); ++i) {
682  // vsq[i] = square(v[i], false);
683  // }
684  // if (fence) world.gop.fence();
685  // return vsq;
686  }
687 
688  /// Sets the threshold in a vector of functions
689 
690  template <typename T, std::size_t NDIM>
691  void set_thresh(World& world,
692  std::vector< Function<T,NDIM> >& v,
693  double thresh,
694  bool fence=true) {
695  for (unsigned int j=0; j<v.size(); ++j) {
696  v[j].set_thresh(thresh,false);
697  }
698  if (fence) world.gop.fence();
699  }
700 
701  /// Returns the complex conjugate of the vector of functions
702 
703  template <typename T, std::size_t NDIM>
704  std::vector< Function<T,NDIM> >
705  conj(World& world,
706  const std::vector< Function<T,NDIM> >& v,
707  bool fence=true){
708  PROFILE_BLOCK(Vconj);
709  std::vector< Function<T,NDIM> > r = copy(world, v); // Currently don't have oop conj
710  for (unsigned int i=0; i<v.size(); ++i) {
711  r[i].conj(false);
712  }
713  if (fence) world.gop.fence();
714  return r;
715  }
716 
717 
718  /// Returns a deep copy of a vector of functions
719 
720  template <typename T, std::size_t NDIM>
721  std::vector< Function<T,NDIM> >
722  copy(World& world,
723  const std::vector< Function<T,NDIM> >& v,
724  bool fence=true) {
725  PROFILE_BLOCK(Vcopy);
726  std::vector< Function<T,NDIM> > r(v.size());
727  for (unsigned int i=0; i<v.size(); ++i) {
728  r[i] = copy(v[i], false);
729  }
730  if (fence) world.gop.fence();
731  return r;
732  }
733 
734  /// Returns a vector of deep copies of of a function
735 
736  template <typename T, std::size_t NDIM>
737  std::vector< Function<T,NDIM> >
738  copy(World& world,
739  const Function<T,NDIM>& v,
740  const unsigned int n,
741  bool fence=true) {
742  PROFILE_BLOCK(Vcopy1);
743  std::vector< Function<T,NDIM> > r(n);
744  for (unsigned int i=0; i<n; ++i) {
745  r[i] = copy(v, false);
746  }
747  if (fence) world.gop.fence();
748  return r;
749  }
750 
751  /// Returns new vector of functions --- q[i] = a[i] + b[i]
752 
753  template <typename T, typename R, std::size_t NDIM>
754  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
755  add(World& world,
756  const std::vector< Function<T,NDIM> >& a,
757  const std::vector< Function<R,NDIM> >& b,
758  bool fence=true,
759  unsigned int blk=1) {
760  PROFILE_BLOCK(Vadd);
761  MADNESS_ASSERT(a.size() == b.size());
762  compress(world, a, blk);
763  compress(world, b, blk);
764 
765  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
766 
767  unsigned int vvsize = a.size();
768  for (unsigned int i=0; i<vvsize; i+= blk) {
769  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
770  r[j] = add(a[j], b[j], false);
771  if (fence && (blk !=1 )) world.gop.fence();
772  }
773  if (fence) world.gop.fence();
774 
775  return r;
776  }
777 
778  /// Returns new vector of functions --- q[i] = a + b[i]
779 
780  template <typename T, typename R, std::size_t NDIM>
781  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
782  add(World& world,
783  const Function<T,NDIM> & a,
784  const std::vector< Function<R,NDIM> >& b,
785  bool fence=true,
786  unsigned int blk=1) {
787 
788  PROFILE_BLOCK(Vadd1);
789  a.compress();
790  compress(world, b, blk);
791 
792  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(b.size());
793 
794  unsigned int vvsize = b.size();
795  for (unsigned int i=0; i<vvsize; i+= blk) {
796  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
797  r[j] = add(a, b[j], false);
798  if (fence && (blk !=1 )) world.gop.fence();
799  }
800  if (fence) world.gop.fence();
801 
802  return r;
803  }
804 
805  template <typename T, typename R, std::size_t NDIM>
806  inline std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
807  add(World& world,
808  const std::vector< Function<R,NDIM> >& b,
809  const Function<T,NDIM> & a,
810  bool fence=true,
811  unsigned int blk=1) {
812  return add(world, a, b, fence, blk);
813  }
814 
815  /// Returns new vector of functions --- q[i] = a[i] - b[i]
816 
817  template <typename T, typename R, std::size_t NDIM>
818  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
819  sub(World& world,
820  const std::vector< Function<T,NDIM> >& a,
821  const std::vector< Function<R,NDIM> >& b,
822  bool fence=true,
823  unsigned int blk=1) {
824  PROFILE_BLOCK(Vsub);
825  MADNESS_ASSERT(a.size() == b.size());
826  compress(world, a, fence, blk);
827  compress(world, b, fence, blk);
828 
829  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
830 
831  unsigned int vvsize = a.size();
832  for (unsigned int i=0; i<vvsize; i+= blk) {
833  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
834  r[j] = sub(a[j], b[j], false);
835  if (fence && (blk !=1 )) world.gop.fence();
836  }
837  if (fence) world.gop.fence();
838  return r;
839  }
840 
841 
842  /// Generalized A*X+Y for vectors of functions ---- a[i] = alpha*a[i] + beta*b[i]
843 
844  template <typename T, typename Q, typename R, std::size_t NDIM>
845  void gaxpy(World& world,
846  Q alpha,
847  std::vector< Function<T,NDIM> >& a,
848  Q beta,
849  const std::vector< Function<R,NDIM> >& b,
850  unsigned int blk=1,
851  bool fence=true) {
852  PROFILE_BLOCK(Vgaxpy);
853  MADNESS_ASSERT(a.size() == b.size());
854  compress(world, a, fence, blk);
855  compress(world, b, fence, blk);
856 
857  unsigned int vvsize = a.size();
858 
859  for (unsigned int i=0; i<vvsize; i+= blk) {
860  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
861  a[j].gaxpy(alpha, b[j], beta, false);
862 
863  if (fence && (blk !=1 )) world.gop.fence();
864 
865  }
866  // for (unsigned int i=0; i<a.size(); ++i) {
867  // a[i].gaxpy(alpha, b[i], beta, false);
868  // }
869  if (fence) world.gop.fence();
870  }
871 
872 
873  /// Applies a vector of operators to a vector of functions --- q[i] = apply(op[i],f[i])
874 
875  template <typename opT, typename R, std::size_t NDIM>
876  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> >
877  apply(World& world,
878  const std::vector< std::shared_ptr<opT> >& op,
879  const std::vector< Function<R,NDIM> > f,
880  const unsigned int blk=1){
881 
882  PROFILE_BLOCK(Vapplyv);
883  MADNESS_ASSERT(f.size()==op.size());
884 
885  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
886 
887  reconstruct(world, f, blk);
888  nonstandard(world, ncf, blk);
889 
890  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> > result(f.size());
891  unsigned int ff = f.size();
892 
893  for (unsigned int i=0; i<ff; ++blk) {
894  for (unsigned int j=i; j<std::min(ff,(i+1)*blk); ++j)
895  result[j] = apply_only(*op[j], f[j], false);
896  if (blk !=1)
897  world.gop.fence();
898  }
899 
900  world.gop.fence();
901 
902  standard(world, ncf, false); // restores promise of logical constness
903  world.gop.fence();
904  reconstruct(world, result, blk);
905 
906  return result;
907  }
908 
909 
910  /// Applies an operator to a vector of functions --- q[i] = apply(op,f[i])
911 
912  template <typename T, typename R, std::size_t NDIM>
913  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
914  apply(World& world,
915  const SeparatedConvolution<T,NDIM>& op,
916  const std::vector< Function<R,NDIM> > f,
917  const unsigned int blk=1) {
918  PROFILE_BLOCK(Vapply);
919 
920  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
921 
922  reconstruct(world, f, blk);
923  nonstandard(world, ncf, blk);
924 
925  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> > result(f.size());
926 
927  unsigned int ff = f.size();
928  for (unsigned int i=0; i<ff; ++blk) {
929  for (unsigned int j=i; j<std::min(ff,(i+1)*blk); ++j)
930  result[j] = apply_only(op, f[j], false);
931  if (blk !=1)
932  world.gop.fence();
933  }
934  world.gop.fence();
935 
936  standard(world, ncf, blk, false); // restores promise of logical constness
937  world.gop.fence();
938  reconstruct(world, result, blk);
939 
940  return result;
941  }
942 
943  /// Normalizes a vector of functions --- v[i] = v[i].scale(1.0/v[i].norm2())
944 
945  template <typename T, std::size_t NDIM>
946  void normalize(World& world,
947  std::vector< Function<T,NDIM> >& v,
948  bool fence=true){
949 
950  PROFILE_BLOCK(Vnormalize);
951  std::vector<double> nn = norm2s(world, v);
952 
953  for (unsigned int i=0; i<v.size(); ++i)
954  v[i].scale(1.0/nn[i],false);
955 
956  if (fence) world.gop.fence();
957  }
958 
959 }
960 #endif // MADNESS_MRA_VMRA_H__INCLUDED
q
double q(double t)
Definition: DKops.h:18
madness::BaseTensor::size
long size() const
Returns the number of elements in the tensor.
Definition: basetensor.h:138
madness::Derivative
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
madness::FunctionFactory
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
madness::Function
A multiresolution adaptive numerical function.
Definition: mra.h:122
madness::PoolTaskInterface::make_id
static std::enable_if< detail::function_traits< fnT >::value||detail::memfunc_traits< fnT >::value >::type make_id(std::pair< void *, unsigned short > &id, fnT fn)
Definition: thread.h:916
madness::TaskInterface
All world tasks must be derived from this public interface.
Definition: taskfn.h:69
madness::TaskInterface::world
volatile World * world
Definition: taskfn.h:72
madness::Tensor
A tensor is a multidimension array.
Definition: tensor.h:317
madness::WorldGopInterface::fence
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
madness::WorldGopInterface::sum
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:870
madness::World
A parallel world class.
Definition: world.h:132
madness::World::gop
WorldGopInterface & gop
Global operations.
Definition: world.h:205
R
static const double R
Definition: csqrt.cc:46
derivative.h
Declaration and initialization of tree traversal functions and generic derivative.
m
const double m
Definition: gfit.cc:199
T
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
beta
const double beta
Definition: gygi_soltion.cc:62
v
static const double v
Definition: hatom_sf_dirac.cc:20
op
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
MADNESS_ASSERT
#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
mra.h
Main include file for MADNESS and defines Function interface.
madness
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
madness::sub
Function< TENSOR_RESULT_TYPE(L, R), NDIM > sub(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator- but with optional fence and no automatic compression.
Definition: mra.h:1971
madness::add
Function< TENSOR_RESULT_TYPE(L, R), NDIM > add(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator+ but with optional fence and no automatic compression.
Definition: mra.h:1927
madness::mul
Function< TENSOR_RESULT_TYPE(Q, T), NDIM > mul(const Q alpha, const Function< T, NDIM > &f, bool fence=true)
Returns new function equal to alpha*f(x) with optional fence.
Definition: mra.h:1701
madness::scale
response_space scale(response_space a, double b)
madness::conj
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
madness::apply
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
madness::truncate
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
madness::norm_tree
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
madness::nonstandard
@ nonstandard
s and d coeffs in internal nodes
Definition: funcdefaults.h:61
madness::zero_functions
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
madness::standard
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
madness::copy
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
madness::compress
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:133
madness::square
Function< T, NDIM > square(const Function< T, NDIM > &f, bool fence=true)
Create a new function that is the square of f - global comm only if not reconstructed.
Definition: mra.h:2681
madness::set_thresh
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
madness::norm2
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
madness::_
static const Slice _(0,-1, 1)
madness::vmulXX
std::vector< Function< TENSOR_RESULT_TYPE(L, R), D > > vmulXX(const Function< L, D > &left, const std::vector< Function< R, D > > &vright, double tol, bool fence=true)
Use the vmra/mul(...) interface instead.
Definition: mra.h:1807
madness::TENSOR_RESULT_TYPE
TENSOR_RESULT_TYPE(T, R) inner(const Function< T
Computes the scalar/inner product between two functions.
Definition: vmra.h:1059
madness::f
NDIM & f
Definition: mra.h:2416
madness::mul_sparse
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition: mra.h:1742
madness::g
NDIM const Function< R, NDIM > & g
Definition: mra.h:2416
madness::inner
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
madness::normalize
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition: vmra.h:1614
madness::apply_only
Function< TENSOR_RESULT_TYPE(typename opT::opT, R), NDIM > apply_only(const opT &op, const Function< R, NDIM > &f, bool fence=true)
Apply operator ONLY in non-standard form - required other steps missing !!
Definition: mra.h:2120
madness::transform
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
madness::matrix_inner
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
madness::norm2s
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
madness::gaxpy
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
madness::reconstruct
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
b
static const double b
Definition: nonlinschro.cc:119
a
static const double a
Definition: nonlinschro.cc:118
Q
double Q(double a)
Definition: relops.cc:20
c
static const double c
Definition: relops.cc:10
L
static const double L
Definition: rk.cc:46
thresh
static const double thresh
Definition: rk.cc:45
D
Definition: test_ar.cc:204
madness::MatrixInnerTask
Definition: vmra1.h:463
madness::MatrixInnerTask::get_id
virtual void get_id(std::pair< void *, unsigned short > &id) const
Get the task id.
Definition: vmra1.h:485
madness::MatrixInnerTask::g
const std::vector< Function< R, NDIM > > & g
Definition: vmra1.h:466
madness::MatrixInnerTask::MatrixInnerTask
MatrixInnerTask(const Tensor< TENSOR_RESULT_TYPE(T, R)> &result, const Function< T, NDIM > &f, const std::vector< Function< R, NDIM > > &g, long jtop)
Definition: vmra1.h:469
madness::MatrixInnerTask::f
const Function< T, NDIM > & f
Definition: vmra1.h:465
madness::MatrixInnerTask::result
Tensor< TENSOR_RESULT_TYPE(T, R)> result
Definition: vmra1.h:464
madness::MatrixInnerTask::run
void run(World &world)
Runs a single-threaded task ... derived classes must implement this.
Definition: vmra1.h:475
madness::MatrixInnerTask::jtop
long jtop
Definition: vmra1.h:467
ncf
Definition: dirac-hatom.cc:108
alpha
const double alpha
Definition: test_coulomb.cc:54
NDIM
static const std::size_t NDIM
Definition: testpdiff.cc:42
PROFILE_BLOCK
#define PROFILE_BLOCK(name)
Definition: worldprofile.h:208
Generated on Sat Oct 26 2024 16:26:22 for MADNESS by doxygen 1.9.1