convolution1d.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id$
*/
#ifndef MADNESS_MRA_CONVOLUTION1D_H__INCLUDED
#define MADNESS_MRA_CONVOLUTION1D_H__INCLUDED

#include <madness/world/array.h>
#include <madness/constants.h>
#include <limits.h>
#include <madness/tensor/tensor.h>
#include <madness/mra/simplecache.h>
#include <madness/mra/adquad.h>
#include <madness/mra/twoscale.h>
#include <madness/tensor/mtxmq.h>
#include <madness/tensor/aligned.h>
#include <madness/tensor/tensor_lapack.h>
#include <algorithm>



namespace madness {

    void aligned_add(long n, double* restrict a, const double* restrict b);
    void aligned_sub(long n, double* restrict a, const double* restrict b);
    void aligned_add(long n, double_complex* restrict a, const double_complex* restrict b);
    void aligned_sub(long n, double_complex* restrict a, const double_complex* restrict b);

    template <typename T>
    static void copy_2d_patch(T* restrict out, long ldout, const T* restrict in, long ldin, long n, long m) {
        for (long i=0; i<n; ++i, out+=ldout, in+=ldin) {
            for (long j=0; j<m; ++j) {
                out[j] = in[j];
            }
        }
    }

    template <typename T>
    inline void fast_transpose(long n, long m, const T* a, T* restrict b) {
        // n will always be k or 2k (k=wavelet order) and m will be anywhere
        // from 2^(NDIM-1) to (2k)^(NDIM-1).

//                  for (long i=0; i<n; ++i)
//                      for (long j=0; j<m; ++j)
//                          b[j*n+i] = a[i*m+j];
//                  return;

        if (n==1 || m==1) {
            long nm=n*m;
            for (long i=0; i<nm; ++i) b[i] = a[i];
            return;
        }

        long n4 = (n>>2)<<2;
        long m4 = m<<2;
        const T* a0 = a;
        for (long i=0; i<n4; i+=4, a0+=m4) {
            const T* a1 = a0+m;
            const T* a2 = a1+m;
            const T* a3 = a2+m;
            T* restrict bi = b+i;
            for (long j=0; j<m; ++j, bi+=n) {
                T tmp0 = a0[j];
                T tmp1 = a1[j];
                T tmp2 = a2[j];
                T tmp3 = a3[j];

                bi[0] = tmp0;
                bi[1] = tmp1;
                bi[2] = tmp2;
                bi[3] = tmp3;
            }
        }

        for (long i=n4; i<n; ++i)
            for (long j=0; j<m; ++j)
                b[j*n+i] = a[i*m+j];

    }


    template <typename T>
    inline T* shrink(long n, long m, long r, const T* a, T* restrict b) {
        T* result = b;
        if (r == 2) {
            for (long i=0; i<n; ++i, a+=m, b+=r) {
                b[0] = a[0];
                b[1] = a[1];
            }
        }
        else if (r == 4) {
            for (long i=0; i<n; ++i, a+=m, b+=r) {
                b[0] = a[0];
                b[1] = a[1];
                b[2] = a[2];
                b[3] = a[3];
            }
        }
        else {
            MADNESS_ASSERT((r&0x1L)==0);
            for (long i=0; i<n; ++i, a+=m, b+=r) {
                for (long j=0; j<r; j+=2) {
                    b[j  ] = a[j  ];
                    b[j+1] = a[j+1];
                }
            }
        }
        return result;
    }


    template <typename Q>
    struct ConvolutionData1D {


        Tensor<Q> R, T;                 
        Tensor<Q> RU, RVT, TU, TVT;     
        Tensor<typename Tensor<Q>::scalar_type> Rs, Ts;     

        // norms for NS form
        double Rnorm, Tnorm, Rnormf, Tnormf, NSnormf;

        // norms for modified NS form
        double N_up, N_diff, N_F;               


        ConvolutionData1D(const Tensor<Q>& R, const Tensor<Q>& T) : R(R), T(T) {
            Rnormf = R.normf();
            // Making the approximations is expensive ... only do it for
            // significant components
            if (Rnormf > 1e-20) {
                Tnormf = T.normf();
                make_approx(T, TU, Ts, TVT, Tnorm);
                make_approx(R, RU, Rs, RVT, Rnorm);
                int k = T.dim(0);

                Tensor<Q> NS = copy(R);
                for (int i=0; i<k; ++i)
                    for (int j=0; j<k; ++j)
                        NS(i,j) = 0.0;
                NSnormf = NS.normf();

            }
            else {
                Rnorm = Tnorm = Rnormf = Tnormf = NSnormf = 0.0;
                N_F = N_up = N_diff = 0.0;
            }
        }

