branch: master commit 05f69a3e91911c0739c03b27ab9360b731aba604 Author: Konstantinos Poulios <logar...@gmail.com> Date: Wed Feb 20 08:58:49 2019 +0100
Code simplifications and whitespace --- src/getfem/bgeot_tensor.h | 20 +- src/getfem/getfem_generic_assembly.h | 14 +- .../getfem_generic_assembly_compile_and_exec.h | 10 +- src/getfem/getfem_model_solvers.h | 238 ++++----- src/getfem_generic_assembly_interpolation.cc | 8 +- src/getfem_generic_assembly_semantic.cc | 3 +- src/getfem_generic_assembly_workspace.cc | 76 ++- src/getfem_mesh.cc | 2 +- src/getfem_model_solvers.cc | 24 +- src/getfem_models.cc | 109 ++-- src/gmm/gmm_vector.h | 570 ++++++++++----------- 11 files changed, 526 insertions(+), 548 deletions(-) diff --git a/src/getfem/bgeot_tensor.h b/src/getfem/bgeot_tensor.h index 43d779d..f424c2d 100644 --- a/src/getfem/bgeot_tensor.h +++ b/src/getfem/bgeot_tensor.h @@ -62,7 +62,7 @@ namespace bgeot { } else resize(1); } - void reset(void) { std::fill(begin(), end(), 0); } + void reset() { std::fill(begin(), end(), 0); } inline bool finished(const multi_index &m) { if (m.size() == 0) @@ -83,7 +83,7 @@ namespace bgeot { : std::vector<size_type>(4) { (*this)[0] = i; (*this)[1] = j; (*this)[2] = k; (*this)[3] = l; } - multi_index(void) {} + multi_index() {} bool is_equal(const multi_index &m) const { if (this->size() != m.size()) return false; @@ -92,7 +92,7 @@ namespace bgeot { return true; } - size_type total_size(void) const { + size_type total_size() const { size_type s = 1; for (size_type k = 0; k < this->size(); ++k) s *= (*this)[k]; return s; @@ -194,10 +194,10 @@ namespace bgeot { return *(this->begin() + d); } - inline size_type size(void) const { return std::vector<T>::size(); } + inline size_type size() const { return std::vector<T>::size(); } inline size_type size(size_type i) const { return sizes_[i]; } - inline const multi_index &sizes(void) const { return sizes_; } - inline size_type order(void) const { return sizes_.size(); } + inline const multi_index &sizes() const { return sizes_; } + inline size_type order() const { return sizes_.size(); } void init(const multi_index &c) { auto it = c.begin(); @@ -236,7 +236,7 @@ namespace bgeot { } inline void adjust_sizes(const multi_index &mi) { init(mi); } - inline void adjust_sizes(void) { init(); } + inline void adjust_sizes() { init(); } inline void adjust_sizes(size_type i) { init(i); } inline void adjust_sizes(size_type i, size_type j) { init(i, j); } inline void adjust_sizes(size_type i, size_type j, size_type k) @@ -294,8 +294,8 @@ namespace bgeot { + sizeof(*this) + sizes_.memsize() + coeff.memsize(); } - std::vector<T> &as_vector(void) { return *this; } - const std::vector<T> &as_vector(void) const { return *this; } + std::vector<T> &as_vector() { return *this; } + const std::vector<T> &as_vector() const { return *this; } tensor<T>& operator +=(const tensor<T>& w) @@ -325,7 +325,7 @@ namespace bgeot { tensor(const multi_index &c) { init(c); } tensor(size_type i, size_type j, size_type k, size_type l) { init(multi_index(i, j, k, l)); } - tensor(void) {} + tensor() {} }; template<class T> void tensor<T>::mat_transp_reduction diff --git a/src/getfem/getfem_generic_assembly.h b/src/getfem/getfem_generic_assembly.h index 5dd045c..5b91109 100644 --- a/src/getfem/getfem_generic_assembly.h +++ b/src/getfem/getfem_generic_assembly.h @@ -108,7 +108,7 @@ namespace getfem { std::map<var_trans_pair, base_tensor> &derivatives, bool compute_derivatives) const = 0; virtual void finalize() const = 0; - virtual std::string expression(void) const { return std::string(); } + virtual std::string expression() const { return std::string(); } virtual ~virtual_interpolate_transformation() {} }; @@ -143,9 +143,9 @@ namespace getfem { public: - const mesh_im &mim(void) const { return mim_; } + const mesh_im &mim() const { return mim_; } virtual const mesh_region &give_region(const mesh &m, - size_type cv, short_type f) const = 0; + size_type cv, short_type f) const = 0; // virtual void init(const ga_workspace &workspace) const = 0; // virtual void finalize() const = 0; @@ -367,7 +367,7 @@ namespace getfem { const mesh_region &rg, const std::string &expr, size_type add_derivative_order, bool scalar_expr, size_type for_interpolation, - const std::string varname_interpolation); + const std::string varname_interpolation); std::shared_ptr<model_real_sparse_matrix> K; @@ -409,7 +409,7 @@ namespace getfem { size_type add_expression(const std::string &expr, const mesh_im &mim, const mesh_region &rg=mesh_region::all_convexes(), size_type add_derivative_order = 2, - const std::string &secondary_dom = ""); + const std::string &secondary_dom = ""); /* Internal use */ void add_function_expression(const std::string &expr); /* Internal use */ @@ -448,9 +448,9 @@ namespace getfem { const model_real_plain_vector &VV); bool used_variables(std::vector<std::string> &vl, - std::vector<std::string> &vl_test1, + std::vector<std::string> &vl_test1, std::vector<std::string> &vl_test2, - std::vector<std::string> &dl, + std::vector<std::string> &dl, size_type order); bool variable_exists(const std::string &name) const; diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h b/src/getfem/getfem_generic_assembly_compile_and_exec.h index e64a1d5..36588f5 100644 --- a/src/getfem/getfem_generic_assembly_compile_and_exec.h +++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h @@ -212,14 +212,14 @@ namespace getfem { void ga_exec(ga_instruction_set &gis, ga_workspace &workspace); void ga_function_exec(ga_instruction_set &gis); void ga_compile(ga_workspace &workspace, ga_instruction_set &gis, - size_type order); + size_type order); void ga_compile_function(ga_workspace &workspace, - ga_instruction_set &gis, bool scalar); + ga_instruction_set &gis, bool scalar); void ga_compile_interpolation(ga_workspace &workspace, - ga_instruction_set &gis); + ga_instruction_set &gis); void ga_interpolation_exec(ga_instruction_set &gis, - ga_workspace &workspace, - ga_interpolation_context &gic); + ga_workspace &workspace, + ga_interpolation_context &gic); void ga_interpolation_single_point_exec (ga_instruction_set &gis, ga_workspace &workspace, const fem_interpolation_context &ctx_x, const base_small_vector &Normal, diff --git a/src/getfem/getfem_model_solvers.h b/src/getfem/getfem_model_solvers.h index 278e7ba..6aa3b5e 100644 --- a/src/getfem/getfem_model_solvers.h +++ b/src/getfem/getfem_model_solvers.h @@ -233,10 +233,10 @@ namespace getfem { // Dummy linesearch for the newton with step control struct newton_search_with_step_control : public abstract_newton_line_search { - + virtual void init_search(double /*r*/, size_t /*git*/, double /*R0*/ = 0.0) { GMM_ASSERT1(false, "Not to be used"); } - + virtual double next_try(void) { GMM_ASSERT1(false, "Not to be used"); } @@ -245,7 +245,7 @@ namespace getfem { newton_search_with_step_control() {} }; - + struct quadratic_newton_line_search : public abstract_newton_line_search { double R0_, R1_; @@ -408,9 +408,11 @@ namespace getfem { /*********************************************************************/ template <typename PB> - void Newton_with_step_control(PB &pb, gmm::iteration &iter, - const abstract_linear_solver<typename PB::MATRIX, - typename PB::VECTOR> &linear_solver) { + void Newton_with_step_control + (PB &pb, gmm::iteration &iter, + const abstract_linear_solver<typename PB::MATRIX, + typename PB::VECTOR> &linear_solver) + { typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T; typedef typename gmm::number_traits<T>::magnitude_type R; gmm::iteration iter_linsolv0 = iter; @@ -436,103 +438,103 @@ namespace getfem { // GMM_ASSERT1(ls, "Internal error"); size_type nit = 0, stagn = 0; R coeff = R(2); - + scalar_type crit = pb.residual_norm() / approx_eln; if (iter.finished(crit)) return; for(;;) { - + crit = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)) - / approx_eln; - if (!iter.converged(crit)) { - gmm::iteration iter_linsolv = iter_linsolv0; - if (iter.get_noisy() > 1) - cout << "starting tangent matrix computation" << endl; - - int is_singular = 1; - while (is_singular) { // Linear system solve - pb.compute_tangent_matrix(); - gmm::clear(dr); - gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b); - gmm::add(gmm::scaled(b0,alpha-R(1)), b); - if (iter.get_noisy() > 1) cout << "starting linear solver" << endl; - iter_linsolv.init(); - linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv); - if (!iter_linsolv.converged()) { - is_singular++; - if (is_singular <= 4) { - if (iter.get_noisy()) - cout << "Singular tangent matrix:" - " perturbation of the state vector." << endl; - pb.perturbation(); - pb.compute_residual(); - } else { - if (iter.get_noisy()) - cout << "Singular tangent matrix: perturbation failed, " - << "aborting." << endl; - return; - } - } - else is_singular = 0; - } - if (iter.get_noisy() > 1) cout << "linear solver done" << endl; - - - gmm::add(dr, pb.state_vector()); - pb.compute_residual(); - R res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); - R dec = R(1)/R(2), coeff2 = coeff * R(1.5); - - while (dec > R(1E-5) && res >= res0 * coeff) { - gmm::add(gmm::scaled(dr, -dec), pb.state_vector()); - pb.compute_residual(); - R res2 = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); - if (res2 < res*R(0.95) || res2 >= res0 * coeff2) { - dec /= R(2); res = res2; coeff2 *= R(1.5); - } else { - gmm::add(gmm::scaled(dr, dec), pb.state_vector()); - break; - } - } - dec *= R(2); - - nit++; - coeff = std::max(R(1.05), coeff*R(0.93)); - bool near_end = (iter.get_iteration() > iter.get_maxiter()/2); - bool cut = (alpha < R(1)) && near_end; - if ((res > minres && nit > 4) || cut) { - stagn++; - - if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) { - alpha = (alpha + alpha0) / R(2); - if (near_end) alpha = R(1); - gmm::copy(Xi, pb.state_vector()); - pb.compute_residual(); - res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); - nit = 0; - stagn = 0; coeff = R(2); - } - } - if (res < minres || (alpha == R(1) && nit < 10)) stagn = 0; - res0 = res; - if (nit < 5) minres = res0; else minres = std::min(minres, res0); - - if (iter.get_noisy()) - cout << "step control [" << std::setw(8) << alpha0 << "," - << std::setw(8) << alpha << "," << std::setw(10) << dec << "]"; - ++iter; - // crit = std::min(res / approx_eln, - // gmm::vect_norm1(dr) / std::max(1E-25,pb.state_norm())); - crit = res / approx_eln; + / approx_eln; + if (!iter.converged(crit)) { + gmm::iteration iter_linsolv = iter_linsolv0; + if (iter.get_noisy() > 1) + cout << "starting tangent matrix computation" << endl; + + int is_singular = 1; + while (is_singular) { // Linear system solve + pb.compute_tangent_matrix(); + gmm::clear(dr); + gmm::copy(gmm::scaled(pb.residual(), pb.scale_residual()), b); + gmm::add(gmm::scaled(b0,alpha-R(1)), b); + if (iter.get_noisy() > 1) cout << "starting linear solver" << endl; + iter_linsolv.init(); + linear_solver(pb.tangent_matrix(), dr, b, iter_linsolv); + if (!iter_linsolv.converged()) { + is_singular++; + if (is_singular <= 4) { + if (iter.get_noisy()) + cout << "Singular tangent matrix:" + " perturbation of the state vector." << endl; + pb.perturbation(); + pb.compute_residual(); + } else { + if (iter.get_noisy()) + cout << "Singular tangent matrix: perturbation failed, " + << "aborting." << endl; + return; + } + } + else is_singular = 0; + } + if (iter.get_noisy() > 1) cout << "linear solver done" << endl; + + + gmm::add(dr, pb.state_vector()); + pb.compute_residual(); + R res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); + R dec = R(1)/R(2), coeff2 = coeff * R(1.5); + + while (dec > R(1E-5) && res >= res0 * coeff) { + gmm::add(gmm::scaled(dr, -dec), pb.state_vector()); + pb.compute_residual(); + R res2 = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); + if (res2 < res*R(0.95) || res2 >= res0 * coeff2) { + dec /= R(2); res = res2; coeff2 *= R(1.5); + } else { + gmm::add(gmm::scaled(dr, dec), pb.state_vector()); + break; + } + } + dec *= R(2); + + nit++; + coeff = std::max(R(1.05), coeff*R(0.93)); + bool near_end = (iter.get_iteration() > iter.get_maxiter()/2); + bool cut = (alpha < R(1)) && near_end; + if ((res > minres && nit > 4) || cut) { + stagn++; + + if ((stagn > 10 && alpha > alpha0 + R(5E-2)) || cut) { + alpha = (alpha + alpha0) / R(2); + if (near_end) alpha = R(1); + gmm::copy(Xi, pb.state_vector()); + pb.compute_residual(); + res = gmm::vect_dist1(pb.residual(), gmm::scaled(b0, R(1)-alpha)); + nit = 0; + stagn = 0; coeff = R(2); + } + } + if (res < minres || (alpha == R(1) && nit < 10)) stagn = 0; + res0 = res; + if (nit < 5) minres = res0; else minres = std::min(minres, res0); + + if (iter.get_noisy()) + cout << "step control [" << std::setw(8) << alpha0 << "," + << std::setw(8) << alpha << "," << std::setw(10) << dec << "]"; + ++iter; + // crit = std::min(res / approx_eln, + // gmm::vect_norm1(dr) / std::max(1E-25,pb.state_norm())); + crit = res / approx_eln; } - + if (iter.finished(crit)) { - if (iter.converged() && alpha < R(1)) { - R a = alpha; - alpha = std::min(R(1), alpha*R(3) - alpha0*R(2)); - alpha0 = a; - gmm::copy(pb.state_vector(), Xi); - nit = 0; stagn = 0; coeff = R(2); - } else return; + if (iter.converged() && alpha < R(1)) { + R a = alpha; + alpha = std::min(R(1), alpha*R(3) - alpha0*R(2)); + alpha0 = a; + gmm::copy(pb.state_vector(), Xi); + nit = 0; stagn = 0; coeff = R(2); + } else return; } } } @@ -544,9 +546,11 @@ namespace getfem { /* ***************************************************************** */ template <typename PB> - void classical_Newton(PB &pb, gmm::iteration &iter, - const abstract_linear_solver<typename PB::MATRIX, - typename PB::VECTOR> &linear_solver) { + void classical_Newton + (PB &pb, gmm::iteration &iter, + const abstract_linear_solver<typename PB::MATRIX, + typename PB::VECTOR> &linear_solver) + { typedef typename gmm::linalg_traits<typename PB::VECTOR>::value_type T; typedef typename gmm::number_traits<T>::magnitude_type R; gmm::iteration iter_linsolv0 = iter; @@ -750,10 +754,10 @@ namespace getfem { #elif GETFEM_PARA_LEVEL > 1 && GETFEM_PARA_SOLVER == MUMPS_PARA_SOLVER if (md.is_symmetric()) return std::make_shared - <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>(); + <linear_solver_distributed_mumps_sym<MATRIX, VECTOR>>(); else - return std::make_shared - <linear_solver_distributed_mumps<MATRIX, VECTOR>>(); + return std::make_shared + <linear_solver_distributed_mumps<MATRIX, VECTOR>>(); #else size_type ndof = md.nb_dof(), max3d = 15000, dim = md.leading_dimension(); # ifdef GMM_USES_MUMPS @@ -762,24 +766,24 @@ namespace getfem { if ((ndof<300000 && dim<=2) || (ndof<max3d && dim<=3) || (ndof<1000)) { # ifdef GMM_USES_MUMPS if (md.is_symmetric()) - return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>(); + return std::make_shared<linear_solver_mumps_sym<MATRIX, VECTOR>>(); else - return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>(); + return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>(); # else return std::make_shared<linear_solver_superlu<MATRIX, VECTOR>>(); # endif } else { if (md.is_coercive()) - return std::make_shared - <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>(); + return std::make_shared + <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>(); else { if (dim <= 2) - return std::make_shared - <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>(); - else - return std::make_shared - <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>(); + return std::make_shared + <linear_solver_gmres_preconditioned_ilut<MATRIX,VECTOR>>(); + else + return std::make_shared + <linear_solver_gmres_preconditioned_ilu<MATRIX,VECTOR>>(); } } #endif @@ -800,7 +804,7 @@ namespace getfem { return std::make_shared<linear_solver_mumps<MATRIX, VECTOR>>(); # else return std::make_shared - <linear_solver_distributed_mumps<MATRIX, VECTOR>>(); + <linear_solver_distributed_mumps<MATRIX, VECTOR>>(); # endif #else GMM_ASSERT1(false, "Mumps is not interfaced"); @@ -808,16 +812,16 @@ namespace getfem { } else if (bgeot::casecmp(name, "cg/ildlt") == 0) return std::make_shared - <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>(); + <linear_solver_cg_preconditioned_ildlt<MATRIX, VECTOR>>(); else if (bgeot::casecmp(name, "gmres/ilu") == 0) return std::make_shared - <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>(); + <linear_solver_gmres_preconditioned_ilu<MATRIX, VECTOR>>(); else if (bgeot::casecmp(name, "gmres/ilut") == 0) return std::make_shared - <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>(); + <linear_solver_gmres_preconditioned_ilut<MATRIX, VECTOR>>(); else if (bgeot::casecmp(name, "gmres/ilutp") == 0) return std::make_shared - <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>(); + <linear_solver_gmres_preconditioned_ilutp<MATRIX, VECTOR>>(); else if (bgeot::casecmp(name, "auto") == 0) return default_linear_solver<MATRIX, VECTOR>(md); else diff --git a/src/getfem_generic_assembly_interpolation.cc b/src/getfem_generic_assembly_interpolation.cc index 08ad1e8..60951cd 100644 --- a/src/getfem_generic_assembly_interpolation.cc +++ b/src/getfem_generic_assembly_interpolation.cc @@ -570,7 +570,7 @@ namespace getfem { var.varname, var.transname, 1); if (tree.root) ga_semantic_analysis(tree, local_workspace, dummy_mesh(), - 1, false, true); + 1, false, true); ga_compile_interpolation(pwi.workspace(), pwi.gis()); } } @@ -815,8 +815,8 @@ namespace getfem { m_x.trans_of_convex(adj_face.cv)); bool converged = true; gic.invert(ctx_x.xreal(), P_ref, converged); - bool is_in = (ctx_x.pgt()->convex_ref()->is_in(P_ref) < 1E-4); - GMM_ASSERT1(is_in && converged, "Geometric transformation inversion " + bool is_in = (ctx_x.pgt()->convex_ref()->is_in(P_ref) < 1E-4); + GMM_ASSERT1(is_in && converged, "Geometric transformation inversion " "has failed in neighbour transformation"); face_num = adj_face.f; cv = adj_face.cv; @@ -962,7 +962,7 @@ namespace getfem { public: virtual const mesh_region &give_region(const mesh &, - size_type, short_type) const + size_type, short_type) const { return region; } // virtual void init(const ga_workspace &workspace) const = 0; // virtual void finalize() const = 0; diff --git a/src/getfem_generic_assembly_semantic.cc b/src/getfem_generic_assembly_semantic.cc index 9f70cc3..905dc62 100644 --- a/src/getfem_generic_assembly_semantic.cc +++ b/src/getfem_generic_assembly_semantic.cc @@ -513,7 +513,8 @@ namespace getfem { pnode->test_function_type = t_type; for (size_type i = 0; i < n; ++i) for (size_type j = 0; j < n; ++j) - pnode->tensor()(i,j) = (i==j) ? scalar_type(1) : scalar_type(0); + pnode->tensor()(i,j) = (i==j) ? scalar_type(1) + : scalar_type(0); } } } diff --git a/src/getfem_generic_assembly_workspace.cc b/src/getfem_generic_assembly_workspace.cc index 506ad16..55a5bcf 100644 --- a/src/getfem_generic_assembly_workspace.cc +++ b/src/getfem_generic_assembly_workspace.cc @@ -300,7 +300,7 @@ namespace getfem { GMM_ASSERT1(false, "A secondary domain with the same " "name already exists"); if (transformations.find(name) != transformations.end()) - GMM_ASSERT1(name.compare("neighbour_elt"), "neighbour_elt is a " + GMM_ASSERT1(name != "neighbour_elt", "neighbour_elt is a " "reserved interpolate transformation name"); transformations[name] = ptrans; } @@ -446,9 +446,11 @@ namespace getfem { } if (!found) { - ind_tree = trees.size(); remain = false; + ind_tree = trees.size(); + remain = false; trees.push_back(tree_description()); - trees.back().mim = &mim; trees.back().m = &m; + trees.back().mim = &mim; + trees.back().m = &m; trees.back().rg = &rg; trees.back().secondary_domain = tree.secondary_domain; trees.back().ptree = new ga_tree; @@ -488,17 +490,6 @@ namespace getfem { } } - ga_workspace::m_tree::~m_tree() { if (ptree) delete ptree; } - ga_workspace::m_tree::m_tree(const m_tree& o) - : ptree(o.ptree), meshdim(o.meshdim), ignore_X(o.ignore_X) - { if (o.ptree) ptree = new ga_tree(*(o.ptree)); } - ga_workspace::m_tree &ga_workspace::m_tree::operator =(const m_tree& o) { - if (ptree) delete ptree; - ptree = o.ptree; meshdim = o.meshdim; ignore_X = o.ignore_X; - if (o.ptree) ptree = new ga_tree(*(o.ptree)); - return *this; - } - size_type ga_workspace::add_expression(const std::string &expr, const mesh_im &mim, const mesh_region &rg_, @@ -543,11 +534,11 @@ namespace getfem { } } - for (size_type i = 0; i < ltrees.size(); ++i) { - if (ltrees[i].root) { - // cout << "adding tree " << ga_tree_to_string(ltrees[i]) << endl; - max_order = std::max(ltrees[i].root->nb_test_functions(), max_order); - add_tree(ltrees[i], mim.linked_mesh(), mim, rg, expr, + for (ga_tree <ree : ltrees) { + if (ltree.root) { + // cout << "adding tree " << ga_tree_to_string(ltree) << endl; + max_order = std::max(ltree.root->nb_test_functions(), max_order); + add_tree(ltree, mim.linked_mesh(), mim, rg, expr, add_derivative_order, true, 0, ""); } } @@ -630,22 +621,18 @@ namespace getfem { size_type order) { bool islin = true; std::set<var_trans_pair> vll, dll; - for (size_type i = 0; i < vl.size(); ++i) - vll.insert(var_trans_pair(vl[i], "")); - for (size_type i = 0; i < dl.size(); ++i) - dll.insert(var_trans_pair(dl[i], "")); + for (const std::string &v : vl) vll.insert(var_trans_pair(v, "")); + for (const std::string &d : dl) dll.insert(var_trans_pair(d, "")); - for (size_type i = 0; i < trees.size(); ++i) { - ga_workspace::tree_description &td = trees[i]; + for (const ga_workspace::tree_description &td : trees) { std::set<var_trans_pair> dllaux; bool fv = ga_extract_variables(td.ptree->root, *this, *(td.m), dllaux, false); - if (td.order == order) { - for (std::set<var_trans_pair>::iterator it = dllaux.begin(); - it!=dllaux.end(); ++it) - dll.insert(*it); - } + if (td.order == order) + for (const auto &t : dllaux) + dll.insert(t); + switch (td.order) { case 0: break; case 1: @@ -659,7 +646,7 @@ namespace getfem { } bool found = false; for (const std::string &t : vl_test1) - if (td.name_test1.compare(t) == 0) + if (td.name_test1 == t) found = true; if (!found) vl_test1.push_back(td.name_test1); @@ -683,8 +670,8 @@ namespace getfem { } bool found = false; for (size_type j = 0; j < vl_test1.size(); ++j) - if ((td.name_test1.compare(vl_test1[j]) == 0) && - (td.name_test2.compare(vl_test2[j]) == 0)) + if ((td.name_test1 == vl_test1[j]) && + (td.name_test2 == vl_test2[j])) found = true; if (!found) { vl_test1.push_back(td.name_test1); @@ -697,11 +684,11 @@ namespace getfem { } vl.clear(); for (const auto &var : vll) - if (vl.size() == 0 || var.varname.compare(vl.back())) + if (vl.size() == 0 || var.varname != vl.back()) vl.push_back(var.varname); dl.clear(); for (const auto &var : dll) - if (dl.size() == 0 || var.varname.compare(dl.back())) + if (dl.size() == 0 || var.varname != dl.back()) dl.push_back(var.varname); return islin; @@ -926,9 +913,7 @@ namespace getfem { std::string ga_workspace::extract_constant_term(const mesh &m) { std::string constant_term; - for (size_type i = 0; i < trees.size(); ++i) { - ga_workspace::tree_description &td = trees[i]; - + for (const ga_workspace::tree_description &td : trees) { if (td.order == 1) { ga_tree local_tree = *(td.ptree); if (local_tree.root) @@ -950,8 +935,7 @@ namespace getfem { std::string ga_workspace::extract_order0_term() { std::string term; - for (size_type i = 0; i < trees.size(); ++i) { - ga_workspace::tree_description &td = trees[i]; + for (const ga_workspace::tree_description &td : trees) { if (td.order == 0) { ga_tree &local_tree = *(td.ptree); if (term.size()) @@ -970,9 +954,8 @@ namespace getfem { std::string ga_workspace::extract_order1_term(const std::string &varname) { std::string term; - for (size_type i = 0; i < trees.size(); ++i) { - ga_workspace::tree_description &td = trees[i]; - if (td.order == 1 && td.name_test1.compare(varname) == 0) { + for (const ga_workspace::tree_description &td : trees) { + if (td.order == 1 && td.name_test1 == varname) { ga_tree &local_tree = *(td.ptree); if (term.size()) term += "+("+ga_tree_to_string(local_tree)+")"; @@ -989,13 +972,12 @@ namespace getfem { std::string ga_workspace::extract_Neumann_term(const std::string &varname) { std::string result; - for (size_type i = 0; i < trees.size(); ++i) { - ga_workspace::tree_description &td = trees[i]; - if (td.order == 1 && td.name_test1.compare(varname) == 0) { + for (const ga_workspace::tree_description &td : trees) { + if (td.order == 1 && td.name_test1 == varname) { ga_tree &local_tree = *(td.ptree); if (local_tree.root) ga_extract_Neumann_term(local_tree, varname, *this, - local_tree.root, result); + local_tree.root, result); } } return result; diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc index edf8d96..e189e6e 100644 --- a/src/getfem_mesh.cc +++ b/src/getfem_mesh.cc @@ -363,7 +363,7 @@ namespace getfem { } bgeot::mesh_structure::sup_convex(ic); if (sup_points) - for (size_type ip = 0; ip < ipt.size(); ++ip) sup_point(ipt[ip]); + for (const size_type &ip : ipt) sup_point(ip); trans_exists.sup(ic); sup_convex_from_regions(ic); if (Bank_info.get()) Bank_sup_convex_from_green(ic); diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc index 00da7d2..2a45398 100644 --- a/src/getfem_model_solvers.cc +++ b/src/getfem_model_solvers.cc @@ -29,7 +29,7 @@ namespace getfem { static rmodel_plsolver_type rdefault_linear_solver(const model &md) { return default_linear_solver<model_real_sparse_matrix, model_real_plain_vector>(md); - } + } static cmodel_plsolver_type cdefault_linear_solver(const model &md) { return default_linear_solver<model_complex_sparse_matrix, @@ -49,7 +49,7 @@ namespace getfem { max_ratio_reached = false; } - double default_newton_line_search::next_try(void) { + double default_newton_line_search::next_try() { alpha_old = alpha; ++it; // alpha *= 0.5; if (alpha >= 0.4) alpha *= 0.5; else alpha *= alpha_mult; @@ -60,7 +60,7 @@ namespace getfem { // cout << "r = " << r << " alpha = " << alpha_old << " count_pat = " << count_pat << endl; if (!max_ratio_reached && r < first_res * alpha_max_ratio) { alpha_max_ratio_reached = alpha_old; r_max_ratio_reached = r; - it_max_ratio_reached = it; max_ratio_reached = true; + it_max_ratio_reached = it; max_ratio_reached = true; } if (max_ratio_reached && r < r_max_ratio_reached * 0.5 && @@ -71,9 +71,9 @@ namespace getfem { if (count == 0 || r < conv_r) { conv_r = r; conv_alpha = alpha_old; count = 1; } if (conv_r < first_res) ++count; - + if (r < first_res * alpha_min_ratio) - { count_pat = 0; return true; } + { count_pat = 0; return true; } if (count>=5 || (alpha < alpha_min && max_ratio_reached) || alpha<1e-15) { if (conv_r < first_res * 0.99) count_pat = 0; if (/*gmm::random() * 50. < -log(conv_alpha)-4.0 ||*/ count_pat >= 3) @@ -91,10 +91,10 @@ namespace getfem { /* ***************************************************************** */ template <typename MATRIX, typename VECTOR, typename PLSOLVER> - void compute_init_values(model &md, gmm::iteration &iter, - PLSOLVER lsolver, - abstract_newton_line_search &ls, const MATRIX &K, - const VECTOR &rhs) { + void compute_init_values(model &md, gmm::iteration &iter, + PLSOLVER lsolver, + abstract_newton_line_search &ls, const MATRIX &K, + const VECTOR &rhs) { VECTOR state(md.nb_dof()); md.from_variables(state); @@ -103,7 +103,7 @@ namespace getfem { scalar_type dt = md.get_time_step(); scalar_type ddt = md.get_init_time_step(); scalar_type t = md.get_time(); - + // Solve for ddt md.set_time_step(ddt); gmm::iteration iter1 = iter; @@ -148,9 +148,9 @@ namespace getfem { else { model_pb<MATRIX, VECTOR> mdpb(md, ls, state, rhs, K); if (dynamic_cast<newton_search_with_step_control *>(&ls)) - Newton_with_step_control(mdpb, iter, *lsolver); + Newton_with_step_control(mdpb, iter, *lsolver); else - classical_Newton(mdpb, iter, *lsolver); + classical_Newton(mdpb, iter, *lsolver); } md.to_variables(state); // copy the state vector into the model variables } diff --git a/src/getfem_models.cc b/src/getfem_models.cc index 6754dd0..e7c9315 100644 --- a/src/getfem_models.cc +++ b/src/getfem_models.cc @@ -705,7 +705,7 @@ namespace getfem { "Cannot explicitly resize a fem variable or data"); GMM_ASSERT1(variables[name].pim_data == 0, "Cannot explicitly resize an im data"); - GMM_ASSERT1(size, "Variable of null size are not allowed"); + GMM_ASSERT1(size, "Variables of null size are not allowed"); variables[name].qdims.resize(1); variables[name].qdims[0] = size; variables[name].set_size(); @@ -742,32 +742,32 @@ namespace getfem { void model::add_initialized_matrix_data(const std::string &name, const base_matrix &M) { - this->add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M), - gmm::mat_ncols(M))); - GMM_ASSERT1(!(this->is_complex()), "Sorry, complex version to be done"); - gmm::copy(M.as_vector(), this->set_real_variable(name)); + add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M), + gmm::mat_ncols(M))); + GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done"); + gmm::copy(M.as_vector(), set_real_variable(name)); } void model::add_initialized_matrix_data(const std::string &name, const base_complex_matrix &M) { - this->add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M), - gmm::mat_ncols(M))); - GMM_ASSERT1(!(this->is_complex()), "Sorry, complex version to be done"); - gmm::copy(M.as_vector(), this->set_complex_variable(name)); + add_fixed_size_data(name, bgeot::multi_index(gmm::mat_nrows(M), + gmm::mat_ncols(M))); + GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done"); + gmm::copy(M.as_vector(), set_complex_variable(name)); } void model::add_initialized_tensor_data(const std::string &name, const base_tensor &t) { - this->add_fixed_size_data(name, t.sizes(), 1); - GMM_ASSERT1(!(this->is_complex()), "Sorry, complex version to be done"); - gmm::copy(t.as_vector(), this->set_real_variable(name)); + add_fixed_size_data(name, t.sizes(), 1); + GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done"); + gmm::copy(t.as_vector(), set_real_variable(name)); } void model::add_initialized_tensor_data(const std::string &name, const base_complex_tensor &t) { - this->add_fixed_size_data(name, t.sizes(), 1); - GMM_ASSERT1(!(this->is_complex()), "Sorry, complex version to be done"); - gmm::copy(t.as_vector(), this->set_complex_variable(name)); + add_fixed_size_data(name, t.sizes(), 1); + GMM_ASSERT1(!(is_complex()), "Sorry, complex version to be done"); + gmm::copy(t.as_vector(), set_complex_variable(name)); } void model::add_im_data(const std::string &name, const im_data &im_data, @@ -1035,13 +1035,13 @@ namespace getfem { is_symmetric_ = is_symmetric_ && pbr->is_symmetric(); is_coercive_ = is_coercive_ && pbr->is_coercive(); - for (size_type i=0; i < varnames.size(); ++i) - GMM_ASSERT1(variables.find(varnames[i]) != variables.end(), - "Undefined model variable " << varnames[i]); + for (const auto &vname : varnames) + GMM_ASSERT1(variables.count(vname), + "Undefined model variable " << vname); // cout << "dl == " << datanames << endl; - for (size_type i=0; i < datanames.size(); ++i) - GMM_ASSERT1(variables.find(datanames[i]) != variables.end(), - "Undefined model data or variable " << datanames[i]); + for (const auto &dname : datanames) + GMM_ASSERT1(variables.count(dname), + "Undefined model data or variable " << dname); return ib; } @@ -1072,25 +1072,23 @@ namespace getfem { GMM_ASSERT1(valid_bricks[ib], "Inexistent brick"); touch_brick(ib); bricks[ib].vlist = vl; - for (size_type i=0; i < vl.size(); ++i) - GMM_ASSERT1(variables.find(vl[i]) != variables.end(), - "Undefined model variable " << vl[i]); + for (const auto &v : vl) + GMM_ASSERT1(variables.count(v), "Undefined model variable " << v); } void model::change_data_of_brick(size_type ib, const varnamelist &dl) { GMM_ASSERT1(valid_bricks[ib], "Inexistent brick"); touch_brick(ib); bricks[ib].dlist = dl; - for (size_type i=0; i < dl.size(); ++i) - GMM_ASSERT1(variables.find(dl[i]) != variables.end(), - "Undefined model variable " << dl[i]); + for (const auto &v : dl) + GMM_ASSERT1(variables.count(v), "Undefined model variable " << v); } void model::change_mims_of_brick(size_type ib, const mimlist &ml) { GMM_ASSERT1(valid_bricks[ib], "Inexistent brick"); touch_brick(ib); bricks[ib].mims = ml; - for (size_type i = 0; i < ml.size(); ++i) add_dependency(*(ml[i])); + for (const auto &mim : ml) add_dependency(*mim); } void model::change_update_flag_of_brick(size_type ib, bool flag) { @@ -1170,11 +1168,11 @@ namespace getfem { if (name_v.size()) { if (is_complex()) { - model_complex_plain_vector v0 = this->complex_variable(name_v); - gmm::copy(v0, this->set_complex_variable(name_previous_v)); + model_complex_plain_vector v0 = complex_variable(name_v); + gmm::copy(v0, set_complex_variable(name_previous_v)); } else { - const model_real_plain_vector &v0 = this->real_variable(name_v); - gmm::copy(v0, this->set_real_variable(name_previous_v)); + const model_real_plain_vector &v0 = real_variable(name_v); + gmm::copy(v0, set_real_variable(name_previous_v)); } } } @@ -1803,16 +1801,14 @@ namespace getfem { void model::first_iter() { context_check(); if (act_size_to_be_done) actualize_sizes(); - for (VAR_SET::iterator it = variables.begin(); it != variables.end(); ++it) - it->second.clear_temporaries(); + for (auto && v : variables) v.second.clear_temporaries(); set_dispatch_coeff(); for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) { brick_description &brick = bricks[ib]; - bool cplx = is_complex() && brick.pbr->is_complex(); if (brick.pdispatch) { - if (cplx) + if (is_complex() && brick.pbr->is_complex()) brick.pdispatch->next_complex_iter(*this, ib, brick.vlist, brick.dlist, brick.cmatlist, brick.cveclist, @@ -1831,9 +1827,8 @@ namespace getfem { for (dal::bv_visitor ib(active_bricks); !ib.finished(); ++ib) { brick_description &brick = bricks[ib]; - bool cplx = is_complex() && brick.pbr->is_complex(); if (brick.pdispatch) { - if (cplx) + if (is_complex() && brick.pbr->is_complex()) brick.pdispatch->next_complex_iter(*this, ib, brick.vlist, brick.dlist, brick.cmatlist, brick.cveclist, @@ -1845,18 +1840,14 @@ namespace getfem { } } - for (VAR_SET::iterator it = variables.begin(); it != variables.end(); - ++it) { - for (size_type i = 1; i < it->second.n_iter; ++i) { + for (auto &&v : variables) + for (size_type i = 1; i < v.second.n_iter; ++i) { if (is_complex()) - gmm::copy(it->second.complex_value[i-1], - it->second.complex_value[i]); + gmm::copy(v.second.complex_value[i-1], v.second.complex_value[i]); else - gmm::copy(it->second.real_value[i-1], - it->second.real_value[i]); - it->second.v_num_data[i] = act_counter(); + gmm::copy(v.second.real_value[i-1], v.second.real_value[i]); + v.second.v_num_data[i] = act_counter(); } - } } bool model::is_var_newer_than_brick(const std::string &varname, @@ -2731,18 +2722,18 @@ namespace getfem { const model_real_plain_vector & model::real_variable(const std::string &name) const { - if (is_old(name)) return real_variable(no_old_prefix_name(name), 1); - else return real_variable(name, size_type(-1)); + return is_old(name) ? real_variable(no_old_prefix_name(name), 1) + : real_variable(name, size_type(-1)); } const model_real_plain_vector & model::real_variable(const std::string &name, size_type niter) const { GMM_ASSERT1(!complex_version, "This model is a complex one"); - GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination with" - " variable version"); + GMM_ASSERT1(!is_old(name), "Please don't use Old_ prefix in combination " + "with variable version"); context_check(); auto it = variables.find(name); - GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name); + GMM_ASSERT1(it != variables.end(), "Undefined variable " << name); if (act_size_to_be_done && it->second.is_fem_dofs) { if (it->second.filter != VDESCRFILTER_NO) actualize_sizes(); @@ -2757,8 +2748,8 @@ namespace getfem { const model_complex_plain_vector & model::complex_variable(const std::string &name) const { - if (is_old(name)) return complex_variable(no_old_prefix_name(name), 1); - else return complex_variable(name, size_type(-1)); + return is_old(name) ? complex_variable(no_old_prefix_name(name), 1) + : complex_variable(name, size_type(-1)); } const model_complex_plain_vector & @@ -2784,8 +2775,8 @@ namespace getfem { model_real_plain_vector & model::set_real_variable(const std::string &name) const { - if (is_old(name)) return set_real_variable(no_old_prefix_name(name), 1); - else return set_real_variable(name, size_type(-1)); + return is_old(name) ? set_real_variable(no_old_prefix_name(name), 1) + : set_real_variable(name, size_type(-1)); } @@ -2811,10 +2802,10 @@ namespace getfem { return it->second.real_value[niter]; } -model_complex_plain_vector & + model_complex_plain_vector & model::set_complex_variable(const std::string &name) const { - if (is_old(name)) return set_complex_variable(no_old_prefix_name(name), 1); - else return set_complex_variable(name, size_type(-1)); + return is_old(name) ? set_complex_variable(no_old_prefix_name(name), 1) + : set_complex_variable(name, size_type(-1)); } model_complex_plain_vector & diff --git a/src/gmm/gmm_vector.h b/src/gmm/gmm_vector.h index e69931d..dbc4029 100644 --- a/src/gmm/gmm_vector.h +++ b/src/gmm/gmm_vector.h @@ -53,7 +53,7 @@ namespace gmm { V *pm; size_type l; - + public : operator T() const { return pm->r(l); } @@ -92,7 +92,7 @@ namespace gmm { V *pm; size_type l; - + public : operator std::complex<T>() const { return pm->r(l); } @@ -139,21 +139,21 @@ namespace gmm { { return std::complex<T>(*this)* v; } std::complex<T> operator /(std::complex<T> v) { return std::complex<T>(*this)/ v; } - }; + }; + - template<typename T, typename V> inline bool operator ==(T v, const ref_elt_vector<T, V> &re) { return (v==T(re)); } template<typename T, typename V> inline bool operator !=(T v, const ref_elt_vector<T, V> &re) { return (v!=T(re)); } template<typename T, typename V> inline - T &operator +=(T &v, const ref_elt_vector<T, V> &re) + T &operator +=(T &v, const ref_elt_vector<T, V> &re) { v += T(re); return v; } template<typename T, typename V> inline T &operator -=(T &v, const ref_elt_vector<T, V> &re) { v -= T(re); return v; } template<typename T, typename V> inline - T &operator *=(T &v, const ref_elt_vector<T, V> &re) + T &operator *=(T &v, const ref_elt_vector<T, V> &re) { v *= T(re); return v; } template<typename T, typename V> inline T &operator /=(T &v, const ref_elt_vector<T, V> &re) @@ -224,7 +224,7 @@ namespace gmm { size_type i; // Current index. T* p; // Pointer to the current position. dsvector<T> *v; // Pointer to the vector. - + typedef T value_type; typedef value_type* pointer; typedef const value_type* const_pointer; @@ -233,20 +233,20 @@ namespace gmm { typedef ptrdiff_t difference_type; typedef std::bidirectional_iterator_tag iterator_category; typedef dsvector_iterator<T> iterator; - + reference operator *() const { return *p; } pointer operator->() const { return &(operator*()); } iterator &operator ++() { for (size_type k = (i & 15); k < 15; ++k) - { ++p; ++i; if (*p != T(0)) return *this; } + { ++p; ++i; if (*p != T(0)) return *this; } v->next_pos(*(const_cast<const_pointer *>(&(p))), i); return *this; } iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; } iterator &operator --() { for (size_type k = (i & 15); k > 0; --k) - { --p; --i; if (*p != T(0)) return *this; } + { --p; --i; if (*p != T(0)) return *this; } v->previous_pos(p, i); return *this; } @@ -256,10 +256,10 @@ namespace gmm { { return (i == it.i && p == it.p && v == it.v); } bool operator !=(const iterator &it) const { return !(it == *this); } - - size_type index(void) const { return i; } - dsvector_iterator(void) : i(size_type(-1)), p(0), v(0) {} + size_type index() const { return i; } + + dsvector_iterator() : i(size_type(-1)), p(0), v(0) {} dsvector_iterator(dsvector<T> &w) : i(size_type(-1)), p(0), v(&w) {}; }; @@ -268,7 +268,7 @@ namespace gmm { size_type i; // Current index. const T* p; // Pointer to the current position. const dsvector<T> *v; // Pointer to the vector. - + typedef T value_type; typedef const value_type* pointer; typedef const value_type& reference; @@ -276,19 +276,19 @@ namespace gmm { typedef ptrdiff_t difference_type; typedef std::bidirectional_iterator_tag iterator_category; typedef dsvector_const_iterator<T> iterator; - + reference operator *() const { return *p; } pointer operator->() const { return &(operator*()); } iterator &operator ++() { for (size_type k = (i & 15); k < 15; ++k) - { ++p; ++i; if (*p != T(0)) return *this; } + { ++p; ++i; if (*p != T(0)) return *this; } v->next_pos(p, i); return *this; } iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; } iterator &operator --() { for (size_type k = (i & 15); k > 0; --k) - { --p; --i; if (*p != T(0)) return *this; } + { --p; --i; if (*p != T(0)) return *this; } v->previous_pos(p, i); return *this; } @@ -298,17 +298,17 @@ namespace gmm { { return (i == it.i && p == it.p && v == it.v); } bool operator !=(const iterator &it) const { return !(it == *this); } - - size_type index(void) const { return i; } - dsvector_const_iterator(void) : i(size_type(-1)), p(0) {} + size_type index() const { return i; } + + dsvector_const_iterator() : i(size_type(-1)), p(0) {} dsvector_const_iterator(const dsvector_iterator<T> &it) : i(it.i), p(it.p), v(it.v) {} dsvector_const_iterator(const dsvector<T> &w) : i(size_type(-1)), p(0), v(&w) {}; }; - + /** Sparse vector built on distribution sort principle. Read and write access have a constant complexity depending only on the @@ -323,7 +323,7 @@ namespace gmm { typedef const T * const_pointer; typedef void * void_pointer; typedef const void * const_void_pointer; - + protected: size_type n; // Potential vector size size_type depth; // Number of row of pointer arrays @@ -337,10 +337,10 @@ namespace gmm { void_pointer p = root_ptr; if (!p) return 0; for (size_type k = 0; k < depth; ++k) { - p = ((void **)(p))[(i & my_mask) >> my_shift]; - if (!p) return 0; - my_mask = (my_mask >> 4); - my_shift -= 4; + p = ((void **)(p))[(i & my_mask) >> my_shift]; + if (!p) return 0; + my_mask = (my_mask >> 4); + my_shift -= 4; } GMM_ASSERT1(my_shift == 0, "internal error"); GMM_ASSERT1(my_mask == 15, "internal error"); @@ -351,32 +351,32 @@ namespace gmm { GMM_ASSERT1(i < n, "index " << i << " out of range (size " << n << ")"); size_type my_mask = mask, my_shift = shift; if (!root_ptr) { - if (depth) { - root_ptr = new void_pointer[16]; - std::memset(root_ptr, 0, 16*sizeof(void_pointer)); - } else { - root_ptr = new T[16]; - for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0); - } + if (depth) { + root_ptr = new void_pointer[16]; + std::memset(root_ptr, 0, 16*sizeof(void_pointer)); + } else { + root_ptr = new T[16]; + for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0); + } } void_pointer p = root_ptr; for (size_type k = 0; k < depth; ++k) { - size_type j = (i & my_mask) >> my_shift; - void_pointer q = ((void_pointer *)(p))[j]; - if (!q) { - if (k+1 != depth) { - q = new void_pointer[16]; - std::memset(q, 0, 16*sizeof(void_pointer)); - } else { - q = new T[16]; - for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0); - } - ((void_pointer *)(p))[j] = q; - } - p = q; - my_mask = (my_mask >> 4); - my_shift -= 4; + size_type j = (i & my_mask) >> my_shift; + void_pointer q = ((void_pointer *)(p))[j]; + if (!q) { + if (k+1 != depth) { + q = new void_pointer[16]; + std::memset(q, 0, 16*sizeof(void_pointer)); + } else { + q = new T[16]; + for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0); + } + ((void_pointer *)(p))[j] = q; + } + p = q; + my_mask = (my_mask >> 4); + my_shift -= 4; } GMM_ASSERT1(my_shift == 0, "internal error"); GMM_ASSERT1(my_mask == 15, "internal error " << my_mask); @@ -392,65 +392,65 @@ namespace gmm { void rec_del(void_pointer p, size_type my_depth) { if (my_depth) { - for (size_type k = 0; k < 16; ++k) - if (((void_pointer *)(p))[k]) - rec_del(((void_pointer *)(p))[k], my_depth-1); - delete[] ((void_pointer *)(p)); + for (size_type k = 0; k < 16; ++k) + if (((void_pointer *)(p))[k]) + rec_del(((void_pointer *)(p))[k], my_depth-1); + delete[] ((void_pointer *)(p)); } else { - delete[] ((T *)(p)); + delete[] ((T *)(p)); } } void rec_clean(void_pointer p, size_type my_depth, double eps) { if (my_depth) { - for (size_type k = 0; k < 16; ++k) - if (((void_pointer *)(p))[k]) - rec_clean(((void_pointer *)(p))[k], my_depth-1, eps); + for (size_type k = 0; k < 16; ++k) + if (((void_pointer *)(p))[k]) + rec_clean(((void_pointer *)(p))[k], my_depth-1, eps); } else { - for (size_type k = 0; k < 16; ++k) - if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0); + for (size_type k = 0; k < 16; ++k) + if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0); } } void rec_clean_i(void_pointer p, size_type my_depth, size_type my_mask, - size_type i, size_type base) { + size_type i, size_type base) { if (my_depth) { - my_mask = (my_mask >> 4); - for (size_type k = 0; k < 16; ++k) - if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i) - rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask, - i, base + k*(my_mask+1)); + my_mask = (my_mask >> 4); + for (size_type k = 0; k < 16; ++k) + if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i) + rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask, + i, base + k*(my_mask+1)); } else { - for (size_type k = 0; k < 16; ++k) - if (base+k > i) ((T *)(p))[k] = T(0); + for (size_type k = 0; k < 16; ++k) + if (base+k > i) ((T *)(p))[k] = T(0); } } - - + + size_type rec_nnz(void_pointer p, size_type my_depth) const { size_type nn = 0; if (my_depth) { - for (size_type k = 0; k < 16; ++k) - if (((void_pointer *)(p))[k]) - nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1); + for (size_type k = 0; k < 16; ++k) + if (((void_pointer *)(p))[k]) + nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1); } else { - for (size_type k = 0; k < 16; ++k) - if (((const T *)(p))[k] != T(0)) nn++; + for (size_type k = 0; k < 16; ++k) + if (((const T *)(p))[k] != T(0)) nn++; } return nn; } void copy_rec(void_pointer &p, const_void_pointer q, size_type my_depth) { if (my_depth) { - p = new void_pointer[16]; - std::memset(p, 0, 16*sizeof(void_pointer)); - for (size_type l = 0; l < 16; ++l) - if (((const const_void_pointer *)(q))[l]) - copy_rec(((void_pointer *)(p))[l], - ((const const_void_pointer *)(q))[l], my_depth-1); + p = new void_pointer[16]; + std::memset(p, 0, 16*sizeof(void_pointer)); + for (size_type l = 0; l < 16; ++l) + if (((const const_void_pointer *)(q))[l]) + copy_rec(((void_pointer *)(p))[l], + ((const const_void_pointer *)(q))[l], my_depth-1); } else { - p = new T[16]; - for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l]; + p = new T[16]; + for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l]; } } @@ -462,75 +462,75 @@ namespace gmm { } void next_pos_rec(void_pointer p, size_type my_depth, size_type my_mask, - const_pointer &pp, size_type &i, size_type base) const { + const_pointer &pp, size_type &i, size_type base) const { size_type ii = i; if (my_depth) { - my_mask = (my_mask >> 4); - for (size_type k = 0; k < 16; ++k) - if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) { - next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask, - pp, i, base + k*(my_mask+1)); - if (i != size_type(-1)) return; else i = ii; - } - i = size_type(-1); pp = 0; + my_mask = (my_mask >> 4); + for (size_type k = 0; k < 16; ++k) + if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) { + next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask, + pp, i, base + k*(my_mask+1)); + if (i != size_type(-1)) return; else i = ii; + } + i = size_type(-1); pp = 0; } else { - for (size_type k = 0; k < 16; ++k) - if (base+k > i && ((const_pointer)(p))[k] != T(0)) - { i = base+k; pp = &(((const_pointer)(p))[k]); return; } - i = size_type(-1); pp = 0; + for (size_type k = 0; k < 16; ++k) + if (base+k > i && ((const_pointer)(p))[k] != T(0)) + { i = base+k; pp = &(((const_pointer)(p))[k]); return; } + i = size_type(-1); pp = 0; } } void previous_pos_rec(void_pointer p, size_type my_depth, size_type my_mask, - const_pointer &pp, size_type &i, - size_type base) const { + const_pointer &pp, size_type &i, + size_type base) const { size_type ii = i; if (my_depth) { - my_mask = (my_mask >> 4); - for (size_type k = 15; k != size_type(-1); --k) - if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) { - previous_pos_rec(((void_pointer *)(p))[k], my_depth-1, - my_mask, pp, i, base + k*(my_mask+1)); - if (i != size_type(-1)) return; else i = ii; - } - i = size_type(-1); pp = 0; + my_mask = (my_mask >> 4); + for (size_type k = 15; k != size_type(-1); --k) + if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) { + previous_pos_rec(((void_pointer *)(p))[k], my_depth-1, + my_mask, pp, i, base + k*(my_mask+1)); + if (i != size_type(-1)) return; else i = ii; + } + i = size_type(-1); pp = 0; } else { - for (size_type k = 15; k != size_type(-1); --k) - if (base+k < i && ((const_pointer)(p))[k] != T(0)) - { i = base+k; pp = &(((const_pointer)(p))[k]); return; } - i = size_type(-1); pp = 0; + for (size_type k = 15; k != size_type(-1); --k) + if (base+k < i && ((const_pointer)(p))[k] != T(0)) + { i = base+k; pp = &(((const_pointer)(p))[k]); return; } + i = size_type(-1); pp = 0; } } - - + + public: void clean(double eps) { if (root_ptr) rec_clean(root_ptr, depth); } void resize(size_type n_) { if (n_ != n) { - n = n_; - if (n_ < n) { // Depth unchanged (a choice) - if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0); - } else { - // may change the depth (add some levels) - size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_; - while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; } - my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth; - if (my_depth > depth || depth == 0) { - if (root_ptr) { - for (size_type k = depth; k < my_depth; ++k) { - void_pointer *q = new void_pointer [16]; - std::memset(q, 0, 16*sizeof(void_pointer)); - q[0] = root_ptr; root_ptr = q; - } - } - mask = my_mask; depth = my_depth; shift = my_shift; - } - } + n = n_; + if (n_ < n) { // Depth unchanged (a choice) + if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0); + } else { + // may change the depth (add some levels) + size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_; + while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; } + my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth; + if (my_depth > depth || depth == 0) { + if (root_ptr) { + for (size_type k = depth; k < my_depth; ++k) { + void_pointer *q = new void_pointer [16]; + std::memset(q, 0, 16*sizeof(void_pointer)); + q[0] = root_ptr; root_ptr = q; + } + } + mask = my_mask; depth = my_depth; shift = my_shift; + } + } } } - - void clear(void) { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; } - + + void clear() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; } + void next_pos(const_pointer &pp, size_type &i) const { if (!root_ptr || i >= n) { pp = 0, i = size_type(-1); return; } next_pos_rec(root_ptr, depth, mask, pp, i, 0); @@ -542,29 +542,29 @@ namespace gmm { previous_pos_rec(root_ptr, depth, mask, pp, i, 0); } - iterator begin(void) { - iterator it(*this); + iterator begin() { + iterator it(*this); if (n && root_ptr) { - it.i = 0; it.p = const_cast<T *>(read_access(0)); - if (!(it.p) || *(it.p) == T(0)) - next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i); + it.i = 0; it.p = const_cast<T *>(read_access(0)); + if (!(it.p) || *(it.p) == T(0)) + next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i); } return it; } - iterator end(void) { return iterator(*this); } + iterator end() { return iterator(*this); } - const_iterator begin(void) const { + const_iterator begin() const { const_iterator it(*this); if (n && root_ptr) { - it.i = 0; it.p = read_access(0); - if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i); + it.i = 0; it.p = read_access(0); + if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i); } return it; } - const_iterator end(void) const { return const_iterator(*this); } - + const_iterator end() const { return const_iterator(*this); } + inline ref_elt_vector<T, dsvector<T> > operator [](size_type c) { return ref_elt_vector<T, dsvector<T> >(this, c); } @@ -580,22 +580,22 @@ namespace gmm { { const T *p = read_access(c); if (p) return *p; else return T(0); } inline T operator [](size_type c) const { return r(c); } - - size_type nnz(void) const + + size_type nnz() const { if (root_ptr) return rec_nnz(root_ptr, depth); else return 0; } - size_type size(void) const { return n; } + size_type size() const { return n; } void swap(dsvector<T> &v) { std::swap(n, v.n); std::swap(root_ptr, v.root_ptr); std::swap(depth, v.depth); std::swap(shift, v.shift); std::swap(mask, v.mask); } - + /* Constructors */ dsvector(const dsvector<T> &v) { init(0); copy(v); } dsvector<T> &operator =(const dsvector<T> &v) { copy(v); return *this; } explicit dsvector(size_type l){ init(l); } - dsvector(void) { init(0); } + dsvector() { init(0); } ~dsvector() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; } }; @@ -621,10 +621,10 @@ namespace gmm { { o->clear(); } static void do_clear(this_type &v) { v.clear(); } static value_type access(const origin_type *o, const const_iterator &, - const const_iterator &, size_type i) + const const_iterator &, size_type i) { return (*o)[i]; } static reference access(origin_type *o, const iterator &, const iterator &, - size_type i) + size_type i) { return (*o)[i]; } static void resize(this_type &v, size_type n) { v.resize(n); } }; @@ -635,12 +635,12 @@ namespace gmm { /******* Optimized operations for dsvector<T> ****************************/ template <typename T> inline void copy(const dsvector<T> &v1, - dsvector<T> &v2) { + dsvector<T> &v2) { GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch"); v2 = v1; } template <typename T> inline void copy(const dsvector<T> &v1, - const dsvector<T> &v2) { + const dsvector<T> &v2) { GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch"); v2 = const_cast<dsvector<T> &>(v1); } @@ -655,23 +655,23 @@ namespace gmm { } template <typename T> inline void copy(const simple_vector_ref<const dsvector<T> *> &v1, - dsvector<T> &v2) + dsvector<T> &v2) { copy(*(v1.origin), v2); } template <typename T> inline void copy(const simple_vector_ref<dsvector<T> *> &v1, dsvector<T> &v2) { copy(*(v1.origin), v2); } template <typename T> inline void copy(const simple_vector_ref<dsvector<T> *> &v1, - const simple_vector_ref<dsvector<T> *> &v2) + const simple_vector_ref<dsvector<T> *> &v2) { copy(*(v1.origin), v2); } template <typename T> inline void copy(const simple_vector_ref<const dsvector<T> *> &v1, - const simple_vector_ref<dsvector<T> *> &v2) + const simple_vector_ref<dsvector<T> *> &v2) { copy(*(v1.origin), v2); } - + template <typename T> inline size_type nnz(const dsvector<T>& l) { return l.nnz(); } - + /*************************************************************************/ /* */ /* Class wsvector: sparse vector optimized for random write operations, */ @@ -679,7 +679,7 @@ namespace gmm { /* Based on std::map */ /* */ /*************************************************************************/ - + template<typename T> struct wsvector_iterator : public std::map<size_type, T>::iterator { typedef typename std::map<size_type, T>::iterator base_it_type; @@ -689,12 +689,12 @@ namespace gmm { // typedef size_t size_type; typedef ptrdiff_t difference_type; typedef std::bidirectional_iterator_tag iterator_category; - + reference operator *() const { return (base_it_type::operator*()).second; } pointer operator->() const { return &(operator*()); } - size_type index(void) const { return (base_it_type::operator*()).first; } + size_type index() const { return (base_it_type::operator*()).first; } - wsvector_iterator(void) {} + wsvector_iterator() {} wsvector_iterator(const base_it_type &it) : base_it_type(it) {} }; @@ -707,12 +707,12 @@ namespace gmm { // typedef size_t size_type; typedef ptrdiff_t difference_type; typedef std::bidirectional_iterator_tag iterator_category; - + reference operator *() const { return (base_it_type::operator*()).second; } pointer operator->() const { return &(operator*()); } - size_type index(void) const { return (base_it_type::operator*()).first; } + size_type index() const { return (base_it_type::operator*()).first; } - wsvector_const_iterator(void) {} + wsvector_const_iterator() {} wsvector_const_iterator(const wsvector_iterator<T> &it) : base_it_type(it) {} wsvector_const_iterator(const base_it_type &it) : base_it_type(it) {} @@ -725,7 +725,7 @@ namespace gmm { */ template<typename T> class wsvector : public std::map<size_type, T> { public: - + typedef typename std::map<int, T>::size_type size_type; typedef std::map<size_type, T> base_type; typedef typename base_type::iterator iterator; @@ -733,11 +733,11 @@ namespace gmm { protected: size_type nbl; - + public: void clean(double eps); void resize(size_type); - + inline ref_elt_vector<T, wsvector<T> > operator [](size_type c) { return ref_elt_vector<T, wsvector<T> >(this, c); } @@ -750,9 +750,9 @@ namespace gmm { inline void wa(size_type c, const T &e) { GMM_ASSERT2(c < nbl, "out of range"); if (e != T(0)) { - iterator it = this->lower_bound(c); - if (it != this->end() && it->first == c) it->second += e; - else base_type::operator [](c) = e; + iterator it = this->lower_bound(c); + if (it != this->end() && it->first == c) it->second += e; + else base_type::operator [](c) = e; } } @@ -764,18 +764,18 @@ namespace gmm { } inline T operator [](size_type c) const { return r(c); } - - size_type nb_stored(void) const { return base_type::size(); } - size_type size(void) const { return nbl; } + + size_type nb_stored() const { return base_type::size(); } + size_type size() const { return nbl; } void swap(wsvector<T> &v) { std::swap(nbl, v.nbl); std::map<size_type, T>::swap(v); } - + /* Constructors */ void init(size_type l) { nbl = l; this->clear(); } explicit wsvector(size_type l){ init(l); } - wsvector(void) { init(0); } + wsvector() { init(0); } }; template<typename T> void wsvector<T>::clean(double eps) { @@ -815,10 +815,10 @@ namespace gmm { { o->clear(); } static void do_clear(this_type &v) { v.clear(); } static value_type access(const origin_type *o, const const_iterator &, - const const_iterator &, size_type i) + const const_iterator &, size_type i) { return (*o)[i]; } static reference access(origin_type *o, const iterator &, const iterator &, - size_type i) + size_type i) { return (*o)[i]; } static void resize(this_type &v, size_type n) { v.resize(n); } }; @@ -829,7 +829,7 @@ namespace gmm { /******* Optimized BLAS for wsvector<T> **********************************/ template <typename T> inline void copy(const wsvector<T> &v1, - wsvector<T> &v2) { + wsvector<T> &v2) { GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch"); v2 = v1; } @@ -844,7 +844,7 @@ namespace gmm { } template <typename T> inline void copy(const simple_vector_ref<const wsvector<T> *> &v1, - wsvector<T> &v2) + wsvector<T> &v2) { copy(*(v1.origin), v2); } template <typename T> inline void copy(const simple_vector_ref<wsvector<T> *> &v1, wsvector<T> &v2) @@ -853,9 +853,9 @@ namespace gmm { template <typename T> inline void clean(wsvector<T> &v, double eps) { typedef typename number_traits<T>::magnitude_type R; typename wsvector<T>::iterator it = v.begin(), ite = v.end(), itc; - while (it != ite) + while (it != ite) if (gmm::abs((*it).second) <= R(eps)) - { itc=it; ++it; v.erase(itc); } else ++it; + { itc=it; ++it; v.erase(itc); } else ++it; } template <typename T> @@ -881,7 +881,7 @@ namespace gmm { size_type c; T e; /* e is initialized by default to avoid some false warnings of valgrind. (from http://valgrind.org/docs/manual/mc-manual.html: - + When memory is read into the CPU's floating point registers, the relevant V bits are read from memory and they are immediately checked. If any are invalid, an uninitialised value error is @@ -890,7 +890,7 @@ namespace gmm { does not have to track the validity status of the floating-point registers. */ - elt_rsvector_(void) : e(0) {} + elt_rsvector_() : e(0) {} elt_rsvector_(size_type cc) : c(cc), e(0) {} elt_rsvector_(size_type cc, const T &ee) : c(cc), e(ee) {} bool operator < (const elt_rsvector_ &a) const { return c < a.c; } @@ -921,8 +921,8 @@ namespace gmm { bool operator ==(const iterator &i) const { return it == i.it; } bool operator !=(const iterator &i) const { return !(i == *this); } - size_type index(void) const { return it->c; } - rsvector_iterator(void) {} + size_type index() const { return it->c; } + rsvector_iterator() {} rsvector_iterator(const IT &i) : it(i) {} }; @@ -940,7 +940,7 @@ namespace gmm { reference operator *() const { return it->e; } pointer operator->() const { return &(operator*()); } - size_type index(void) const { return it->c; } + size_type index() const { return it->c; } iterator &operator ++() { ++it; return *this; } iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; } @@ -950,18 +950,18 @@ namespace gmm { bool operator ==(const iterator &i) const { return it == i.it; } bool operator !=(const iterator &i) const { return !(i == *this); } - rsvector_const_iterator(void) {} + rsvector_const_iterator() {} rsvector_const_iterator(const rsvector_iterator<T> &i) : it(i.it) {} rsvector_const_iterator(const IT &i) : it(i) {} }; /** sparse vector built upon std::vector. Read access is fast, - but insertion is O(n) + but insertion is O(n) */ template<typename T> class rsvector : public std::vector<elt_rsvector_<T> > { public: - + typedef std::vector<elt_rsvector_<T> > base_type_; typedef typename base_type_::iterator iterator; typedef typename base_type_::const_iterator const_iterator; @@ -969,14 +969,14 @@ namespace gmm { typedef T value_type; protected: - size_type nbl; /* size of the vector. */ - + size_type nbl; // size of the vector + public: void sup(size_type j); void base_resize(size_type n) { base_type_::resize(n); } void resize(size_type); - + ref_elt_vector<T, rsvector<T> > operator [](size_type c) { return ref_elt_vector<T, rsvector<T> >(this, c); } @@ -986,16 +986,16 @@ namespace gmm { void swap_indices(size_type i, size_type j); inline T operator [](size_type c) const { return r(c); } - - size_type nb_stored(void) const { return base_type_::size(); } - size_type size(void) const { return nbl; } - void clear(void) { base_type_::resize(0); } + + size_type nb_stored() const { return base_type_::size(); } + size_type size() const { return nbl; } + void clear() { base_type_::resize(0); } void swap(rsvector<T> &v) { std::swap(nbl, v.nbl); std::vector<elt_rsvector_<T> >::swap(v); } /* Constructeurs */ explicit rsvector(size_type l) : nbl(l) { } - rsvector(void) : nbl(0) { } + rsvector() : nbl(0) { } }; template <typename T> @@ -1012,18 +1012,18 @@ namespace gmm { switch (situation) { case 1 : a = *iti; a.c = j; it = iti; ++it; ite = this->end(); - for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it; - *iti = a; - break; + for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it; + *iti = a; + break; case 2 : a = *itj; a.c = i; it = itj; ite = this->begin(); - if (it != ite) { - --it; - while (it->c >= i) { *itj = *it; --itj; if (it==ite) break; --it; } - } - *itj = a; - break; + if (it != ite) { + --it; + while (it->c >= i) { *itj = *it; --itj; if (it==ite) break; --it; } + } + *itj = a; + break; case 3 : std::swap(iti->e, itj->e); - break; + break; } } } @@ -1033,8 +1033,8 @@ namespace gmm { elt_rsvector_<T> ev(j); iterator it = std::lower_bound(this->begin(), this->end(), ev); if (it != this->end() && it->c == j) { - for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1); - base_resize(nb_stored()-1); + for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1); + base_resize(nb_stored()-1); } } } @@ -1042,7 +1042,7 @@ namespace gmm { template<typename T> void rsvector<T>::resize(size_type n) { if (n < nbl) { for (size_type i = 0; i < nb_stored(); ++i) - if (base_type_::operator[](i).c >= n) { base_resize(i); break; } + if (base_type_::operator[](i).c >= n) { base_resize(i); break; } } nbl = n; } @@ -1053,24 +1053,24 @@ namespace gmm { else { elt_rsvector_<T> ev(c, e); if (nb_stored() == 0) { - base_type_::push_back(ev); + base_type_::push_back(ev); } else { - iterator it = std::lower_bound(this->begin(), this->end(), ev); - if (it != this->end() && it->c == c) it->e = e; - else { - size_type ind = it - this->begin(), nb = this->nb_stored(); + iterator it = std::lower_bound(this->begin(), this->end(), ev); + if (it != this->end() && it->c == c) it->e = e; + else { + size_type ind = it - this->begin(), nb = this->nb_stored(); if (nb - ind > 1100) GMM_WARNING2("Inefficient addition of element in rsvector with " << this->nb_stored() - ind << " non-zero entries"); - base_type_::push_back(ev); - if (ind != nb) { - it = this->begin() + ind; - iterator ite = this->end(); --ite; iterator itee = ite; - for (; ite != it; --ite) { --itee; *ite = *itee; } - *it = ev; - } - } + base_type_::push_back(ev); + if (ind != nb) { + it = this->begin() + ind; + iterator ite = this->end(); --ite; iterator itee = ite; + for (; ite != it; --ite) { --itee; *ite = *itee; } + *it = ev; + } + } } } } @@ -1080,31 +1080,31 @@ namespace gmm { if (e != T(0)) { elt_rsvector_<T> ev(c, e); if (nb_stored() == 0) { - base_type_::push_back(ev); + base_type_::push_back(ev); } else { - iterator it = std::lower_bound(this->begin(), this->end(), ev); - if (it != this->end() && it->c == c) it->e += e; - else { - size_type ind = it - this->begin(), nb = this->nb_stored(); + iterator it = std::lower_bound(this->begin(), this->end(), ev); + if (it != this->end() && it->c == c) it->e += e; + else { + size_type ind = it - this->begin(), nb = this->nb_stored(); if (nb - ind > 1100) GMM_WARNING2("Inefficient addition of element in rsvector with " << this->nb_stored() - ind << " non-zero entries"); - base_type_::push_back(ev); - if (ind != nb) { - it = this->begin() + ind; - iterator ite = this->end(); --ite; iterator itee = ite; - for (; ite != it; --ite) { --itee; *ite = *itee; } - *it = ev; - } - } + base_type_::push_back(ev); + if (ind != nb) { + it = this->begin() + ind; + iterator ite = this->end(); --ite; iterator itee = ite; + for (; ite != it; --ite) { --itee; *ite = *itee; } + *it = ev; + } + } } } } - + template <typename T> T rsvector<T>::r(size_type c) const { - GMM_ASSERT2(c < nbl, "out of range. Index " << c - << " for a length of " << nbl); + GMM_ASSERT2(c < nbl, "out of range. Index " << c + << " for a length of " << nbl); if (nb_stored() != 0) { elt_rsvector_<T> ev(c); const_iterator it = std::lower_bound(this->begin(), this->end(), ev); @@ -1137,10 +1137,10 @@ namespace gmm { { o->clear(); } static void do_clear(this_type &v) { v.clear(); } static value_type access(const origin_type *o, const const_iterator &, - const const_iterator &, size_type i) + const const_iterator &, size_type i) { return (*o)[i]; } static reference access(origin_type *o, const iterator &, const iterator &, - size_type i) + size_type i) { return (*o)[i]; } static void resize(this_type &v, size_type n) { v.resize(n); } }; @@ -1151,7 +1151,7 @@ namespace gmm { /******* Optimized operations for rsvector<T> ****************************/ template <typename T> inline void copy(const rsvector<T> &v1, - rsvector<T> &v2) { + rsvector<T> &v2) { GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch"); v2 = v1; } @@ -1166,39 +1166,39 @@ namespace gmm { } template <typename T> inline void copy(const simple_vector_ref<const rsvector<T> *> &v1, - rsvector<T> &v2) + rsvector<T> &v2) { copy(*(v1.origin), v2); } template <typename T> inline void copy(const simple_vector_ref<rsvector<T> *> &v1, rsvector<T> &v2) { copy(*(v1.origin), v2); } template <typename V, typename T> inline void add(const V &v1, - rsvector<T> &v2) { + rsvector<T> &v2) { if ((const void *)(&v1) != (const void *)(&v2)) { GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch"); - add_rsvector(v1, v2, typename linalg_traits<V>::storage_type()); + add_rsvector(v1, v2, typename linalg_traits<V>::storage_type()); } } - template <typename V, typename T> + template <typename V, typename T> inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_dense) { add(v1, v2, abstract_dense(), abstract_sparse()); } - template <typename V, typename T> + template <typename V, typename T> inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline) { add(v1, v2, abstract_skyline(), abstract_sparse()); } - template <typename V, typename T> + template <typename V, typename T> void add_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) { add_rsvector(v1, v2, typename linalg_traits<V>::index_sorted()); } - template <typename V, typename T> + template <typename V, typename T> void add_rsvector(const V &v1, rsvector<T> &v2, linalg_false) { add(v1, v2, abstract_sparse(), abstract_sparse()); } - template <typename V, typename T> + template <typename V, typename T> void add_rsvector(const V &v1, rsvector<T> &v2, linalg_true) { typename linalg_traits<V>::const_iterator it1 = vect_const_begin(v1), ite1 = vect_const_end(v1); @@ -1227,16 +1227,16 @@ namespace gmm { if ((const void *)(&v1) != (const void *)(&v2)) { GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch"); if (same_origin(v1, v2)) - GMM_WARNING2("a conflict is possible in vector copy\n"); + GMM_WARNING2("a conflict is possible in vector copy\n"); copy_rsvector(v1, v2, typename linalg_traits<V>::storage_type()); } } - template <typename V, typename T> + template <typename V, typename T> void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_dense) { copy_vect(v1, v2, abstract_dense(), abstract_sparse()); } - template <typename V, typename T> + template <typename V, typename T> void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline) { copy_vect(v1, v2, abstract_skyline(), abstract_sparse()); } @@ -1244,7 +1244,7 @@ namespace gmm { void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) { copy_rsvector(v1, v2, typename linalg_traits<V>::index_sorted()); } - + template <typename V, typename T2> void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_true) { typedef typename linalg_traits<V>::value_type T1; @@ -1271,7 +1271,7 @@ namespace gmm { v2.base_resize(nn); std::sort(v2.begin(), v2.end()); } - + template <typename T> inline void clean(rsvector<T> &v, double eps) { typedef typename number_traits<T>::magnitude_type R; typename rsvector<T>::iterator it = v.begin(), ite = v.end(); @@ -1280,7 +1280,7 @@ namespace gmm { typename rsvector<T>::iterator itc = it; size_type erased = 1; for (++it; it != ite; ++it) - { *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; } + { *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; } v.base_resize(v.nb_stored() - erased); } } @@ -1294,7 +1294,7 @@ namespace gmm { clean(*pv, eps); svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv); } - + template <typename T> inline size_type nnz(const rsvector<T>& l) { return l.nb_stored(); } @@ -1316,8 +1316,8 @@ namespace gmm { base_iterator it; size_type shift; - - + + iterator &operator ++() { ++it; ++shift; return *this; } iterator &operator --() @@ -1336,21 +1336,21 @@ namespace gmm { { iterator tmp = *this; return (tmp -= i); } difference_type operator -(const iterator &i) const { return it - i.it; } - + reference operator *() const { return *it; } reference operator [](int ii) { return *(it + ii); } - + bool operator ==(const iterator &i) const { return it == i.it; } bool operator !=(const iterator &i) const { return !(i == *this); } bool operator < (const iterator &i) const { return it < i.it; } - size_type index(void) const { return shift; } + size_type index() const { return shift; } - slvector_iterator(void) {} + slvector_iterator() {} slvector_iterator(const base_iterator &iter, size_type s) : it(iter), shift(s) {} }; @@ -1367,8 +1367,8 @@ namespace gmm { base_iterator it; size_type shift; - - + + iterator &operator ++() { ++it; ++shift; return *this; } iterator &operator --() @@ -1387,21 +1387,21 @@ namespace gmm { { iterator tmp = *this; return (tmp -= i); } difference_type operator -(const iterator &i) const { return it - i.it; } - + value_type operator *() const { return *it; } value_type operator [](int ii) { return *(it + ii); } - + bool operator ==(const iterator &i) const { return it == i.it; } bool operator !=(const iterator &i) const { return !(i == *this); } bool operator < (const iterator &i) const { return it < i.it; } - size_type index(void) const { return shift; } + size_type index() const { return shift; } - slvector_const_iterator(void) {} + slvector_const_iterator() {} slvector_const_iterator(const slvector_iterator<T>& iter) : it(iter.it), shift(iter.shift) {} slvector_const_iterator(const base_iterator &iter, size_type s) @@ -1412,7 +1412,7 @@ namespace gmm { /** skyline vector. */ template <typename T> class slvector { - + public : typedef slvector_iterator<T> iterators; typedef slvector_const_iterator<T> const_iterators; @@ -1427,17 +1427,17 @@ namespace gmm { public : - size_type size(void) const { return size_; } - size_type first(void) const { return shift; } - size_type last(void) const { return shift + data.size(); } + size_type size() const { return size_; } + size_type first() const { return shift; } + size_type last() const { return shift + data.size(); } ref_elt_vector<T, slvector<T> > operator [](size_type c) { return ref_elt_vector<T, slvector<T> >(this, c); } - typename std::vector<T>::iterator data_begin(void) { return data.begin(); } - typename std::vector<T>::iterator data_end(void) { return data.end(); } - typename std::vector<T>::const_iterator data_begin(void) const + typename std::vector<T>::iterator data_begin() { return data.begin(); } + typename std::vector<T>::iterator data_end() { return data.end(); } + typename std::vector<T>::const_iterator data_begin() const { return data.begin(); } - typename std::vector<T>::const_iterator data_end(void) const + typename std::vector<T>::const_iterator data_end() const { return data.end(); } void w(size_type c, const T &e); @@ -1450,7 +1450,7 @@ namespace gmm { inline T operator [](size_type c) const { return r(c); } void resize(size_type); - void clear(void) { data.resize(0); shift = 0; } + void clear() { data.resize(0); shift = 0; } void swap(slvector<T> &v) { std::swap(data, v.data); std::swap(shift, v.shift); @@ -1458,7 +1458,7 @@ namespace gmm { } - slvector(void) : data(0), shift(0), size_(0) {} + slvector() : data(0), shift(0), size_(0) {} explicit slvector(size_type l) : data(0), shift(0), size_(l) {} slvector(size_type l, size_type d, size_type s) : data(d), shift(s), size_(l) {} @@ -1477,7 +1477,7 @@ namespace gmm { size_type s = data.size(); if (!s) { data.resize(1); shift = c; } else if (c < shift) { - data.resize(s + shift - c); + data.resize(s + shift - c); typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1; typename std::vector<T>::iterator it3 = it2 - shift + c; for (; it3 >= it; --it3, --it2) *it2 = *it3; @@ -1496,7 +1496,7 @@ namespace gmm { size_type s = data.size(); if (!s) { data.resize(1, e); shift = c; return; } else if (c < shift) { - data.resize(s + shift - c); + data.resize(s + shift - c); typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1; typename std::vector<T>::iterator it3 = it2 - shift + c; for (; it3 >= it; --it3, --it2) *it2 = *it3; @@ -1513,8 +1513,8 @@ namespace gmm { } data[c - shift] += e; } - - + + template <typename T> struct linalg_traits<slvector<T> > { typedef slvector<T> this_type; typedef this_type origin_type; @@ -1541,10 +1541,10 @@ namespace gmm { { o->clear(); } static void do_clear(this_type &v) { v.clear(); } static value_type access(const origin_type *o, const const_iterator &, - const const_iterator &, size_type i) + const const_iterator &, size_type i) { return (*o)[i]; } static reference access(origin_type *o, const iterator &, const iterator &, - size_type i) + size_type i) { return (*o)[i]; } static void resize(this_type &v, size_type n) { v.resize(n); } };