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 $Id$
33 */
34
35 /// \file mra/testgconv.cc
36 /// \brief Test convolution with Gaussian * polyn
37
38 #include <madness/mra/mra.h>
39 #include <madness/mra/operator.h>
40 #include <madness/constants.h>
41
42 using namespace madness;
43
44 static const int k = 10;
45 static const double thresh = 1e-6;
46 static const double L = 17;
47
48 // exp(-r^2) / sqrt(pi) = normalized gaussian
g(const coord_1d & r)49 double g(const coord_1d& r) {
50 static const double fac = 1.0/sqrt(constants::pi);
51 return exp(-r[0]*r[0]) * fac;
52 }
53
54
55 // g' = exp(-r^2) / sqrt(pi) = derivative of a normalized gaussian
gprime(const coord_1d & r)56 double gprime(const coord_1d& r) {
57 static const double fac = 1.0/sqrt(constants::pi);
58 return -2.* r[0] * exp(-r[0]*r[0]) * fac;
59 }
60
61 // conv g()
convg(const coord_1d & r)62 double convg(const coord_1d& r) {
63 static const double fac = 1.0/sqrt(2.0*constants::pi);
64 return exp(-r[0]*r[0]*0.5) * fac;
65 }
66
67 // sqrt(8)*x*exp(-x^2)
h(const coord_1d & r)68 double h(const coord_1d& r) {
69 static const double fac = sqrt(8.0);
70 return exp(-r[0]*r[0]) * fac * r[0];
71 }
72
73 // g() conv h() == h() conv g()
gconvh(const coord_1d & r)74 double gconvh(const coord_1d& r) {
75 return exp(-0.5*r[0]*r[0]) * r[0];
76 }
77
78 // D(conv) (g())
conv_prime_g(const coord_1d & r)79 double conv_prime_g(const coord_1d& r) {
80 static const double fac = 1.0/sqrt(2.0*constants::pi);
81 return -exp(-0.5*r[0]*r[0]) * r[0] *fac;
82 }
83
84
test_gconv(World & world)85 int test_gconv(World& world) {
86 coord_1d origin(0.0), lo(-L), hi(L);
87 double width = 2.0*L;
88 int success=0;
89
90 if (world.rank() == 0) print("Test gconv operation");
91
92 const real_function_1d f = real_factory_1d(world).f(g);
93 double error=f.trace()-1.0;
94 print("error in integral(g) ", error);
95 if (error>FunctionDefaults<1>::get_thresh()) success++;
96 print("success 0 ", success);
97
98 // convolve with a normalized Gaussian kernel
99 std::vector< std::shared_ptr< Convolution1D<double> > > ops(1);
100 ops[0].reset(new GaussianConvolution1D<double>(k, width/sqrt(constants::pi),
101 width*width, 0, false));
102 real_convolution_1d op(world, ops);
103
104 real_function_1d opf = op(f);
105 error=opf.trace()-1.0;
106 print("error in integral(op(g)) ", error);
107 if (error>FunctionDefaults<1>::get_thresh()) success++;
108 print("success 1 ", success);
109
110 real_function_1d exact = real_factory_1d(world).f(convg);
111 print("norm2(g conv g - exact)", (opf-exact).norm2());
112 error=(opf-exact).norm2();
113 if (error>FunctionDefaults<1>::get_thresh()) success++;
114 print("success 2 ", success);
115
116 real_function_1d q = real_factory_1d(world).f(h);
117 error=q.trace();
118 print("error in integral(h) ", error);
119 if (error>FunctionDefaults<1>::get_thresh()) success++;
120 print("success 3 ", success);
121
122 error=q.norm2() - sqrt(sqrt(2.0*constants::pi));
123 print("error in norm2(h)", error);
124 if (error>FunctionDefaults<1>::get_thresh()) success++;
125 print("success 4 ", success);
126
127 real_function_1d opq = op(q);
128 exact = real_factory_1d(world).f(gconvh);
129 error=(opq-exact).norm2();
130 print("norm2(g conv h - exact)", error);
131 if (error>FunctionDefaults<1>::get_thresh()) success++;
132 print("success 5 ", success);
133
134
135
136 // test the convolution with a derivative Gaussian:
137 // result(y) = \int g'(x-y) f(x) = \int g(x-y) f'(x)
138 // where we use
139 // f(x) = exp(-x^2) / sqrt(pi)
140 // f'(x) = -2 x exp(-x^2) / sqrt(pi)
141 // result(y) = -y exp(-y/2) / sqrt(2 pi)
142
143 // the derivative Gaussian convolution kernel:
144 // note the scaling of the coeffs because the derivative operator brings
145 // down the scaling factor of the exponent
146 ops[0].reset(new GaussianConvolution1D<double>(k, 1.0/sqrt(constants::pi),
147 width*width, 1, false));
148
149 real_convolution_1d oph(world, ops);
150
151 // this is the result hardwired
152 const real_function_1d convpg=real_factory_1d(world).f(conv_prime_g);
153
154 // apply the derivative Gaussian on f
155 opq = oph(f);
156
157 // apply the Gaussian on the derivative of f
158 const real_function_1d fp=real_factory_1d(world).f(gprime);
159 real_function_1d opfp=op(fp);
160
161 // the error
162 const double error1=(opq-convpg).norm2();
163 const double error2=(opfp-convpg).norm2();
164 print("norm2(conv' g - exact)", error1);
165 print("norm2(conv g' - exact)", error2);
166 if (error1>FunctionDefaults<1>::get_thresh()) success++;
167 print("success 6a ", success);
168 if (error2>FunctionDefaults<1>::get_thresh()) success++;
169 print("success 6b ", success);
170
171 plot_line("opf.dat", 1001, lo, hi, q, opq, convpg);
172
173 world.gop.fence();
174 return success;
175 }
176
177
main(int argc,char ** argv)178 int main(int argc, char**argv) {
179 initialize(argc,argv);
180 World world(SafeMPI::COMM_WORLD);
181
182 int success=0;
183 try {
184 startup(world,argc,argv);
185
186 FunctionDefaults<1>::set_cubic_cell(-L,L);
187 FunctionDefaults<1>::set_k(k);
188 FunctionDefaults<1>::set_thresh(thresh);
189 FunctionDefaults<1>::set_truncate_mode(1);
190 FunctionDefaults<1>::set_initial_level(5);
191 if (world.rank()==0) {
192 print(" threshold ", thresh);
193 print(" polynomial ", k,"\n");
194 }
195 success+=test_gconv(world);
196
197 }
198 catch (const SafeMPI::Exception& e) {
199 print(e);
200 error("caught an MPI exception");
201 }
202 catch (const madness::MadnessException& e) {
203 print(e);
204 error("caught a MADNESS exception");
205 }
206 catch (const madness::TensorException& e) {
207 print(e);
208 error("caught a Tensor exception");
209 }
210 catch (char* s) {
211 print(s);
212 error("caught a c-string exception");
213 }
214 catch (const char* s) {
215 print(s);
216 error("caught a c-string exception");
217 }
218 catch (const std::string& s) {
219 print(s);
220 error("caught a string (class) exception");
221 }
222 catch (const std::exception& e) {
223 print(e.what());
224 error("caught an STL exception");
225 }
226 catch (...) {
227 error("caught unhandled exception");
228 }
229
230 world.gop.fence();
231 finalize();
232
233 return success;
234 }
235
236