        ConvolutionData1D(const Tensor<Q>& R, const Tensor<Q>& T, const bool modified) : R(R), T(T) {

            // note that R can be small, but T still be large

            Rnormf = R.normf();
            Tnormf = T.normf();
            // Making the approximations is expensive ... only do it for
            // significant components
            if (Rnormf > 1e-20) make_approx(R, RU, Rs, RVT, Rnorm);
            if (Tnormf > 1e-20) make_approx(T, TU, Ts, TVT, Tnorm);

            // norms for modified NS form: follow Beylkin, 2008, Eq. (21) ff
            N_F=Rnormf;
            N_up=Tnormf;
            N_diff=(R-T).normf();
        }



        void make_approx(const Tensor<Q>& R,
                         Tensor<Q>& RU, Tensor<typename Tensor<Q>::scalar_type>& Rs, Tensor<Q>& RVT, double& norm) {
            int n = R.dim(0);
            svd(R, RU, Rs, RVT);
            for (int i=0; i<n; ++i) {
                for (int j=0; j<n; ++j) {
                    RVT(i,j) *= Rs[i];
                }
            }
            for (int i=n-1; i>1; --i) { // Form cumulative sum of norms
                Rs[i-1] += Rs[i];
            }

            norm = Rs[0];
            if (Rs[0]>0.0) { // Turn into relative errors
                double rnorm = 1.0/norm;
                for (int i=0; i<n; ++i) {
                    Rs[i] *= rnorm;
                }
            }
        }
    };


    template <typename Q>
    class Convolution1D {
    public:
        typedef Q opT;  
        int k;          
        int npt;        
        int maxR;       
        Tensor<double> quad_x;
        Tensor<double> quad_w;
        Tensor<double> c;
        Tensor<double> hgT, hg;
        Tensor<double> hgT2k;
        double arg;

        mutable SimpleCache<Tensor<Q>, 1> rnlp_cache;
        mutable SimpleCache<Tensor<Q>, 1> rnlij_cache;
        mutable SimpleCache<ConvolutionData1D<Q>, 1> ns_cache;
        mutable SimpleCache<ConvolutionData1D<Q>, 2> mod_ns_cache;

        virtual ~Convolution1D() {};

        Convolution1D(int k, int npt, int maxR, double arg = 0.0)
                : k(k)
                , npt(npt)
                , maxR(maxR)
                , quad_x(npt)
                , quad_w(npt)
                , arg(arg)
        {

            MADNESS_ASSERT(autoc(k,&c));

            gauss_legendre(npt,0.0,1.0,quad_x.ptr(),quad_w.ptr());
            MADNESS_ASSERT(two_scale_hg(k,&hg));
            hgT = transpose(hg);
            MADNESS_ASSERT(two_scale_hg(2*k,&hgT2k));
            hgT2k = transpose(hgT2k);

            // Cannot construct the coefficients here since the
            // derived class is not yet constructed so cannot call
            // (even indirectly) a virtual method
        }

        virtual Tensor<Q> rnlp(Level n, Translation lx) const = 0;

        virtual bool issmall(Level n, Translation lx) const = 0;

        bool get_issmall(Level n, Translation lx) const {
            if (maxR == 0) {
                return issmall(n, lx);
            }
            else {
                Translation twon = Translation(1)<<n;
                for (int R=-maxR; R<=maxR; ++R) {
                    if (!issmall(n, R*twon+lx)) return false;
                }
                return true;
            }
        }

        //virtual Level natural_level() const {
        //    return 13;
        //}
        virtual Level natural_level() const {return 13;}


        const Tensor<Q>& rnlij(Level n, Translation lx, bool do_transpose=false) const {
            const Tensor<Q>* p=rnlij_cache.getptr(n,lx);
            if (p) return *p;

            PROFILE_MEMBER_FUNC(Convolution1D);

            long twok = 2*k;
            Tensor<Q> R(2*twok);
            R(Slice(0,twok-1)) = get_rnlp(n,lx-1);
            R(Slice(twok,2*twok-1)) = get_rnlp(n,lx);

            R.scale(pow(0.5,0.5*n));
            R = inner(c,R);
            if (do_transpose) R = transpose(R);
            rnlij_cache.set(n,lx,R);
            return *rnlij_cache.getptr(n,lx);
        };



