Permalink
Browse files

Merge pull request #1233 from dgasmith/grid_cache

DFT Grid Cache and VV10 Gradients
  • Loading branch information...
loriab committed Oct 10, 2018
2 parents e20a281 + be14b44 commit cb341495898ff07b54cbe0f076ec5ccfe4f3918f
@@ -72,14 +72,14 @@ def scf_compute_energy(self):
try: try:
self.iterations() self.iterations()
except SCFConvergenceError: except SCFConvergenceError:
self.finalize() self.finalize()
raise SCFConvergenceError("""SCF DF preiterations""", self.iteration_, self, 0, 0) raise SCFConvergenceError("""SCF DF preiterations""", self.iteration_, self, 0, 0)
core.print_out("\n DF guess converged.\n\n") core.print_out("\n DF guess converged.\n\n")
# reset the DIIS & JK objects in prep for DIRECT # reset the DIIS & JK objects in prep for DIRECT
if self.initialized_diis_manager_: if self.initialized_diis_manager_:
self.diis_manager().reset_subspace() self.diis_manager().reset_subspace()
self.integrals() self.initialize_jk(self.memory_jk_)
else: else:
self.initialize() self.initialize()
@@ -104,13 +104,58 @@ def scf_compute_energy(self):
def scf_initialize(self): def scf_initialize(self):
"""Specialized initialization, compute integrals and does everything to prepare for iterations""" """Specialized initialization, compute integrals and does everything to prepare for iterations"""
# Figure out memory distributions
# Get memory in terms of doubles
total_memory = (core.get_memory() / 8) * core.get_global_option("SCF_MEM_SAFETY_FACTOR")
# Figure out how large the DFT collocation matrices are
vbase = self.V_potential()
if vbase:
collocation_size = vbase.grid().collocation_size()
if vbase.functional().ansatz() == 1:
collocation_size *= 4 # First derivs
elif vbase.functional().ansatz() == 2:
collocation_size *= 10 # Second derivs
else:
collocation_size = 0
# Change allocation for collocation matrices based on DFT type
scf_type = core.get_global_option('SCF_TYPE').upper()
nbf = self.get_basisset("ORBITAL").nbf()
naux = self.get_basisset("DF_BASIS_SCF").nbf()
if "DIRECT" == scf_type:
jk_size = total_memory * 0.1
elif scf_type.endswith('DF'):
jk_size = naux * nbf * nbf
else:
jk_size = nbf ** 4
# Give remaining to collocation
if total_memory > jk_size:
collocation_memory = total_memory - jk_size
# Give up to 10% to collocation
elif (total_memory * 0.1) > collocation_size:
collocation_memory = collocation_size
else:
collocation_memory = total_memory * 0.1
if collocation_memory > collocation_size:
collocation_memory = collocation_size
# Set constants
self.iteration_ = 0 self.iteration_ = 0
efp_enabled = hasattr(self.molecule(), 'EFP') self.memory_jk_ = int(total_memory - collocation_memory)
self.memory_collocation_ = int(collocation_memory)
# Print out initial docc/socc/etc data
if core.get_option('SCF', "PRINT") > 0: if core.get_option('SCF', "PRINT") > 0:
core.print_out(" ==> Pre-Iterations <==\n\n") core.print_out(" ==> Pre-Iterations <==\n\n")
self.print_preiterations() self.print_preiterations()
# Initialize EFP
efp_enabled = hasattr(self.molecule(), 'EFP')
if efp_enabled: if efp_enabled:
# EFP: Set QM system, options, and callback. Display efp geom in [A] # EFP: Set QM system, options, and callback. Display efp geom in [A]
efpobj = self.molecule().EFP efpobj = self.molecule().EFP
@@ -123,13 +168,16 @@ def scf_initialize(self):
efpobj.set_electron_density_field_fn(field_fn) efpobj.set_electron_density_field_fn(field_fn)
# Initilize all integratals and perform the first guess
if self.attempt_number_ == 1: if self.attempt_number_ == 1:
mints = core.MintsHelper(self.basisset()) mints = core.MintsHelper(self.basisset())
if core.get_global_option('RELATIVISTIC') in ['X2C', 'DKH']: if core.get_global_option('RELATIVISTIC') in ['X2C', 'DKH']:
mints.set_rel_basisset(self.get_basisset('BASIS_RELATIVISTIC')) mints.set_rel_basisset(self.get_basisset('BASIS_RELATIVISTIC'))
mints.one_electron_integrals() mints.one_electron_integrals()
self.integrals() self.initialize_jk(self.memory_jk_)
if self.V_potential():
self.V_potential().build_collocation_cache(self.memory_collocation_)
core.timer_on("HF: Form core H") core.timer_on("HF: Form core H")
self.form_H() self.form_H()
@@ -507,6 +555,8 @@ def scf_finalize_energy(self):
core.set_variable(k, v) core.set_variable(k, v)
self.finalize() self.finalize()
if self.V_potential():
self.V_potential().clear_collocation_cache()
core.print_out("\nComputation Completed\n") core.print_out("\nComputation Completed\n")
@@ -556,7 +606,7 @@ def scf_print_energies(self):
self.set_variable('SCF ITERATIONS', self.iteration_) self.set_variable('SCF ITERATIONS', self.iteration_)
# Bind functions to core.HF class
core.HF.initialize = scf_initialize core.HF.initialize = scf_initialize
core.HF.iterations = scf_iterate core.HF.iterations = scf_iterate
core.HF.compute_energy = scf_compute_energy core.HF.compute_energy = scf_compute_energy
@@ -163,7 +163,8 @@ void export_functional(py::module &m) {
.def("get_block", &VBase::get_block, "Returns the requested BlockOPoints.") .def("get_block", &VBase::get_block, "Returns the requested BlockOPoints.")
.def("nblocks", &VBase::nblocks, "Total number of blocks.") .def("nblocks", &VBase::nblocks, "Total number of blocks.")
.def("quadrature_values", &VBase::quadrature_values, "Returns the quadrature values.") .def("quadrature_values", &VBase::quadrature_values, "Returns the quadrature values.")
.def("build_collocation_cache", &VBase::build_collocation_cache, "Constructs a collocation cache to prevent recomputation.")
.def("clear_collocation_cache", &VBase::clear_collocation_cache, "Clears the collocation cache.")
.def("set_D", &VBase::set_D, "Sets the internal density.") .def("set_D", &VBase::set_D, "Sets the internal density.")
.def("Dao", &VBase::set_D, "Returns internal AO density.") .def("Dao", &VBase::set_D, "Returns internal AO density.")
.def("compute_V", &VBase::compute_V, "doctsring") .def("compute_V", &VBase::compute_V, "doctsring")
@@ -260,6 +261,8 @@ void export_functional(py::module &m) {
"Returns the maximum number of points in a block.") "Returns the maximum number of points in a block.")
.def("max_functions", &MolecularGrid::max_functions, .def("max_functions", &MolecularGrid::max_functions,
"Returns the maximum number of functions in a block.") "Returns the maximum number of functions in a block.")
.def("collocation_size", &MolecularGrid::collocation_size,
"Returns the total collocation size of all blocks.")
.def("blocks", &MolecularGrid::blocks, "Returns a list of blocks."); .def("blocks", &MolecularGrid::blocks, "Returns a list of blocks.");
py::class_<DFTGrid, std::shared_ptr<DFTGrid>, MolecularGrid>(m, "DFTGrid", "docstring") py::class_<DFTGrid, std::shared_ptr<DFTGrid>, MolecularGrid>(m, "DFTGrid", "docstring")
@@ -185,7 +185,7 @@ void export_wavefunction(py::module& m) {
.def("form_H", &scf::HF::form_H, "Forms the core Hamiltonian") .def("form_H", &scf::HF::form_H, "Forms the core Hamiltonian")
.def("form_Shalf", &scf::HF::form_Shalf, "Forms the S^1/2 matrix") .def("form_Shalf", &scf::HF::form_Shalf, "Forms the S^1/2 matrix")
.def("guess", &scf::HF::guess, "Forms the guess (guarantees C, D, and E)") .def("guess", &scf::HF::guess, "Forms the guess (guarantees C, D, and E)")
.def("integrals", &scf::HF::integrals, "Sets up the JK object") .def("initialize_jk", &scf::HF::initialize_jk, "Sets up the JK object")
.def("onel_Hx", &scf::HF::onel_Hx, "One-electron Hessian-vector products.") .def("onel_Hx", &scf::HF::onel_Hx, "One-electron Hessian-vector products.")
.def("twoel_Hx", &scf::HF::twoel_Hx, "Two-electron Hessian-vector products") .def("twoel_Hx", &scf::HF::twoel_Hx, "Two-electron Hessian-vector products")
.def("cphf_Hx", &scf::HF::cphf_Hx, "CPHF Hessian-vector prodcuts (4 * J - K - K.T).") .def("cphf_Hx", &scf::HF::cphf_Hx, "CPHF Hessian-vector prodcuts (4 * J - K - K.T).")
@@ -113,7 +113,7 @@ void DFCorrGrad::print_header() const {
outfile->Printf(" OpenMP threads: %11d\n", nthreads_); outfile->Printf(" OpenMP threads: %11d\n", nthreads_);
outfile->Printf(" Integrals threads: %11d\n", df_ints_num_threads_); outfile->Printf(" Integrals threads: %11d\n", df_ints_num_threads_);
outfile->Printf(" Memory (MB): %11ld\n", (memory_ * 8L) / (1024L * 1024L)); outfile->Printf(" Memory [GiB]: %11.3f\n", ((double) memory_ * 8.0) / (1024.0 * 1024.0 * 1024.0));
outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_); outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_);
outfile->Printf(" Fitting Condition: %11.0E\n\n", condition_); outfile->Printf(" Fitting Condition: %11.0E\n\n", condition_);
@@ -232,14 +232,13 @@ void DFHelper::AO_core() {
required += 3 * nao_ * nao_ * Qshell_max_; required += 3 * nao_ * nao_ * Qshell_max_;
// a fraction of memory to use, do we want it as an option? // a fraction of memory to use, do we want it as an option?
double fraction_of_memory = 0.8; if (memory_ < required) AO_core_ = false;
if (memory_ * fraction_of_memory < required) AO_core_ = false;
if (print_lvl_ > 0) { if (print_lvl_ > 0) {
outfile->Printf(" DFHelper Memory: AOs need %.3f [GiB]; user supplied %.3f [GiB]. ", outfile->Printf(" DFHelper Memory: AOs need %.3f GiB; user supplied %.3f GiB. ",
(required / fraction_of_memory * 8 / (1024 * 1024 * 1024.0)), (required * 8 / (1024 * 1024 * 1024.0)),
(memory_ * 8 / (1024 * 1024 * 1024.0))); (memory_ * 8 / (1024 * 1024 * 1024.0)));
outfile->Printf("%s in-core AOs.\n\n", (memory_ * fraction_of_memory < required) ? "Turning off" : "Using"); outfile->Printf("%s in-core AOs.\n\n", (memory_ < required) ? "Turning off" : "Using");
} }
} }
void DFHelper::print_header() { void DFHelper::print_header() {
@@ -190,7 +190,7 @@ void CubicScalarGrid::populate_grid() {
} }
} }
} }
blocks_.push_back(std::make_shared<BlockOPoints>(block_size, xp, yp, zp, wp, extents_)); blocks_.push_back(std::make_shared<BlockOPoints>(0, block_size, xp, yp, zp, wp, extents_));
} }
} }
} }
@@ -158,7 +158,7 @@ void CDJK::print_header() const
} }
outfile->Printf( " OpenMP threads: %11d\n", omp_nthread_); outfile->Printf( " OpenMP threads: %11d\n", omp_nthread_);
outfile->Printf( " Integrals threads: %11d\n", df_ints_num_threads_); outfile->Printf( " Integrals threads: %11d\n", df_ints_num_threads_);
outfile->Printf( " Memory (MB): %11ld\n", (memory_ *8L) / (1024L * 1024L)); outfile->Printf( " Memory [MiB]: %11ld\n", (memory_ *8L) / (1024L * 1024L));
outfile->Printf( " Algorithm: %11s\n", (is_core_ ? "Core" : "Disk")); outfile->Printf( " Algorithm: %11s\n", (is_core_ ? "Core" : "Disk"));
outfile->Printf( " Integral Cache: %11s\n", df_ints_io_.c_str()); outfile->Printf( " Integral Cache: %11s\n", df_ints_io_.c_str());
outfile->Printf( " Schwarz Cutoff: %11.0E\n", cutoff_); outfile->Printf( " Schwarz Cutoff: %11.0E\n", cutoff_);
@@ -84,7 +84,7 @@ void DirectJK::print_header() const
if (do_wK_) if (do_wK_)
outfile->Printf( " Omega: %11.3E\n", omega_); outfile->Printf( " Omega: %11.3E\n", omega_);
outfile->Printf( " Integrals threads: %11d\n", df_ints_num_threads_); outfile->Printf( " Integrals threads: %11d\n", df_ints_num_threads_);
//outfile->Printf( " Memory (MB): %11ld\n", (memory_ *8L) / (1024L * 1024L)); //outfile->Printf( " Memory [MiB]: %11ld\n", (memory_ *8L) / (1024L * 1024L));
outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_); outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_);
} }
} }
@@ -612,7 +612,7 @@ void DirectJK::print_header() const
if (do_wK_) if (do_wK_)
outfile->Printf( " Omega: %11.3E\n", omega_); outfile->Printf( " Omega: %11.3E\n", omega_);
outfile->Printf( " OpenMP threads: %11d\n", omp_nthread_); outfile->Printf( " OpenMP threads: %11d\n", omp_nthread_);
outfile->Printf( " Memory (MB): %11ld\n", (memory_ *8L) / (1024L * 1024L)); outfile->Printf( " Memory [MiB]: %11ld\n", (memory_ *8L) / (1024L * 1024L));
outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_); outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_);
} }
} }
@@ -257,7 +257,7 @@ void DiskDFJK::print_header() const {
if (do_wK_) outfile->Printf(" Omega: %11.3E\n", omega_); if (do_wK_) outfile->Printf(" Omega: %11.3E\n", omega_);
outfile->Printf(" OpenMP threads: %11d\n", omp_nthread_); outfile->Printf(" OpenMP threads: %11d\n", omp_nthread_);
outfile->Printf(" Integrals threads: %11d\n", df_ints_num_threads_); outfile->Printf(" Integrals threads: %11d\n", df_ints_num_threads_);
outfile->Printf(" Memory (MB): %11ld\n", (memory_ * 8L) / (1024L * 1024L)); outfile->Printf(" Memory [MiB]: %11ld\n", (memory_ * 8L) / (1024L * 1024L));
outfile->Printf(" Algorithm: %11s\n", (is_core_ ? "Core" : "Disk")); outfile->Printf(" Algorithm: %11s\n", (is_core_ ? "Core" : "Disk"));
outfile->Printf(" Integral Cache: %11s\n", df_ints_io_.c_str()); outfile->Printf(" Integral Cache: %11s\n", df_ints_io_.c_str());
outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_); outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_);
@@ -70,7 +70,7 @@ void DiskJK::print_header() const
outfile->Printf( " J tasked: %11s\n", (do_J_ ? "Yes" : "No")); outfile->Printf( " J tasked: %11s\n", (do_J_ ? "Yes" : "No"));
outfile->Printf( " K tasked: %11s\n", (do_K_ ? "Yes" : "No")); outfile->Printf( " K tasked: %11s\n", (do_K_ ? "Yes" : "No"));
outfile->Printf( " wK tasked: %11s\n", (do_wK_ ? "Yes" : "No")); outfile->Printf( " wK tasked: %11s\n", (do_wK_ ? "Yes" : "No"));
outfile->Printf( " Memory (MB): %11ld\n", (memory_ *8L) / (1024L * 1024L)); outfile->Printf( " Memory [MiB]: %11ld\n", (memory_ *8L) / (1024L * 1024L));
if (do_wK_) if (do_wK_)
outfile->Printf( " Omega: %11.3E\n", omega_); outfile->Printf( " Omega: %11.3E\n", omega_);
outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_); outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_);
@@ -130,7 +130,7 @@ void MemDFJK::print_header() const {
outfile->Printf(" wK tasked: %11s\n", (do_wK_ ? "Yes" : "No")); outfile->Printf(" wK tasked: %11s\n", (do_wK_ ? "Yes" : "No"));
if (do_wK_) outfile->Printf(" Omega: %11.3E\n", omega_); if (do_wK_) outfile->Printf(" Omega: %11.3E\n", omega_);
outfile->Printf(" OpenMP threads: %11d\n", omp_nthread_); outfile->Printf(" OpenMP threads: %11d\n", omp_nthread_);
outfile->Printf(" Memory (MB): %11ld\n", (memory_ * 8L) / (1024L * 1024L)); outfile->Printf(" Memory [MiB]: %11ld\n", (memory_ * 8L) / (1024L * 1024L));
outfile->Printf(" Algorithm: %11s\n", (dfh_->get_AO_core() ? "Core" : "Disk")); outfile->Printf(" Algorithm: %11s\n", (dfh_->get_AO_core() ? "Core" : "Disk"));
outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_); outfile->Printf(" Schwarz Cutoff: %11.0E\n", cutoff_);
outfile->Printf(" Mask sparsity (%%): %11.4f\n", 100. * dfh_->ao_sparsity()); outfile->Printf(" Mask sparsity (%%): %11.4f\n", 100. * dfh_->ao_sparsity());
@@ -87,7 +87,7 @@ void PKJK::print_header() const
outfile->Printf( " wK tasked: %11s\n", (do_wK_ ? "Yes" : "No")); outfile->Printf( " wK tasked: %11s\n", (do_wK_ ? "Yes" : "No"));
if (do_wK_) if (do_wK_)
outfile->Printf( " Omega: %11.3E\n", omega_); outfile->Printf( " Omega: %11.3E\n", omega_);
outfile->Printf( " Memory (MB): %11ld\n", (memory_ *8L) / (1024L * 1024L)); outfile->Printf( " Memory [MiB]: %11ld\n", (memory_ *8L) / (1024L * 1024L));
outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_); outfile->Printf( " Schwarz Cutoff: %11.0E\n\n", cutoff_);
outfile->Printf( " OpenMP threads: %11d\n\n", nthreads_); outfile->Printf( " OpenMP threads: %11d\n\n", nthreads_);
} }
@@ -3839,7 +3839,8 @@ void BasisExtents::print(std::string out)
} }
BlockOPoints::BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w, BlockOPoints::BlockOPoints(SharedVector x, SharedVector y, SharedVector z, SharedVector w,
std::shared_ptr<BasisExtents> extents) std::shared_ptr<BasisExtents> extents)
: npoints_(x->dimpi().sum()), : index_(0),
npoints_(x->dimpi().sum()),
xvec_(x), xvec_(x),
yvec_(y), yvec_(y),
zvec_(z), zvec_(z),
@@ -3852,9 +3853,9 @@ BlockOPoints::BlockOPoints(SharedVector x, SharedVector y, SharedVector z, Share
bound(); bound();
populate(); populate();
} }
BlockOPoints::BlockOPoints(int npoints, double *x, double *y, double *z, double *w, BlockOPoints::BlockOPoints(size_t index, size_t npoints, double *x, double *y, double *z, double *w,
std::shared_ptr<BasisExtents> extents) std::shared_ptr<BasisExtents> extents)
: npoints_(npoints), x_(x), y_(y), z_(z), w_(w), extents_(extents) { : index_(index), npoints_(npoints), x_(x), y_(y), z_(z), w_(w), extents_(extents) {
bound(); bound();
populate(); populate();
} }
@@ -3926,6 +3927,7 @@ void BlockOPoints::populate()
} }
} }
} }
local_nbf_ = functions_local_to_global_.size();
} }
void BlockOPoints::print(std::string out, int print) void BlockOPoints::print(std::string out, int print)
{ {
@@ -4129,6 +4131,7 @@ void MolecularGrid::block(int max_points, int min_points, double max_radius)
npoints_ = blocker->npoints(); npoints_ = blocker->npoints();
max_points_ = blocker->max_points(); max_points_ = blocker->max_points();
max_functions_ = blocker->max_functions(); max_functions_ = blocker->max_functions();
collocation_size_ = blocker->collocation_size();
const std::vector<std::shared_ptr<BlockOPoints> >& block = blocker->blocks(); const std::vector<std::shared_ptr<BlockOPoints> >& block = blocker->blocks();
for (size_t i = 0; i < block.size(); i++) { for (size_t i = 0; i < block.size(); i++) {
@@ -4190,18 +4193,19 @@ void MolecularGrid::print(std::string out, int /*print*/) const
std::shared_ptr<psi::PsiOutStream> printer=(out=="outfile"?outfile: std::shared_ptr<psi::PsiOutStream> printer=(out=="outfile"?outfile:
std::make_shared<PsiOutStream>(out)); std::make_shared<PsiOutStream>(out));
printer->Printf(" => Molecular Quadrature <=\n\n"); printer->Printf(" => Molecular Quadrature <=\n\n");
printer->Printf(" Radial Scheme = %14s\n" , RadialGridMgr::SchemeName(options_.radscheme)); printer->Printf(" Radial Scheme = %14s\n" , RadialGridMgr::SchemeName(options_.radscheme));
printer->Printf(" Pruning Scheme = %14s\n" , RadialPruneMgr::SchemeName(options_.prunescheme)); printer->Printf(" Pruning Scheme = %14s\n" , RadialPruneMgr::SchemeName(options_.prunescheme));
printer->Printf(" Nuclear Scheme = %14s\n", NuclearWeightMgr::SchemeName(options_.nucscheme)); printer->Printf(" Nuclear Scheme = %14s\n", NuclearWeightMgr::SchemeName(options_.nucscheme));
printer->Printf("\n"); printer->Printf("\n");
printer->Printf(" BS radius alpha = %14g\n", options_.bs_radius_alpha); printer->Printf(" BS radius alpha = %14g\n", options_.bs_radius_alpha);
printer->Printf(" Pruning alpha = %14g\n", options_.pruning_alpha); printer->Printf(" Pruning alpha = %14g\n", options_.pruning_alpha);
printer->Printf(" Radial Points = %14d\n", options_.nradpts); printer->Printf(" Radial Points = %14d\n", options_.nradpts);
printer->Printf(" Spherical Points = %14d\n", options_.nangpts); printer->Printf(" Spherical Points = %14d\n", options_.nangpts);
printer->Printf(" Total Points = %14d\n", npoints_); printer->Printf(" Total Points = %14d\n", npoints_);
printer->Printf(" Total Blocks = %14zu\n", blocks_.size()); printer->Printf(" Total Blocks = %14zu\n", blocks_.size());
printer->Printf(" Max Points = %14d\n", max_points_); printer->Printf(" Max Points = %14d\n", max_points_);
printer->Printf(" Max Functions = %14d\n", max_functions_); printer->Printf(" Max Functions = %14d\n", max_functions_);
// printer->Printf(" Collocation Size [MiB] = %14d\n", (int)((8.0 * collocation_size_) / (1024.0 * 1024.0)));
printer->Printf("\n"); printer->Printf("\n");
} }
@@ -4250,6 +4254,7 @@ void NaiveGridBlocker::block()
npoints_ = npoints_ref_; npoints_ = npoints_ref_;
max_points_ = tol_max_points_; max_points_ = tol_max_points_;
max_functions_ = extents_->basis()->nbf(); max_functions_ = extents_->basis()->nbf();
collocation_size_ = max_points_ * max_functions_;
x_ = new double[npoints_]; x_ = new double[npoints_];
y_ = new double[npoints_]; y_ = new double[npoints_];
@@ -4264,9 +4269,9 @@ void NaiveGridBlocker::block()
::memcpy((void*)index_,(void*)index_ref_, sizeof(int)*npoints_); ::memcpy((void*)index_,(void*)index_ref_, sizeof(int)*npoints_);
blocks_.clear(); blocks_.clear();
for (int Q = 0; Q < npoints_; Q += max_points_) { for (size_t Q = 0; Q < npoints_; Q += max_points_) {
int n = (Q + max_points_ >= npoints_ ? npoints_ - Q : max_points_); size_t n = (Q + max_points_ >= npoints_ ? npoints_ - Q : max_points_);
blocks_.push_back(std::make_shared<BlockOPoints>(n,&x_[Q],&y_[Q],&z_[Q],&w_[Q], extents_)); blocks_.push_back(std::make_shared<BlockOPoints>(Q, n,&x_[Q],&y_[Q],&z_[Q],&w_[Q], extents_));
} }
} }
OctreeGridBlocker::OctreeGridBlocker(const int npoints_ref, double const* x_ref, double const* y_ref, double const* z_ref, OctreeGridBlocker::OctreeGridBlocker(const int npoints_ref, double const* x_ref, double const* y_ref, double const* z_ref,
@@ -4467,24 +4472,28 @@ void OctreeGridBlocker::block()
if (block.size()) unique_block++; if (block.size()) unique_block++;
} }
blocks_.clear(); blocks_.clear();
index = 0; index = 0;
max_points_ = 0; max_points_ = 0;
for (size_t A = 0; A < completed_tree.size(); A++) { for (size_t A = 0; A < completed_tree.size(); A++) {
std::vector<int> block = completed_tree[A]; std::vector<int> block = completed_tree[A];
if (!block.size()) continue; if (!block.size()) continue;
blocks_.push_back(std::make_shared<BlockOPoints>(block.size(),&x_[index],&y_[index],&z_[index],&w_[index],extents_)); blocks_.push_back(std::make_shared<BlockOPoints>(A, block.size(), &x_[index],
if ((size_t)max_points_ < block.size()) { &y_[index], &z_[index],
max_points_ = block.size(); &w_[index], extents_));
} if ((size_t)max_points_ < block.size()) {
index += block.size(); max_points_ = block.size();
}
index += block.size();
} }
max_functions_ = 0; max_functions_ = 0;
collocation_size_ = 0;
for (size_t A = 0; A < blocks_.size(); A++) { for (size_t A = 0; A < blocks_.size(); A++) {
if ((size_t)max_functions_ < blocks_[A]->functions_local_to_global().size()) collocation_size_ += blocks_[A]->local_nbf() * blocks_[A]->npoints();
max_functions_ = blocks_[A]->functions_local_to_global().size(); if ((size_t)max_functions_ < blocks_[A]->local_nbf()){
max_functions_ = blocks_[A]->local_nbf();
}
} }
if (bench_) { if (bench_) {
Oops, something went wrong.

0 comments on commit cb34149

Please sign in to comment.