From cbbc6b73c857126bdd4cd6f95a35f33b020ec99f Mon Sep 17 00:00:00 2001 From: fbischoff Date: Mon, 18 Sep 2023 13:37:20 +0200 Subject: [PATCH] low-rank function seems to work for exchange operator --- src/madness/chem/lowrankfunction.h | 239 ++++++++++++++++----- src/madness/chem/test_low_rank_function.cc | 205 ++++++++++++++---- src/madness/mra/vmra.h | 18 +- src/madness/world/timing_utilities.h | 7 +- 4 files changed, 369 insertions(+), 100 deletions(-) diff --git a/src/madness/chem/lowrankfunction.h b/src/madness/chem/lowrankfunction.h index bb6ee4bf237..8486d04352b 100644 --- a/src/madness/chem/lowrankfunction.h +++ b/src/madness/chem/lowrankfunction.h @@ -15,6 +15,88 @@ namespace madness { + struct LowRankFunctionParameters : QCCalculationParametersBase { + + LowRankFunctionParameters() : QCCalculationParametersBase() { + + // initialize with: key, value, comment (optional), allowed values (optional) + initialize("radius",2.0,"the radius"); + initialize("gamma",1.0,"the exponent of the correlation factor"); + initialize("volume_element",0.1,"volume covered by each grid point"); + initialize("hard_shell",true,"radius is hard"); + initialize("tol",1.e-8,"rank-reduced cholesky tolerance"); + initialize("f12type","Slater","correlation factor",{"Slater","SlaterF12"}); + initialize("orthomethod","cholesky","orthonormalization",{"cholesky","canonical","symmetric"}); + initialize("transpose","slater2","transpose of the matrix",{"slater1","slater2"}); + initialize("gridtype","random","the grid type",{"random","cartesian","spherical"}); + initialize("rhsfunctiontype","exponential","the type of function",{"delta","exponential"}); + initialize("optimize",1,"number of optimization iterations"); + } + + void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) { + read_input_and_commandline_options(world,parser,tag); + } + + double radius() const {return get("radius");} + double gamma() const {return get("gamma");} + double volume_element() const {return get("volume_element");} + double tol() const {return get("tol");} + bool hard_shell() const {return get("hard_shell");} + int optimize() const {return get("optimize");} + std::string gridtype() const {return get("gridtype");} + std::string orthomethod() const {return get("orthomethod");} + std::string rhsfunctiontype() const {return get("rhsfunctiontype");} + std::string f12type() const {return get("f12type");} + }; + + + + template + class cgrid { + public: + double volume_element=0.1; + double radius=3; + bool hard_shell=true; + + cgrid(const double volume_element, const double radius, bool hard_shell) + :volume_element(volume_element), radius(radius), hard_shell(hard_shell) { + }; + + std::vector> get_coordinates() const { + // 1D grid + double volume_element_1d=std::pow(volume_element,1./NDIM); + long ngrid=std::ceil(radius/volume_element_1d); + double stepsize=radius/ngrid; + double scale=1.0; + if (not hard_shell) scale=std::pow(2.0,1.0/(ngrid+1)); + print("scale",scale); + + std::vector coord1d; + print("volume element, stepsize, ngrid" ,volume_element, std::pow(stepsize,NDIM),stepsize,ngrid); + for (int i=0; i> result; + for (int i=0; i c{coord1d[i],coord1d[j],coord1d[k]}; + double cutoff = hard_shell ? radius : 2.0*radius; + if (c.normf() struct cartesian_grid { @@ -25,14 +107,19 @@ namespace madness { long total_n; Vector increment; + cartesian_grid(const double volume_per_gridpoint, const double radius) { + double length_per_gridpoint=std::pow(volume_per_gridpoint,1.0/NDIM); + n_per_dim=ceil(2.0*radius/length_per_gridpoint); + print("length per gridpoint, n_per_dim",length_per_gridpoint,n_per_dim); + print("volume_per_gridpoint",std::pow(length_per_gridpoint,NDIM)); + initialize(-radius,radius); + print("increment",increment); + } + + cartesian_grid(const long n_per_dim, const double lo, const double hi) : n_per_dim(n_per_dim) { - lovec.fill(lo); - hivec.fill(hi); - increment=(hivec-lovec)*(1.0/double(n_per_dim-1)); - stride=std::vector(NDIM,1l); - total_n=std::pow(n_per_dim,NDIM); - for (long i=NDIM-2; i>=0; --i) stride[i]=n_per_dim*stride[i+1]; + initialize(lo,hi); } cartesian_grid(const cartesian_grid& other) : lovec(other.lovec), @@ -46,6 +133,16 @@ namespace madness { return *this; } + void initialize(const double lo, const double hi) { + lovec.fill(lo); + hivec.fill(hi); + increment=(hivec-lovec)*(1.0/double(n_per_dim-1)); + stride=std::vector(NDIM,1l); + total_n=std::pow(n_per_dim,NDIM); + for (long i=NDIM-2; i>=0; --i) stride[i]=n_per_dim*stride[i+1]; + + } + double volume_per_gridpoint() const{ double volume=1.0; for (int i=0; i> g,h; LRFunctor lrfunctor; particle p1=particle::particle1(); @@ -372,36 +469,26 @@ namespace madness { return LowRankFunction(g * a, h); } - void project(const double volume_per_point, const double radius, const std::string gridtype, std::string rhsfunctiontype, - double tol1) { - rank_revealing_tol=tol1; - long rank=0; - if (gridtype=="random") { - // number of points within radius = variance: 0.67 * #total points = 0.67*rank - // volume of sphere 4/3 pi *r^3 - // volume per point=volume/#points in radius = volume / (0.67 * rank) - // rank= volume/(0.67/volume_per_point) - double volume=4.0/3.0*constants::pi *std::pow(radius,3.0); - rank = lround(volume/(0.67*volume_per_point )); - } - project(rank,radius,gridtype,rhsfunctiontype); - } - /// following Halko /// ||A - Q Q^T A||< epsilon /// Y = A Omega && Q = QR(Y) /// || f(1,2) - \sum_i g_i(1) h_i(2) || < epsilon /// Y_i(1) = \int f(1,2) Omega_i(2) d2 && g_i(1) = QR(Y_i(1)) && h_i(2) = \int g_i^*(1) f(1,2) d1 - void project(const long rank, const double radius, const std::string gridtype, std::string rhsfunctiontype) { + void project(const LowRankFunctionParameters& params) { timer t1(world); t1.do_print=do_print; + orthomethod=params.orthomethod(); + rank_revealing_tol=params.tol(); + // make grid std::vector> grid; - if (gridtype=="random") { + if (params.gridtype()=="random") { + double volume=4.0/3.0*constants::pi *std::pow(params.radius(),3.0); + long rank = lround(volume/(0.67*params.volume_element() )); for (int i=0; i::gaussian_random_distribution(0,radius); + auto tmp=randomgaussian::gaussian_random_distribution(0,params.radius()); auto cell=FunctionDefaults::get_cell(); auto is_in_cell = [&cell](const Vector& r) { for (int d=0; dcell(d,1)) return false; @@ -411,39 +498,57 @@ namespace madness { grid.push_back(tmp); } if (world.rank()==0 and do_print) { - double volume = 4.0/3.0*constants::pi * std::pow(radius,3.0); + double volume = 4.0/3.0*constants::pi * std::pow(params.radius(),3.0); print("volume element in random grid",volume/(0.67*rank)); } - } else if (gridtype=="cartesian") { - long nperdim=std::lround(std::pow(double(rank),1.0/3.0)); - cartesian_grid cg(nperdim,-radius,radius); - for (cg.index=0; cg(); ++cg) grid.push_back(cg.get_coordinates()); - if (world.rank()==0 and do_print) print("volume element in cartesian grid",cg.volume_per_gridpoint()); +// } else if (gridtype=="cartesian") { +// long nperdim=std::lround(std::pow(double(rank),1.0/3.0)); +// cartesian_grid cg(nperdim,-radius,radius); +// for (cg.index=0; cg(); ++cg) grid.push_back(cg.get_coordinates()); +// if (world.rank()==0 and do_print) print("volume element in cartesian grid",cg.volume_per_gridpoint()); + } else if (params.gridtype()=="spherical") { + cgrid cg(params.volume_element(),params.radius(),params.hard_shell()); + grid=cg.get_coordinates(); +// cartesian_grid cg(volume_per_point,radius); +// for (cg.index=0; cg(); ++cg) grid.push_back(cg.get_coordinates()); +// if (world.rank()==0 and do_print) print("volume element in cartesian grid",cg.volume_per_gridpoint()); } else { MADNESS_EXCEPTION("unknown grid type in project",1); } - if (world.rank()==0 and do_print) { - print("grid is",gridtype,"with radius",radius,"and",grid.size(),"gridpoints"); - print("rhs functions are",rhsfunctiontype); - } - auto Y=Yformer(grid,rhsfunctiontype); + auto Y=Yformer(grid,params.rhsfunctiontype()); t1.tag("Yforming"); + auto norms=norm2s(world,Y); + for (int i=0; irank_revealing_tol) g.push_back(Y[i]); + normalize(world,g); + print("g.size after normalization",g.size()); std::ostringstream oss; oss << std::scientific << std::setprecision(1) << rank_revealing_tol; std::string scientificString = oss.str(); double err=1.0; - double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1; - g=Y; - while (err>1.e-10) { - g=truncate(orthonormalize_rrcd(g,rank_revealing_tol),std::min(1.e-3,tight_thresh*100.0)); + // need only approximate orthonormality here, as we will re-orthonormalize in the optimization step + while (err>1.e-1) { + double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1; + err=check_orthonormality(g); + t1.tag("check orthonormality"); + auto ovlp=matrix_inner(world,g,g); + t1.tag("compute ovlp"); +// g=truncate(orthonormalize_rrcd(g,ovlp,rank_revealing_tol),tight_thresh); + g=truncate(orthonormalize_rrcd(g,ovlp,rank_revealing_tol)); + t1.tag("rrcd/truncate/thresh"); +// ovlp=matrix_inner(world,g,g); +// t1.tag("compute ovlp"); +// g=truncate(madness::orthonormalize(g)); +// t1.tag("ortho/symmetric/truncate"); + auto sz=get_size(world,g); + print("gsize",sz); + t1.tag("Y orthonormalizing rrcd "); err=check_orthonormality(g); - print("error of non-orthonormality",err); - t1.tag("Y orthonormalizing with rank_revealing_tol loose_thresh "+scientificString); + t1.tag("Y orthonormalizing with rank_revealing_tol normal thresh "+scientificString); } g=truncate(g); @@ -467,12 +572,13 @@ namespace madness { } else if (rhsfunctiontype=="exponential") { std::vector> omega; + double coeff=std::pow(2.0*exponent/constants::pi,0.25*LDIM); for (const auto& point : grid) { omega.push_back(FunctionFactory(world) - .functor([&point,&exponent](const Vector& r) + .functor([&point,&exponent,&coeff](const Vector& r) { auto r_rel=r-point; - return exp(-exponent*madness::inner(r_rel,r_rel)); + return coeff*exp(-exponent*madness::inner(r_rel,r_rel)); })); } Y=inner(lrfunctor,omega,p2,p1); @@ -490,31 +596,58 @@ namespace madness { return fapprox; } + + std::vector> orthonormalize(const std::vector>& g) const { + + double tol=rank_revealing_tol; + print("orthonormalizing with method/tol",orthomethod,tol); + std::vector> g2; + auto ovlp=matrix_inner(world,g,g); + if (orthomethod=="symmetric") g2=orthonormalize_symmetric(g,ovlp); + else if (orthomethod=="canonical") g2=orthonormalize_canonical(g,ovlp,tol); + else if (orthomethod=="cholesky") g2=orthonormalize_rrcd(g,ovlp,tol); + else { + MADNESS_EXCEPTION("no such orthomethod",1); + } + double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1; + return truncate(g2,tight_thresh); + } + + /// optimize using Cholesky decomposition - /// if stable_power_iteration is true, orthonormalize in between applications of the kernel (Alg. 4.4 in Halko) /// @param[in] nopt number of iterations (wrt to Alg. 4.3 in Halko) void optimize(const long nopt=1) { timer t(world); t.do_print=do_print; - double tight_thresh=FunctionDefaults<3>::get_thresh()*0.1; + double tight_thresh=FunctionDefaults<3>::get_thresh(); + print("tight_thresh",tight_thresh); for (int i=0; i>& v) const { + timer t(world); Tensor ovlp=matrix_inner(world,v,v); for (int i=0; i("radius",2.0,"the radius"); - initialize("gamma",1.0,"the exponent of the correlation factor"); - initialize("volume_element",0.1,"volume covered by each grid point"); - initialize("rank",500,"the number of grid points in random grids"); - initialize("stable_power_iteration",true,"use stable power iteration algorithm (orthonormalize)"); - initialize("tol",1.e-12,"rank-reduced cholesky tolerance"); - initialize("f12type","Slater","correlation factor",{"Slater","SlaterF12"}); - initialize("transpose","slater2","transpose of the matrix",{"slater1","slater2"}); - initialize("gridtype","random","the grid type",{"random","cartesian"}); - initialize("rhsfunctiontype","exponential","the type of function",{"delta","exponential"}); - initialize("optimize",1,"number of optimization iterations"); - } - - void read_and_set_derived_values(World& world, const commandlineparser& parser, std::string tag) { - read_input_and_commandline_options(world,parser,tag); - } - - double radius() const {return get("radius");} - double gamma() const {return get("gamma");} - double volume_element() const {return get("volume_element");} - double tol() const {return get("tol");} - long rank() const {return get("rank");} - bool stable_power_iteration() const {return get("stable_power_iteration");} - int optimize() const {return get("optimize");} - std::string gridtype() const {return get("gridtype");} - std::string rhsfunctiontype() const {return get("rhsfunctiontype");} -}; - int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { @@ -61,19 +28,19 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { print("eps, k, NDIM, id",FunctionDefaults::get_thresh(),FunctionDefaults::get_k(),NDIM,id); parameters.print("grid","end"); - std::string f12type=parameters.get("f12type"); + std::string f12type=parameters.f12type(); std::string transpose=parameters.get("transpose"); json j; std::string jsonfilename="test_low_rank_function."+id+".json"; j["radius"]=parameters.radius(); - j["f12type"]=parameters.get("f12type"); + j["f12type"]=parameters.f12type(); j["gamma"]=parameters.gamma(); j["volume_element"]=parameters.volume_element(); j["tol"]=parameters.tol(); - j["rank"]=parameters.rank(); + j["hard_shell"]=parameters.hard_shell(); j["transpose"]=transpose; - j["stable_power_iteration"]=parameters.stable_power_iteration(); + j["orthomethod"]=parameters.orthomethod(); j["gridtype"]=parameters.gridtype(); j["rhsfunctiontype"]=parameters.rhsfunctiontype(); j["optimize"]=parameters.optimize(); @@ -140,14 +107,13 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { output(0.0,0.0,0.0,0.0,0.0,0.0); LowRankFunction lrf(f12, copy(phi1), copy(phi2)); - lrf.stable_power_iteration=parameters.stable_power_iteration(); plot_plane<6>(world,lrf.lrfunctor,"plot_original."+id,PlotParameters(world).set_plane({"x1","x4"})); double cpu0=cpu_time(); - lrf.project(parameters.volume_element(),parameters.radius(),parameters.gridtype(),parameters.rhsfunctiontype(),parameters.tol()); + lrf.project(parameters); double cpu1=cpu_time(); double error1=lrf.l2error(); print("l2error projection",error1); - plot_plane<6>(world,lrf,"plot_lrf_projection."+id,PlotParameters(world).set_plane({"x1","x4"})); +// plot_plane<6>(world,lrf,"plot_lrf_projection."+id,PlotParameters(world).set_plane({"x1","x4"})); // compare // \phi(1) \bar \phi(1) = \int phi(1) \phi(2) f(1,2) d2 @@ -156,8 +122,8 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { real_function_3d result=compute_result(lrf); double projection_error=compute_relative_error(reference,result,lrf); auto diff=reference-result; - plot_plane<3>(world,diff,"plot_diff_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); - plot_plane<3>(world,result,"plot_lrf_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); +// plot_plane<3>(world,diff,"plot_diff_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); +// plot_plane<3>(world,result,"plot_lrf_int_projection."+id,PlotParameters(world).set_plane({"x1","x2"})); output(projection_error,lrf.rank(),cpu1-cpu0,0.0,0.0,0.0); j["projection_error_integrated"]=projection_error; j["projection_error_l2"]=error1; @@ -191,6 +157,154 @@ int test_lowrank_function(World& world, LowRankFunctionParameters& parameters) { } +/// test the K commutator of the He atom + +/// < ij | K f12 | ij > +int test_Kcommutator(World& world, LowRankFunctionParameters& parameters) { + test_output t1("CCPairFunction::low exchange commutator"); + t1.set_cout_to_terminal(); + madness::default_random_generator.setstate(int(cpu_time())%4149); + std::string id=unique_fileid(); + + constexpr std::size_t LDIM=3; + constexpr std::size_t NDIM=2*LDIM; + print("eps, k, NDIM, id",FunctionDefaults::get_thresh(),FunctionDefaults::get_k(),NDIM,id); + parameters.print("grid"); + + real_convolution_3d g12=(CoulombOperator(world,1.e-6,FunctionDefaults::get_thresh())); + std::shared_ptr f12ptr; + std::string f12type=parameters.f12type(); + if (f12type=="slaterf12") f12ptr.reset(SlaterF12OperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else if (f12type=="slater") f12ptr.reset(SlaterOperatorPtr(world,parameters.gamma(),1.e-6,FunctionDefaults::get_thresh())); + else { + MADNESS_EXCEPTION(std::string("unknown f12type"+f12type).c_str(),1); + } + real_convolution_3d& f12=*f12ptr; + + real_function_3d phi=real_factory_3d(world).f([](const coord_3d& r){return exp(-r.normf());}); + double n=phi.norm2(); + phi.scale(1/n); + real_function_3d phi_k=phi; // lookys silly, helps reading. + + + // reference term ( < ij | K(1) ) = = < Ki(1) j(2) | f12 | i(1) j(2) > = ",reference); + + json j; + std::string jsonfilename="test_kcommuntator."+id+".json"; + j["radius"]=parameters.radius(); + j["f12type"]=parameters.f12type(); + j["gamma"]=parameters.gamma(); + j["thresh"]=FunctionDefaults<3>::get_thresh(); + j["volume_element"]=parameters.volume_element(); + j["tol"]=parameters.tol(); + j["hard_shell"]=parameters.hard_shell(); + j["orthomethod"]=parameters.orthomethod(); + j["gridtype"]=parameters.gridtype(); + j["rhsfunctiontype"]=parameters.rhsfunctiontype(); + j["optimize"]=parameters.optimize(); + j["reference"]=reference; + + auto json2file= [](const json& j, const std::string& jsonfilename) { + std::ofstream of(jsonfilename, std::ios::out); + of << j; + of.close(); + }; + + json2file(j,jsonfilename); + timer t(world); + + { + // lowrankfunction left phi: lrf(1',2) = f12(1',2) i(1') + // K f12 ij = \sum_k k(1) \int g(1,1') f12(1'2) i(1') j(2) k(1') d1' + // = \sum_kr k(1) j(2) \int g(1,1') g_r(1') h_r(2) k(1') d1' + // = \sum_r j(2) h_r(2) \sum_k k(1) \int g(1,1') g_r(1') k(1') d1' + real_function_3d one = real_factory_3d(world).f([](const coord_3d& r) { return 1.0; }); + LowRankFunction fi_one(f12ptr, copy(phi), copy(one)); + fi_one.project(parameters); + j["left_project_time"]=t.tag("left_project_time"); + json2file(j,jsonfilename); + + { + auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 + auto hj = fi_one.h * phi; // function of 2 + Tensor j_hj = inner(world, phi, hj); + Tensor i_gk = inner(world, phi, gk); + double result_left = j_hj.trace(i_gk); + print("result_left, project only", result_left); + j["left_project"]=result_left-reference; + j["left_project_rank"]=fi_one.rank(); + j["left_project_compute_time"]=t.tag("left_project_compute_time"); + } + json2file(j,jsonfilename); + + fi_one.optimize(); + j["left_optimize_time"]=t.tag("left_optimize_time"); + json2file(j,jsonfilename); + { + auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 + auto hj = fi_one.h * phi; // function of 2 + Tensor j_hj = inner(world, phi, hj); + Tensor i_gk = inner(world, phi, gk); + double result_left = j_hj.trace(i_gk); + print("result_left, optimize", result_left); + j["left_optimize"]=result_left-reference; + j["left_optimize_rank"]=fi_one.rank(); + j["left_optimize_compute_time"]=t.tag("left_optimize_compute_time"); + } + json2file(j,jsonfilename); + } + + // lowrankfunction right phi: lrf(1',2) = f12(1',2) i(1') + { + real_function_3d one = real_factory_3d(world).f([](const coord_3d &r) { return 1.0; }); + LowRankFunction fi_one(f12ptr, copy(one), copy(phi)); + fi_one.project(parameters); + std::swap(fi_one.g,fi_one.h); + j["right_project_time"]=t.tag("right_project_time"); + json2file(j,jsonfilename); + + + { + auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 + auto hj = fi_one.h * phi; // function of 2 + Tensor j_hj = inner(world, phi, hj); + Tensor i_gk = inner(world, phi, gk); + double result_right = j_hj.trace(i_gk); + print("result_right, project only", result_right); + j["right_project"]=result_right-reference; + j["right_project_rank"]=fi_one.rank(); + j["left_optimize_compute_time"]=t.tag("left_optimize_compute_time"); + j["right_project_compute_time"]=t.tag("right_project_compute_time"); + } + json2file(j,jsonfilename); + std::swap(fi_one.g,fi_one.h); + fi_one.optimize(); + std::swap(fi_one.g,fi_one.h); + { + auto gk = mul(world, phi_k, g12(fi_one.g * phi_k)); // function of 1 + auto hj = fi_one.h * phi; // function of 2 + Tensor j_hj = inner(world, phi, hj); + Tensor i_gk = inner(world, phi, gk); + double result_right = j_hj.trace(i_gk); + print("result_right, optimize", result_right); + j["right_optimize"]=result_right-reference; + j["right_optimize_rank"]=fi_one.rank(); + j["right_optimize_compute_time"]=t.tag("right_optimize_compute_time"); + } + json2file(j,jsonfilename); + + } + + + return 0; + +} + + int main(int argc, char **argv) { madness::World& world = madness::initialize(argc, argv); @@ -221,19 +335,24 @@ int main(int argc, char **argv) { FunctionDefaults<5>::set_cubic_cell(-10.,10.); FunctionDefaults<6>::set_cubic_cell(-10.,10.); + FunctionDefaults<2>::set_tensor_type(TT_FULL); print("numerical parameters: k, eps(3D), eps(6D)", FunctionDefaults<3>::get_k(), FunctionDefaults<3>::get_thresh(), FunctionDefaults<6>::get_thresh()); LowRankFunctionParameters parameters; parameters.read_and_set_derived_values(world,parser,"grid"); + parameters.print("grid"); int isuccess=0; #ifdef USE_GENTENSOR + try { parser.set_keyval("geometry", "he"); parser.print_map(); - isuccess+=test_lowrank_function(world,parameters); + +// isuccess+=test_lowrank_function(world,parameters); + isuccess+=test_Kcommutator(world,parameters); } catch (std::exception& e) { madness::print("an error occured"); madness::print(e.what()); diff --git a/src/madness/mra/vmra.h b/src/madness/mra/vmra.h index bf8568c0371..de50bce735d 100644 --- a/src/madness/mra/vmra.h +++ b/src/madness/mra/vmra.h @@ -518,9 +518,12 @@ namespace madness { return v; } + auto sv=copy(ovlp); rr_cholesky(ovlp,tol,piv,rank); // destroys ovlp and gives back Upper ∆ Matrix from CCD + // ovlp zeroed such that input = inner(transpose(output),output). - // rearrange and truncate the functions according to the pivoting of the rr_cholesky + + // rearrange and truncate the functions according to the pivoting of the rr_cholesky std::vector > pv(rank); for(integer i=0;i Linv = inverse(L); Tensor U = transpose(Linv); +// // L L^T = ovlp +// // Linv ovlp LT inv = 1 +// Tensor LTinv=inverse(ovlp); +// Tensor test=inner(LTinv, inner(sv,Linv)); +// print("LTinv, inner(sv,Linv)"); +// print(test); +// +// print("ovlp - LLt"); +// print(sv-inner(L,L,1,0)); +// print("Linv ovlp Ltinv",inner(Linv, inner(sv,LTinv)).normf()); +// +// +// throw; World& world=v.front().world(); return transform(world, pv, U,tol); } diff --git a/src/madness/world/timing_utilities.h b/src/madness/world/timing_utilities.h index 8b2a80266fd..35b443b12ed 100644 --- a/src/madness/world/timing_utilities.h +++ b/src/madness/world/timing_utilities.h @@ -17,7 +17,7 @@ struct timer { sss = cpu_time(); } - void tag(const std::string msg) { + double tag(const std::string msg) { world.gop.fence(); double tt1 = wall_time() - ttt; double ss1 = cpu_time() - sss; @@ -29,10 +29,11 @@ struct timer { } ttt = wall_time(); sss = cpu_time(); + return ss1; } - void end(const std::string msg) { - tag(msg); + double end(const std::string msg) { + return tag(msg); // world.gop.fence(); // double tt1 = wall_time() - ttt; // double ss1 = cpu_time() - sss;