From d8c34e0e66db350a044dde91bc51e54f93eb46d3 Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Tue, 2 Apr 2024 07:26:51 -0400 Subject: [PATCH 1/6] initial version of writing function coeffs --- src/examples/CMakeLists.txt | 2 +- src/examples/writecoeff.cc | 82 +++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+), 1 deletion(-) create mode 100644 src/examples/writecoeff.cc diff --git a/src/examples/CMakeLists.txt b/src/examples/CMakeLists.txt index 7a6e447c9e3..8a70b74992e 100644 --- a/src/examples/CMakeLists.txt +++ b/src/examples/CMakeLists.txt @@ -6,7 +6,7 @@ set(EXAMPLE_SOURCES dataloadbal hatom_1d binaryop dielectric hehf 3dharmonic testsolver testspectralprop dielectric_external_field tiny h2dynamic newsolver testcomplexfunctionsolver helium_exact density_smoothing siam_example ac_corr dirac-hatom - derivatives array_worldobject) + derivatives array_worldobject writecoeff) if(LIBXC_FOUND) list(APPEND EXAMPLE_SOURCES hefxc) diff --git a/src/examples/writecoeff.cc b/src/examples/writecoeff.cc new file mode 100644 index 00000000000..d57e95245be --- /dev/null +++ b/src/examples/writecoeff.cc @@ -0,0 +1,82 @@ +#include +#include +#include +#include + +using namespace madness; + +static const size_t D = 2; +typedef Vector coordT; +typedef std::shared_ptr< FunctionFunctorInterface,D> > functorT; +typedef Function,D> cfunctionT; +typedef FunctionFactory,D> factoryT; +typedef SeparatedConvolution,D> operatorT; + +static const double R = 1.4; // bond length +static const double L = 32.0*R; // box size +static const long k = 3; // wavelet order +static const double thresh = 1e-3; // precision + +static std::complex f(const coordT& r) +{ + return std::complex(0.0,2.0); +} + +template +class WriteCoeffImpl : public Mutex { + const int k; + std::ostream& out; +public: + WriteCoeffImpl(const int k, std::ostream& out) + : k(k) + , out(out) + {} + + void operator()(const Key& key, const Tensor< T >& t) const { + ScopedMutex obolus(*this); + out << key << " " << t << std::endl; + } +}; + +template +class WriteCoeff { + std::shared_ptr> impl; + +public: + WriteCoeff(const int k, std::ostream& out) + : impl(new WriteCoeffImpl(k, out)) + {} + + void operator()(const Key& key, const Tensor< T >& t) const { + (*impl)(key, t); + } +}; + + + +int main(int argc, char** argv) +{ + initialize(argc, argv); + World world(SafeMPI::COMM_WORLD); + + startup(world,argc,argv); + std::cout.precision(6); + + FunctionDefaults::set_k(k); + FunctionDefaults::set_thresh(thresh); + FunctionDefaults::set_refine(true); + FunctionDefaults::set_initial_level(2); + FunctionDefaults::set_truncate_mode(0); + FunctionDefaults::set_cubic_cell(-L/2, L/2); + + cfunctionT fun = factoryT(world).f(f); + fun.truncate(); + + cfunctionT sqrt_of_fun = copy(fun); + auto op = WriteCoeff,D>(k, std::cout); + fun.unaryop(op); + world.gop.fence(); + + finalize(); + return 0; +} From a09e349be4a268a0fb7e4972342fbe8685ded74e Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Tue, 2 Apr 2024 08:16:00 -0400 Subject: [PATCH 2/6] next version of writing coeffs that can run distributed --- src/examples/CMakeLists.txt | 2 +- src/examples/writecoeff2.cc | 78 +++++++++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 1 deletion(-) create mode 100644 src/examples/writecoeff2.cc diff --git a/src/examples/CMakeLists.txt b/src/examples/CMakeLists.txt index 8a70b74992e..69744a8e6b7 100644 --- a/src/examples/CMakeLists.txt +++ b/src/examples/CMakeLists.txt @@ -6,7 +6,7 @@ set(EXAMPLE_SOURCES dataloadbal hatom_1d binaryop dielectric hehf 3dharmonic testsolver testspectralprop dielectric_external_field tiny h2dynamic newsolver testcomplexfunctionsolver helium_exact density_smoothing siam_example ac_corr dirac-hatom - derivatives array_worldobject writecoeff) + derivatives array_worldobject writecoeff writecoeff2) if(LIBXC_FOUND) list(APPEND EXAMPLE_SOURCES hefxc) diff --git a/src/examples/writecoeff2.cc b/src/examples/writecoeff2.cc new file mode 100644 index 00000000000..08e66c89e70 --- /dev/null +++ b/src/examples/writecoeff2.cc @@ -0,0 +1,78 @@ +#include +#include +#include + +using namespace madness; + +static const size_t D = 2; +typedef Vector coordT; +typedef Key keyT; +typedef std::shared_ptr< FunctionFunctorInterface,D> > functorT; +typedef Function,D> cfunctionT; +typedef FunctionFactory,D> factoryT; +typedef SeparatedConvolution,D> operatorT; + +static const double R = 1.4; // bond length +static const double L = 32.0*R; // box size +static const long k = 3; // wavelet order +static const double thresh = 1e-3; // precision + +static std::complex f(const coordT& r) +{ + return std::complex(0.0,2.0); +} + +template +void write_function_coeffs(const typename FunctionImpl::dcT& coeffs, std::ostream& out, const Key& key) { + auto it = coeffs.find(key).get(); + if (it == coeffs.end()) { + for (int i=0; i " << coeffs.owner(key) << "\n"; + } + else { + const auto& node = it->second; + for (int i=0; i " << "\n"; + if (node.has_children()) { + for (KeyChildIterator kit(key); kit; ++kit) { + write_function_coeffs(coeffs, out, kit.key()); + } + } + } +} + +template +void write_function(const Function& f, std::ostream& out) { + if (f.get_impl()->world.rank() == 0) { + out << NDIM << " " << f.k() << " " << FunctionDefaults::get_cell() << std::endl; + + write_function_coeffs(f.get_impl()->get_coeffs(), out, Key(0)); + } + f.get_impl()->world.gop.fence(); +} + +void test(World& world) { + cfunctionT fun = factoryT(world).f(f); + fun.truncate(); + write_function(fun,std::cout); +} + +int main(int argc, char** argv) +{ + World& world = initialize(argc, argv); + startup(world,argc,argv); + std::cout.precision(6); + + FunctionDefaults::set_k(k); + FunctionDefaults::set_thresh(thresh); + FunctionDefaults::set_refine(true); + FunctionDefaults::set_initial_level(2); + FunctionDefaults::set_truncate_mode(0); + FunctionDefaults::set_cubic_cell(-L/2, L/2); + + test(world); + + world.gop.fence(); + finalize(); + return 0; +} From 7c3c50f6fae394a1171b9f7f111641ed1cacefcb Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Tue, 2 Apr 2024 08:52:04 -0400 Subject: [PATCH 3/6] correctly print values --- src/examples/writecoeff2.cc | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/src/examples/writecoeff2.cc b/src/examples/writecoeff2.cc index 08e66c89e70..08cb05cb29f 100644 --- a/src/examples/writecoeff2.cc +++ b/src/examples/writecoeff2.cc @@ -7,23 +7,25 @@ using namespace madness; static const size_t D = 2; typedef Vector coordT; typedef Key keyT; -typedef std::shared_ptr< FunctionFunctorInterface,D> > functorT; -typedef Function,D> cfunctionT; -typedef FunctionFactory,D> factoryT; -typedef SeparatedConvolution,D> operatorT; +typedef double dataT; // was std::complex +typedef std::shared_ptr< FunctionFunctorInterface > functorT; +typedef Function cfunctionT; +typedef FunctionFactory factoryT; +typedef SeparatedConvolution operatorT; -static const double R = 1.4; // bond length -static const double L = 32.0*R; // box size -static const long k = 3; // wavelet order +static const double L = 4.0; +static const long k = 5; // wavelet order static const double thresh = 1e-3; // precision -static std::complex f(const coordT& r) +static dataT f(const coordT& r) { - return std::complex(0.0,2.0); + double R = r.normf(); + return std::exp(-R*R); } template -void write_function_coeffs(const typename FunctionImpl::dcT& coeffs, std::ostream& out, const Key& key) { +void write_function_coeffs(const Function& f, std::ostream& out, const Key& key) { + const auto& coeffs = f.get_impl()->get_coeffs(); auto it = coeffs.find(key).get(); if (it == coeffs.end()) { for (int i=0; i::dcT& coeffs, std } else { const auto& node = it->second; - for (int i=0; i " << "\n"; + if (node.has_coeff()) { + auto values = f.get_impl()->coeffs2values(key, node.coeff()); + for (int i=0; i kit(key); kit; ++kit) { - write_function_coeffs(coeffs, out, kit.key()); + write_function_coeffs(f, out, kit.key()); } } } @@ -46,7 +53,7 @@ void write_function(const Function& f, std::ostream& out) { if (f.get_impl()->world.rank() == 0) { out << NDIM << " " << f.k() << " " << FunctionDefaults::get_cell() << std::endl; - write_function_coeffs(f.get_impl()->get_coeffs(), out, Key(0)); + write_function_coeffs(f, out, Key(0)); } f.get_impl()->world.gop.fence(); } From 3cf8130de511858dcefb77805cf7a178b735f0f0 Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Tue, 2 Apr 2024 10:30:13 -0400 Subject: [PATCH 4/6] working (?) example --- src/examples/writecoeff2.cc | 100 +++++++++++++++++++++++++++++++++--- 1 file changed, 94 insertions(+), 6 deletions(-) diff --git a/src/examples/writecoeff2.cc b/src/examples/writecoeff2.cc index 08cb05cb29f..69dd18da1e2 100644 --- a/src/examples/writecoeff2.cc +++ b/src/examples/writecoeff2.cc @@ -1,3 +1,5 @@ +#include +#include #include #include #include @@ -9,7 +11,7 @@ typedef Vector coordT; typedef Key keyT; typedef double dataT; // was std::complex typedef std::shared_ptr< FunctionFunctorInterface > functorT; -typedef Function cfunctionT; +typedef Function functionT; typedef FunctionFactory factoryT; typedef SeparatedConvolution operatorT; @@ -36,8 +38,10 @@ void write_function_coeffs(const Function& f, std::ostream& out, const K if (node.has_coeff()) { auto values = f.get_impl()->coeffs2values(key, node.coeff()); for (int i=0; i& f, std::ostream& out, const K template void write_function(const Function& f, std::ostream& out) { + auto flags = out.flags(); + auto precision = out.precision(); + out << std::setprecision(17); + out << std::scientific; + if (f.get_impl()->world.rank() == 0) { - out << NDIM << " " << f.k() << " " << FunctionDefaults::get_cell() << std::endl; + out << NDIM << std::endl; + const auto& cell = FunctionDefaults::get_cell(); + for (int d=0; d(0)); } f.get_impl()->world.gop.fence(); + + out << std::setprecision(precision); + out.setf(flags); +} + +template +void read_function_coeffs(Function& f, std::istream& in, const Key& key) { + auto& coeffs = f.get_impl()->get_coeffs(); + + while (1) { + Level n; + Vector l; + long dims[NDIM]; + in >> n; + if (in.eof()) break; + + for (int i=0; i> l[i]; + dims[i] = f.k(); + } + Key key(n,l); + + Tensor values(NDIM,dims); + for (size_t i=0; i< values.size(); i++) in >> values.ptr()[i]; + auto t = f.get_impl()->values2coeffs(key, values); + + //f.get_impl()->accumulate2(t, coeffs, key); + coeffs.task(key, &FunctionNode::accumulate2, t, coeffs, key); + } +} + +template +Function read_function(World& world, std::istream& in) { + size_t ndim; + in >> ndim; + MADNESS_CHECK(ndim == NDIM); + + Tensor cell(NDIM,2); + for (int d=0; d> cell(d,i); + } + FunctionDefaults::set_cell(cell); + + int k; + in >> k; + FunctionFactory factory(world); + Function f(factory.k(k).empty()); + world.gop.fence(); + + read_function_coeffs(f, in, Key(0)); + + f.verify_tree(); + + return f; } void test(World& world) { - cfunctionT fun = factoryT(world).f(f); + functionT fun = factoryT(world).f(f); fun.truncate(); - write_function(fun,std::cout); + + { + double norm = fun.norm2(); + if (world.rank() == 0) std::cout << "norm = " << norm << std::endl; + std::ofstream out("fun.dat", std::ios::out); + write_function(fun,out); + out.close(); + //fun.print_tree(); + } + + { + std::ifstream in("fun.dat", std::ios::in); + functionT fun2 = read_function(world, in); + double norm = fun2.norm2(); + if (world.rank() == 0) std::cout << "norm = " << norm << std::endl; + //write_function(fun2,std::cout); + //fun2.print_tree(); + double err = (fun - fun2).norm2(); + if (world.rank() == 0) std::cout << "error = " << err << std::endl; + } } int main(int argc, char** argv) From a781b6ae79f1c454f1d9c282bb2ac43a6d015b5b Mon Sep 17 00:00:00 2001 From: "Robert J. Harrison" Date: Wed, 3 Apr 2024 05:09:43 -0400 Subject: [PATCH 5/6] count leaf nodes --- src/examples/writecoeff2.cc | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/examples/writecoeff2.cc b/src/examples/writecoeff2.cc index 69dd18da1e2..a2685662670 100644 --- a/src/examples/writecoeff2.cc +++ b/src/examples/writecoeff2.cc @@ -52,8 +52,26 @@ void write_function_coeffs(const Function& f, std::ostream& out, const K } } +template +size_t count_leaf_nodes(const Function& f) { + const auto& coeffs = f.get_impl()->get_coeffs(); + size_t count = 0; + for (auto it = coeffs.begin(); it != coeffs.end(); ++it) { + const auto& key = it->first; + const auto& node = it->second; + if (node.has_coeff()) { + count++; + } + } + f.get_impl()->world.gop.sum(count); + return count; +} + template void write_function(const Function& f, std::ostream& out) { + f.reconstruct(); + std::cout << "NUMBER OF LEAF NODES: " << count_leaf_nodes(f) << std::endl; + auto flags = out.flags(); auto precision = out.precision(); out << std::setprecision(17); From 31b59ace9e5d239c2605d7d5c4cc1281dc7cfd76 Mon Sep 17 00:00:00 2001 From: Adrian Hurtado Date: Wed, 3 Apr 2024 14:06:38 +0200 Subject: [PATCH 6/6] output number of leaves in function file --- src/examples/writecoeff2.cc | 150 +++++++++++++++++------------------- 1 file changed, 72 insertions(+), 78 deletions(-) diff --git a/src/examples/writecoeff2.cc b/src/examples/writecoeff2.cc index a2685662670..68b3e2c7203 100644 --- a/src/examples/writecoeff2.cc +++ b/src/examples/writecoeff2.cc @@ -1,77 +1,71 @@ -#include -#include #include -#include +#include +#include #include +#include using namespace madness; static const size_t D = 2; -typedef Vector coordT; +typedef Vector coordT; typedef Key keyT; -typedef double dataT; // was std::complex -typedef std::shared_ptr< FunctionFunctorInterface > functorT; -typedef Function functionT; -typedef FunctionFactory factoryT; -typedef SeparatedConvolution operatorT; +typedef double dataT;// was std::complex +typedef std::shared_ptr> functorT; +typedef Function functionT; +typedef FunctionFactory factoryT; +typedef SeparatedConvolution operatorT; static const double L = 4.0; -static const long k = 5; // wavelet order -static const double thresh = 1e-3; // precision +static const long k = 5; // wavelet order +static const double thresh = 1e-3;// precision -static dataT f(const coordT& r) -{ +static dataT f(const coordT &r) { double R = r.normf(); - return std::exp(-R*R); + return std::exp(-R * R); } -template -void write_function_coeffs(const Function& f, std::ostream& out, const Key& key) { - const auto& coeffs = f.get_impl()->get_coeffs(); +template +void write_function_coeffs(const Function &f, std::ostream &out, const Key &key) { + const auto &coeffs = f.get_impl()->get_coeffs(); auto it = coeffs.find(key).get(); if (it == coeffs.end()) { - for (int i=0; i " << coeffs.owner(key) << "\n"; - } - else { - const auto& node = it->second; + } else { + const auto &node = it->second; if (node.has_coeff()) { auto values = f.get_impl()->coeffs2values(key, node.coeff()); - for (int i=0; i kit(key); kit; ++kit) { - write_function_coeffs(f, out, kit.key()); - } + for (KeyChildIterator kit(key); kit; ++kit) { write_function_coeffs(f, out, kit.key()); } } } } -template -size_t count_leaf_nodes(const Function& f) { - const auto& coeffs = f.get_impl()->get_coeffs(); +template +size_t count_leaf_nodes(const Function &f) { + const auto &coeffs = f.get_impl()->get_coeffs(); size_t count = 0; for (auto it = coeffs.begin(); it != coeffs.end(); ++it) { - const auto& key = it->first; - const auto& node = it->second; - if (node.has_coeff()) { - count++; - } + const auto &key = it->first; + const auto &node = it->second; + if (node.has_coeff()) { count++; } } - f.get_impl()->world.gop.sum(count); + f.get_impl()->world.gop.sum(count); return count; } -template -void write_function(const Function& f, std::ostream& out) { +template +void write_function(const Function &f, std::ostream &out) { f.reconstruct(); std::cout << "NUMBER OF LEAF NODES: " << count_leaf_nodes(f) << std::endl; - + auto flags = out.flags(); auto precision = out.precision(); out << std::setprecision(17); @@ -79,13 +73,14 @@ void write_function(const Function& f, std::ostream& out) { if (f.get_impl()->world.rank() == 0) { out << NDIM << std::endl; - const auto& cell = FunctionDefaults::get_cell(); - for (int d=0; d::get_cell(); + for (int d = 0; d < NDIM; ++d) { + for (int i = 0; i < 2; ++i) out << cell(d, i) << " "; out << std::endl; } out << f.k() << std::endl; - + out << count_leaf_nodes(f) << std::endl; + write_function_coeffs(f, out, Key(0)); } f.get_impl()->world.gop.fence(); @@ -94,58 +89,58 @@ void write_function(const Function& f, std::ostream& out) { out.setf(flags); } -template -void read_function_coeffs(Function& f, std::istream& in, const Key& key) { - auto& coeffs = f.get_impl()->get_coeffs(); +template +void read_function_coeffs(Function &f, std::istream &in) { + auto &coeffs = f.get_impl()->get_coeffs(); - while (1) { + while (true) { Level n; - Vector l; + Vector l; long dims[NDIM]; in >> n; if (in.eof()) break; - - for (int i=0; i> l[i]; dims[i] = f.k(); } - Key key(n,l); - - Tensor values(NDIM,dims); - for (size_t i=0; i< values.size(); i++) in >> values.ptr()[i]; + Key key(n, l); + + Tensor values(NDIM, dims); + for (size_t i = 0; i < values.size(); i++) in >> values.ptr()[i]; auto t = f.get_impl()->values2coeffs(key, values); - //f.get_impl()->accumulate2(t, coeffs, key); - coeffs.task(key, &FunctionNode::accumulate2, t, coeffs, key); + // f.get_impl()->accumulate2(t, coeffs, key); + coeffs.task(key, &FunctionNode::accumulate2, t, coeffs, key); } } -template -Function read_function(World& world, std::istream& in) { +template +Function read_function(World &world, std::istream &in) { size_t ndim; in >> ndim; MADNESS_CHECK(ndim == NDIM); - - Tensor cell(NDIM,2); - for (int d=0; d> cell(d,i); + + Tensor cell(NDIM, 2); + for (int d = 0; d < NDIM; ++d) { + for (int i = 0; i < 2; ++i) in >> cell(d, i); } FunctionDefaults::set_cell(cell); int k; in >> k; - FunctionFactory factory(world); - Function f(factory.k(k).empty()); + FunctionFactory factory(world); + Function f(factory.k(k).empty()); world.gop.fence(); - - read_function_coeffs(f, in, Key(0)); + + read_function_coeffs(f, in); f.verify_tree(); - + return f; } -void test(World& world) { +void test(World &world) { functionT fun = factoryT(world).f(f); fun.truncate(); @@ -153,27 +148,26 @@ void test(World& world) { double norm = fun.norm2(); if (world.rank() == 0) std::cout << "norm = " << norm << std::endl; std::ofstream out("fun.dat", std::ios::out); - write_function(fun,out); + write_function(fun, out); out.close(); - //fun.print_tree(); + // fun.print_tree(); } { std::ifstream in("fun.dat", std::ios::in); - functionT fun2 = read_function(world, in); + functionT fun2 = read_function(world, in); double norm = fun2.norm2(); if (world.rank() == 0) std::cout << "norm = " << norm << std::endl; - //write_function(fun2,std::cout); - //fun2.print_tree(); + // write_function(fun2,std::cout); + // fun2.print_tree(); double err = (fun - fun2).norm2(); if (world.rank() == 0) std::cout << "error = " << err << std::endl; } } -int main(int argc, char** argv) -{ - World& world = initialize(argc, argv); - startup(world,argc,argv); +int main(int argc, char **argv) { + World &world = initialize(argc, argv); + startup(world, argc, argv); std::cout.precision(6); FunctionDefaults::set_k(k); @@ -181,7 +175,7 @@ int main(int argc, char** argv) FunctionDefaults::set_refine(true); FunctionDefaults::set_initial_level(2); FunctionDefaults::set_truncate_mode(0); - FunctionDefaults::set_cubic_cell(-L/2, L/2); + FunctionDefaults::set_cubic_cell(-L / 2, L / 2); test(world);