branch: master commit f374fd75eabd2596d3a2fd766ec08932c3307e40 Author: Konstantinos Poulios <logar...@gmail.com> Date: Tue Feb 19 23:15:04 2019 +0100
Whitespace, typos, file encodings --- interface/src/gf_asm.cc | 80 +++---- src/getfem/getfem_assembling.h | 8 +- src/getfem/getfem_generic_assembly_semantic.h | 34 +-- src/getfem/getfem_generic_assembly_tree.h | 16 +- src/getfem/getfem_models.h | 16 +- src/getfem_generic_assembly_tree.cc | 2 +- src/getfem_generic_assembly_workspace.cc | 14 +- src/gmm/gmm_matrix.h | 330 +++++++++++++------------- 8 files changed, 248 insertions(+), 252 deletions(-) diff --git a/interface/src/gf_asm.cc b/interface/src/gf_asm.cc index 42416af..30c2010 100644 --- a/interface/src/gf_asm.cc +++ b/interface/src/gf_asm.cc @@ -70,7 +70,7 @@ public: } const bgeot::multi_index &sizes(getfem::size_type) const { return sizes_; } virtual void compute(getfem::fem_interpolation_context& ctx, - bgeot::base_tensor &t) { + bgeot::base_tensor &t) { bgeot::size_type cv = ctx.convex_num(); coeff.resize(mf.nb_basic_dof_of_element(cv)); gmm::copy @@ -95,7 +95,7 @@ void asm_lsneuman_matrix nterm(ls.get_mesh_fem(), ls.values()); getfem::generic_assembly assem("t=comp(Base(#2).Grad(#1).NonLin(#3));" - "M(#2, #1)+= t(:,:,i,i)"); + "M(#2, #1)+= t(:,:,i,i)"); assem.push_mi(mim); assem.push_mf(mf); assem.push_mf(mf_mult); @@ -107,7 +107,7 @@ void asm_lsneuman_matrix /** - generic normal grad level set matrix (on the whole mesh level set or + generic normal grad level set matrix (on the whole mesh level set or on the specified element set level set or boundary level set) */ @@ -123,7 +123,7 @@ template<typename MAT> void asm_nlsgrad_matrix getfem::generic_assembly assem("t=comp(Grad(#1).NonLin(#3).Grad(#2).NonLin(#3));" - "M(#1, #2)+= sym(t(:,i,i,:,j,j))"); + "M(#1, #2)+= sym(t(:,i,i,:,j,j))"); assem.push_mi(mim); assem.push_mf(mf1); assem.push_mf(mf2); @@ -228,11 +228,11 @@ void asm_stabilization_patch_matrix int options[5] = {0,0,0,0,0}; // float ubvec[1] = {1.03f}; //METIS_mCPartGraphKway(&ne, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), &(adjwgt[0]), &wgtflag, - // &numflag, &nparts, &(ubvec[0]), options, &edgecut, &(part[0])); + // &numflag, &nparts, &(ubvec[0]), options, &edgecut, &(part[0])); //METIS_mCPartGraphRecursive(&ne, &ncon, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), &(adjwgt[0]), &wgtflag, - // &numflag, &nparts, options, &edgecut, &(part[0])); + // &numflag, &nparts, options, &edgecut, &(part[0])); //METIS_PartGraphKway(&ne, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), &(adjwgt[0]), &wgtflag, - // &numflag, &nparts, options, &edgecut, &(part[0])); + // &numflag, &nparts, options, &edgecut, &(part[0])); METIS_PartGraphRecursive(&ne, &(xadj[0]), &(adjncy[0]), &(vwgt[0]), &(adjwgt[0]), &wgtflag, &numflag, &nparts, options, &edgecut, &(part[0])); # else @@ -280,7 +280,7 @@ void asm_stabilization_patch_matrix bgeot::size_type cv = mf_P0.first_convex_of_basic_dof(r); int p=part[indelt[cv]]; gmm::copy(gmm::scaled(gmm::mat_row(MAT_aux, p), 1./size_patch[p]), - gmm::mat_row(MAT_proj, r)); + gmm::mat_row(MAT_proj, r)); } gmm::mult(M0, MAT_proj, M1); @@ -446,16 +446,16 @@ static void do_high_level_generic_assembly(mexargs_in& in, mexargs_out& out) { while (in.remaining()) { std::string varname = in.pop().to_string(); if (varname.compare("select_output") == 0 || - varname.compare("select output") == 0) { + varname.compare("select output") == 0) { GMM_ASSERT1(order > 0, "select_output option is for order 1 or 2" "assemblies only"); with_select_output = true; select_var1 = in.pop().to_string(); if (order == 2) select_var2 = in.pop().to_string(); } else if (varname.compare("Secondary_domain") == 0 || - varname.compare("Secondary_Domain") == 0) { + varname.compare("Secondary_Domain") == 0) { GMM_ASSERT1(!with_secondary, - "Only one secondary domain can be specified"); + "Only one secondary domain can be specified"); secondary_domain = in.pop().to_string(); with_secondary = true; } else { @@ -463,37 +463,37 @@ static void do_high_level_generic_assembly(mexargs_in& in, mexargs_out& out) { const getfem::mesh_fem *mf(0); const getfem::im_data *mimd(0); if (is_meshfem_object(in.front())) - mf = to_meshfem_object(in.pop()); + mf = to_meshfem_object(in.pop()); else if (is_meshimdata_object(in.front())) - mimd = to_meshimdata_object(in.pop()); + mimd = to_meshimdata_object(in.pop()); darray U = in.pop().to_darray(); - GMM_ASSERT1(vectors.find(varname) == vectors.end(), - "The same variable/constant name is repeated twice: " - << varname); + GMM_ASSERT1(vectors.count(varname) == 0, + "The same variable/constant name is repeated twice: " + << varname); GMM_ASSERT1(!with_model || !md.variable_exists(varname), - "The same variable/constant name is already defined in " - "the model: " << varname); + "The same variable/constant name is already defined in " + "the model: " << varname); gmm::resize(vectors[varname], U.size()); gmm::copy(U, vectors[varname]); if (is_cte) { - if (mf) - workspace.add_fem_constant(varname, *mf, vectors[varname]); - else if (mimd) - workspace.add_im_data(varname, *mimd, vectors[varname]); - else - workspace.add_fixed_size_constant(varname, vectors[varname]); + if (mf) + workspace.add_fem_constant(varname, *mf, vectors[varname]); + else if (mimd) + workspace.add_im_data(varname, *mimd, vectors[varname]); + else + workspace.add_fixed_size_constant(varname, vectors[varname]); } else { - if (mf) { - gmm::sub_interval I(nbdof, mf->nb_dof()); - nbdof += mf->nb_dof(); - workspace.add_fem_variable(varname, *mf, I, vectors[varname]); - } else if (mimd) { - THROW_BADARG("Data defined on integration points can not be a variable"); - } else { - gmm::sub_interval I(nbdof, U.size()); - nbdof += U.size(); - workspace.add_fixed_size_variable(varname, I, vectors[varname]); - } + if (mf) { + gmm::sub_interval I(nbdof, mf->nb_dof()); + nbdof += mf->nb_dof(); + workspace.add_fem_variable(varname, *mf, I, vectors[varname]); + } else if (mimd) { + THROW_BADARG("Data defined on integration points can not be a variable"); + } else { + gmm::sub_interval I(nbdof, U.size()); + nbdof += U.size(); + workspace.add_fixed_size_variable(varname, I, vectors[varname]); + } } } } @@ -748,7 +748,7 @@ template <typename T> static inline void dummy_func(T &) {} getfemint::mexargs_out& out) \ { dummy_func(in); dummy_func(out); code } \ }; \ - psub_command psubc = std::make_shared<subc>(); \ + psub_command psubc = std::make_shared<subc>(); \ psubc->arg_in_min = arginmin; psubc->arg_in_max = arginmax; \ psubc->arg_out_min = argoutmin; psubc->arg_out_max = argoutmax; \ subc_tab[cmd_normalize(name)] = psubc; \ @@ -799,9 +799,9 @@ void gf_asm(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { functions are only available for variables, not for constants. `select_output` is an optional parameter which allows to reduce the - output vecotr (for `order` equal to 1) or the matrix (for `order` + output vector (for `order` equal to 1) or the matrix (for `order` equal to 2) to the degrees of freedom of the specified variables. - One variable has to be specified for a vector ouptut and two for a + One variable has to be specified for a vector output and two for a matrix output. Note that if several variables are given, the assembly of the @@ -963,7 +963,7 @@ void gf_asm(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { getfem::asm_nonlinear_elasticity_rhs(B, *mim, *mf_u, U, mf_d, param,*law, rg); } else if (cmd_strmatch(what, "incompressible tangent matrix")) { - const getfem::mesh_fem *mf_p = to_meshfem_object(in.pop()); + const getfem::mesh_fem *mf_p = to_meshfem_object(in.pop()); darray P = in.pop().to_darray(int(mf_p->nb_dof())); gf_real_sparse_by_col K(mf_u->nb_dof(), mf_u->nb_dof()); gf_real_sparse_by_col B(mf_u->nb_dof(), mf_p->nb_dof()); @@ -972,7 +972,7 @@ void gf_asm(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { out.pop().from_sparse(K); out.pop().from_sparse(B); } else if (cmd_strmatch(what, "incompressible rhs")) { - const getfem::mesh_fem *mf_p = to_meshfem_object(in.pop()); + const getfem::mesh_fem *mf_p = to_meshfem_object(in.pop()); darray P = in.pop().to_darray(int(mf_p->nb_dof())); darray RU = out.pop().create_darray_v(unsigned(mf_u->nb_dof())); darray RB = out.pop().create_darray_v(unsigned(mf_p->nb_dof())); diff --git a/src/getfem/getfem_assembling.h b/src/getfem/getfem_assembling.h index d958318..885cbdf 100644 --- a/src/getfem/getfem_assembling.h +++ b/src/getfem/getfem_assembling.h @@ -944,7 +944,7 @@ namespace getfem { } /** - Stiffness matrix for linear elasticity, with Lam� coefficients + Stiffness matrix for linear elasticity, with Lamé coefficients @ingroup asm */ template<class MAT, class VECT> @@ -986,7 +986,7 @@ namespace getfem { /** - Stiffness matrix for linear elasticity, with constant Lam� coefficients + Stiffness matrix for linear elasticity, with constant Lamé coefficients @ingroup asm */ template<class MAT, class VECT> @@ -1589,7 +1589,7 @@ namespace getfem { -> the constraint is simplified: we replace \int{(H_j.psi_j)*phi_i}=\int{R_j.psi_j} (sum over j) with H_j*phi_i = R_j - --> Le principe peut �tre faux : non identique � la projection + --> Le principe peut être faux : non identique à la projection L^2 et peut entrer en conccurence avec les autres ddl -> a revoir */ if (tdof_u == tdof_rh && @@ -1684,7 +1684,7 @@ namespace getfem { gmm::mult(H, e, aux); for (size_type j = 0; j < nb_bimg; ++j) { T c = gmm::vect_sp(aux, base_img[j]); - // if (gmm::abs(c > 1.0E-6) { // � scaler sur l'ensemble de H ... + // if (gmm::abs(c > 1.0E-6) { // à scaler sur l'ensemble de H ... if (c != T(0)) { gmm::add(gmm::scaled(base_img[j], -c), aux); gmm::add(gmm::scaled(base_img_inv[j], -c), f); diff --git a/src/getfem/getfem_generic_assembly_semantic.h b/src/getfem/getfem_generic_assembly_semantic.h index 98ae7b4..f129065 100644 --- a/src/getfem/getfem_generic_assembly_semantic.h +++ b/src/getfem/getfem_generic_assembly_semantic.h @@ -50,33 +50,33 @@ namespace getfem { option = 0 : strict analysis, 1 : do not complain about incompatible test functions but store them, - 2 : cut incompatible test function branches with respect to the + 2 : cut incompatible test function branches with respect to the one in workspace.selected_pair 3 : do not complain about incompatible test functions neither store them. */ void ga_semantic_analysis(ga_tree &tree, - const ga_workspace &workspace, - const mesh &m, - size_type ref_elt_dim, - bool eval_fixed_size, - bool ignore_X, int option = 0); + const ga_workspace &workspace, + const mesh &m, + size_type ref_elt_dim, + bool eval_fixed_size, + bool ignore_X, int option = 0); /* Extract the variables used in a sub-tree, ignoring or not the data. Variable groups are taken into account. Return if at least one variable as been detected or not */ bool ga_extract_variables(const pga_tree_node pnode, - const ga_workspace &workspace, - const mesh &m, - std::set<var_trans_pair> &vars, - bool ignore_data); + const ga_workspace &workspace, + const mesh &m, + std::set<var_trans_pair> &vars, + bool ignore_data); /* Extract a sub tree which consists of the corresponding node and of the terms multiplying this term, but not the term in addition. The aim is to expand an expression is a sum of elementary factors. Complains if a nonlinear term is encountered. */ void ga_extract_factor(ga_tree &result_tree, pga_tree_node pnode, - pga_tree_node &new_pnode); + pga_tree_node &new_pnode); /* Extract the constant term of degree 1 expressions. */ bool ga_node_extract_constant_term @@ -93,16 +93,16 @@ namespace getfem { The tree is modified and should be copied first and passed to ga_semantic_analysis after for enrichment. */ void ga_derivative(ga_tree &tree, const ga_workspace &workspace, - const mesh &m, const std::string &varname, - const std::string &interpolatename, - size_type order); + const mesh &m, const std::string &varname, + const std::string &interpolatename, + size_type order); std::string ga_derivative_scalar_function(const std::string &expr, - const std::string &var); + const std::string &var); bool ga_is_affine(const ga_tree &tree, const ga_workspace &workspace, - const std::string &varname, - const std::string &interpolatename); + const std::string &varname, + const std::string &interpolatename); // Function of internal use inline size_type ref_elt_dim_of_mesh(const mesh &m) { diff --git a/src/getfem/getfem_generic_assembly_tree.h b/src/getfem/getfem_generic_assembly_tree.h index 84b0dcb..f7628d6 100644 --- a/src/getfem/getfem_generic_assembly_tree.h +++ b/src/getfem/getfem_generic_assembly_tree.h @@ -201,11 +201,11 @@ namespace getfem { typedef std::shared_ptr<std::string> pstring; // Print error message indicating the position in the assembly string void ga_throw_error_msg(pstring expr, size_type pos, - const std::string &msg); + const std::string &msg); # define ga_throw_error(expr, pos, msg) \ - { std::stringstream ss; ss << msg; \ - ga_throw_error_msg(expr, pos, ss.str()); \ + { std::stringstream ss; ss << msg; \ + ga_throw_error_msg(expr, pos, ss.str()); \ GMM_ASSERT1(false, "Error in assembly string" ); \ } @@ -391,7 +391,7 @@ namespace getfem { : node_type(GA_NODE_VOID), test_function_type(-1), qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0), pos(0), expr(0), der1(0), der2(0), symmetric_op(false), hash_value(0) {} - ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_) + ga_tree_node(GA_NODE_TYPE ty, size_type p, pstring expr_) : node_type(ty), test_function_type(-1), qdim1(0), qdim2(0), nbc1(0), nbc2(0), nbc3(0), pos(p), expr(expr_), der1(0), der2(0), symmetric_op(false), @@ -421,7 +421,7 @@ namespace getfem { void add_scalar(scalar_type val, size_type pos, pstring expr); void add_allindices(size_type pos, pstring expr); void add_name(const char *name, size_type length, size_type pos, - pstring expr); + pstring expr); void add_sub_tree(ga_tree &sub_tree); void add_params(size_type pos, pstring expr); void add_matrix(size_type pos, pstring expr); @@ -446,9 +446,9 @@ namespace getfem { std::swap(secondary_domain, tree.secondary_domain); } - ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {} + ga_tree() : root(nullptr), current_node(nullptr), secondary_domain() {} - ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr), + ga_tree(const ga_tree &tree) : root(nullptr), current_node(nullptr), secondary_domain(tree.secondary_domain) { if (tree.root) copy_node(tree.root, nullptr, root); } @@ -473,7 +473,7 @@ namespace getfem { // Transform the expression of a node and its sub-nodes in the equivalent // assembly string sent to ostream str void ga_print_node(const pga_tree_node pnode, - std::ostream &str); + std::ostream &str); // The same for the whole tree, the result is a std::string std::string ga_tree_to_string(const ga_tree &tree); diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h index 3dd3b4c..b8b697c 100644 --- a/src/getfem/getfem_models.h +++ b/src/getfem/getfem_models.h @@ -185,7 +185,7 @@ namespace getfem { model_complex_plain_vector affine_complex_value; scalar_type alpha; // Factor for the affine dependent variables std::string org_name; // Name of the original variable for affine - // dependent variables + // dependent variables // im data description const im_data *pim_data; @@ -239,7 +239,7 @@ namespace getfem { { return is_complex ? complex_value[0].size() : real_value[0].size(); } void set_size(); - }; + }; // struct var_description public: @@ -506,21 +506,21 @@ namespace getfem { active_bricks.add(ib); } - /** Disable a variable (and its attached mutlipliers). */ + /** Disable a variable (and its attached mutlipliers). */ void disable_variable(const std::string &name); - /** Enable a variable (and its attached mutlipliers). */ + /** Enable a variable (and its attached mutlipliers). */ void enable_variable(const std::string &name); - /** Says if a name corresponds to a declared variable. */ + /** States if a name corresponds to a declared variable. */ bool variable_exists(const std::string &name) const; bool is_disabled_variable(const std::string &name) const; - /** Says if a name corresponds to a declared data or disabled variable. */ + /** States if a name corresponds to a declared data or disabled variable. */ bool is_data(const std::string &name) const; - /** Says if a name corresponds to a declared data. */ + /** States if a name corresponds to a declared data. */ bool is_true_data(const std::string &name) const; bool is_affine_dependent_variable(const std::string &name) const; @@ -882,7 +882,7 @@ namespace getfem { return cTM; } - /** Gives the access to the right hand side of the tangent linear system. + /** Gives access to the right hand side of the tangent linear system. For the real version. An assembly of the rhs has to be done first. */ const model_real_plain_vector &real_rhs() const { GMM_ASSERT1(!complex_version, "This model is a complex one"); diff --git a/src/getfem_generic_assembly_tree.cc b/src/getfem_generic_assembly_tree.cc index c553596..2d673c7 100644 --- a/src/getfem_generic_assembly_tree.cc +++ b/src/getfem_generic_assembly_tree.cc @@ -1559,7 +1559,7 @@ namespace getfem { case GA_MINUS: // unary - tree.add_op(GA_UNARY_MINUS, token_pos, expr); state = 1; break; - + case GA_PLUS: // unary + state = 1; break; diff --git a/src/getfem_generic_assembly_workspace.cc b/src/getfem_generic_assembly_workspace.cc index 15cb98e..506ad16 100644 --- a/src/getfem_generic_assembly_workspace.cc +++ b/src/getfem_generic_assembly_workspace.cc @@ -59,7 +59,7 @@ namespace getfem { (const std::string &name, const mesh_fem &mf, const model_real_plain_vector &VV) { GMM_ASSERT1(mf.nb_dof(), "The provided mesh_fem of variable" << name - << "has zero degrees of freedom."); + << "has zero degrees of freedom"); size_type Q = gmm::vect_size(VV)/mf.nb_dof(); if (Q == 0) Q = size_type(1); variables[name] = var_description(false, true, &mf, @@ -374,9 +374,6 @@ namespace getfem { } - - - const mesh_region & ga_workspace::register_region(const mesh &m, const mesh_region ®ion) { if (&m == &dummy_mesh()) @@ -389,7 +386,6 @@ namespace getfem { return lmr.back(); } - void ga_workspace::add_tree(ga_tree &tree, const mesh &m, const mesh_im &mim, const mesh_region &rg, const std::string &expr, @@ -454,7 +450,7 @@ namespace getfem { trees.push_back(tree_description()); trees.back().mim = &mim; trees.back().m = &m; trees.back().rg = &rg; - trees.back().secondary_domain = tree.secondary_domain; + trees.back().secondary_domain = tree.secondary_domain; trees.back().ptree = new ga_tree; trees.back().ptree->swap(tree); pga_tree_node root = trees.back().ptree->root; @@ -507,7 +503,7 @@ namespace getfem { const mesh_im &mim, const mesh_region &rg_, size_type add_derivative_order, - const std::string &secondary_dom) { + const std::string &secondary_dom) { const mesh_region &rg = register_region(mim.linked_mesh(), rg_); // cout << "adding expression " << expr << endl; GA_TIC; @@ -516,7 +512,7 @@ namespace getfem { ga_read_string(expr, ltrees[0], macro_dictionary()); if (secondary_dom.size()) { GMM_ASSERT1(secondary_domain_exists(secondary_dom), - "Unknow secondary domain " << secondary_dom); + "Unknown secondary domain " << secondary_dom); ltrees[0].secondary_domain = secondary_dom; } // cout << "read : " << ga_tree_to_string(ltrees[0]) << endl; @@ -910,7 +906,7 @@ namespace getfem { { if (ptree) delete ptree; ptree = 0; } ga_workspace::ga_workspace(const getfem::model &md_, - bool enable_all_variables) + bool enable_all_variables) : md(&md_), parent_workspace(0), enable_all_md_variables(enable_all_variables), macro_dict(md_.macro_dictionary()) diff --git a/src/gmm/gmm_matrix.h b/src/gmm/gmm_matrix.h index 23fb9d2..3b133bb 100644 --- a/src/gmm/gmm_matrix.h +++ b/src/gmm/gmm_matrix.h @@ -48,9 +48,9 @@ namespace gmm { /* ******************************************************************** */ - /* */ - /* Identity matrix */ - /* */ + /* */ + /* Identity matrix */ + /* */ /* ******************************************************************** */ struct identity_matrix { @@ -71,7 +71,7 @@ namespace gmm void mult(const identity_matrix&, const V1 &v1, V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline - void mult(const identity_matrix&, const V1 &v1, const V2 &v2) + void mult(const identity_matrix&, const V1 &v1, const V2 &v2) { copy(v1, v2); } template <typename V1, typename V2, typename V3> inline void mult(const identity_matrix&, const V1 &v1, const V2 &v2, V3 &v3) @@ -83,25 +83,25 @@ namespace gmm void left_mult(const identity_matrix&, const V1 &v1, V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline - void left_mult(const identity_matrix&, const V1 &v1, const V2 &v2) + void left_mult(const identity_matrix&, const V1 &v1, const V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline void right_mult(const identity_matrix&, const V1 &v1, V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline - void right_mult(const identity_matrix&, const V1 &v1, const V2 &v2) + void right_mult(const identity_matrix&, const V1 &v1, const V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline void transposed_left_mult(const identity_matrix&, const V1 &v1, V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline - void transposed_left_mult(const identity_matrix&, const V1 &v1,const V2 &v2) + void transposed_left_mult(const identity_matrix&, const V1 &v1,const V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline void transposed_right_mult(const identity_matrix&, const V1 &v1, V2 &v2) { copy(v1, v2); } template <typename V1, typename V2> inline - void transposed_right_mult(const identity_matrix&,const V1 &v1,const V2 &v2) + void transposed_right_mult(const identity_matrix&,const V1 &v1,const V2 &v2) { copy(v1, v2); } template <typename M> void copy_ident(const identity_matrix&, M &m) { size_type i = 0, n = std::min(mat_nrows(m), mat_ncols(m)); @@ -109,7 +109,7 @@ namespace gmm for (; i < n; ++i) m(i,i) = typename linalg_traits<M>::value_type(1); } template <typename M> inline void copy(const identity_matrix&, M &m) - { copy_ident(identity_matrix(), m); } + { copy_ident(identity_matrix(), m); } template <typename M> inline void copy(const identity_matrix &, const M &m) { copy_ident(identity_matrix(), linalg_const_cast(m)); } template <typename V1, typename V2> inline @@ -124,24 +124,24 @@ namespace gmm inline bool is_identity(const identity_matrix&) { return true; } /* ******************************************************************** */ - /* */ - /* Row matrix */ - /* */ + /* */ + /* Row matrix */ + /* */ /* ******************************************************************** */ template<typename V> class row_matrix { protected : std::vector<V> li; /* array of rows. */ size_type nc; - + public : - + typedef typename linalg_traits<V>::reference reference; typedef typename linalg_traits<V>::value_type value_type; - + row_matrix(size_type r, size_type c) : li(r, V(c)), nc(c) {} row_matrix(void) : nc(0) {} - reference operator ()(size_type l, size_type c) + reference operator ()(size_type l, size_type c) { return li[l][c]; } value_type operator ()(size_type l, size_type c) const { return li[l][c]; } @@ -151,19 +151,19 @@ namespace gmm typename std::vector<V>::iterator begin(void) { return li.begin(); } - typename std::vector<V>::iterator end(void) + typename std::vector<V>::iterator end(void) { return li.end(); } typename std::vector<V>::const_iterator begin(void) const { return li.begin(); } typename std::vector<V>::const_iterator end(void) const { return li.end(); } - - + + V& row(size_type i) { return li[i]; } const V& row(size_type i) const { return li[i]; } V& operator[](size_type i) { return li[i]; } const V& operator[](size_type i) const { return li[i]; } - + inline size_type nrows(void) const { return li.size(); } inline size_type ncols(void) const { return nc; } @@ -176,7 +176,7 @@ namespace gmm li.resize(m); for (size_type i=nr; i < m; ++i) gmm::resize(li[i], n); if (n != nc) { - for (size_type i=0; i < nr; ++i) gmm::resize(li[i], n); + for (size_type i=0; i < nr; ++i) gmm::resize(li[i], n); nc = n; } } @@ -213,7 +213,7 @@ namespace gmm { return m.end(); } static const_sub_row_type row(const const_row_iterator &it) { return const_sub_row_type(*it); } - static sub_row_type row(const row_iterator &it) + static sub_row_type row(const row_iterator &it) { return sub_row_type(*it); } static origin_type* origin(this_type &m) { return &m; } static const origin_type* origin(const this_type &m) { return &m; } @@ -232,21 +232,21 @@ namespace gmm (std::ostream &o, const row_matrix<V>& m) { gmm::write(o,m); return o; } /* ******************************************************************** */ - /* */ - /* Column matrix */ - /* */ + /* */ + /* Column matrix */ + /* */ /* ******************************************************************** */ template<typename V> class col_matrix { protected : std::vector<V> li; /* array of columns. */ size_type nr; - + public : - + typedef typename linalg_traits<V>::reference reference; typedef typename linalg_traits<V>::value_type value_type; - + col_matrix(size_type r, size_type c) : li(c, V(r)), nr(r) { } col_matrix(void) : nr(0) {} reference operator ()(size_type l, size_type c) @@ -264,13 +264,13 @@ namespace gmm typename std::vector<V>::iterator begin(void) { return li.begin(); } - typename std::vector<V>::iterator end(void) + typename std::vector<V>::iterator end(void) { return li.end(); } typename std::vector<V>::const_iterator begin(void) const { return li.begin(); } typename std::vector<V>::const_iterator end(void) const { return li.end(); } - + inline size_type ncols(void) const { return li.size(); } inline size_type nrows(void) const { return nr; } @@ -283,7 +283,7 @@ namespace gmm li.resize(n); for (size_type i=nc; i < n; ++i) gmm::resize(li[i], m); if (m != nr) { - for (size_type i=0; i < nc; ++i) gmm::resize(li[i], m); + for (size_type i=0; i < nc; ++i) gmm::resize(li[i], m); nr = m; } } @@ -319,7 +319,7 @@ namespace gmm { return m.end(); } static const_sub_col_type col(const const_col_iterator &it) { return *it; } - static sub_col_type col(const col_iterator &it) + static sub_col_type col(const col_iterator &it) { return *it; } static origin_type* origin(this_type &m) { return &m; } static const origin_type* origin(const this_type &m) { return &m; } @@ -338,9 +338,9 @@ namespace gmm (std::ostream &o, const col_matrix<V>& m) { gmm::write(o,m); return o; } /* ******************************************************************** */ - /* */ - /* Dense matrix */ - /* */ + /* */ + /* Dense matrix */ + /* */ /* ******************************************************************** */ template<typename T> class dense_matrix : public std::vector<T> { @@ -350,12 +350,12 @@ namespace gmm typedef typename std::vector<T>::const_iterator const_iterator; typedef typename std::vector<T>::reference reference; typedef typename std::vector<T>::const_reference const_reference; - + protected: size_type nbc, nbl; - + public: - + inline const_reference operator ()(size_type l, size_type c) const { GMM_ASSERT2(l < nbl && c < nbc, "out of range"); return *(this->begin() + c*nbl+l); @@ -371,13 +371,13 @@ namespace gmm void resize(size_type, size_type); void base_resize(size_type, size_type); void reshape(size_type, size_type); - + void fill(T a, T b = T(0)); inline size_type nrows(void) const { return nbl; } inline size_type ncols(void) const { return nbc; } void swap(dense_matrix<T> &m) { std::vector<T>::swap(m); std::swap(nbc, m.nbc); std::swap(nbl, m.nbl); } - + dense_matrix(size_type l, size_type c) : std::vector<T>(c*l), nbc(c), nbl(l) {} dense_matrix(void) { nbl = nbc = 0; } @@ -389,33 +389,33 @@ namespace gmm } template<typename T> void dense_matrix<T>::base_resize(size_type m, - size_type n) + size_type n) { std::vector<T>::resize(n*m); nbl = m; nbc = n; } - + template<typename T> void dense_matrix<T>::resize(size_type m, size_type n) { if (n*m > nbc*nbl) std::vector<T>::resize(n*m); if (m < nbl) { for (size_type i = 1; i < std::min(nbc, n); ++i) - std::copy(this->begin()+i*nbl, this->begin()+(i*nbl+m), - this->begin()+i*m); + std::copy(this->begin()+i*nbl, this->begin()+(i*nbl+m), + this->begin()+i*m); for (size_type i = std::min(nbc, n); i < n; ++i) - std::fill(this->begin()+(i*m), this->begin()+(i+1)*m, T(0)); + std::fill(this->begin()+(i*m), this->begin()+(i+1)*m, T(0)); } else if (m > nbl) { /* do nothing when the nb of rows does not change */ for (size_type i = std::min(nbc, n); i > 1; --i) - std::copy(this->begin()+(i-1)*nbl, this->begin()+i*nbl, - this->begin()+(i-1)*m); + std::copy(this->begin()+(i-1)*nbl, this->begin()+i*nbl, + this->begin()+(i-1)*m); for (size_type i = 0; i < std::min(nbc, n); ++i) - std::fill(this->begin()+(i*m+nbl), this->begin()+(i+1)*m, T(0)); + std::fill(this->begin()+(i*m+nbl), this->begin()+(i+1)*m, T(0)); } if (n*m < nbc*nbl) std::vector<T>::resize(n*m); nbl = m; nbc = n; } - + template<typename T> void dense_matrix<T>::fill(T a, T b) { std::fill(this->begin(), this->end(), b); size_type n = std::min(nbl, nbc); - if (a != b) for (size_type i = 0; i < n; ++i) (*this)(i,i) = a; + if (a != b) for (size_type i = 0; i < n; ++i) (*this)(i,i) = a; } template <typename T> struct linalg_traits<dense_matrix<T> > { @@ -427,25 +427,25 @@ namespace gmm typedef T& reference; typedef abstract_dense storage_type; typedef tab_ref_reg_spaced_with_origin<typename this_type::iterator, - this_type> sub_row_type; + this_type> sub_row_type; typedef tab_ref_reg_spaced_with_origin<typename this_type::const_iterator, - this_type> const_sub_row_type; + this_type> const_sub_row_type; typedef dense_compressed_iterator<typename this_type::iterator, - typename this_type::iterator, - this_type *> row_iterator; + typename this_type::iterator, + this_type *> row_iterator; typedef dense_compressed_iterator<typename this_type::const_iterator, - typename this_type::iterator, - const this_type *> const_row_iterator; - typedef tab_ref_with_origin<typename this_type::iterator, - this_type> sub_col_type; + typename this_type::iterator, + const this_type *> const_row_iterator; + typedef tab_ref_with_origin<typename this_type::iterator, + this_type> sub_col_type; typedef tab_ref_with_origin<typename this_type::const_iterator, - this_type> const_sub_col_type; + this_type> const_sub_col_type; typedef dense_compressed_iterator<typename this_type::iterator, - typename this_type::iterator, - this_type *> col_iterator; + typename this_type::iterator, + this_type *> col_iterator; typedef dense_compressed_iterator<typename this_type::const_iterator, - typename this_type::iterator, - const this_type *> const_col_iterator; + typename this_type::iterator, + const this_type *> const_col_iterator; typedef col_and_row sub_orientation; typedef linalg_true index_sorted; static size_type nrows(const this_type &m) { return m.nrows(); } @@ -493,7 +493,7 @@ namespace gmm /* ******************************************************************** */ /* */ - /* Read only compressed sparse column matrix */ + /* Read only compressed sparse column matrix */ /* */ /* ******************************************************************** */ @@ -518,7 +518,7 @@ namespace gmm template <typename PT1, typename PT2, typename PT3, int cshift> void init_with(const csc_matrix_ref<PT1,PT2,PT3,cshift>& B) { init_with_good_format(B); } - template <typename U, int cshift> + template <typename U, int cshift> void init_with(const csc_matrix<U, cshift>& B) { init_with_good_format(B); } @@ -529,9 +529,9 @@ namespace gmm size_type nrows(void) const { return nr; } size_type ncols(void) const { return nc; } - void swap(csc_matrix<T, shift> &m) { - std::swap(pr, m.pr); - std::swap(ir, m.ir); std::swap(jc, m.jc); + void swap(csc_matrix<T, shift> &m) { + std::swap(pr, m.pr); + std::swap(ir, m.ir); std::swap(jc, m.jc); std::swap(nc, m.nc); std::swap(nr, m.nr); } value_type operator()(size_type i, size_type j) const @@ -552,30 +552,30 @@ namespace gmm for (size_type j = 0; j < nc; ++j) { col_type col = mat_const_col(B, j); typename linalg_traits<typename org_type<col_type>::t>::const_iterator - it = vect_const_begin(col), ite = vect_const_end(col); + it = vect_const_begin(col), ite = vect_const_end(col); for (size_type k = 0; it != ite; ++it, ++k) { - pr[jc[j]-shift+k] = *it; - ir[jc[j]-shift+k] = IND_TYPE(it.index() + shift); + pr[jc[j]-shift+k] = *it; + ir[jc[j]-shift+k] = IND_TYPE(it.index() + shift); } } } - + template <typename T, int shift> template <typename Matrix> void csc_matrix<T, shift>::init_with(const Matrix &A) { col_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A)); copy(A, B); init_with_good_format(B); } - + template <typename T, int shift> void csc_matrix<T, shift>::init_with_identity(size_type n) { - nc = nr = n; + nc = nr = n; pr.resize(nc); ir.resize(nc); jc.resize(nc+1); for (size_type j = 0; j < nc; ++j) { ir[j] = jc[j] = shift + j; pr[j] = T(1); } jc[nc] = shift + nc; } - + template <typename T, int shift> csc_matrix<T, shift>::csc_matrix(size_type nnr, size_type nnc) : nc(nnc), nr(nnr) { @@ -601,7 +601,7 @@ namespace gmm typedef cs_vector_ref<const T *, const IND_TYPE *, shift> const_sub_col_type; typedef sparse_compressed_iterator<const T *, const IND_TYPE *, - const IND_TYPE *, shift> + const IND_TYPE *, shift> const_col_iterator; typedef abstract_null_type col_iterator; typedef col_major sub_orientation; @@ -612,12 +612,12 @@ namespace gmm { return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0], m.nr, &m.pr[0]); } static const_col_iterator col_end(const this_type &m) { return const_col_iterator(&m.pr[0],&m.ir[0],&m.jc[0]+m.nc, - m.nr,&m.pr[0]); + m.nr,&m.pr[0]); } static const_sub_col_type col(const const_col_iterator &it) { return const_sub_col_type(it.pr + *(it.jc) - shift, - it.ir + *(it.jc) - shift, - *(it.jc + 1) - *(it.jc), it.n); + it.ir + *(it.jc) - shift, + *(it.jc + 1) - *(it.jc), it.n); } static const origin_type* origin(const this_type &m) { return &m.pr[0]; } static void do_clear(this_type &m) { m.do_clear(); } @@ -629,7 +629,7 @@ namespace gmm std::ostream &operator << (std::ostream &o, const csc_matrix<T, shift>& m) { gmm::write(o,m); return o; } - + template <typename T, int shift> inline void copy(const identity_matrix &, csc_matrix<T, shift>& M) { M.init_with_identity(mat_nrows(M)); } @@ -640,7 +640,7 @@ namespace gmm /* ******************************************************************** */ /* */ - /* Read only compressed sparse row matrix */ + /* Read only compressed sparse row matrix */ /* */ /* ******************************************************************** */ @@ -678,16 +678,16 @@ namespace gmm size_type nrows(void) const { return nr; } size_type ncols(void) const { return nc; } - void swap(csr_matrix<T, shift> &m) { - std::swap(pr, m.pr); - std::swap(ir,m.ir); std::swap(jc, m.jc); + void swap(csr_matrix<T, shift> &m) { + std::swap(pr, m.pr); + std::swap(ir,m.ir); std::swap(jc, m.jc); std::swap(nc, m.nc); std::swap(nr,m.nr); } - + value_type operator()(size_type i, size_type j) const { return mat_row(*this, i)[j]; } }; - + template <typename T, int shift> template <typename Matrix> void csr_matrix<T, shift>::init_with_good_format(const Matrix &B) { typedef typename linalg_traits<Matrix>::const_sub_row_type row_type; @@ -702,24 +702,24 @@ namespace gmm for (size_type j = 0; j < nr; ++j) { row_type row = mat_const_row(B, j); typename linalg_traits<typename org_type<row_type>::t>::const_iterator - it = vect_const_begin(row), ite = vect_const_end(row); + it = vect_const_begin(row), ite = vect_const_end(row); for (size_type k = 0; it != ite; ++it, ++k) { - pr[jc[j]-shift+k] = *it; - ir[jc[j]-shift+k] = IND_TYPE(it.index()+shift); + pr[jc[j]-shift+k] = *it; + ir[jc[j]-shift+k] = IND_TYPE(it.index()+shift); } } } - template <typename T, int shift> template <typename Matrix> - void csr_matrix<T, shift>::init_with(const Matrix &A) { - row_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A)); - copy(A, B); + template <typename T, int shift> template <typename Matrix> + void csr_matrix<T, shift>::init_with(const Matrix &A) { + row_matrix<wsvector<T> > B(mat_nrows(A), mat_ncols(A)); + copy(A, B); init_with_good_format(B); } - template <typename T, int shift> + template <typename T, int shift> void csr_matrix<T, shift>::init_with_identity(size_type n) { - nc = nr = n; + nc = nr = n; pr.resize(nr); ir.resize(nr); jc.resize(nr+1); for (size_type j = 0; j < nr; ++j) { ir[j] = jc[j] = shift + j; pr[j] = T(1); } @@ -753,7 +753,7 @@ namespace gmm typedef cs_vector_ref<const T *, const IND_TYPE *, shift> const_sub_row_type; typedef sparse_compressed_iterator<const T *, const IND_TYPE *, - const IND_TYPE *, shift> + const IND_TYPE *, shift> const_row_iterator; typedef abstract_null_type row_iterator; typedef row_major sub_orientation; @@ -766,8 +766,8 @@ namespace gmm { return const_row_iterator(&m.pr[0], &m.ir[0], &m.jc[0] + m.nr, m.nc, &m.pr[0]); } static const_sub_row_type row(const const_row_iterator &it) { return const_sub_row_type(it.pr + *(it.jc) - shift, - it.ir + *(it.jc) - shift, - *(it.jc + 1) - *(it.jc), it.n); + it.ir + *(it.jc) - shift, + *(it.jc + 1) - *(it.jc), it.n); } static const origin_type* origin(const this_type &m) { return &m.pr[0]; } static void do_clear(this_type &m) { m.do_clear(); } @@ -779,7 +779,7 @@ namespace gmm std::ostream &operator << (std::ostream &o, const csr_matrix<T, shift>& m) { gmm::write(o,m); return o; } - + template <typename T, int shift> inline void copy(const identity_matrix &, csr_matrix<T, shift>& M) { M.init_with_identity(mat_nrows(M)); } @@ -789,9 +789,9 @@ namespace gmm { M.init_with(A); } /* ******************************************************************** */ - /* */ - /* Block matrix */ - /* */ + /* */ + /* Block matrix */ + /* */ /* ******************************************************************** */ template <typename MAT> class block_matrix { @@ -811,7 +811,7 @@ namespace gmm size_type ncolblocks(void) const { return ncolblocks_; } const sub_interval &subrowinterval(size_type i) const { return introw[i]; } const sub_interval &subcolinterval(size_type i) const { return intcol[i]; } - const MAT &block(size_type i, size_type j) const + const MAT &block(size_type i, size_type j) const { return blocks[j*ncolblocks_+i]; } MAT &block(size_type i, size_type j) { return blocks[j*ncolblocks_+i]; } @@ -820,20 +820,20 @@ namespace gmm value_type operator() (size_type i, size_type j) const { size_type k, l; for (k = 0; k < nrowblocks_; ++k) - if (i >= introw[k].min && i < introw[k].max) break; + if (i >= introw[k].min && i < introw[k].max) break; for (l = 0; l < nrowblocks_; ++l) - if (j >= introw[l].min && j < introw[l].max) break; + if (j >= introw[l].min && j < introw[l].max) break; return (block(k, l))(i - introw[k].min, j - introw[l].min); } reference operator() (size_type i, size_type j) { size_type k, l; for (k = 0; k < nrowblocks_; ++k) - if (i >= introw[k].min && i < introw[k].max) break; + if (i >= introw[k].min && i < introw[k].max) break; for (l = 0; l < nrowblocks_; ++l) - if (j >= introw[l].min && j < introw[l].max) break; + if (j >= introw[l].min && j < introw[l].max) break; return (block(k, l))(i - introw[k].min, j - introw[l].min); } - + template <typename CONT> void resize(const CONT &c1, const CONT &c2); template <typename CONT> block_matrix(const CONT &c1, const CONT &c2) { resize(c1, c2); } @@ -864,17 +864,17 @@ namespace gmm static origin_type* origin(this_type &m) { return &m; } static const origin_type* origin(const this_type &m) { return &m; } static void do_clear(this_type &m) { m.do_clear(); } - // access to be done ... + // access to be done ... static void resize(this_type &, size_type , size_type) { GMM_ASSERT1(false, "Sorry, to be done"); } static void reshape(this_type &, size_type , size_type) { GMM_ASSERT1(false, "Sorry, to be done"); } }; - template <typename MAT> void block_matrix<MAT>::do_clear(void) { + template <typename MAT> void block_matrix<MAT>::do_clear(void) { for (size_type j = 0, l = 0; j < ncolblocks_; ++j) for (size_type i = 0, k = 0; i < nrowblocks_; ++i) - clear(block(i,j)); + clear(block(i,j)); } template <typename MAT> template <typename CONT> @@ -886,8 +886,8 @@ namespace gmm for (size_type j = 0, l = 0; j < ncolblocks_; ++j) { intcol[j] = sub_interval(l, c2[j]); l += c2[j]; for (size_type i = 0, k = 0; i < nrowblocks_; ++i) { - if (j == 0) { introw[i] = sub_interval(k, c1[i]); k += c1[i]; } - block(i, j) = MAT(c1[i], c2[j]); + if (j == 0) { introw[i] = sub_interval(k, c1[i]); k += c1[i]; } + block(i, j) = MAT(c1[i], c2[j]); } } } @@ -896,14 +896,14 @@ namespace gmm void copy(const block_matrix<M1> &m1, M2 &m2) { for (size_type j = 0; j < m1.ncolblocks(); ++j) for (size_type i = 0; i < m1.nrowblocks(); ++i) - copy(m1.block(i,j), sub_matrix(m2, m1.subrowinterval(i), - m1.subcolinterval(j))); + copy(m1.block(i,j), sub_matrix(m2, m1.subrowinterval(i), + m1.subcolinterval(j))); } template <typename M1, typename M2> void copy(const block_matrix<M1> &m1, const M2 &m2) { copy(m1, linalg_const_cast(m2)); } - + template <typename MAT, typename V1, typename V2> void mult(const block_matrix<MAT> &m, const V1 &v1, V2 &v2) { @@ -911,9 +911,9 @@ namespace gmm typename sub_vector_type<V2 *, sub_interval>::vector_type sv; for (size_type i = 0; i < m.nrowblocks() ; ++i) for (size_type j = 0; j < m.ncolblocks() ; ++j) { - sv = sub_vector(v2, m.subrowinterval(i)); - mult(m.block(i,j), - sub_vector(v1, m.subcolinterval(j)), sv, sv); + sv = sub_vector(v2, m.subrowinterval(i)); + mult(m.block(i,j), + sub_vector(v1, m.subcolinterval(j)), sv, sv); } } @@ -922,16 +922,16 @@ namespace gmm typename sub_vector_type<V3 *, sub_interval>::vector_type sv; for (size_type i = 0; i < m.nrowblocks() ; ++i) for (size_type j = 0; j < m.ncolblocks() ; ++j) { - sv = sub_vector(v3, m.subrowinterval(i)); - if (j == 0) - mult(m.block(i,j), - sub_vector(v1, m.subcolinterval(j)), - sub_vector(v2, m.subrowinterval(i)), sv); - else - mult(m.block(i,j), - sub_vector(v1, m.subcolinterval(j)), sv, sv); + sv = sub_vector(v3, m.subrowinterval(i)); + if (j == 0) + mult(m.block(i,j), + sub_vector(v1, m.subcolinterval(j)), + sub_vector(v2, m.subrowinterval(i)), sv); + else + mult(m.block(i,j), + sub_vector(v1, m.subcolinterval(j)), sv, sv); } - + } template <typename MAT, typename V1, typename V2> @@ -939,15 +939,15 @@ namespace gmm { mult(m, v1, linalg_const_cast(v2)); } template <typename MAT, typename V1, typename V2, typename V3> - void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2, - const V3 &v3) + void mult(const block_matrix<MAT> &m, const V1 &v1, const V2 &v2, + const V3 &v3) { mult_const(m, v1, v2, linalg_const_cast(v3)); } } /* ******************************************************************** */ - /* */ - /* Distributed matrices */ - /* */ + /* */ + /* Distributed matrices */ + /* */ /* ******************************************************************** */ #ifdef GMM_USES_MPI @@ -955,8 +955,8 @@ namespace gmm namespace gmm { - - + + template <typename T> inline MPI_Datatype mpi_type(T) { GMM_ASSERT1(false, "Sorry unsupported type"); return MPI_FLOAT; } inline MPI_Datatype mpi_type(double) { return MPI_DOUBLE; } @@ -980,7 +980,7 @@ namespace gmm { const MAT &local_matrix(void) const { return M; } MAT &local_matrix(void) { return M; } }; - + template <typename MAT> inline MAT &eff_matrix(MAT &m) { return m; } template <typename MAT> inline const MAT &eff_matrix(const MAT &m) { return m; } @@ -988,29 +988,29 @@ namespace gmm { MAT &eff_matrix(mpi_distributed_matrix<MAT> &m) { return m.M; } template <typename MAT> inline const MAT &eff_matrix(const mpi_distributed_matrix<MAT> &m) { return m.M; } - + template <typename MAT1, typename MAT2> inline void copy(const mpi_distributed_matrix<MAT1> &m1, - mpi_distributed_matrix<MAT2> &m2) + mpi_distributed_matrix<MAT2> &m2) { copy(eff_matrix(m1), eff_matrix(m2)); } template <typename MAT1, typename MAT2> inline void copy(const mpi_distributed_matrix<MAT1> &m1, - const mpi_distributed_matrix<MAT2> &m2) + const mpi_distributed_matrix<MAT2> &m2) { copy(m1.M, m2.M); } - + template <typename MAT1, typename MAT2> inline void copy(const mpi_distributed_matrix<MAT1> &m1, MAT2 &m2) { copy(m1.M, m2); } template <typename MAT1, typename MAT2> inline void copy(const mpi_distributed_matrix<MAT1> &m1, const MAT2 &m2) { copy(m1.M, m2); } - + template <typename MATSP, typename V1, typename V2> inline typename strongest_value_type3<V1,V2,MATSP>::value_type vect_sp(const mpi_distributed_matrix<MATSP> &ps, const V1 &v1, - const V2 &v2) { + const V2 &v2) { typedef typename strongest_value_type3<V1,V2,MATSP>::value_type T; T res = vect_sp(ps.M, v1, v2), rest; MPI_Allreduce(&res, &rest, 1, mpi_type(T()), MPI_SUM,MPI_COMM_WORLD); @@ -1019,7 +1019,7 @@ namespace gmm { template <typename MAT, typename V1, typename V2> inline void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - V2 &v2) { + V2 &v2) { typedef typename linalg_traits<V2>::value_type T; std::vector<T> v3(vect_size(v2)), v4(vect_size(v2)); static double tmult_tot = 0.0; @@ -1029,7 +1029,7 @@ namespace gmm { if (is_sparse(v2)) GMM_WARNING2("Using a plain temporary, here."); double t_ref2 = MPI_Wtime(); MPI_Allreduce(&(v3[0]), &(v4[0]),gmm::vect_size(v2), mpi_type(T()), - MPI_SUM,MPI_COMM_WORLD); + MPI_SUM,MPI_COMM_WORLD); tmult_tot2 = MPI_Wtime()-t_ref2; cout << "reduce mult mpi = " << tmult_tot2 << endl; gmm::add(v4, v2); @@ -1039,59 +1039,59 @@ namespace gmm { template <typename MAT, typename V1, typename V2> void mult_add(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - const V2 &v2_) + const V2 &v2_) { mult_add(m, v1, const_cast<V2 &>(v2_)); } template <typename MAT, typename V1, typename V2> inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - const V2 &v2_) + const V2 &v2_) { V2 &v2 = const_cast<V2 &>(v2_); clear(v2); mult_add(m, v1, v2); } template <typename MAT, typename V1, typename V2> inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - V2 &v2) + V2 &v2) { clear(v2); mult_add(m, v1, v2); } template <typename MAT, typename V1, typename V2, typename V3> inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - const V2 &v2, const V3 &v3_) + const V2 &v2, const V3 &v3_) { V3 &v3 = const_cast<V3 &>(v3_); gmm::copy(v2, v3); mult_add(m, v1, v3); } template <typename MAT, typename V1, typename V2, typename V3> inline void mult(const mpi_distributed_matrix<MAT> &m, const V1 &v1, - const V2 &v2, V3 &v3) + const V2 &v2, V3 &v3) { gmm::copy(v2, v3); mult_add(m, v1, v3); } - + template <typename MAT> inline - size_type mat_nrows(const mpi_distributed_matrix<MAT> &M) + size_type mat_nrows(const mpi_distributed_matrix<MAT> &M) { return mat_nrows(M.M); } template <typename MAT> inline - size_type mat_ncols(const mpi_distributed_matrix<MAT> &M) + size_type mat_ncols(const mpi_distributed_matrix<MAT> &M) { return mat_nrows(M.M); } template <typename MAT> inline void resize(mpi_distributed_matrix<MAT> &M, size_type m, size_type n) { resize(M.M, m, n); } template <typename MAT> inline void clear(mpi_distributed_matrix<MAT> &M) { clear(M.M); } - + // For compute reduced system template <typename MAT1, typename MAT2> inline void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2, - mpi_distributed_matrix<MAT2> &M3) + mpi_distributed_matrix<MAT2> &M3) { mult(M1, M2.M, M3.M); } template <typename MAT1, typename MAT2> inline void mult(const mpi_distributed_matrix<MAT2> &M2, - const MAT1 &M1, mpi_distributed_matrix<MAT2> &M3) + const MAT1 &M1, mpi_distributed_matrix<MAT2> &M3) { mult(M2.M, M1, M3.M); } template <typename MAT1, typename MAT2, typename MAT3> inline void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2, - MAT3 &M3) + MAT3 &M3) { mult(M1, M2.M, M3); } template <typename MAT1, typename MAT2, typename MAT3> inline void mult(const MAT1 &M1, const mpi_distributed_matrix<MAT2> &M2, - const MAT3 &M3) + const MAT3 &M3) { mult(M1, M2.M, M3); } template <typename M, typename SUBI1, typename SUBI2> @@ -1112,16 +1112,16 @@ namespace gmm { template <typename MAT, typename SUBI1, typename SUBI2> inline typename select_return<typename sub_matrix_type<const MAT *, SUBI1, SUBI2> ::matrix_type, typename sub_matrix_type<MAT *, SUBI1, SUBI2>::matrix_type, - const MAT *>::return_type + const MAT *>::return_type sub_matrix(const mpi_distributed_matrix<MAT> &m, const SUBI1 &si1, - const SUBI2 &si2) + const SUBI2 &si2) { return sub_matrix(m.M, si1, si2); } template <typename M, typename SUBI1> inline typename select_return<typename sub_matrix_type<const M *, SUBI1, SUBI1> ::matrix_type, typename sub_matrix_type<M *, SUBI1, SUBI1>::matrix_type, M *>::return_type - sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1) + sub_matrix(mpi_distributed_matrix<M> &m, const SUBI1 &si1) { return sub_matrix(m.M, si1, si1); } template <typename M, typename SUBI1> inline @@ -1132,11 +1132,11 @@ namespace gmm { { return sub_matrix(m.M, si1, si1); } - template <typename L> struct transposed_return<const mpi_distributed_matrix<L> *> + template <typename L> struct transposed_return<const mpi_distributed_matrix<L> *> { typedef abstract_null_type return_type; }; - template <typename L> struct transposed_return<mpi_distributed_matrix<L> *> + template <typename L> struct transposed_return<mpi_distributed_matrix<L> *> { typedef abstract_null_type return_type; }; - + template <typename L> inline typename transposed_return<const L *>::return_type transposed(const mpi_distributed_matrix<L> &l) { return transposed(l.M); } @@ -1185,10 +1185,10 @@ namespace std { template <typename T> void swap(gmm::dense_matrix<T> &m1, gmm::dense_matrix<T> &m2) { m1.swap(m2); } - template <typename T, int shift> void + template <typename T, int shift> void swap(gmm::csc_matrix<T,shift> &m1, gmm::csc_matrix<T,shift> &m2) { m1.swap(m2); } - template <typename T, int shift> void + template <typename T, int shift> void swap(gmm::csr_matrix<T,shift> &m1, gmm::csr_matrix<T,shift> &m2) { m1.swap(m2); } }