        const ConvolutionData1D<Q>* mod_nonstandard(const Key<2>& op_key) const {

            const Level& n=op_key.level();
            const Translation& sx=op_key.translation()[0];      // source translation
            const Translation& tx=op_key.translation()[1];      // target translation
            const Translation  lx=tx-sx;                        // displacement
            const Translation  s_off=sx%2;
            const Translation  t_off=tx%2;

            // we cache translation and source offset
            const Key<2> cache_key(n,Vector<Translation,2>(vec(lx,s_off)));
            const ConvolutionData1D<Q>* p = mod_ns_cache.getptr(cache_key);
            if (p) return p;

            // for paranoid me
            MADNESS_ASSERT(sx>=0 and tx>=0);


            Tensor<Q> R, T, Rm;
//            if (!get_issmall(n, lx)) {
//                print("no issmall", lx, source, n);

                const Translation lx_half = tx/2 - sx/2;
                const Slice s0(0,k-1), s1(k,2*k-1);
//                print("sx, tx",lx,lx_half,sx, tx,"off",s_off,t_off);

                // this is the operator matrix in its actual level
                R = rnlij(n,lx);

                // this is the upsampled operator matrix
                Rm = Tensor<Q>(2*k,2*k);
                if (n>0) Rm(s0,s0)=rnlij(n-1,lx_half);
                {
                    PROFILE_BLOCK(Convolution1D_nstran);
                    Rm = transform(Rm,hg);
                }

                {
                    PROFILE_BLOCK(Convolution1D_nscopy);
                    T=Tensor<Q>(k,k);
                    if (t_off==0 and s_off==0) T=copy(Rm(s0,s0));
                    if (t_off==0 and s_off==1) T=copy(Rm(s0,s1));
                    if (t_off==1 and s_off==0) T=copy(Rm(s1,s0));
                    if (t_off==1 and s_off==1) T=copy(Rm(s1,s1));
//                    if (t_off==0 and s_off==0) T=copy(Rm(s0,s0));
//                    if (t_off==1 and s_off==0) T=copy(Rm(s0,s1));
//                    if (t_off==0 and s_off==1) T=copy(Rm(s1,s0));
//                    if (t_off==1 and s_off==1) T=copy(Rm(s1,s1));
                }

                {
                    PROFILE_BLOCK(Convolution1D_trans);

                    Tensor<Q> RT(k,k), TT(k,k);
                    fast_transpose(k,k,R.ptr(), RT.ptr());
                    fast_transpose(k,k,T.ptr(), TT.ptr());
                    R = RT;
                    T = TT;
                }

//            }

            mod_ns_cache.set(cache_key,ConvolutionData1D<Q>(R,T,true));
            return mod_ns_cache.getptr(cache_key);
        }

        const ConvolutionData1D<Q>* nonstandard(Level n, Translation lx) const {
            const ConvolutionData1D<Q>* p = ns_cache.getptr(n,lx);
            if (p) return p;

            PROFILE_MEMBER_FUNC(Convolution1D);

            Tensor<Q> R, T;
            if (!get_issmall(n, lx)) {
                Translation lx2 = lx*2;
#if 0 // UNUSED VARIABLES
                Slice s0(0,k-1), s1(k,2*k-1);
#endif
                const Tensor<Q> r0 = rnlij(n+1,lx2);
                const Tensor<Q> rp = rnlij(n+1,lx2+1);
                const Tensor<Q> rm = rnlij(n+1,lx2-1);

                R = Tensor<Q>(2*k,2*k);

//                 R(s0,s0) = r0;
//                 R(s1,s1) = r0;
//                 R(s1,s0) = rp;
//                 R(s0,s1) = rm;

                {
                    PROFILE_BLOCK(Convolution1D_nscopy);
                    copy_2d_patch(R.ptr(),           2*k, r0.ptr(), k, k, k);
                    copy_2d_patch(R.ptr()+2*k*k + k, 2*k, r0.ptr(), k, k, k);
                    copy_2d_patch(R.ptr()+2*k*k,     2*k, rp.ptr(), k, k, k);
                    copy_2d_patch(R.ptr()       + k, 2*k, rm.ptr(), k, k, k);
                }

                //print("R ", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());


                {
                    PROFILE_BLOCK(Convolution1D_nstran);
                    R = transform(R,hgT);
                }

                //print("RX", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());

                {
                    PROFILE_BLOCK(Convolution1D_trans);

                    Tensor<Q> RT(2*k,2*k);
                    fast_transpose(2*k, 2*k, R.ptr(), RT.ptr());
                    R = RT;

                    //print("RT", n, lx, R.normf(), r0.normf(), rp.normf(), rm.normf());

                    //T = copy(R(s0,s0));
                    T = Tensor<Q>(k,k);
                    copy_2d_patch(T.ptr(), k, R.ptr(), 2*k, k, k);
                }

                //print("NS", n, lx, R.normf(), T.normf());
            }

            ns_cache.set(n,lx,ConvolutionData1D<Q>(R,T));

            return ns_cache.getptr(n,lx);
        };

