MADNESS 0.10.1
convolution1d.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
32#ifndef MADNESS_MRA_CONVOLUTION1D_H__INCLUDED
33#define MADNESS_MRA_CONVOLUTION1D_H__INCLUDED
34
36#include <madness/constants.h>
37#include <limits.h>
40#include <madness/mra/adquad.h>
45#include <algorithm>
46
47/// \file mra/convolution1d.h
48/// \brief Computes most matrix elements over 1D operators (including Gaussians)
49
50/// \ingroup function
51
52namespace madness {
53
54 void aligned_add(long n, double* MADNESS_RESTRICT a, const double* MADNESS_RESTRICT b);
55 void aligned_sub(long n, double* MADNESS_RESTRICT a, const double* MADNESS_RESTRICT b);
58
59 template <typename T>
60 static void copy_2d_patch(T* MADNESS_RESTRICT out, long ldout, const T* MADNESS_RESTRICT in, long ldin, long n, long m) {
61 for (long i=0; i<n; ++i, out+=ldout, in+=ldin) {
62 for (long j=0; j<m; ++j) {
63 out[j] = in[j];
64 }
65 }
66 }
67
68 /// a(n,m) --> b(m,n) ... optimized for smallish matrices
69 template <typename T>
70 inline void fast_transpose(long n, long m, const T* a, T* MADNESS_RESTRICT b) {
71 // n will always be k or 2k (k=wavelet order) and m will be anywhere
72 // from 2^(NDIM-1) to (2k)^(NDIM-1).
73
74// for (long i=0; i<n; ++i)
75// for (long j=0; j<m; ++j)
76// b[j*n+i] = a[i*m+j];
77// return;
78
79 if (n==1 || m==1) {
80 long nm=n*m;
81 for (long i=0; i<nm; ++i) b[i] = a[i];
82 return;
83 }
84
85 long n4 = (n>>2)<<2;
86 long m4 = m<<2;
87 const T* a0 = a;
88 for (long i=0; i<n4; i+=4, a0+=m4) {
89 const T* a1 = a0+m;
90 const T* a2 = a1+m;
91 const T* a3 = a2+m;
93 for (long j=0; j<m; ++j, bi+=n) {
94 T tmp0 = a0[j];
95 T tmp1 = a1[j];
96 T tmp2 = a2[j];
97 T tmp3 = a3[j];
98
99 bi[0] = tmp0;
100 bi[1] = tmp1;
101 bi[2] = tmp2;
102 bi[3] = tmp3;
103 }
104 }
105
106 for (long i=n4; i<n; ++i)
107 for (long j=0; j<m; ++j)
108 b[j*n+i] = a[i*m+j];
109
110 }
111
112 // /// a(i,j) --> b(i,j) for i=0..n-1 and j=0..r-1 noting dimensions are a(n,m) and b(n,r).
113
114 // /// returns b
115 // template <typename T>
116 // inline T* shrink(long n, long m, long r, const T* a, T* MADNESS_RESTRICT b) {
117 // T* result = b;
118 // if (r == 0) {
119 // ;
120 // }
121 // else if (r == 1) {
122 // for (long i=0; i<n; ++i) {
123 // b[i] = a[i];
124 // }
125 // }
126 // else if (r == 2) {
127 // for (long i=0; i<n; ++i, a+=m, b+=2) {
128 // b[0] = a[0];
129 // b[1] = a[1];
130 // }
131 // }
132 // else if (r == 3) {
133 // for (long i=0; i<n; ++i, a+=m, b+=3) {
134 // b[0] = a[0];
135 // b[1] = a[1];
136 // b[2] = a[2];
137 // }
138 // }
139 // else if (r == 4) {
140 // for (long i=0; i<n; ++i, a+=m, b+=4) {
141 // b[0] = a[0];
142 // b[1] = a[1];
143 // b[2] = a[2];
144 // b[3] = a[3];
145 // }
146 // }
147 // else {
148 // for (long i=0; i<n; ++i, a+=m, b+=r) {
149 // for (long j=0; j<r; j++) {
150 // b[j] = a[j];
151 // }
152 // }
153 // }
154 // return result;
155 // }
156
157 /// actual data for 1 dimension and for 1 term and for 1 displacement for a convolution operator
158 /// here we keep the transformation matrices
159
160 /// !!! Note that if Rnormf is zero then ***ALL*** of the tensors are empty
161 template <typename Q>
163
164
165 Tensor<Q> R, T; ///< if NS: R=ns, T=T part of ns; if modified NS: T=\uparrow r^(n-1)
166 Tensor<Q> RU, RVT, TU, TVT; ///< SVD approximations to R and T
167 Tensor<typename Tensor<Q>::scalar_type> Rs, Ts; ///< hold relative errors, NOT the singular values..
168
169 // norms for NS form
171
172 // norms for modified NS form
173 double N_up, N_diff, N_F; ///< the norms according to Beylkin 2008, Eq. (21) ff
174
175
176 /// ctor for NS form
177 /// make the operator matrices r^n and \uparrow r^(n-1)
178 /// @param[in] R operator matrix of the requested level; NS: unfilter(r^(n+1)); modified NS: r^n
179 /// @param[in] T upsampled operator matrix from level n-1; NS: r^n; modified NS: filter( r^(n-1) )
180 ConvolutionData1D(const Tensor<Q>& R, const Tensor<Q>& T) : R(R), T(T) {
181 Rnormf = R.normf();
182 // Making the approximations is expensive ... only do it for
183 // significant components
184 if (Rnormf > 1e-20) {
185 Tnormf = T.normf();
186 make_approx(T, TU, Ts, TVT, Tnorm);
187 make_approx(R, RU, Rs, RVT, Rnorm);
188 int k = T.dim(0);
189
190 Tensor<Q> NS = copy(R);
191 for (int i=0; i<k; ++i)
192 for (int j=0; j<k; ++j)
193 NS(i,j) = 0.0;
194 NSnormf = NS.normf();
195
196 }
197 else {
198 Rnorm = Tnorm = Rnormf = Tnormf = NSnormf = 0.0;
199 N_F = N_up = N_diff = 0.0;
200 }
201 }
202
203 /// ctor for modified NS form
204 /// make the operator matrices r^n and \uparrow r^(n-1)
205 /// @param[in] R operator matrix of the requested level; NS: unfilter(r^(n+1)); modified NS: r^n
206 /// @param[in] T upsampled operator matrix from level n-1; NS: r^n; modified NS: filter( r^(n-1) )
207 /// @param[in] modified use (un) modified NS form
208 ConvolutionData1D(const Tensor<Q>& R, const Tensor<Q>& T, const bool modified) : R(R), T(T) {
209
210 // note that R can be small, but T still be large
211
212 Rnormf = R.normf();
213 Tnormf = T.normf();
214 // Making the approximations is expensive ... only do it for
215 // significant components
216 if (Rnormf > 1e-20) make_approx(R, RU, Rs, RVT, Rnorm);
217 if (Tnormf > 1e-20) make_approx(T, TU, Ts, TVT, Tnorm);
218
219 // norms for modified NS form: follow Beylkin, 2008, Eq. (21) ff
220 N_F=Rnormf;
221 N_up=Tnormf;
222 N_diff=(R-T).normf();
223 }
224
225
226
227 /// approximate the operator matrices using SVD, and abuse Rs to hold the error instead of
228 /// the singular values (seriously, who named this??)
230 Tensor<Q>& RU, Tensor<typename Tensor<Q>::scalar_type>& Rs, Tensor<Q>& RVT, double& norm) {
231 int n = R.dim(0);
232 svd(R, RU, Rs, RVT);
233 for (int i=0; i<n; ++i) {
234 for (int j=0; j<n; ++j) {
235 RVT(i,j) *= Rs[i];
236 }
237 }
238 for (int i=n-1; i>1; --i) { // Form cumulative sum of norms
239 Rs[i-1] += Rs[i];
240 }
241
242 norm = Rs[0];
243 if (Rs[0]>0.0) { // Turn into relative errors
244 double rnorm = 1.0/norm;
245 for (int i=0; i<n; ++i) {
246 Rs[i] *= rnorm;
247 }
248 }
249 }
250 };
251
252 /// Provides the common functionality/interface of all 1D convolutions
253
254 /// interface for 1 term and for 1 dimension;
255 /// the actual data are kept in ConvolutionData1D
256 /// Derived classes must implement rnlp, issmall, natural_level
257 template <typename Q>
259 public:
260 typedef Q opT; ///< The apply function uses this to infer resultT=opT*inputT
261 int k; ///< Wavelet order
262 int npt; ///< Number of quadrature points (is this used?)
263 int maxR; ///< Number of lattice translations for sum
264 double bloch_k; ///< k in exp(i k R) Bloch phase factor folded into lattice sum
265 KernelRange range; ///< if range is nonnull, kernel range limited to to range (in simulation cell units), useful for finite-range convolutions with periodic functions
271
276
277 bool lattice_summed() const { return maxR != 0; }
278 bool range_restricted() const { return range.finite(); }
279
280 virtual ~Convolution1D() {};
281
282 Convolution1D(int k, int npt, int maxR,
283 double bloch_k = 0.0,
284 KernelRange rng = {})
285 : k(k)
286 , npt(npt)
287 , maxR(maxR)
288 , quad_x(npt)
289 , quad_w(npt)
291 , range(rng)
292 {
293 auto success = autoc(k,&c);
295
296 gauss_legendre(npt,0.0,1.0,quad_x.ptr(),quad_w.ptr());
299 hgT = transpose(hg);
303
304 // Cannot construct the coefficients here since the
305 // derived class is not yet constructed so cannot call
306 // (even indirectly) a virtual method
307 }
308
309 /// Compute the projection of the operator onto the double order polynomials
310 virtual Tensor<Q> rnlp(Level n, Translation lx) const = 0;
311
312 /// @return true if the block of [r^n_l]_ij is expected to be small
313 virtual bool issmall(Level n, Translation lx) const = 0;
314
315 /// @return true if the block of [r^n_l]_ij is expected to be small
316 /// @note unlike issmall(), this handles periodicity and range restriction
317 bool get_issmall(Level n, Translation lx) const {
318 // issmall modufulated by range restriction
319 auto is_small = [this,n](Translation l) {
320 if (!range_restricted())
321 return issmall(n, l);
322 else {
323 // [r^n_l]_ij = superposition of [r^n_l]_p and [r^n_l-1]_p
324 // so rnlij is out of the range only if both rnlp contributions are
325 const auto closest_l = l<=0 ? l : l-1;
326 return rnlp_is_zero(n, closest_l) || issmall(n, l);
327 }
328 };
329
330 // handle lattice summation, if needed
331 if (lattice_summed()) {
332 const Translation twon = Translation(1) << n;
333 for (int R = -maxR; R <= maxR; ++R) {
334 if (!is_small(R * twon + lx))
335 return false;
336 }
337 return true;
338 } else { // !lattice_summed
339 return is_small(lx);
340 }
341 }
342
343 /// @return true if `[r^n_l]` is zero due to range restriction
344 bool rnlp_is_zero(Level n, Translation l) const {
345 bool result = false;
346 if (range_restricted()) {
347 if (n == 0) {
348 // result = l > 0 || l < -1;
349 if (l >= 0)
350 result = Translation(range.iextent_x2()) <= 2*l;
351 else
352 result = Translation(range.iextent_x2()) <= 2*(-l-1);
353 } else { // n > 0
354 if (l >= 0)
355 result = (1 << (n - 1)) * Translation(range.iextent_x2()) <= l;
356 else
357 result = ((1 << (n - 1)) * Translation(range.iextent_x2())) <= (-l-1);
358 }
359 }
360 return result;
361 }
362
363 /// Returns the level for projection
364 virtual Level natural_level() const {return 13;}
365
366 /// Computes the transition matrix elements for the convolution for n,l
367
368 /// Returns the tensor
369 /// \code
370 /// r(i,j) = int(K(x-y) phi[n0](x) phi[nl](y), x=0..1, y=0..1)
371 /// \endcode
372 /// This is computed from the matrix elements over the correlation
373 /// function which in turn are computed from the matrix elements
374 /// over the double order legendre polynomials.
375 /// \note if `this->range_restricted()==true`, `θ(D/2 - |x-y|) K(x-y)` is used as the kernel
376 const Tensor<Q>& rnlij(Level n, Translation lx, bool do_transpose=false) const {
377 const Tensor<Q>* p=rnlij_cache.getptr(n,lx);
378 if (p) return *p;
379
380 // PROFILE_MEMBER_FUNC(Convolution1D); // Too fine grain for routine profiling
381
382 long twok = 2*k;
383 Tensor<Q> R(2*twok);
384 R(Slice(0,twok-1)) = get_rnlp(n,lx-1);
385 R(Slice(twok,2*twok-1)) = get_rnlp(n,lx);
386
387 R.scale(pow(0.5,0.5*n));
388 R = inner(c,R);
389 if (do_transpose) R = transpose(R);
390 rnlij_cache.set(n,lx,R);
391 return *rnlij_cache.getptr(n,lx);
392 };
393
394
395 /// Returns a pointer to the cached modified make_nonstandard form of the operator
396
397 /// @param[in] op_key holds the scale and the source and target translations
398 /// @return a pointer to the cached modified make_nonstandard form of the operator
400
401 const Level& n=op_key.level();
402 const Translation& sx=op_key.translation()[0]; // source translation
403 const Translation& tx=op_key.translation()[1]; // target translation
404 const Translation lx=tx-sx; // displacement
405 const Translation s_off=sx%2;
406 const Translation t_off=tx%2;
407
408 // we cache translation and source offset
411 if (p) return p;
412
413 // for paranoid me
414 MADNESS_ASSERT(sx>=0 and tx>=0);
415
416
417 Tensor<Q> R, T, Rm;
418// if (!get_issmall(n, lx)) {
419// print("no issmall", lx, source, n);
420
421 const Translation lx_half = tx/2 - sx/2;
422 const Slice s0(0,k-1), s1(k,2*k-1);
423// print("sx, tx",lx,lx_half,sx, tx,"off",s_off,t_off);
424
425 // this is the operator matrix in its actual level
426 R = rnlij(n,lx);
427
428 // this is the upsampled operator matrix
429 Rm = Tensor<Q>(2*k,2*k);
430 if (n>0) Rm(s0,s0)=rnlij(n-1,lx_half);
431 {
432 // PROFILE_BLOCK(Convolution1D_nstran); // Too fine grain for routine profiling
433 Rm = transform(Rm,hg);
434 }
435
436 {
437 // PROFILE_BLOCK(Convolution1D_nscopy); // Too fine grain for routine profiling
438 T=Tensor<Q>(k,k);
439 if (t_off==0 and s_off==0) T=copy(Rm(s0,s0));
440 if (t_off==0 and s_off==1) T=copy(Rm(s0,s1));
441 if (t_off==1 and s_off==0) T=copy(Rm(s1,s0));
442 if (t_off==1 and s_off==1) T=copy(Rm(s1,s1));
443// if (t_off==0 and s_off==0) T=copy(Rm(s0,s0));
444// if (t_off==1 and s_off==0) T=copy(Rm(s0,s1));
445// if (t_off==0 and s_off==1) T=copy(Rm(s1,s0));
446// if (t_off==1 and s_off==1) T=copy(Rm(s1,s1));
447 }
448
449 {
450 // PROFILE_BLOCK(Convolution1D_trans); // Too fine grain for routine profiling
451
452 Tensor<Q> RT(k,k), TT(k,k);
453 fast_transpose(k,k,R.ptr(), RT.ptr());
454 fast_transpose(k,k,T.ptr(), TT.ptr());
455 R = RT;
456 T = TT;
457 }
458
459// }
460
462 return mod_ns_cache.getptr(cache_key);
463 }
464
465 /// Returns a pointer to the cached make_nonstandard form of the operator
467 const ConvolutionData1D<Q>* p = ns_cache.getptr(n,lx);
468 if (p) return p;
469
470 // PROFILE_MEMBER_FUNC(Convolution1D); // Too fine grain for routine profiling
471
472 Tensor<Q> R, T;
473 if (!get_issmall(n, lx)) {
474 Translation lx2 = lx*2;
475#if 0 // UNUSED VARIABLES
476 Slice s0(0,k-1), s1(k,2*k-1);
477#endif
478 const Tensor<Q> r0 = rnlij(n+1,lx2);
479 const Tensor<Q> rp = rnlij(n+1,lx2+1);
480 const Tensor<Q> rm = rnlij(n+1,lx2-1);
481
482 R = Tensor<Q>(2*k,2*k);
483
484 // this does not match Eq 18 of the MRAQC appendix .. parity off??
485 // either rnlij dox swap bra/ket or hgT include extra parity phases
486// R(s0,s0) = r0;
487// R(s1,s1) = r0;
488// R(s1,s0) = rp;
489// R(s0,s1) = rm;
490
491 {
492 // PROFILE_BLOCK(Convolution1D_nscopy); // Too fine grain for routine profiling
493 copy_2d_patch(R.ptr(), 2*k, r0.ptr(), k, k, k);
494 copy_2d_patch(R.ptr()+2*k*k + k, 2*k, r0.ptr(), k, k, k);
495 copy_2d_patch(R.ptr()+2*k*k, 2*k, rp.ptr(), k, k, k);
496 copy_2d_patch(R.ptr() + k, 2*k, rm.ptr(), k, k, k);
497 }
498
499 //print("R ", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());
500
501
502 {
503 // PROFILE_BLOCK(Convolution1D_nstran); // Too fine grain for routine profiling
504 R = transform(R,hgT);
505 }
506
507 //print("RX", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());
508
509 {
510 // PROFILE_BLOCK(Convolution1D_trans); // Too fine grain for routine profiling
511
512 Tensor<Q> RT(2*k,2*k);
513 fast_transpose(2*k, 2*k, R.ptr(), RT.ptr());
514 R = RT;
515
516 //print("RT", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());
517
518 //T = copy(R(s0,s0));
519 T = Tensor<Q>(k,k);
520 copy_2d_patch(T.ptr(), k, R.ptr(), 2*k, k, k);
521 }
522
523 //print("NS", n, lx, R.normf(), T.normf());
524 }
525
526 ns_cache.set(n,lx,ConvolutionData1D<Q>(R,T));
527
528 return ns_cache.getptr(n,lx);
529 };
530
531 Q phase(double R) const {
532 if constexpr (std::is_arithmetic_v<Q>)
533 return 1;
534 else
535 return exp(Q(0.0,bloch_k*R));
536 }
537
538
539 const Tensor<Q>& get_rnlp(Level n, Translation lx) const {
540 const Tensor<Q>* p=rnlp_cache.getptr(n,lx);
541 if (p) return *p;
542
543 // PROFILE_MEMBER_FUNC(Convolution1D); // Too fine grain for routine profiling
544
545 long twok = 2*k;
546 Tensor<Q> r;
547
548 if (get_issmall(n, lx)) {
549 r = Tensor<Q>(twok);
550 }
551 else if (n < natural_level()) {
552 Tensor<Q> R(2*twok);
553 R(Slice(0,twok-1)) = get_rnlp(n+1,2*lx);
554 R(Slice(twok,2*twok-1)) = get_rnlp(n+1,2*lx+1);
555
556 R = transform(R, hgT2k);
557 r = copy(R(Slice(0,twok-1)));
558 }
559 else {
560 // PROFILE_BLOCK(Convolution1Drnlp); // Too fine grain for routine profiling
561
562 if (lattice_summed()) {
563 Translation twon = Translation(1)<<n;
564 r = Tensor<Q>(2*k);
565 for (int R=-maxR; R<=maxR; ++R) {
566 r.gaxpy(1.0, rnlp(n,R*twon+lx), phase(R));
567 }
568 }
569 else {
570 r = rnlp(n, lx);
571 }
572 }
573
574 rnlp_cache.set(n, lx, r);
575 //print(" SET rnlp", n, lx, r);
576 return *rnlp_cache.getptr(n,lx);
577 }
578 };
579
580 /// Array of 1D convolutions (one / dimension)
581
582 /// data for 1 term and all dimensions
583 template <typename Q, int NDIM>
585 std::array<std::shared_ptr<Convolution1D<Q> >, NDIM> ops;
587
588 public:
589 ConvolutionND() : fac(1.0) {}
590
591 ConvolutionND(const ConvolutionND& other) : fac(other.fac)
592 {
593 ops = other.ops;
594 }
595
596 ConvolutionND(std::shared_ptr<Convolution1D<Q> > op, Q fac=1.0) : fac(fac)
597 {
598 std::fill(ops.begin(), ops.end(), op);
599 }
600
601 void setop(int dim, const std::shared_ptr<Convolution1D<Q> >& op) {
602 ops[dim] = op;
603 }
604
605 std::shared_ptr<Convolution1D<Q> > getop(int dim) const {
606 return ops[dim];
607 }
608
609 void setfac(Q value) {
610 fac = value;
611 }
612
613 Q getfac() const {
614 return fac;
615 }
616
617 /// @return whether lattice sum is performed along each axis
619 array_of_bools<NDIM> result(false);
620 for (int d = 0; d != NDIM; ++d) {
622 result[d] = ops[d]->lattice_summed();
623 }
624 return result;
625 }
626 };
627
628 // To test generic convolution by comparing with GaussianConvolution1D
629 template <typename Q>
631 private:
633 double exponent;
634 int m;
636
637 public:
638 // coeff * exp(-exponent*x^2) * x^m
641 natlev(Level(0.5*log(exponent)/log(2.0)+1)) {}
642
643 Q operator()(double x) const {
644 Q ee = coeff*exp(-exponent*x*x);
645 for (int mm=0; mm<m; ++mm) ee *= x;
646 return ee;
647 }
648 Level natural_level() const {return natlev;}
649 };
650
651
652 /// Generic 1D convolution using brute force (i.e., slow) adaptive quadrature for rnlp
653
654 /// Calls op(x) with x in *simulation coordinates* to evaluate the function.
655 template <typename Q, typename opT>
657 private:
659 long maxl; ///< At natural level is l beyond which operator is zero
660 public:
661
663
664 GenericConvolution1D(int k, const opT& op, int maxR, double bloch_k = 0.0)
665 : Convolution1D<Q>(k, 20, maxR, bloch_k), op(op), maxl(LONG_MAX-1) {
666 // PROFILE_MEMBER_FUNC(GenericConvolution1D); // Too fine grain for routine profiling
667
668 // For efficiency carefully compute outwards at the "natural" level
669 // until several successive boxes are determined to be zero. This
670 // then defines the future range of the operator and also serves
671 // to precompute the values used in the rnlp cache.
672
674 int nzero = 0;
675 for (Translation lx=0; lx<(1L<<natl); ++lx) {
676 const Tensor<Q>& rp = this->get_rnlp(natl, lx);
677 const Tensor<Q>& rm = this->get_rnlp(natl,-lx);
678 if (rp.normf()<1e-12 && rm.normf()<1e-12) ++nzero;
679 if (nzero == 3) {
680 maxl = lx-2;
681 break;
682 }
683 }
684 }
685
686 virtual Level natural_level() const final {return op.natural_level();}
687
688 struct Shmoo {
693
696
697 returnT operator()(double x) const {
698 int twok = q.k*2;
699 double fac = std::pow(0.5,n);
700 double phix[twok];
702 Q f = q.op(fac*x)*sqrt(fac);
703 returnT v(twok);
704 for (long p=0; p<twok; ++p) v(p) += f*phix[p];
705 return v;
706 }
707 };
708
709 Tensor<Q> rnlp(Level n, Translation lx) const final {
710 return adq1(lx, lx+1, Shmoo(n, lx, this), 1e-12,
711 this->npt, this->quad_x.ptr(), this->quad_w.ptr(), 0);
712 }
713
714 bool issmall(Level n, Translation lx) const final {
715 if (lx < 0) lx = 1 - lx;
716 // Always compute contributions to nearest neighbor coupling
717 // ... we are two levels below so 0,1 --> 0,1,2,3 --> 0,...,7
718 if (lx <= 7) return false;
719
720 n = natural_level()-n;
721 if (n >= 0) lx = lx << n;
722 else lx = lx >> n;
723
724 return lx >= maxl;
725 }
726 };
727
728
729 /// 1D convolution with (derivative) Gaussian; coeff and expnt given in *simulation* coordinates [0,1]
730
731 /// Note that the derivative is computed in *simulation* coordinates so
732 /// you must be careful to scale the coefficients correctly.
733 template <typename Q>
735 // Returns range of Gaussian for periodic lattice sum in simulation coords
736 // N.B. for range-restricted kernels lattice summation range may or may not be limited by the kernel range
737 static int maxR(bool periodic, double expnt, const KernelRange& rng = {}) {
738 if (periodic) {
739 // kernel is 1e-16 past this many simulation cells due to range-restriction
740 const int maxR_rng = rng.finite() ? (rng.iextent_x2(1e-16) + 1)/2 : std::numeric_limits<int>::max();
741 // kernel is 1e-16 past this many simulation cells due to decay of Gaussian kernel
742 const int maxR_G = std::max(1, int(sqrt(16.0 * 2.3 / expnt) + 1));
743 return std::min(maxR_rng,maxR_G);
744 }
745 else {
746 return 0;
747 }
748 }
749 public:
750 const Q coeff; ///< Coefficient
751 const double expnt; ///< Exponent
752 const Level natlev; ///< Level to evaluate
753 const int m; ///< Order of derivative (0, 1, or 2 only)
754
755 explicit GaussianConvolution1D(int k, Q coeff, double expnt,
756 int m, bool periodic, double bloch_k = 0.0,
757 KernelRange rng = {})
758 : Convolution1D<Q>(k,k+11,maxR(periodic,expnt,rng),bloch_k, rng)
759 , coeff(coeff)
760 , expnt(expnt)
761 , natlev(Level(0.5*log(expnt)/log(2.0)+1))
762 , m(m)
763 {
764 MADNESS_ASSERT(m>=0 && m<=2);
765 // std::cout << "GC expnt=" << expnt << " coeff=" << coeff << " natlev=" << natlev << " maxR=" << maxR(periodic,expnt) << std::endl;
766 // for (Level n=0; n<5; n++) {
767 // for (Translation l=0; l<(1<<n); l++) {
768 // std::cout << "RNLP " << n << " " << l << " " << this->get_rnlp(n,l).normf() << std::endl;
769 // }
770 // std::cout << std::endl;
771 // }
772 }
773
775
776 virtual Level natural_level() const final {
777 return natlev;
778 }
779
780 /// Compute the projection of the operator onto the double order polynomials
781
782 /// The returned reference is to a cached tensor ... if you want to
783 /// modify it, take a copy first.
784 ///
785 /// Return in \c v[p] \c p=0..2*k-1
786 /// \code
787 /// r(n,l,p) = 2^(-n) * int(K(2^(-n)*(z+l)) * phi(p,z), z=0..1)
788 /// \endcode
789 /// The kernel is coeff*exp(-expnt*z^2)*z^m (with m>0). This is equivalent to
790 /// \code
791 /// r(n,l,p) = 2^(-n*(m+1))*coeff * int( ((d/dz)^m exp(-beta*z^2)) * phi(p,z-l), z=l..l+1)
792 /// \endcode
793 /// where
794 /// \code
795 /// beta = alpha * 2^(-2*n)
796 /// \endcode
797 Tensor<Q> rnlp(Level n, const Translation lx) const final {
798 int twok = 2*this->k;
799 Tensor<Q> v(twok); // Can optimize this away by passing in
801 constexpr bool use_kahan = false; // change to true to use Kahan accumulator
802
803 // integration range is [0,1] ...
804 std::pair<double, double> integration_limits{0,1};
805 // ... unless using hard (step-wise) range restriction
806 if (this->range.finite_hard()) {
807 const auto two_to_nm1 = (1ul << n) * 0.5;
808 if (lx < 0) {
809 integration_limits = std::make_pair(
810 std::min(std::max(-two_to_nm1 * this->range.iextent_x2() - lx, 0.), 1.), 1.);
811 } else {
812 integration_limits = std::make_pair(
813 0., std::max(std::min(two_to_nm1 * this->range.iextent_x2() - lx, 1.), 0.));
814 }
815 // early return if empty integration range (this indicates that
816 // the range restriction makes the kernel zero everywhere in the box)
817 if (integration_limits.first == integration_limits.second) {
818 MADNESS_ASSERT(this->rnlp_is_zero(n, lx));
819 return v;
820 }
821 else {
822 MADNESS_ASSERT(!this->rnlp_is_zero(n, lx));
823 }
824 }
825 // integration range lower bound, upper bound, length
826 const auto x0 = integration_limits.first;
827 const auto x1 = integration_limits.second;
828 const auto L = x1 - x0;
829
830 /* Apply high-order Gauss Legendre onto subintervals
831
832 coeff*int(exp(-beta(x+l)**2) * z^m * phi[p](x),x=0..1);
833
834 The translations internally considered are all +ve, so
835 significant pieces will be on the left. Finish after things
836 become insignificant.
837
838 The resulting coefficients are accurate to about 1e-20.
839 */
840
841 // Rescale expnt & coeff onto level n so integration range
842 // is [l,l+1]
843 Q scaledcoeff = coeff*pow(0.5,0.5*n*(2*m+1));
844
845 // Subdivide interval into nbox boxes of length h
846 // ... estimate appropriate size from the exponent. A
847 // Gaussian with real-part of the exponent beta falls in
848 // magnitude by a factor of 1/e at x=1/sqrt(beta), and by
849 // a factor of e^-49 ~ 5e-22 at x=7/sqrt(beta) (and with
850 // polyn of z^2 it is 1e-20). So, if we use a box of size
851 // 1/sqrt(beta) we will need at most 7 boxes. Incorporate
852 // the coefficient into the screening since it may be
853 // large. We can represent exp(-x^2) over [l,l+1] with a
854 // polynomial of order 21 to a relative
855 // precision of better than machine precision for
856 // l=0,1,2,3 and for l>3 the absolute error is less than
857 // 1e-23. We want to compute matrix elements with
858 // polynomials of order 2*k-1+m, so the total order is
859 // 2*k+20+m, which can be integrated with a quadrature rule
860 // of npt=k+11+(m+1)/2. npt is set in the constructor.
861
862 double fourn = std::pow(4.0,double(n));
863 double beta = expnt * pow(0.25,double(n));
864 double h = 1.0/sqrt(beta); // 2.0*sqrt(0.5/beta);
865 long nbox = long(1.0/h);
866 if (nbox < 1) nbox = 1;
867
868 // corner case: soft range restriction, range boundary within in interval
869 // since the integrand changes rapidly at x=0, 1/2, or 1, need finer
870 // integration
871 // N.B. this should have almost no impact on performance since this
872 // will produce a large number of boxes only when 2^n * sigma << 1
873 if (this->range.finite_soft()) {
874 // range boundary within in interval if 2^{n-1}*range.N() \in [l,l+1] (l>=0) or [-l,-l+1] (l<0)
875 const auto range_edge_in_interval = [&]() {
876 const auto two_to_nm1 = (1ul << n) * 0.5;
877 const auto range_boundary = two_to_nm1 * this->range.N();
878 if (lx >= 0) {
879 return range_boundary >= lx && range_boundary <= lx+1;
880 }
881 else { // lx < 0
882 return range_boundary <= -lx && range_boundary >= -lx-1;
883 }
884 };
885
886 // ensure that box is at least 1/sigma in size
888 nbox = std::max(nbox, static_cast<long>(ceil(1./((1ul << n) * this->range.sigma()))));
889 }
890 }
891
892 h = L/nbox;
893
894 // Find argmax such that h*scaledcoeff*exp(-argmax)=1e-22 ... if
895 // beta*xlo*xlo is already greater than argmax we can neglect this
896 // and subsequent boxes.
897
898 // The derivatives add a factor of expnt^m to the size of
899 // the function at long range.
900 double sch = std::abs(scaledcoeff*h);
901 if (m == 1) sch *= expnt;
902 else if (m == 2) sch *= expnt*expnt;
903 double argmax = std::abs(log(1e-22/sch)); // perhaps should be -log(1e-22/sch) ?
904
905 // to screen need to iterate over boxes in the order of decreasing kernel values
906 const bool left_to_right = lx >= 0;
907 // if going left-to-right, start at left, else at right
908 const double xstartedge = left_to_right ? x0+lx : lx + 1;
909
910 // with oscillatory integrands the heuristic for reducing roundoff
911 // is to sum from large to small, i.e. proceed in same direction as the order of boxes
912 // WARNING: the grid points in quad_{x,w} are in order of decreasing x!
913 // hence decrement grid point indices for left_to_right, increment otherwise
914 const long first_pt = left_to_right ? this->npt-1: 0;
915 const long sentinel_pt = left_to_right ? -1 : this->npt;
916 const auto next_pt = [lx, left_to_right](auto i) { return left_to_right ? i-1 : i+1; };
917
919 double xhi;
920 for (long box=0; box!=nbox; ++box, xlo = (left_to_right ? xhi : xlo-h)) {
921
922 // can ignore this and rest of boxes if the Gaussian has decayed enough at the side of the box closest to the origin
923 xhi=xlo+h;
924 const auto xabs_min = std::min(std::abs(xhi),std::abs(xlo));
925 if (beta*xabs_min*xabs_min > argmax) break;
926
927 for (long i=first_pt; i!=sentinel_pt; i=next_pt(i)) {
928#ifdef IBMXLC
929 double phix[80];
930#else
931 double phix[twok];
932#endif
933 double xx = xlo + h*this->quad_x(i);
934 Q ee = scaledcoeff*exp(-beta*xx*xx)*this->quad_w(i)*h;
935 if (this->range && this->range.finite()) {
936 const auto x = xx * pow(0.5,double(n));
937 ee *= this->range.value(x);
938 }
939
940 // Differentiate as necessary
941 if (m == 1) {
942 ee *= -2.0*expnt*xx;
943 }
944 else if (m == 2) {
945 ee *= (4.0*xx*xx*expnt*expnt - 2.0*expnt*fourn);
946 }
947
949 for (long p=0; p<twok; ++p) {
950 if constexpr (use_kahan)
951 v_accumulator[p] += ee * phix[p];
952 else
953 v(p) += ee * phix[p];
954 }
955 }
956 }
957
958 if constexpr (use_kahan) {
959 for (long p = 0; p < twok; ++p)
960 v(p) = static_cast<Q>(v_accumulator[p]);
961 }
962
963 return v;
964 }
965
966 /// @return true if the block of [r^n_l]_ij is expected to be small
967 bool issmall(Level n, Translation lx) const final {
968 // [r^n_l]_ij = superposition of [r^n_l]_p and [r^n_l-1]_ij
969 // lx>0? the nearest box is lx-1 -> the edge closest to the origin is lx - 1
970 // lx<0? the nearest box is lx -> the edge closest to the origin is lx + 1
971 // lx==0? interactions within same box are never small
972 if (lx == 0)
973 return false;
974 else {
975 const double beta = expnt * pow(0.25,double(n));
976 const double overly_large_beta_r2 = 49.0; // 49 -> 5e-22 69 -> 1e-30
977 if (lx > 0) {
978 const auto ll = lx - 1;
979 return beta * ll * ll > overly_large_beta_r2;
980 } else {
981 const auto ll = lx + 1;
982 return beta * ll * ll > overly_large_beta_r2;
983 }
984 }
985 };
986 };
987
988
989 template <typename Q>
994
995 static std::shared_ptr< GaussianConvolution1D<Q> > get(int k, double expnt, int m, bool periodic,
996 double bloch_k = 0.0,
997 const KernelRange& range = {}) {
998 hashT key = hash_value(expnt);
999 hash_combine(key, k);
1000 hash_combine(key, m);
1001 hash_combine(key, int(periodic));
1002 hash_combine(key, bloch_k);
1003 if (range) hash_combine(key, range);
1004
1006 MADNESS_PRAGMA_CLANG(diagnostic ignored "-Wundefined-var-template")
1007
1008 iterator it = map.find(key);
1009 if (it == map.end()) {
1010 [[maybe_unused]] auto&& [tmpit, inserted] = map.insert(datumT(key, std::make_shared< GaussianConvolution1D<Q> >(k,
1011 Q(sqrt(expnt/constants::pi)),
1012 expnt,
1013 m,
1014 periodic,
1015 bloch_k,
1016 range
1017 )));
1019 it = map.find(key);
1020 //printf("conv1d: making %d %.8e\n",k,expnt);
1021 }
1022 else {
1023 //printf("conv1d: reusing %d %.8e\n",k,expnt);
1024 }
1025 auto& result = it->second;
1026 MADNESS_ASSERT(result->expnt == expnt &&
1027 result->k == k &&
1028 result->m == m &&
1029 result->lattice_summed() == periodic &&
1030 result->range == range &&
1031 result->bloch_k == bloch_k);
1032 return result;
1033
1035
1036 }
1037 };
1038
1039 // instantiated in mra1.cc
1040 template <>
1043
1044 // instantiated in mra1.cc
1045 template <>
1048
1049}
1050
1051#endif // MADNESS_MRA_CONVOLUTION1D_H__INCLUDED
Provides routines for internal use optimized for aligned data.
std::complex< double > double_complex
Definition cfft.h:14
long dim(int i) const
Returns the size of dimension i.
Definition basetensor.h:147
Definition worldhashmap.h:396
Provides the common functionality/interface of all 1D convolutions.
Definition convolution1d.h:258
SimpleCache< Tensor< Q >, 1 > rnlp_cache
Definition convolution1d.h:272
Tensor< double > hgT2k
Definition convolution1d.h:270
virtual Level natural_level() const
Returns the level for projection.
Definition convolution1d.h:364
bool lattice_summed() const
Definition convolution1d.h:277
int maxR
Number of lattice translations for sum.
Definition convolution1d.h:263
Tensor< double > hg
Definition convolution1d.h:269
int k
Wavelet order.
Definition convolution1d.h:261
SimpleCache< ConvolutionData1D< Q >, 1 > ns_cache
Definition convolution1d.h:274
Tensor< double > hgT
Definition convolution1d.h:269
Tensor< double > c
Definition convolution1d.h:268
double bloch_k
k in exp(i k R) Bloch phase factor folded into lattice sum
Definition convolution1d.h:264
const Tensor< Q > & get_rnlp(Level n, Translation lx) const
Definition convolution1d.h:539
bool rnlp_is_zero(Level n, Translation l) const
Definition convolution1d.h:344
SimpleCache< Tensor< Q >, 1 > rnlij_cache
Definition convolution1d.h:273
virtual Tensor< Q > rnlp(Level n, Translation lx) const =0
Compute the projection of the operator onto the double order polynomials.
Convolution1D(int k, int npt, int maxR, double bloch_k=0.0, KernelRange rng={})
Definition convolution1d.h:282
bool range_restricted() const
Definition convolution1d.h:278
int npt
Number of quadrature points (is this used?)
Definition convolution1d.h:262
KernelRange range
if range is nonnull, kernel range limited to to range (in simulation cell units), useful for finite-r...
Definition convolution1d.h:265
virtual ~Convolution1D()
Definition convolution1d.h:280
Tensor< double > quad_w
Definition convolution1d.h:267
Q phase(double R) const
Definition convolution1d.h:531
SimpleCache< ConvolutionData1D< Q >, 2 > mod_ns_cache
Definition convolution1d.h:275
bool get_issmall(Level n, Translation lx) const
Definition convolution1d.h:317
Tensor< double > quad_x
Definition convolution1d.h:266
const Tensor< Q > & rnlij(Level n, Translation lx, bool do_transpose=false) const
Computes the transition matrix elements for the convolution for n,l.
Definition convolution1d.h:376
virtual bool issmall(Level n, Translation lx) const =0
const ConvolutionData1D< Q > * mod_nonstandard(const Key< 2 > &op_key) const
Returns a pointer to the cached modified make_nonstandard form of the operator.
Definition convolution1d.h:399
Q opT
The apply function uses this to infer resultT=opT*inputT.
Definition convolution1d.h:260
const ConvolutionData1D< Q > * nonstandard(Level n, Translation lx) const
Returns a pointer to the cached make_nonstandard form of the operator.
Definition convolution1d.h:466
Array of 1D convolutions (one / dimension)
Definition convolution1d.h:584
std::shared_ptr< Convolution1D< Q > > getop(int dim) const
Definition convolution1d.h:605
void setop(int dim, const std::shared_ptr< Convolution1D< Q > > &op)
Definition convolution1d.h:601
ConvolutionND(const ConvolutionND &other)
Definition convolution1d.h:591
ConvolutionND()
Definition convolution1d.h:589
std::array< std::shared_ptr< Convolution1D< Q > >, NDIM > ops
Definition convolution1d.h:585
Q fac
Definition convolution1d.h:586
ConvolutionND(std::shared_ptr< Convolution1D< Q > > op, Q fac=1.0)
Definition convolution1d.h:596
Q getfac() const
Definition convolution1d.h:613
array_of_bools< NDIM > lattice_summed() const
Definition convolution1d.h:618
void setfac(Q value)
Definition convolution1d.h:609
1D convolution with (derivative) Gaussian; coeff and expnt given in simulation coordinates [0,...
Definition convolution1d.h:734
const int m
Order of derivative (0, 1, or 2 only)
Definition convolution1d.h:753
Tensor< Q > rnlp(Level n, const Translation lx) const final
Compute the projection of the operator onto the double order polynomials.
Definition convolution1d.h:797
GaussianConvolution1D(int k, Q coeff, double expnt, int m, bool periodic, double bloch_k=0.0, KernelRange rng={})
Definition convolution1d.h:755
virtual ~GaussianConvolution1D()
Definition convolution1d.h:774
const double expnt
Exponent.
Definition convolution1d.h:751
const Level natlev
Level to evaluate.
Definition convolution1d.h:752
static int maxR(bool periodic, double expnt, const KernelRange &rng={})
Definition convolution1d.h:737
const Q coeff
Coefficient.
Definition convolution1d.h:750
bool issmall(Level n, Translation lx) const final
Definition convolution1d.h:967
virtual Level natural_level() const final
Returns the level for projection.
Definition convolution1d.h:776
Definition convolution1d.h:630
Q coeff
Definition convolution1d.h:632
int m
Definition convolution1d.h:634
Level natural_level() const
Definition convolution1d.h:648
Level natlev
Definition convolution1d.h:635
GaussianGenericFunctor(Q coeff, double exponent, int m=0)
Definition convolution1d.h:639
Q operator()(double x) const
Definition convolution1d.h:643
double exponent
Definition convolution1d.h:633
Generic 1D convolution using brute force (i.e., slow) adaptive quadrature for rnlp.
Definition convolution1d.h:656
bool issmall(Level n, Translation lx) const final
Definition convolution1d.h:714
GenericConvolution1D(int k, const opT &op, int maxR, double bloch_k=0.0)
Definition convolution1d.h:664
opT op
Definition convolution1d.h:658
Tensor< Q > rnlp(Level n, Translation lx) const final
Compute the projection of the operator onto the double order polynomials.
Definition convolution1d.h:709
virtual Level natural_level() const final
Returns the level for projection.
Definition convolution1d.h:686
long maxl
At natural level is l beyond which operator is zero.
Definition convolution1d.h:659
GenericConvolution1D()
Definition convolution1d.h:662
Definition kernelrange.h:19
bool finite_soft() const
Definition kernelrange.h:118
double sigma() const
Definition kernelrange.h:111
int iextent_x2(double epsilon=extent_default_epsilon) const
Definition kernelrange.h:157
bool finite() const
Definition kernelrange.h:114
unsigned int N() const
Definition kernelrange.h:109
double value(double r) const
Definition kernelrange.h:135
bool finite_hard() const
Definition kernelrange.h:116
Key is the index for a node of the 2^NDIM-tree.
Definition key.h:68
Simplified interface around hash_map to cache stuff for 1D.
Definition simplecache.h:46
A slice defines a sub-range or patch of a dimension.
Definition slice.h:103
A tensor is a multidimensional array.
Definition tensor.h:317
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition tensor.h:1726
TensorTypeData< T >::scalar_type scalar_type
C++ typename of the real type associated with a complex type.
Definition tensor.h:409
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1824
Tensor< T > & gaxpy(T alpha, const Tensor< T > &t, T beta)
Inplace generalized saxpy ... this = this*alpha + other*beta.
Definition tensor.h:1804
A simple, fixed dimension vector.
Definition vector.h:64
syntactic sugar for std::array<bool, N>
Definition array_of_bools.h:19
Defines common mathematical and physical constants.
static const double R
Definition csqrt.cc:46
char * p(char *buf, const char *name, int k, int initial_level, double thresh, int order)
Definition derivatives.cc:72
auto T(World &world, response_space &f) -> response_space
Definition global_functions.cc:34
const double beta
Definition gygi_soltion.cc:62
static const double v
Definition hatom_sf_dirac.cc:20
Tensor< double > op(const Tensor< double > &x)
Definition kain.cc:508
static double pow(const double *a, const double *b)
Definition lda.h:74
#define max(a, b)
Definition lda.h:51
#define final(a, b, c)
Definition lookup3.c:153
#define MADNESS_PRAGMA_CLANG(x)
Definition madness_config.h:200
#define MADNESS_CHECK(condition)
Check a condition — even in a release build the condition is always evaluated so it can have side eff...
Definition madness_exception.h:182
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition madness_exception.h:134
constexpr double pi
Mathematical constant .
Definition constants.h:48
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
bool two_scale_hg(int k, Tensor< double > *hg)
Definition twoscale.cc:151
void aligned_add(long n, double *MADNESS_RESTRICT a, const double *MADNESS_RESTRICT b)
void fast_transpose(long n, long m, const T *a, T *MADNESS_RESTRICT b)
a(n,m) --> b(m,n) ... optimized for smallish matrices
Definition convolution1d.h:70
void legendre_scaling_functions(double x, long k, double *p)
Evaluate the first k Legendre scaling functions.
Definition legendre.cc:85
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
response_space transpose(response_space &f)
Definition basic_operators.cc:10
void hash_combine(hashT &seed, const T &v)
Combine hash values.
Definition worldhash.h:260
int64_t Translation
Definition key.h:56
int Level
Definition key.h:57
static void copy_2d_patch(T *MADNESS_RESTRICT out, long ldout, const T *MADNESS_RESTRICT in, long ldin, long n, long m)
Definition convolution1d.h:60
bool autoc(int k, Tensor< double > *c)
Return the autocorrelation coefficients for scaling functions of given order.
Definition twoscale.cc:234
bool gauss_legendre(int n, double xlo, double xhi, double *x, double *w)
Definition legendre.cc:226
static double pop(std::vector< double > &v)
Definition SCF.cc:113
void svd(const Tensor< T > &a, Tensor< T > &U, Tensor< typename Tensor< T >::scalar_type > &s, Tensor< T > &VT)
Compute the singluar value decomposition of an n-by-m matrix using *gesvd.
Definition lapack.cc:739
NDIM & f
Definition mra.h:2448
std::size_t hashT
The hash value type.
Definition worldhash.h:145
double inner(response_space &a, response_space &b)
Definition response_functions.h:442
void aligned_sub(long n, double *MADNESS_RESTRICT a, const double *MADNESS_RESTRICT b)
funcT::returnT adq1(double lo, double hi, const funcT &func, double thresh, int n, const double *x, const double *w, int level)
Definition adquad.h:64
static XNonlinearSolver< std::vector< Function< T, NDIM > >, T, vector_function_allocator< T, NDIM > > nonlinear_vector_solver(World &world, const long nvec)
Definition nonlinsol.h:284
madness::hashT hash_value(const std::array< T, N > &a)
Hash std::array with madness hash.
Definition array_addons.h:78
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition mra.h:2034
Definition mraimpl.h:50
static long abs(long a)
Definition tensor.h:218
static const double b
Definition nonlinschro.cc:119
static const double d
Definition nonlinschro.cc:121
static const double a
Definition nonlinschro.cc:118
double Q(double a)
Definition relops.cc:20
static const double m
Definition relops.cc:9
static const double L
Definition rk.cc:46
static const long k
Definition rk.cc:44
!!! Note that if Rnormf is zero then ALL of the tensors are empty
Definition convolution1d.h:162
Tensor< Q > T
if NS: R=ns, T=T part of ns; if modified NS: T=\uparrow r^(n-1)
Definition convolution1d.h:165
double Rnorm
Definition convolution1d.h:170
Tensor< typename Tensor< Q >::scalar_type > Rs
Definition convolution1d.h:167
void make_approx(const Tensor< Q > &R, Tensor< Q > &RU, Tensor< typename Tensor< Q >::scalar_type > &Rs, Tensor< Q > &RVT, double &norm)
Definition convolution1d.h:229
double N_up
Definition convolution1d.h:173
double NSnormf
Definition convolution1d.h:170
double Rnormf
Definition convolution1d.h:170
Tensor< Q > TU
Definition convolution1d.h:166
ConvolutionData1D(const Tensor< Q > &R, const Tensor< Q > &T)
Definition convolution1d.h:180
double N_F
the norms according to Beylkin 2008, Eq. (21) ff
Definition convolution1d.h:173
double N_diff
Definition convolution1d.h:173
Tensor< Q > RU
Definition convolution1d.h:166
double Tnorm
Definition convolution1d.h:170
Tensor< typename Tensor< Q >::scalar_type > Ts
hold relative errors, NOT the singular values..
Definition convolution1d.h:167
ConvolutionData1D(const Tensor< Q > &R, const Tensor< Q > &T, const bool modified)
Definition convolution1d.h:208
double Tnormf
Definition convolution1d.h:170
Tensor< Q > R
Definition convolution1d.h:165
Tensor< Q > TVT
SVD approximations to R and T.
Definition convolution1d.h:166
Tensor< Q > RVT
Definition convolution1d.h:166
Definition convolution1d.h:990
static std::shared_ptr< GaussianConvolution1D< Q > > get(int k, double expnt, int m, bool periodic, double bloch_k=0.0, const KernelRange &range={})
Definition convolution1d.h:995
ConcurrentHashMap< hashT, std::shared_ptr< GaussianConvolution1D< Q > > >::iterator iterator
Definition convolution1d.h:992
static ConcurrentHashMap< hashT, std::shared_ptr< GaussianConvolution1D< Q > > > map
Definition convolution1d.h:991
ConcurrentHashMap< hashT, std::shared_ptr< GaussianConvolution1D< Q > > >::datumT datumT
Definition convolution1d.h:993
Definition convolution1d.h:688
returnT operator()(double x) const
Definition convolution1d.h:697
Tensor< Q > returnT
Definition convolution1d.h:689
Translation lx
Definition convolution1d.h:691
Level n
Definition convolution1d.h:690
Shmoo(Level n, Translation lx, const GenericConvolution1D< Q, opT > *q)
Definition convolution1d.h:694
const GenericConvolution1D< Q, opT > & q
Definition convolution1d.h:692
Definition kahan_accumulator.h:14
static const double x0
Definition tdse1d.cc:145
static const double s0
Definition tdse4.cc:83
Defines and implements most of Tensor.
Prototypes for a partial interface from Tensor to LAPACK.
bool is_small(const double &val, const double &eps)
Definition test6.cc:56
double norm(const T i1)
Definition test_cloud.cc:72
void e()
Definition test_sig.cc:75
constexpr std::size_t NDIM
Definition testgconv.cc:54
double h(const coord_1d &r)
Definition testgconv.cc:175
Implement the madness:Vector class, an extension of std::array that supports some mathematical operat...
const double a2
Definition vnucso.cc:86
const double a1
Definition vnucso.cc:85