[GiNaC-devel] Speeding up simplify_indexed(): 2nd try
Chris Dams
Chris.Dams at mi.infn.it
Mon Oct 3 12:10:51 CEST 2005
Dear all,
Here is a new attempt at the patch. The problem with the previous one was
that if one symmetrizes only amongst indices of precisely the same type,
also renaming of indices should be done taking these types into account.
Another advantage of this is that if a user uses for instance mu and nu
for indices with a variance, these will no longer end up being used as
indices without variance. I think this makes output more readable.
Best,
Chris
-------------- next part --------------
Index: ginac/indexed.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/indexed.cpp,v
retrieving revision 1.96
diff -c -r1.96 indexed.cpp
*** ginac/indexed.cpp 19 May 2005 14:10:40 -0000 1.96
--- ginac/indexed.cpp 3 Oct 2005 09:33:21 -0000
***************
*** 532,537 ****
--- 532,546 ----
return f.get_free_indices();
}
+ template<class T> size_t number_of_type(const exvector&v)
+ {
+ size_t number = 0;
+ for(exvector::const_iterator i=v.begin(); i!=v.end(); ++i)
+ if(is_exactly_a<T>(*i))
+ ++number;
+ return number;
+ }
+
/** Rename dummy indices in an expression.
*
* @param e Expression to work on
***************
*** 540,549 ****
* @param global_dummy_indices The set of dummy indices that have appeared
* before and which we would like to use in "e", too. This gets updated
* by the function */
! static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
{
! size_t global_size = global_dummy_indices.size(),
! local_size = local_dummy_indices.size();
// Any local dummy indices at all?
if (local_size == 0)
--- 549,558 ----
* @param global_dummy_indices The set of dummy indices that have appeared
* before and which we would like to use in "e", too. This gets updated
* by the function */
! template<class T> static ex rename_dummy_indices(const ex & e, exvector & global_dummy_indices, exvector & local_dummy_indices)
{
! size_t global_size = number_of_type<T>(global_dummy_indices),
! local_size = number_of_type<T>(local_dummy_indices);
// Any local dummy indices at all?
if (local_size == 0)
***************
*** 557,563 ****
int remaining = local_size - global_size;
exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
while (it != itend && remaining > 0) {
! if (find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(op0_is_equal(), *it)) == global_dummy_indices.end()) {
global_dummy_indices.push_back(*it);
global_size++;
remaining--;
--- 566,572 ----
int remaining = local_size - global_size;
exvector::const_iterator it = local_dummy_indices.begin(), itend = local_dummy_indices.end();
while (it != itend && remaining > 0) {
! if (is_a<T>(*it) && find_if(global_dummy_indices.begin(), global_dummy_indices.end(), bind2nd(idx_is_equal_ignore_dim(), *it)) == global_dummy_indices.end()) {
global_dummy_indices.push_back(*it);
global_size++;
remaining--;
***************
*** 575,585 ****
exvector local_syms, global_syms;
local_syms.reserve(local_size);
global_syms.reserve(local_size);
! for (size_t i=0; i<local_size; i++)
! local_syms.push_back(local_dummy_indices[i].op(0));
shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
! for (size_t i=0; i<local_size; i++) // don't use more global symbols than necessary
! global_syms.push_back(global_dummy_indices[i].op(0));
shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
// Remove common indices
--- 584,596 ----
exvector local_syms, global_syms;
local_syms.reserve(local_size);
global_syms.reserve(local_size);
! for (size_t i=0; local_syms.size()!=local_size; i++)
! if(is_exactly_a<T>(local_dummy_indices[i]))
! local_syms.push_back(local_dummy_indices[i].op(0));
shaker_sort(local_syms.begin(), local_syms.end(), ex_is_less(), ex_swap());
! for (size_t i=0; global_syms.size()!=local_size; i++) // don't use more global symbols than necessary
! if(is_exactly_a<T>(global_dummy_indices[i]))
! global_syms.push_back(global_dummy_indices[i].op(0));
shaker_sort(global_syms.begin(), global_syms.end(), ex_is_less(), ex_swap());
// Remove common indices
***************
*** 704,709 ****
--- 715,732 ----
}
}
+ template<class T> ex idx_symmetrization(const ex& r,const exvector& local_dummy_indices)
+ { exvector dummy_syms;
+ dummy_syms.reserve(r.nops());
+ for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
+ if(is_exactly_a<T>(*it))
+ dummy_syms.push_back(it->op(0));
+ if(dummy_syms.size() < 2)
+ return r;
+ ex q=symmetrize(r, dummy_syms);
+ return q;
+ }
+
/** Simplify product of indexed expressions (commutative, noncommutative and
* simple squares), return list of free indices. */
ex simplify_indexed_product(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
***************
*** 864,882 ****
// The result should be symmetric with respect to exchange of dummy
// indices, so if the symmetrization vanishes, the whole expression is
// zero. This detects things like eps.i.j.k * p.j * p.k = 0.
! if (local_dummy_indices.size() >= 2) {
! exvector dummy_syms;
! dummy_syms.reserve(local_dummy_indices.size());
! for (exvector::const_iterator it = local_dummy_indices.begin(); it != local_dummy_indices.end(); ++it)
! dummy_syms.push_back(it->op(0));
! if (symmetrize(r, dummy_syms).is_zero()) {
! free_indices.clear();
! return _ex0;
! }
}
// Dummy index renaming
! r = rename_dummy_indices(r, dummy_indices, local_dummy_indices);
// Product of indexed object with a scalar?
if (is_exactly_a<mul>(r) && r.nops() == 2
--- 887,912 ----
// The result should be symmetric with respect to exchange of dummy
// indices, so if the symmetrization vanishes, the whole expression is
// zero. This detects things like eps.i.j.k * p.j * p.k = 0.
! ex q = idx_symmetrization<idx>(r, local_dummy_indices);
! if (q.is_zero()) {
! free_indices.clear();
! return _ex0;
! }
! q = idx_symmetrization<varidx>(q, local_dummy_indices);
! if (q.is_zero()) {
! free_indices.clear();
! return _ex0;
! }
! q = idx_symmetrization<spinidx>(q, local_dummy_indices);
! if (q.is_zero()) {
! free_indices.clear();
! return _ex0;
}
// Dummy index renaming
! r = rename_dummy_indices<idx>(r, dummy_indices, local_dummy_indices);
! r = rename_dummy_indices<varidx>(r, dummy_indices, local_dummy_indices);
! r = rename_dummy_indices<spinidx>(r, dummy_indices, local_dummy_indices);
// Product of indexed object with a scalar?
if (is_exactly_a<mul>(r) && r.nops() == 2
***************
*** 943,948 ****
--- 973,989 ----
}
};
+ bool hasindex(const ex &x, const ex &sym)
+ {
+ if(is_a<idx>(x) && x.op(0)==sym)
+ return true;
+ else
+ for(size_t i=0; i<x.nops(); ++i)
+ if(hasindex(x.op(i), sym))
+ return true;
+ return false;
+ }
+
/** Simplify indexed expression, return list of free indices. */
ex simplify_indexed(const ex & e, exvector & free_indices, exvector & dummy_indices, const scalar_products & sp)
{
***************
*** 971,977 ****
}
// Rename the dummy indices
! return rename_dummy_indices(e_expanded, dummy_indices, local_dummy_indices);
}
// Simplification of sum = sum of simplifications, check consistency of
--- 1012,1021 ----
}
// Rename the dummy indices
! e_expanded = rename_dummy_indices<idx>(e_expanded, dummy_indices, local_dummy_indices);
! e_expanded = rename_dummy_indices<varidx>(e_expanded, dummy_indices, local_dummy_indices);
! e_expanded = rename_dummy_indices<spinidx>(e_expanded, dummy_indices, local_dummy_indices);
! return e_expanded;
}
// Simplification of sum = sum of simplifications, check consistency of
***************
*** 1015,1032 ****
if (num_terms_orig < 2 || dummy_indices.size() < 2)
return sum;
- // Yes, construct vector of all dummy index symbols
- exvector dummy_syms;
- dummy_syms.reserve(dummy_indices.size());
- for (exvector::const_iterator it = dummy_indices.begin(); it != dummy_indices.end(); ++it)
- dummy_syms.push_back(it->op(0));
-
// Chop the sum into terms and symmetrize each one over the dummy
// indices
std::vector<terminfo> terms;
for (size_t i=0; i<sum.nops(); i++) {
const ex & term = sum.op(i);
! ex term_symm = symmetrize(term, dummy_syms);
if (term_symm.is_zero())
continue;
terms.push_back(terminfo(term, term_symm));
--- 1059,1077 ----
if (num_terms_orig < 2 || dummy_indices.size() < 2)
return sum;
// Chop the sum into terms and symmetrize each one over the dummy
// indices
std::vector<terminfo> terms;
for (size_t i=0; i<sum.nops(); i++) {
const ex & term = sum.op(i);
! exvector dummy_indices_of_term;
! dummy_indices_of_term.reserve(dummy_indices.size());
! for(exvector::iterator i=dummy_indices.begin(); i!=dummy_indices.end(); ++i)
! if(hasindex(term,i->op(0)))
! dummy_indices_of_term.push_back(*i);
! ex term_symm = idx_symmetrization<idx>(term, dummy_indices_of_term);
! term_symm = idx_symmetrization<varidx>(term_symm, dummy_indices_of_term);
! term_symm = idx_symmetrization<spinidx>(term_symm, dummy_indices_of_term);
if (term_symm.is_zero())
continue;
terms.push_back(terminfo(term, term_symm));
More information about the GiNaC-devel
mailing list