        Q phase(double R) const {
                return 1.0;
        }

        Q phase(double_complex R) const {
                return exp(double_complex(0.0,arg)*R);
        }


        const Tensor<Q>& get_rnlp(Level n, Translation lx) const {
            const Tensor<Q>* p=rnlp_cache.getptr(n,lx);
            if (p) return *p;

            PROFILE_MEMBER_FUNC(Convolution1D);

            long twok = 2*k;
            Tensor<Q> r;

            if (get_issmall(n, lx)) {
                r = Tensor<Q>(twok);
            }
            else if (n < natural_level()) {
                Tensor<Q>  R(2*twok);
                R(Slice(0,twok-1)) = get_rnlp(n+1,2*lx);
                R(Slice(twok,2*twok-1)) = get_rnlp(n+1,2*lx+1);

                R = transform(R, hgT2k);
                r = copy(R(Slice(0,twok-1)));
            }
            else {
                PROFILE_BLOCK(Convolution1Drnlp);

                if (maxR > 0) {
                    Translation twon = Translation(1)<<n;
                    r = Tensor<Q>(2*k);
                    for (int R=-maxR; R<=maxR; ++R) {
                        r.gaxpy(1.0, rnlp(n,R*twon+lx), phase(Q(R)));
                    }
                }
                else {
                    r = rnlp(n, lx);
                }
            }

            rnlp_cache.set(n, lx, r);
            //print("   SET rnlp", n, lx, r);
            return *rnlp_cache.getptr(n,lx);
        }
    };


    template <typename Q, int NDIM>
    class ConvolutionND {
        std::array<std::shared_ptr<Convolution1D<Q> >, NDIM> ops;
        Q fac;

    public:
        ConvolutionND() : fac(1.0) {}

        ConvolutionND(const ConvolutionND& other) : fac(other.fac)
        {
          ops = other.ops;
        }

        ConvolutionND(std::shared_ptr<Convolution1D<Q> > op, Q fac=1.0) : fac(fac)
        {
            std::fill(ops.begin(), ops.end(), op);
        }

        void setop(int dim, const std::shared_ptr<Convolution1D<Q> >& op)  {
            ops[dim] = op;
        }

        std::shared_ptr<Convolution1D<Q> > getop(int dim) const  {
            return ops[dim];
        }

        void setfac(Q value) {
            fac = value;
        }

        Q getfac() const {
            return fac;
        }
    };

    // To test generic convolution by comparing with GaussianConvolution1D
    template <typename Q>
    class GaussianGenericFunctor {
    private:
        Q coeff;
        double exponent;
        int m;
        Level natlev;

    public:
        // coeff * exp(-exponent*x^2) * x^m
        GaussianGenericFunctor(Q coeff, double exponent, int m=0)
            : coeff(coeff), exponent(exponent), m(m),
              natlev(Level(0.5*log(exponent)/log(2.0)+1)) {}

        Q operator()(double x) const {
            Q ee = coeff*exp(-exponent*x*x);
            for (int mm=0; mm<m; ++mm) ee *= x;
            return ee;
        }
        Level natural_level() const {return natlev;}
    };



    template <typename Q, typename opT>
    class GenericConvolution1D : public Convolution1D<Q> {
    private:
        opT op;
        long maxl;    //< At natural level is l beyond which operator is zero
    public:

        GenericConvolution1D() {}

