MADNESS  version 0.9
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 
123 #include <madness/mra/mra.h>
124 #include <madness/mra/derivative.h>
125 #include <cstdio>
126 
127 namespace madness {
128 
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 
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 
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].nonstandard(false,false);
187  }
188  if ( blk!=1 && fence) world.gop.fence();
189  }
190  if (fence) world.gop.fence();
191  } //nonstand
192 
193 
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 
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 
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 
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 
275 
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<resultT,NDIM>(world, m);
298  compress(world, v, blk);
299  compress(world, vc, blk); // what to do ?
300 
301  for (unsigned int i=0; i<m; i+= blki) {
302  for (unsigned int ii=i; ii<std::min(m,(i+1)*blki); ii++) {
303  for (unsigned int j=0; j<n; j+= blkj) {
304  for (unsigned int jj=j; jj<std::min(n, (j+1)*blkj); jj++)
305  if (c(jj,ii) != R(0.0)) vc[ii].gaxpy(1.0,v[jj],c(jj,ii),false);
306  if (fence && (blkj!=1)) world.gop.fence();
307  }
308  }
309  if (fence && (blki!=1)) world.gop.fence(); // a bit conservative
310  }
311 
312 
313  // for (unsigned int i=0; i<m; ++i) {
314  // for (unsigned int j=0; j<n; ++j) {
315  // if (c(j,i) != R(0.0)) vc[i].gaxpy(1.0,v[j],c(j,i),false);
316  // }
317  // }
318 
319  if (fence) world.gop.fence();
320  return vc;
321  }
322 
323 
324  template <typename L, typename R, std::size_t NDIM>
325  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> >
326  transform(World& world,
327  const std::vector< Function<L,NDIM> >& v,
328  const Tensor<R>& c,
329  const double tol,
330  const unsigned int blki=1,
331  const bool fence)
332  {
333  PROFILE_BLOCK(Vtransform);
334  MADNESS_ASSERT(v.size() == (unsigned int)(c.dim(0)));
335 
336  std::vector< Function<TENSOR_RESULT_TYPE(L,R),NDIM> > vresult(c.dim(1));
337  unsigned int m=c.dim(1);
338 
339  for (unsigned int i=0; i<m; i+= blki) {
340  for (unsigned int ii=i; ii<std::min(m,(i+1)*blki); ii++) {
341  vresult[ii] = Function<TENSOR_RESULT_TYPE(L,R),NDIM>(FunctionFactory<TENSOR_RESULT_TYPE(L,R),NDIM>(world));
342  }
343  if (fence && (blki!=1)) world.gop.fence(); // a bit conservative
344  }
345 
346 
347  // for (unsigned int i=0; i<c.dim(1); ++i) {
348  // vresult[i] = Function<TENSOR_RESULT_TYPE(L,R),NDIM>(FunctionFactory<TENSOR_RESULT_TYPE(L,R),NDIM>(world));
349  // }
350  compress(world, v, blki, false);
351  compress(world, vresult, blki, false);
352  world.gop.fence();
353  vresult[0].vtransform(v, c, vresult, tol, fence);
354  return vresult;
355  }
356 
358 
359  template <typename T, typename Q, std::size_t NDIM>
360  void scale(World& world,
361  std::vector< Function<T,NDIM> >& v,
362  const std::vector<Q>& factors,
363  const unsigned int blk=1,
364  const bool fence=true)
365  {
366  PROFILE_BLOCK(Vscale);
367 
368  unsigned int vvsize = v.size();
369  for (unsigned int i=0; i<vvsize; i+= blk) {
370  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
371  v[j].scale(factors[j],false);
372  }
373  if (fence && blk!=1 ) world.gop.fence();
374  }
375  if (fence) world.gop.fence();
376  }
377 
379 
380  template <typename T, typename Q, std::size_t NDIM>
381  void scale(World& world,
382  std::vector< Function<T,NDIM> >& v,
383  const Q factor,
384  const unsigned int blk=1,
385  const bool fence=true){
386 
387  PROFILE_BLOCK(Vscale); // shouldn't need blocking since it is local
388 
389  unsigned int vvsize = v.size();
390  for (unsigned int i=0; i<vvsize; i+= blk) {
391  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
392  v[j].scale(factor,false);
393  }
394  if (fence && blk!=1 ) world.gop.fence();
395  }
396  if (fence) world.gop.fence();
397  }
398 
400  template <typename T, std::size_t NDIM>
401  std::vector<double> norm2s(World& world,
402  const std::vector< Function<T,NDIM> >& v,
403  const unsigned int blk=1,
404  const bool fence=true){
405 
406  PROFILE_BLOCK(Vnorm2);
407  unsigned int vvsize = v.size();
408  std::vector<double> norms(vvsize);
409 
410  for (unsigned int i=0; i<vvsize; i+= blk) {
411  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j) {
412  norms[j] = v[j].norm2sq_local();
413  }
414  if (fence && (blk!=1)) world.gop.fence();
415  }
416  if (fence ) world.gop.fence();
417 
418  world.gop.sum(&norms[0], norms.size());
419 
420  for (unsigned int i=0; i<vvsize; i+= blk) {
421  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
422  norms[j] = sqrt(norms[j]);
423  if (fence && (blk!=1)) world.gop.fence();
424  }
425 
426  world.gop.fence();
427  return norms;
428  }
429 
431  // should be local; norms[0] contains the result
432 
433  template <typename T, std::size_t NDIM>
434  double norm2(World& world,
435  const std::vector< Function<T,NDIM> >& v)
436  {
437 
438  PROFILE_BLOCK(Vnorm2);
439 
440  std::vector<double> norms(v.size());
441 
442  for (unsigned int i=0; i<v.size(); ++i)
443  norms[i] = v[i].norm2sq_local();
444 
445  world.gop.sum(&norms[0], norms.size());
446 
447  for (unsigned int i=1; i<v.size(); ++i)
448  norms[0] += norms[i];
449 
450  world.gop.fence();
451  return sqrt(norms[0]);
452  }
453 
454  inline double conj(double x) {
455  return x;
456  }
457 
458  inline double conj(float x) {
459  return x;
460  }
461 
462 
463  template <typename T, typename R, std::size_t NDIM>
464  struct MatrixInnerTask : public TaskInterface {
465  Tensor<TENSOR_RESULT_TYPE(T,R)> result; // Must be a copy
466  const Function<T,NDIM>& f;
467  const std::vector< Function<R,NDIM> >& g;
468  long jtop;
469 
470  MatrixInnerTask(const Tensor<TENSOR_RESULT_TYPE(T,R)>& result,
471  const Function<T,NDIM>& f,
472  const std::vector< Function<R,NDIM> >& g,
473  long jtop)
474  : result(result), f(f), g(g), jtop(jtop) {}
475 
476  void run(World& world) {
477  for (long j=0; j<jtop; ++j) {
478  result(j) = f.inner_local(g[j]);
479  }
480  }
481 
482  private:
484 
486  virtual void get_id(std::pair<void*,unsigned short>& id) const {
487  PoolTaskInterface::make_id(id, *this);
488  }
489  }; // struct MatrixInnerTask
490 
492 
496 
497  template <typename T, typename R, std::size_t NDIM>
498  Tensor< TENSOR_RESULT_TYPE(T,R) > matrix_inner(World& world,
499  const std::vector< Function<T,NDIM> >& f,
500  const std::vector< Function<R,NDIM> >& g,
501  bool sym=false) {
502  PROFILE_BLOCK(Vmatrix_inner);
503  unsigned int n=f.size(), m=g.size();
504  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n,m);
505  if (sym) MADNESS_ASSERT(n==m);
506 
507  world.gop.fence();
508  compress(world, f);
509  if (&f != &g) compress(world, g);
510 
511  // for (long i=0; i<n; ++i) {
512  // long jtop = m;
513  // if (sym) jtop = i+1;
514  // for (long j=0; j<jtop; ++j) {
515  // r(i,j) = f[i].inner_local(g[j]);
516  // if (sym) r(j,i) = conj(r(i,j));
517  // }
518  // }
519 
520  for (unsigned int i=n-1; i>=0; --i) {
521  unsigned int jtop = m;
522  if (sym) jtop = i+1;
523  world.taskq.add(new MatrixInnerTask<T,R,NDIM>(r(i,_), f[i], g, jtop));
524  }
525  world.gop.fence();
526  world.gop.sum(r.ptr(),n*m);
527 
528  if (sym) {
529  for (unsigned int i=0; i<n; ++i) {
530  for (unsigned int j=0; j<i; ++j) {
531  r(j,i) = conj(r(i,j));
532  }
533  }
534  }
535  return r;
536  }
537 
539 
540  template <typename T, typename R, std::size_t NDIM>
541  Tensor< TENSOR_RESULT_TYPE(T,R) > inner(World& world,
542  const std::vector< Function<T,NDIM> >& f,
543  const std::vector< Function<R,NDIM> >& g) {
544  PROFILE_BLOCK(Vinnervv);
545  long n=f.size(), m=g.size();
546  MADNESS_ASSERT(n==m);
547  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
548 
549  compress(world, f);
550  compress(world, g);
551 
552  for (long i=0; i<n; ++i) {
553  r(i) = f[i].inner_local(g[i]);
554  }
555 
556  world.taskq.fence();
557  world.gop.sum(r.ptr(),n);
558  world.gop.fence();
559  return r;
560  }
561 
562 
564 
565  template <typename T, typename R, std::size_t NDIM>
566  Tensor< TENSOR_RESULT_TYPE(T,R) > inner(World& world,
567  const Function<T,NDIM>& f,
568  const std::vector< Function<R,NDIM> >& g) {
569  PROFILE_BLOCK(Vinner);
570  long n=g.size();
571  Tensor< TENSOR_RESULT_TYPE(T,R) > r(n);
572 
573  f.compress();
574  compress(world, g);
575 
576  for (long i=0; i<n; ++i) {
577  r(i) = f.inner_local(g[i]);
578  }
579 
580  world.taskq.fence();
581  world.gop.sum(r.ptr(),n);
582  world.gop.fence();
583  return r;
584  }
585 
586 
588 
589  template <typename T, typename R, std::size_t NDIM>
590  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
591  mul(World& world,
592  const Function<T,NDIM>& a,
593  const std::vector< Function<R,NDIM> >& v,
594  const unsigned int blk=1,
595  const bool fence=true) {
596 
597  PROFILE_BLOCK(Vmul);
598  a.reconstruct(false);
599  reconstruct(world, v, blk, false);
600  world.gop.fence();
601  return vmulXX(a, v, 0.0, fence);
602  }
603 
605  template <typename T, typename R, std::size_t NDIM>
606  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
607  mul_sparse(World& world,
608  const Function<T,NDIM>& a,
609  const std::vector< Function<R,NDIM> >& v,
610  const double tol,
611  const bool fence=true,
612  const unsigned int blk=1)
613  {
614  PROFILE_BLOCK(Vmulsp);
615  a.reconstruct(false);
616  reconstruct(world, v, blk, false);
617  world.gop.fence();
618 
619  unsigned int vvsize = v.size();
620  for (unsigned int i=0; i<vvsize; i+= blk) {
621  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
622  v[j].norm_tree(false);
623  if ( fence && (blk == 1)) world.gop.fence();
624  }
625  a.norm_tree();
626  return vmulXX(a, v, tol, fence);
627  }
628 
630  template <typename T, std::size_t NDIM>
631  void norm_tree(World& world,
632  const std::vector< Function<T,NDIM> >& v,
633  bool fence=true,
634  unsigned int blk=1){
635  PROFILE_BLOCK(Vnorm_tree);
636 
637  unsigned int vvsize = v.size();
638  for (unsigned int i=0; i<vvsize; i+= blk) {
639  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
640  v[j].norm_tree(false);
641  if (fence && blk!=1 ) world.gop.fence();
642  }
643  if (fence) world.gop.fence();
644  }
645 
647 
648  template <typename T, typename R, std::size_t NDIM>
649  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
650  mul(World& world,
651  const std::vector< Function<T,NDIM> >& a,
652  const std::vector< Function<R,NDIM> >& b,
653  bool fence=true,
654  unsigned int blk=1){
655  PROFILE_BLOCK(Vmulvv);
656  reconstruct(world, a, blk, false);
657  if (&a != &b) reconstruct(world, b, blk, false);
658  world.gop.fence();
659 
660  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > q(a.size());
661 
662  unsigned int vvsize = a.size();
663  for (unsigned int i=0; i<vvsize; i+= blk) {
664  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
665  q[j] = mul(a[j], b[j], false);
666  if (fence && (blk !=1 )) world.gop.fence();
667  }
668  if (fence) world.gop.fence();
669  return q;
670  }
671 
672 
674 
675  template <typename T, std::size_t NDIM>
676  std::vector< Function<T,NDIM> >
677  square(World& world,
678  const std::vector< Function<T,NDIM> >& v,
679  bool fence=true) {
680  return mul<T,T,NDIM>(world, v, v, fence);
681  // std::vector< Function<T,NDIM> > vsq(v.size());
682  // for (unsigned int i=0; i<v.size(); ++i) {
683  // vsq[i] = square(v[i], false);
684  // }
685  // if (fence) world.gop.fence();
686  // return vsq;
687  }
688 
690 
691  template <typename T, std::size_t NDIM>
692  void set_thresh(World& world,
693  std::vector< Function<T,NDIM> >& v,
694  double thresh,
695  bool fence=true) {
696  for (unsigned int j=0; j<v.size(); ++j) {
697  v[j].set_thresh(thresh,false);
698  }
699  if (fence) world.gop.fence();
700  }
701 
703 
704  template <typename T, std::size_t NDIM>
705  std::vector< Function<T,NDIM> >
706  conj(World& world,
707  const std::vector< Function<T,NDIM> >& v,
708  bool fence=true){
709  PROFILE_BLOCK(Vconj);
710  std::vector< Function<T,NDIM> > r = copy(world, v); // Currently don't have oop conj
711  for (unsigned int i=0; i<v.size(); ++i) {
712  r[i].conj(false);
713  }
714  if (fence) world.gop.fence();
715  return r;
716  }
717 
718 
720 
721  template <typename T, std::size_t NDIM>
722  std::vector< Function<T,NDIM> >
723  copy(World& world,
724  const std::vector< Function<T,NDIM> >& v,
725  bool fence=true) {
726  PROFILE_BLOCK(Vcopy);
727  std::vector< Function<T,NDIM> > r(v.size());
728  for (unsigned int i=0; i<v.size(); ++i) {
729  r[i] = copy(v[i], false);
730  }
731  if (fence) world.gop.fence();
732  return r;
733  }
734 
736 
737  template <typename T, std::size_t NDIM>
738  std::vector< Function<T,NDIM> >
739  copy(World& world,
740  const Function<T,NDIM>& v,
741  const unsigned int n,
742  bool fence=true) {
743  PROFILE_BLOCK(Vcopy1);
744  std::vector< Function<T,NDIM> > r(n);
745  for (unsigned int i=0; i<n; ++i) {
746  r[i] = copy(v, false);
747  }
748  if (fence) world.gop.fence();
749  return r;
750  }
751 
753 
754  template <typename T, typename R, std::size_t NDIM>
755  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
756  add(World& world,
757  const std::vector< Function<T,NDIM> >& a,
758  const std::vector< Function<R,NDIM> >& b,
759  bool fence=true,
760  unsigned int blk=1) {
761  PROFILE_BLOCK(Vadd);
762  MADNESS_ASSERT(a.size() == b.size());
763  compress(world, a, blk);
764  compress(world, b, blk);
765 
766  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
767 
768  unsigned int vvsize = a.size();
769  for (unsigned int i=0; i<vvsize; i+= blk) {
770  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
771  r[j] = add(a[j], b[j], false);
772  if (fence && (blk !=1 )) world.gop.fence();
773  }
774  if (fence) world.gop.fence();
775 
776  return r;
777  }
778 
780 
781  template <typename T, typename R, std::size_t NDIM>
782  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
783  add(World& world,
784  const Function<T,NDIM> & a,
785  const std::vector< Function<R,NDIM> >& b,
786  bool fence=true,
787  unsigned int blk=1) {
788 
789  PROFILE_BLOCK(Vadd1);
790  a.compress();
791  compress(world, b, blk);
792 
793  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(b.size());
794 
795  unsigned int vvsize = b.size();
796  for (unsigned int i=0; i<vvsize; i+= blk) {
797  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
798  r[j] = add(a, b[j], false);
799  if (fence && (blk !=1 )) world.gop.fence();
800  }
801  if (fence) world.gop.fence();
802 
803  return r;
804  }
805 
806  template <typename T, typename R, std::size_t NDIM>
807  inline std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
808  add(World& world,
809  const std::vector< Function<R,NDIM> >& b,
810  const Function<T,NDIM> & a,
811  bool fence=true,
812  unsigned int blk=1) {
813  return add(world, a, b, fence, blk);
814  }
815 
817 
818  template <typename T, typename R, std::size_t NDIM>
819  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
820  sub(World& world,
821  const std::vector< Function<T,NDIM> >& a,
822  const std::vector< Function<R,NDIM> >& b,
823  bool fence=true,
824  unsigned int blk=1) {
825  PROFILE_BLOCK(Vsub);
826  MADNESS_ASSERT(a.size() == b.size());
827  compress(world, a, fence, blk);
828  compress(world, b, fence, blk);
829 
830  std::vector< Function<TENSOR_RESULT_TYPE(T,R),NDIM> > r(a.size());
831 
832  unsigned int vvsize = a.size();
833  for (unsigned int i=0; i<vvsize; i+= blk) {
834  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
835  r[j] = sub(a[j], b[j], false);
836  if (fence && (blk !=1 )) world.gop.fence();
837  }
838  if (fence) world.gop.fence();
839  return r;
840  }
841 
842 
844 
845  template <typename T, typename Q, typename R, std::size_t NDIM>
846  void gaxpy(World& world,
847  Q alpha,
848  std::vector< Function<T,NDIM> >& a,
849  Q beta,
850  const std::vector< Function<R,NDIM> >& b,
851  unsigned int blk=1,
852  bool fence=true) {
853  PROFILE_BLOCK(Vgaxpy);
854  MADNESS_ASSERT(a.size() == b.size());
855  compress(world, a, fence, blk);
856  compress(world, b, fence, blk);
857 
858  unsigned int vvsize = a.size();
859 
860  for (unsigned int i=0; i<vvsize; i+= blk) {
861  for (unsigned int j=i; j<std::min(vvsize,(i+1)*blk); ++j)
862  a[j].gaxpy(alpha, b[j], beta, false);
863 
864  if (fence && (blk !=1 )) world.gop.fence();
865 
866  }
867  // for (unsigned int i=0; i<a.size(); ++i) {
868  // a[i].gaxpy(alpha, b[i], beta, false);
869  // }
870  if (fence) world.gop.fence();
871  }
872 
873 
875 
876  template <typename opT, typename R, std::size_t NDIM>
877  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> >
878  apply(World& world,
879  const std::vector< std::shared_ptr<opT> >& op,
880  const std::vector< Function<R,NDIM> > f,
881  const unsigned int blk=1){
882 
883  PROFILE_BLOCK(Vapplyv);
884  MADNESS_ASSERT(f.size()==op.size());
885 
886  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
887 
888  reconstruct(world, f, blk);
889  nonstandard(world, ncf, blk);
890 
891  std::vector< Function<TENSOR_RESULT_TYPE(typename opT::opT,R), NDIM> > result(f.size());
892  unsigned int ff = f.size();
893 
894  for (unsigned int i=0; i<ff; ++blk) {
895  for (unsigned int j=i; j<std::min(ff,(i+1)*blk); ++j)
896  result[j] = apply_only(*op[j], f[j], false);
897  if (blk !=1)
898  world.gop.fence();
899  }
900 
901  world.gop.fence();
902 
903  standard(world, ncf, false); // restores promise of logical constness
904  world.gop.fence();
905  reconstruct(world, result, blk);
906 
907  return result;
908  }
909 
910 
912 
913  template <typename T, typename R, std::size_t NDIM>
914  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> >
915  apply(World& world,
916  const SeparatedConvolution<T,NDIM>& op,
917  const std::vector< Function<R,NDIM> > f,
918  const unsigned int blk=1) {
919  PROFILE_BLOCK(Vapply);
920 
921  std::vector< Function<R,NDIM> >& ncf = *const_cast< std::vector< Function<R,NDIM> >* >(&f);
922 
923  reconstruct(world, f, blk);
924  nonstandard(world, ncf, blk);
925 
926  std::vector< Function<TENSOR_RESULT_TYPE(T,R), NDIM> > result(f.size());
927 
928  unsigned int ff = f.size();
929  for (unsigned int i=0; i<ff; ++blk) {
930  for (unsigned int j=i; j<std::min(ff,(i+1)*blk); ++j)
931  result[j] = apply_only(op, f[j], false);
932  if (blk !=1)
933  world.gop.fence();
934  }
935  world.gop.fence();
936 
937  standard(world, ncf, blk, false); // restores promise of logical constness
938  world.gop.fence();
939  reconstruct(world, result, blk);
940 
941  return result;
942  }
943 
945 
946  template <typename T, std::size_t NDIM>
947  void normalize(World& world,
948  std::vector< Function<T,NDIM> >& v,
949  bool fence=true){
950 
951  PROFILE_BLOCK(Vnormalize);
952  std::vector<double> nn = norm2s(world, v);
953 
954  for (unsigned int i=0; i<v.size(); ++i)
955  v[i].scale(1.0/nn[i],false);
956 
957  if (fence) world.gop.fence();
958  }
959 
960 }
961 #endif // MADNESS_MRA_VMRA_H__INCLUDED
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: chem/distpm.cc:38
madness::World::gop
WorldGopInterface & gop
Global operations.
Definition: worldfwd.h:462
madness::MatrixInnerTask::jtop
long jtop
Definition: vmra1.h:468
madness::tr1::shptr::shared_ptr
Definition: shared_ptr_bits.h:38
madness::reconstruct
void reconstruct(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Reconstruct a vector of functions.
Definition: vmra.h:149
R
const double R
Definition: dielectric.cc:191
mra.h
Main include file for MADNESS and defines Function interface.
madness::gaxpy
void gaxpy(World &world, Q alpha, std::vector< Function< T, NDIM > > &a, Q beta, const std::vector< Function< R, NDIM > > &b, bool fence=true)
Generalized A*X+Y for vectors of functions -— a[i] = alpha*a[i] + beta*b[i].
Definition: vmra.h:680
madness::MatrixInnerTask::run
void run(World &world)
Runs a single-threaded task ... derived classes must implement this.
Definition: vmra1.h:476
madness::Function::compress
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition: mra.h:683
madness::Derivative
Implements derivatives operators with variety of boundary conditions on simulation domain...
Definition: derivative.h:272
madness::truncate
void truncate(World &world, std::vector< Function< T, NDIM > > &v, double tol=0.0, bool fence=true)
Truncates a vector of functions.
Definition: vmra.h:194
madness::nonstandard
void nonstandard(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Generates non-standard form of a vector of functions.
Definition: vmra.h:167
madness::f
NDIM & f
Definition: mra.h:2179
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:1942
derivative.h
Declaration and initialization of tree traversal functions and generic derivative.
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:1734
madness::TENSOR_RESULT_TYPE
TENSOR_RESULT_TYPE(T, R) inner(const Function<T
Computes the scalar/inner product between two functions.
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:506
beta
const double beta
Definition: gygi_soltion.cc:63
madness::inner
T inner(const vecfunc< T, NDIM > &a, const vecfunc< T, NDIM > &b)
the non-linear solver requires an inner product
Definition: nemo.h:112
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:312
madness::MatrixInnerTask
Definition: vmra1.h:464
PROFILE_BLOCK
#define PROFILE_BLOCK(name)
Definition: worldprofile.h:197
std::tr1::T
const T1 &f1 return GTEST_2_TUPLE_() T(f0, f1)
madness::scale
void scale(World &world, std::vector< Function< T, NDIM > > &v, const std::vector< Q > &factors, bool fence=true)
Scales inplace a vector of functions by distinct values.
Definition: vmra.h:290
madness::Function::reconstruct
void reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm...
Definition: mra.h:737
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:1835
mpfr::min
const mpreal min(const mpreal &x, const mpreal &y)
Definition: mpreal.h:2675
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:1778
a
FLOAT a(int j, FLOAT z)
Definition: y1.cc:86
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:1634
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:470
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:556
madness::WorldGopInterface::sum
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:767
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:181
sqrt
tensorT sqrt(const tensorT &s, double tol=1e-8)
Computes matrix square root (not used any more?)
Definition: DFcode/moldft.cc:446
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:1806
madness::World
A parallel world with full functionality wrapping an MPI communicator.
Definition: worldfwd.h:416
madness::WorldGopInterface::fence
void fence()
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:52
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:325
madness::Function
A multiresolution adaptive numerical function.
Definition: derivative.h:61
madness::MatrixInnerTask::g
const std::vector< Function< R, NDIM > > & g
Definition: vmra1.h:467
m
const double m
Definition: gfit.cc:199
madness::Function::fence
impl world gop fence()
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 alread...
Definition: mra.h:1569
madness::PoolTaskInterface::make_id
static enable_if_c< detail::function_traits< fnT >::value||detail::memfunc_traits< fnT >::value >::type make_id(std::pair< void *, unsigned short > &id, fnT fn)
Definition: worldthread.h:680
madness::Function::size
std::size_t size() const
Returns the number of coefficients in the function ... collective global sum.
Definition: mra.h:515
op
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
madness::MatrixInnerTask::f
const Function< T, NDIM > & f
Definition: vmra1.h:466
madness::Function::norm_tree
void norm_tree(bool fence=true) const
Initializes information about the function norm at all length scales.
Definition: mra.h:663
madness::TaskInterface
All world tasks must be derived from this public interface.
Definition: taskfn.h:68
madness::FunctionFactory
FunctionFactory implements the named-parameter idiom for Function.
Definition: funcimpl.h:70
madness::MatrixInnerTask::result
Tensor< TENSOR_RESULT_TYPE(T, R)> result
Definition: vmra1.h:465
madness
Holds machinery to set up Functions/FuncImpls using various Factories and Interfaces.
Definition: chem/atomutil.cc:45
madness::apply
Function< T, NDIM > apply(const Derivative< T, NDIM > &D, const Function< T, NDIM > &f, bool fence=true)
Applies derivative operator to function (for syntactic equivalence to integral operator apply) ...
Definition: derivative.h:613
c
const double c
Definition: gfit.cc:200
madness::transform
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const DistributedMatrix< R > &c, bool fence=true)
Definition: chem/SCF.cc:86
b
FLOAT b(int j, FLOAT z)
Definition: y1.cc:79
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:765
madness::zero_functions
std::vector< Function< T, NDIM > > zero_functions(World &world, int n)
Generates a vector of zero functions.
Definition: vmra.h:229
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:1528
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:1879
madness::compress
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:130
Generated on Fri Sep 26 2014 10:38:59 for MADNESS by   doxygen 1.8.8