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
127namespace 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::SeparatedConvolution
Convolutions in separated form (including Gaussian)
Definition operator.h:136
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.
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
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
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::scale
response_space scale(response_space a, double b)
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::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::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::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::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::nonstandard
@ nonstandard
s and d coeffs in internal nodes
Definition funcdefaults.h:61
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::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::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::reconstruct
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
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::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::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::f
NDIM & f
Definition mra.h:2416
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::g
NDIM const Function< R, NDIM > & g
Definition mra.h:2416
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::inner
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
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::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:1671
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::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::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::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
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
m
static const double m
Definition relops.cc:9
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
static const double alpha
Definition testcosine.cc:10
NDIM
static const std::size_t NDIM
Definition testpdiff.cc:42
TENSOR_RESULT_TYPE
#define TENSOR_RESULT_TYPE(L, R)
This macro simplifies access to TensorResultType.
Definition type_data.h:205
PROFILE_BLOCK
#define PROFILE_BLOCK(name)
Definition worldprofile.h:208
Generated on Thu Dec 19 2024 21:22:53 for MADNESS by doxygen 1.9.8