        GenericConvolution1D(int k, const opT& op, int maxR, double arg = 0.0)
            : Convolution1D<Q>(k, 20, maxR, arg), op(op), maxl(LONG_MAX-1) {
            PROFILE_MEMBER_FUNC(GenericConvolution1D);

            // For efficiency carefully compute outwards at the "natural" level
            // until several successive boxes are determined to be zero.  This
            // then defines the future range of the operator and also serves
            // to precompute the values used in the rnlp cache.

            Level natl = natural_level();
            int nzero = 0;
            for (Translation lx=0; lx<(1L<<natl); ++lx) {
                const Tensor<Q>& rp = this->get_rnlp(natl, lx);
                const Tensor<Q>& rm = this->get_rnlp(natl,-lx);
                if (rp.normf()<1e-12 && rm.normf()<1e-12) ++nzero;
                if (nzero == 3) {
                    maxl = lx-2;
                    break;
                }
            }
        }

        virtual Level natural_level() const {return op.natural_level();}

        struct Shmoo {
            typedef Tensor<Q> returnT;
            Level n;
            Translation lx;
            const GenericConvolution1D<Q,opT>& q;

            Shmoo(Level n, Translation lx, const GenericConvolution1D<Q,opT>* q)
                    : n(n), lx(lx), q(*q) {}

            returnT operator()(double x) const {
                int twok = q.k*2;
                double fac = std::pow(0.5,n);
                double phix[twok];
                legendre_scaling_functions(x-lx,twok,phix);
                Q f = q.op(fac*x)*sqrt(fac);
                returnT v(twok);
                for (long p=0; p<twok; ++p) v(p) += f*phix[p];
                return v;
            }
        };

        Tensor<Q> rnlp(Level n, Translation lx) const {
            return adq1(lx, lx+1, Shmoo(n, lx, this), 1e-12,
                        this->npt, this->quad_x.ptr(), this->quad_w.ptr(), 0);
        }

        bool issmall(Level n, Translation lx) const {
            if (lx < 0) lx = 1 - lx;
            // Always compute contributions to nearest neighbor coupling
            // ... we are two levels below so 0,1 --> 0,1,2,3 --> 0,...,7
            if (lx <= 7) return false;

            n = natural_level()-n;
            if (n >= 0) lx = lx << n;
            else lx = lx >> n;

            return lx >= maxl;
        }
    };



    template <typename Q>
    class GaussianConvolution1D : public Convolution1D<Q> {
        // Returns range of Gaussian for periodic lattice sum in simulation coords
        static int maxR(bool periodic, double expnt) {
            if (periodic) {
                return std::max(1,int(sqrt(16.0*2.3/expnt)));
            }
            else {
                return 0;
            }
        }
    public:
        const Q coeff;          
        const double expnt;     
        const Level natlev;     
        const int m;            

        explicit GaussianConvolution1D(int k, Q coeff, double expnt,
                        int m, bool periodic, double arg = 0.0)
            : Convolution1D<Q>(k,k+11,maxR(periodic,expnt),arg)
            , coeff(coeff)
            , expnt(expnt)
            , natlev(Level(0.5*log(expnt)/log(2.0)+1))
            , m(m)
        {
            MADNESS_ASSERT(m>=0 && m<=2);
            // std::cout << "GC expnt=" << expnt << " coeff="  << coeff << " natlev=" << natlev << " maxR=" << maxR(periodic,expnt) << std::endl;
            // for (Level n=0; n<5; n++) {
            //     for (Translation l=0; l<(1<<n); l++) {
            //         std::cout << "RNLP " << n << " " << l << " " << this->get_rnlp(n,l).normf() << std::endl;
            //     }
            //     std::cout << std::endl;
            // }
        }

        virtual ~GaussianConvolution1D() {}

        virtual Level natural_level() const {
            return natlev;
        }


