[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Wed, 20 Feb 2019 03:33:51 -0500 (EST) |
branch: master
commit 05f69a3e91911c0739c03b27ab9360b731aba604
Author: Konstantinos Poulios <address@hidden>
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); }
};