1 /*
2  * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
3  * All rights reserved. Use is subject to license terms. See LICENSE.TXT
4  * For use in publications, see ACKNOWLEDGE.TXT
5  */
6 
7 #include "gtest/gtest.h"
8 #include "alps/gf/mesh.hpp"
9 #include "alps/gf/grid.hpp"
10 #include "alps/gf/piecewise_polynomial.hpp"
11 #include "gf_test.hpp"
12 
13 #include <alps/testing/unique_file.hpp>
14 
TEST(Mesh,RealFrequencyLoadSave)15 TEST(Mesh, RealFrequencyLoadSave) {
16     alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
17     const std::string&  filename = ufile.name();
18     namespace g=alps::gf;
19     g::grid::linear_real_frequency_grid grid(-3,3,100);
20     g::real_frequency_mesh mesh1(grid);
21     g::real_frequency_mesh mesh2;
22     {
23         alps::hdf5::archive oar(filename,"w");
24         mesh1.save(oar,"/gf");
25     }
26     {
27         alps::hdf5::archive iar(filename);
28         g::real_frequency_mesh mesh;
29         mesh2.load(iar,"/gf");
30     }
31     EXPECT_EQ(mesh1, mesh2);
32 }
33 
TEST(Mesh,RealFrequencyLoadSaveStream)34 TEST(Mesh, RealFrequencyLoadSaveStream) {
35     alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
36     const std::string&  filename = ufile.name();
37     namespace g=alps::gf;
38     g::grid::linear_real_frequency_grid grid(-3,3,100);
39     g::real_frequency_mesh mesh1(grid);
40     g::real_frequency_mesh mesh2;
41     {
42         alps::hdf5::archive oar(filename,"w");
43         oar["/gf"] << mesh1;
44     }
45     {
46         alps::hdf5::archive iar(filename);
47         g::real_frequency_mesh mesh;
48         iar["/gf"] >> mesh2;
49     }
50     EXPECT_EQ(mesh1, mesh2);
51 }
52 
TEST(Mesh,RealFrequencyMeshQuadric)53 TEST(Mesh, RealFrequencyMeshQuadric) {
54     double spread = 5.0;
55     int nfreq = 41;
56     alps::gf::grid::quadratic_real_frequency_grid grid(spread, nfreq);
57     alps::gf::real_frequency_mesh mesh1(grid);
58     EXPECT_EQ(mesh1.extent(), nfreq);
59 }
60 
TEST(Mesh,RealFrequencyMeshLogarithmic)61 TEST(Mesh, RealFrequencyMeshLogarithmic) {
62     double tmax = 5, tmin = 0.001;
63     int nfreq = 41;
64     alps::gf::grid::logarithmic_real_frequency_grid grid(tmin, tmax, nfreq);
65     alps::gf::real_frequency_mesh mesh1(grid);
66     EXPECT_EQ(mesh1.extent(), nfreq);
67 }
68 
TEST(Mesh,RealFrequencyMeshLinear)69 TEST(Mesh, RealFrequencyMeshLinear) {
70     double Emin = -5;
71     double Emax = 5;
72     int nfreq = 20;
73     alps::gf::grid::linear_real_frequency_grid grid(Emin, Emax, nfreq);
74     alps::gf::real_frequency_mesh mesh1(grid);
75     EXPECT_EQ(mesh1.extent(), nfreq);
76 }
77 
TEST(Mesh,BosonicMatsubara)78 TEST(Mesh, BosonicMatsubara) {
79     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh1(5.0, 20, alps::gf::statistics::BOSONIC);
80     EXPECT_EQ(mesh1.statistics(), alps::gf::statistics::BOSONIC);
81 }
82 
TEST(Mesh,SwapMatsubara)83 TEST(Mesh,SwapMatsubara) {
84     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh1(5.0, 20);
85     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh2(7.0, 40);
86     mesh1.swap(mesh2);
87 
88     EXPECT_EQ(mesh1.beta(), 7.0);
89     EXPECT_EQ(mesh2.beta(), 5.0);
90 }
91 
TEST(Mesh,SwapMatsubaraCheckStatistics)92 TEST(Mesh,SwapMatsubaraCheckStatistics) {
93     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh1(5.0, 20, alps::gf::statistics::BOSONIC);
94     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh2(7.0, 40, alps::gf::statistics::FERMIONIC);
95 
96     EXPECT_ANY_THROW(mesh1.swap(mesh2));
97 }
98 
TEST(Mesh,CompareMatsubara)99 TEST(Mesh,CompareMatsubara) {
100   alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh1(5.0, 20);
101   alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh2(5.0, 20);
102   alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh3(4.0, 20);
103 
104   EXPECT_TRUE(mesh1==mesh2);
105   EXPECT_TRUE(mesh1!=mesh3);
106 
107   EXPECT_FALSE(mesh1==mesh3);
108   EXPECT_FALSE(mesh1!=mesh2);
109 
110 }
111 
TEST(Mesh,MatsubaraPositiveLoadSave)112 TEST(Mesh, MatsubaraPositiveLoadSave) {
113   alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
114   const std::string&  filename = ufile.name();
115   namespace g=alps::gf;
116   double beta = 100;
117   int nfr = 10;
118   typedef g::matsubara_positive_mesh::index_type index;
119   g::matsubara_positive_mesh mesh1(beta, nfr);
120   g::matsubara_positive_mesh mesh2;
121   {
122     alps::hdf5::archive oar(filename,"w");
123     mesh1.save(oar,"/gf");
124   }
125   {
126     alps::hdf5::archive iar(filename);
127     mesh2.load(iar,"/gf");
128   }
129   EXPECT_EQ(mesh1, mesh2);
130   for(index i(0); i<mesh2.extent(); ++i) {
131     EXPECT_EQ(mesh1(i), mesh2(i));
132   }
133 }
134 
TEST(Mesh,MatsubaraPNLoadSave)135 TEST(Mesh, MatsubaraPNLoadSave) {
136   alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
137   const std::string&  filename = ufile.name();
138   namespace g=alps::gf;
139   typedef g::matsubara_pn_mesh::index_type index;
140   double beta = 100;
141   int nfr = 10;
142   g::matsubara_pn_mesh mesh1(beta, nfr);
143   g::matsubara_pn_mesh mesh2;
144   {
145     alps::hdf5::archive oar(filename,"w");
146     mesh1.save(oar,"/gf");
147   }
148   {
149     alps::hdf5::archive iar(filename);
150     mesh2.load(iar,"/gf");
151   }
152   EXPECT_EQ(mesh1, mesh2);
153   for(index i(0); i<mesh2.extent(); ++i) {
154     EXPECT_EQ(mesh1(i), mesh2(i));
155   }
156 }
157 
158 
TEST(Mesh,CompareITime)159 TEST(Mesh,CompareITime) {
160   alps::gf::itime_mesh mesh1(5.0, 20);
161   alps::gf::itime_mesh mesh2(5.0, 20);
162   alps::gf::itime_mesh mesh3(4.0, 20);
163 
164   EXPECT_TRUE(mesh1==mesh2);
165   EXPECT_TRUE(mesh1!=mesh3);
166 
167   EXPECT_FALSE(mesh1==mesh3);
168   EXPECT_FALSE(mesh1!=mesh2);
169 }
170 
TEST(Mesh,PowerMeshHasRightSize)171 TEST(Mesh,PowerMeshHasRightSize) {
172   alps::gf::power_mesh mesh1(20, 12, 16);
173   EXPECT_EQ(417u, mesh1.points().size());
174   alps::gf::power_mesh mesh2(20, 12, 256);
175   EXPECT_EQ(6657u, mesh2.points().size());
176 
177 }
TEST(Mesh,ComparePower)178 TEST(Mesh,ComparePower) {
179   alps::gf::power_mesh mesh1(12, 16, 20);
180   alps::gf::power_mesh mesh2(12, 16, 20);
181   alps::gf::power_mesh mesh3(11, 16, 20);
182   alps::gf::power_mesh mesh4(12, 12, 20);
183   alps::gf::power_mesh mesh5(12, 16, 22);
184 
185   EXPECT_TRUE(mesh1==mesh2);
186   EXPECT_TRUE(mesh1!=mesh3);
187   EXPECT_TRUE(mesh1!=mesh4);
188   EXPECT_TRUE(mesh1!=mesh5);
189 
190   EXPECT_FALSE(mesh1==mesh3);
191   EXPECT_FALSE(mesh1!=mesh2);
192 }
193 
TEST(Mesh,CompareMomentum)194 TEST(Mesh,CompareMomentum) {
195   alps::gf::momentum_index_mesh::container_type points1(boost::extents[20][3]);
196   alps::gf::momentum_index_mesh::container_type points2(boost::extents[20][3]);
197   alps::gf::momentum_index_mesh::container_type points3(boost::extents[20][3]);
198   alps::gf::momentum_index_mesh::container_type points4(boost::extents[3][20]);
199   for (std::size_t i=0; i<points1.num_elements(); ++i) {
200     *(points1.origin()+i)=i;
201     *(points2.origin()+i)=i;
202     *(points3.origin()+i)=i+1;
203     *(points4.origin()+i)=i;
204   }
205 
206   alps::gf::momentum_index_mesh mesh1(points1);
207   alps::gf::momentum_index_mesh mesh2(points2);
208   alps::gf::momentum_index_mesh mesh3(points3);
209   alps::gf::momentum_index_mesh mesh4(points4);
210 
211   EXPECT_TRUE(mesh1==mesh2);
212   EXPECT_TRUE(mesh1!=mesh3);
213   EXPECT_TRUE(mesh1!=mesh4);
214 
215   EXPECT_FALSE(mesh1==mesh3);
216   EXPECT_FALSE(mesh1!=mesh2);
217   EXPECT_FALSE(mesh1==mesh4);
218 }
219 
220 
TEST(Mesh,CompareRealSpace)221 TEST(Mesh,CompareRealSpace) {
222   alps::gf::real_space_index_mesh::container_type points1(boost::extents[20][3]);
223   alps::gf::real_space_index_mesh::container_type points2(boost::extents[20][3]);
224   alps::gf::real_space_index_mesh::container_type points3(boost::extents[20][3]);
225   alps::gf::real_space_index_mesh::container_type points4(boost::extents[3][20]);
226   for (std::size_t i=0; i<points1.num_elements(); ++i) {
227     *(points1.origin()+i)=i;
228     *(points2.origin()+i)=i;
229     *(points3.origin()+i)=i+1;
230     *(points4.origin()+i)=i;
231   }
232 
233   alps::gf::real_space_index_mesh mesh1(points1);
234   alps::gf::real_space_index_mesh mesh2(points2);
235   alps::gf::real_space_index_mesh mesh3(points3);
236   alps::gf::real_space_index_mesh mesh4(points4);
237 
238   EXPECT_TRUE(mesh1==mesh2);
239   EXPECT_TRUE(mesh1!=mesh3);
240   EXPECT_TRUE(mesh1!=mesh4);
241 
242   EXPECT_FALSE(mesh1==mesh3);
243   EXPECT_FALSE(mesh1!=mesh2);
244   EXPECT_FALSE(mesh1==mesh4);
245 }
246 
TEST(Mesh,CompareIndex)247 TEST(Mesh,CompareIndex) {
248   alps::gf::index_mesh mesh1(20);
249   alps::gf::index_mesh mesh2(20);
250   alps::gf::index_mesh mesh3(19);
251 
252   EXPECT_TRUE(mesh1==mesh2);
253   EXPECT_TRUE(mesh1!=mesh3);
254 
255   EXPECT_FALSE(mesh1==mesh3);
256   EXPECT_FALSE(mesh1!=mesh2);
257 }
258 
TEST(Mesh,Legendre)259 TEST(Mesh,Legendre) {
260     alps::gf::legendre_mesh mesh_f(5.0, 20, alps::gf::statistics::FERMIONIC);
261     EXPECT_EQ(mesh_f.statistics(), alps::gf::statistics::FERMIONIC);
262 
263     alps::gf::legendre_mesh mesh_b(5.0, 20, alps::gf::statistics::BOSONIC);
264     EXPECT_EQ(mesh_b.statistics(), alps::gf::statistics::BOSONIC);
265 }
266 
TEST(Mesh,CompareLegendre)267 TEST(Mesh,CompareLegendre) {
268     alps::gf::legendre_mesh mesh_f1(5.0, 20, alps::gf::statistics::FERMIONIC);
269     alps::gf::legendre_mesh mesh_f2(1.0, 20, alps::gf::statistics::FERMIONIC);
270     alps::gf::legendre_mesh mesh_f3(5.0, 20, alps::gf::statistics::FERMIONIC);
271 
272     alps::gf::legendre_mesh mesh_b1(5.0, 20, alps::gf::statistics::BOSONIC);
273     alps::gf::legendre_mesh mesh_b2(1.0, 20, alps::gf::statistics::BOSONIC);
274     alps::gf::legendre_mesh mesh_b3(5.0, 20, alps::gf::statistics::BOSONIC);
275 
276     EXPECT_FALSE(mesh_f1==mesh_f2);
277     EXPECT_TRUE(mesh_f1==mesh_f3);
278 
279     EXPECT_FALSE(mesh_b1==mesh_b2);
280     EXPECT_TRUE(mesh_b1==mesh_b3);
281 
282     EXPECT_FALSE(mesh_f1==mesh_b1);
283 }
284 
TEST(Mesh,LegendreLoadSave)285 TEST(Mesh, LegendreLoadSave) {
286     namespace g=alps::gf;
287     alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
288     const std::string&  filename = ufile.name();
289 
290     g::legendre_mesh mesh1(5.0, 20, alps::gf::statistics::FERMIONIC);
291     g::legendre_mesh mesh2;
292     {
293         alps::hdf5::archive oar(filename,"w");
294         mesh1.save(oar,"/gf");
295     }
296     {
297         alps::hdf5::archive iar(filename);
298         g::legendre_mesh mesh;
299         mesh2.load(iar,"/gf");
300     }
301     EXPECT_EQ(mesh1, mesh2);
302 }
303 
TEST(Mesh,LegendreLoadSaveStream)304 TEST(Mesh, LegendreLoadSaveStream) {
305     namespace g=alps::gf;
306     alps::testing::unique_file ufile("gf.h5.", alps::testing::unique_file::REMOVE_NOW);
307     const std::string&  filename = ufile.name();
308 
309     g::legendre_mesh mesh1(5.0, 20, alps::gf::statistics::FERMIONIC);
310     g::legendre_mesh mesh2;
311     {
312         alps::hdf5::archive oar(filename,"w");
313         oar["/gf"]<<mesh1;
314     }
315     {
316         alps::hdf5::archive iar(filename);
317         g::legendre_mesh mesh;
318         iar["/gf"]>>mesh2;
319     }
320     EXPECT_EQ(mesh1, mesh2);
321 }
322 
TEST(Mesh,SwapLegendre)323 TEST(Mesh,SwapLegendre) {
324     alps::gf::legendre_mesh mesh_1(5.0, 20, alps::gf::statistics::FERMIONIC);
325     alps::gf::legendre_mesh mesh_1r(mesh_1);
326     alps::gf::legendre_mesh mesh_2(10.0, 40, alps::gf::statistics::FERMIONIC);
327     alps::gf::legendre_mesh mesh_2r(mesh_2);
328 
329     mesh_1.swap(mesh_2);
330     EXPECT_EQ(mesh_1, mesh_2r);
331     EXPECT_EQ(mesh_2, mesh_1r);
332 }
333 
TEST(Mesh,SwapLegendreDifferentStatistics)334 TEST(Mesh,SwapLegendreDifferentStatistics) {
335     alps::gf::legendre_mesh mesh_1(5.0, 20, alps::gf::statistics::BOSONIC);
336     alps::gf::legendre_mesh mesh_1r(mesh_1);
337     alps::gf::legendre_mesh mesh_2(10.0, 40, alps::gf::statistics::FERMIONIC);
338     alps::gf::legendre_mesh mesh_2r(mesh_2);
339 
340     EXPECT_THROW(mesh_1.swap(mesh_2), std::runtime_error);
341 }
342 
TEST(Mesh,PrintMatsubaraMeshHeader)343 TEST(Mesh,PrintMatsubaraMeshHeader) {
344   double beta=5.;
345   int n=20;
346   {
347     std::stringstream header_line;
348     header_line << "# MATSUBARA mesh: N: "<<n<<" beta: "<<beta<<" statistics: "<<"FERMIONIC"<<" POSITIVE_ONLY"<<std::endl;
349     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_ONLY> mesh1(beta, n);
350     std::stringstream header_line_from_mesh;
351     header_line_from_mesh << mesh1;
352     EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
353   }
354   {
355     std::stringstream header_line;
356     header_line << "# MATSUBARA mesh: N: "<<n<<" beta: "<<beta<<" statistics: "<<"FERMIONIC"<<" POSITIVE_NEGATIVE"<<std::endl;
357     alps::gf::matsubara_mesh<alps::gf::mesh::POSITIVE_NEGATIVE> mesh1(beta, n);
358     std::stringstream header_line_from_mesh;
359     header_line_from_mesh << mesh1;
360     EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
361   }
362 
363 }
364 
TEST(Mesh,PrintImagTimeMeshHeader)365 TEST(Mesh,PrintImagTimeMeshHeader) {
366   double beta=5.;
367   int ntau=200;
368   std::stringstream header_line;
369   header_line << "# IMAGINARY_TIME mesh: N: "<<ntau<<" beta: "<<beta<<" statistics: "<<"FERMIONIC"<<std::endl;
370   alps::gf::itime_mesh mesh1(beta, ntau);
371   std::stringstream header_line_from_mesh;
372   header_line_from_mesh << mesh1;
373   EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
374 }
TEST(Mesh,PrintPowerMeshHeader)375 TEST(Mesh,PrintPowerMeshHeader) {
376   double beta=5.;
377   int power=12;
378   int uniform=16;
379   int ntau=417;
380   std::stringstream header_line;
381   header_line << "# POWER mesh: power: "<<power<<" uniform: "<<uniform<<" N: "<<ntau<<" beta: "<<beta<<" statistics: "<<"FERMIONIC"<<std::endl;
382   alps::gf::power_mesh mesh1(beta, power, uniform);
383   std::stringstream header_line_from_mesh;
384   header_line_from_mesh << mesh1;
385   EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
386 }
387 
TEST(Mesh,PrintMomentumMeshHeader)388 TEST(Mesh,PrintMomentumMeshHeader) {
389   alps::gf::momentum_index_mesh::container_type data=get_data_for_momentum_mesh();
390   std::stringstream header_line;
391   header_line << "# MOMENTUM_INDEX mesh: N: "<<data.shape()[0]<<" dimension: "<<data.shape()[1]<<" points: ";
392   for(std::size_t i=0;i<data.shape()[0];++i){
393     header_line<<"(";
394     for(std::size_t d=0;d<data.shape()[1]-1;++d){ header_line<<data[i][d]<<","; } header_line<<data[i][data.shape()[1]-1]<<") ";
395   }
396   header_line<<std::endl;
397   alps::gf::momentum_index_mesh mesh1(data);
398   std::stringstream header_line_from_mesh;
399   header_line_from_mesh << mesh1;
400   EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
401 }
TEST(Mesh,PrintRealSpaceMeshHeader)402 TEST(Mesh,PrintRealSpaceMeshHeader) {
403   alps::gf::real_space_index_mesh::container_type data=get_data_for_real_space_mesh();
404   std::stringstream header_line;
405   header_line << "# REAL_SPACE_INDEX mesh: N: "<<data.shape()[0]<<" dimension: "<<data.shape()[1]<<" points: ";
406   for(std::size_t i=0;i<data.shape()[0];++i){
407     header_line<<"(";
408     for(std::size_t d=0;d<data.shape()[1]-1;++d){ header_line<<data[i][d]<<","; } header_line<<data[i][data.shape()[1]-1]<<") ";
409   }
410   header_line<<std::endl;
411   alps::gf::real_space_index_mesh mesh1(data);
412   std::stringstream header_line_from_mesh;
413   header_line_from_mesh << mesh1;
414   EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
415 }
TEST(Mesh,PrintIndexMeshHeader)416 TEST(Mesh,PrintIndexMeshHeader) {
417   int n=2;
418   std::stringstream header_line;
419   header_line << "# INDEX mesh: N: "<<n<<std::endl;
420   alps::gf::index_mesh mesh1(2);
421   std::stringstream header_line_from_mesh;
422   header_line_from_mesh << mesh1;
423   EXPECT_EQ(header_line.str(), header_line_from_mesh.str());
424 }
TEST(Mesh,PowerWeightsAddUpToOne)425 TEST(Mesh,PowerWeightsAddUpToOne) {
426   alps::gf::power_mesh mesh1(20, 12, 16);
427   //check that the integration weights add up to 1:
428   double sum=0;
429   for(int i=0;i<mesh1.extent();++i) sum+=mesh1.weights()[i];
430   EXPECT_NEAR(1, sum, 1.e-10);
431 }
432 
TEST(Mesh,SwapNumericalMesh)433 TEST(Mesh,SwapNumericalMesh) {
434     const int n_section = 2, k = 3;
435     const double beta = 100.0;
436     typedef double Scalar;
437     typedef alps::gf::piecewise_polynomial<Scalar> pp_type;
438 
439     std::vector<double> section_edges(n_section+1);
440     section_edges[0] = -1.0;
441     section_edges[1] =  0.0;
442     section_edges[2] =  1.0;
443 
444     boost::multi_array<Scalar,2> coeff(boost::extents[n_section][k+1]);
445     std::fill(coeff.origin(), coeff.origin()+coeff.num_elements(), 0.0);
446 
447     pp_type p(n_section, section_edges, coeff);
448 
449     std::vector<pp_type> basis_functions;
450     basis_functions.push_back(p);
451     basis_functions.push_back(p);
452 
453     std::vector<pp_type> basis_functions2;
454     basis_functions2.push_back(p);
455 
456     alps::gf::numerical_mesh<double> mesh1(beta, basis_functions, alps::gf::statistics::FERMIONIC);
457     alps::gf::numerical_mesh<double> mesh2(beta, basis_functions2, alps::gf::statistics::FERMIONIC);
458     alps::gf::numerical_mesh<double> mesh3(beta, basis_functions2, alps::gf::statistics::BOSONIC);
459 
460     mesh1.swap(mesh2);
461     ASSERT_TRUE(mesh1.extent()==1);
462     ASSERT_TRUE(mesh2.extent()==2);
463 
464     ASSERT_THROW(mesh1.swap(mesh3), std::runtime_error);
465 }
466 
TEST(Mesh,NumericalMeshSave)467 TEST(Mesh,NumericalMeshSave) {
468     const int n_section = 2, k = 3;
469     const double beta = 100.0;
470     typedef double Scalar;
471     typedef alps::gf::piecewise_polynomial<Scalar> pp_type;
472 
473     alps::testing::unique_file ufile("nm.h5.", alps::testing::unique_file::REMOVE_NOW);
474     const std::string&  filename = ufile.name();
475 
476     std::vector<double> section_edges(n_section+1);
477     section_edges[0] = -1.0;
478     section_edges[1] =  0.0;
479     section_edges[2] =  1.0;
480 
481     boost::multi_array<Scalar,2> coeff(boost::extents[n_section][k+1]);
482     std::fill(coeff.origin(), coeff.origin()+coeff.num_elements(), 0.0);
483 
484     pp_type p(n_section, section_edges, coeff);
485     std::vector<pp_type> basis_functions;
486     basis_functions.push_back(p);
487     basis_functions.push_back(p);
488 
489     alps::gf::numerical_mesh<double> mesh1(beta, basis_functions, alps::gf::statistics::FERMIONIC);
490     {
491         alps::hdf5::archive oar(filename,"w");
492         mesh1.save(oar,"/nm");
493 
494     }
495 
496     alps::gf::numerical_mesh<double> mesh2;
497     {
498         alps::hdf5::archive iar(filename);
499         mesh2.load(iar,"/nm");
500     }
501 }
502 
TEST(Mesh,NumericalMeshSaveStream)503 TEST(Mesh,NumericalMeshSaveStream) {
504     const int n_section = 2, k = 3;
505     const double beta = 100.0;
506     typedef double Scalar;
507     typedef alps::gf::piecewise_polynomial<Scalar> pp_type;
508 
509     alps::testing::unique_file ufile("nm.h5.", alps::testing::unique_file::REMOVE_NOW);
510     const std::string&  filename = ufile.name();
511 
512     std::vector<double> section_edges(n_section+1);
513     section_edges[0] = -1.0;
514     section_edges[1] =  0.0;
515     section_edges[2] =  1.0;
516 
517     boost::multi_array<Scalar,2> coeff(boost::extents[n_section][k+1]);
518     std::fill(coeff.origin(), coeff.origin()+coeff.num_elements(), 0.0);
519 
520     pp_type p(n_section, section_edges, coeff);
521     std::vector<pp_type> basis_functions;
522     basis_functions.push_back(p);
523     basis_functions.push_back(p);
524 
525     alps::gf::numerical_mesh<double> mesh1(beta, basis_functions, alps::gf::statistics::FERMIONIC);
526     {
527         alps::hdf5::archive oar(filename,"w");
528         oar["/nm"] << mesh1;
529 
530     }
531 
532     alps::gf::numerical_mesh<double> mesh2;
533     {
534         alps::hdf5::archive iar(filename);
535         iar["/nm"] >> mesh2;
536     }
537 }
538 
TEST(Mesh,DefaultConstructive)539 TEST(Mesh,DefaultConstructive) {
540     alps::testing::unique_file ufile("m.h5.", alps::testing::unique_file::REMOVE_NOW);
541     const std::string&  filename = ufile.name();
542     alps::hdf5::archive oar(filename, "w");
543 
544     //FIXME: Can we use TYPED_TEST to cover all mesh types?
545     {
546         alps::gf::matsubara_pn_mesh m;
547         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
548     }
549 
550     {
551         alps::gf::power_mesh m;
552         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
553     }
554 
555     {
556         alps::gf::matsubara_positive_mesh m;
557         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
558     }
559 
560     {
561         alps::gf::momentum_index_mesh m;
562         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
563     }
564 
565     {
566         alps::gf::itime_mesh m;
567         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
568     }
569 
570     {
571         alps::gf::real_frequency_mesh m;
572         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
573     }
574 
575     {
576         alps::gf::index_mesh m;
577         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
578     }
579 
580     {
581         alps::gf::legendre_mesh m;
582         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
583     }
584 
585     {
586         alps::gf::numerical_mesh<double> m;
587         EXPECT_THROW(oar["/mesh"] << m, std::runtime_error);
588     }
589 }
590