[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: |
Tue, 19 Feb 2019 17:15:14 -0500 (EST) |
branch: master
commit f374fd75eabd2596d3a2fd766ec08932c3307e40
Author: Konstantinos Poulios <address@hidden>
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); }
}