        Tensor<Q> rnlp(Level n, Translation lx) const {
            int twok = 2*this->k;
            Tensor<Q> v(twok);       // Can optimize this away by passing in

            Translation lkeep = lx;
            if (lx<0) lx = -lx-1;

            /* Apply high-order Gauss Legendre onto subintervals

               coeff*int(exp(-beta(x+l)**2) * z^m * phi[p](x),x=0..1);

               The translations internally considered are all +ve, so
               signficant pieces will be on the left.  Finish after things
               become insignificant.

               The resulting coefficients are accurate to about 1e-20.
            */

            // Rescale expnt & coeff onto level n so integration range
            // is [l,l+1]
            Q scaledcoeff = coeff*pow(0.5,0.5*n*(2*m+1));

            // Subdivide interval into nbox boxes of length h
            // ... estimate appropriate size from the exponent.  A
            // Gaussian with real-part of the exponent beta falls in
            // magnitude by a factor of 1/e at x=1/sqrt(beta), and by
            // a factor of e^-49 ~ 5e-22 at x=7/sqrt(beta) (and with
            // polyn of z^2 it is 1e-20).  So, if we use a box of size
            // 1/sqrt(beta) we will need at most 7 boxes.  Incorporate
            // the coefficient into the screening since it may be
            // large.  We can represent exp(-x^2) over [l,l+1] with a
            // polynomial of order 21 to a relative
            // precision of better than machine precision for
            // l=0,1,2,3 and for l>3 the absolute error is less than
            // 1e-23.  We want to compute matrix elements with
            // polynomials of order 2*k-1+m, so the total order is
            // 2*k+20+m, which can be integrated with a quadrature rule
            // of npt=k+11+(m+1)/2.  npt is set in the constructor.

            double fourn = std::pow(4.0,double(n));
            double beta = expnt * pow(0.25,double(n));
            double h = 1.0/sqrt(beta);  // 2.0*sqrt(0.5/beta);
            long nbox = long(1.0/h);
            if (nbox < 1) nbox = 1;
            h = 1.0/nbox;

            // Find argmax such that h*scaledcoeff*exp(-argmax)=1e-22 ... if
            // beta*xlo*xlo is already greater than argmax we can neglect this
            // and subsequent boxes.

            // The derivatives add a factor of expnt^m to the size of
            // the function at long range.
            double sch = std::abs(scaledcoeff*h);
            if (m == 1) sch *= expnt;
            else if (m == 2) sch *= expnt*expnt;
            double argmax = std::abs(log(1e-22/sch)); // perhaps should be -log(1e-22/sch) ?

            for (long box=0; box<nbox; ++box) {
                double xlo = box*h + lx;
                if (beta*xlo*xlo > argmax) break;
                for (long i=0; i<this->npt; ++i) {
#ifdef IBMXLC
                    double phix[80];
#else
                    double phix[twok];
#endif
                    double xx = xlo + h*this->quad_x(i);
                    Q ee = scaledcoeff*exp(-beta*xx*xx)*this->quad_w(i)*h;

                    // Differentiate as necessary
                    if (m == 1) {
                        ee *= -2.0*expnt*xx;
                    }
                    else if (m == 2) {
                        ee *= (4.0*xx*xx*expnt*expnt - 2.0*expnt*fourn);
                    }

                    legendre_scaling_functions(xx-lx,twok,phix);
                    for (long p=0; p<twok; ++p) v(p) += ee*phix[p];
                }
            }

            if (lkeep < 0) {
                /* phi[p](1-z) = (-1)^p phi[p](z) */
                if (m == 1)
                    for (long p=0; p<twok; ++p) v(p) = -v(p);
                for (long p=1; p<twok; p+=2) v(p) = -v(p);
            }

            return v;
        };

        bool issmall(Level n, Translation lx) const {
            double beta = expnt * pow(0.25,double(n));
            Translation ll;
            if (lx > 0)
                ll = lx - 1;
            else if (lx < 0)
                ll = -1 - lx;
            else
                ll = 0;

            return (beta*ll*ll > 49.0);      // 49 -> 5e-22     69 -> 1e-30
        };
    };


    template <typename Q>
    struct GaussianConvolution1DCache {
        static ConcurrentHashMap<hashT, std::shared_ptr< GaussianConvolution1D<Q> > > map;
        typedef typename ConcurrentHashMap<hashT, std::shared_ptr< GaussianConvolution1D<Q> > >::iterator iterator;
        typedef typename ConcurrentHashMap<hashT, std::shared_ptr< GaussianConvolution1D<Q> > >::datumT datumT;

        static std::shared_ptr< GaussianConvolution1D<Q> > get(int k, double expnt, int m, bool periodic) {
            hashT key = hash_value(expnt);
            hash_combine(key, k);
            hash_combine(key, m);
            hash_combine(key, int(periodic));
            iterator it = map.find(key);
            if (it == map.end()) {
                map.insert(datumT(key, std::shared_ptr< GaussianConvolution1D<Q> >(new GaussianConvolution1D<Q>(k,
                                                                                                          Q(sqrt(expnt/constants::pi)),
                                                                                                          expnt,
                                                                                                          m,
                                                                                                          periodic
                                                                                                          ))));
                it = map.find(key);
                //printf("conv1d: making  %d %.8e\n",k,expnt);
            }
            else {
                //printf("conv1d: reusing %d %.8e\n",k,expnt);
            }
            return it->second;
        }
    };
}

#endif // MADNESS_MRA_CONVOLUTION1D_H__INCLUDED