From 9d8d0f4924171b0acb7ec6a333897fe1ec545707 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Thu, 22 Mar 2001 18:52:08 +0000 Subject: [PATCH] sums of indexed matrices are now possible --- doc/tutorial/ginac.texi | 15 ++++++++------- ginac/basic.cpp | 19 +++++++++++++++++-- ginac/basic.h | 1 + ginac/indexed.cpp | 14 ++++++++------ ginac/matrix.cpp | 35 +++++++++++++++++++++++++++++++++++ ginac/matrix.h | 1 + 6 files changed, 70 insertions(+), 15 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index fe09793d..ce36a1b5 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -1724,7 +1724,8 @@ tensor). @subsection Linear algebra The @code{matrix} class can be used with indices to do some simple linear -algebra (products of vectors and matrices, traces and scalar products): +algebra (sums and products of vectors and matrices, traces and scalar +products): @example @{ @@ -1743,20 +1744,20 @@ algebra (products of vectors and matrices, traces and scalar products): cout << e.simplify_indexed() << endl; // -> [[ [[2*y+x]], [[4*y+3*x]] ]].i - e = indexed(A, i, j) * indexed(X, i); + e = indexed(A, i, j) * indexed(X, i) + indexed(X, j); cout << e.simplify_indexed() << endl; - // -> [[ [[3*y+x,4*y+2*x]] ]].j + // -> [[ [[3*y+2*x,5*y+2*x]] ]].j @} @end example -You can of course obtain the same results with the @code{matrix::mul()} -and @code{matrix::trace()} methods but with indices you don't have to -worry about transposing matrices. +You can of course obtain the same results with the @code{matrix::add()}, +@code{matrix::mul()} and @code{matrix::trace()} methods but with indices you +don't have to worry about transposing matrices. Matrix indices always start at 0 and their dimension must match the number of rows/columns of the matrix. Matrices with one row or one column are vectors and can have one or two indices (it doesn't matter whether it's a -row or a columnt vector). Other matrices must have two indices. +row or a column vector). Other matrices must have two indices. You should be careful when using indices with variance on matrices. GiNaC doesn't look at the variance and doesn't know that @samp{F~mu~nu} and diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 9fc4fa30..6415f33f 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -297,16 +297,31 @@ ex basic::eval_indexed(const basic & i) const return i.hold(); } +/** Add two indexed expressions. They are guaranteed to be of class indexed + * (or a subclass) and their indices are compatible. This function is used + * internally by simplify_indexed(). + * + * @param self First indexed expression; it's base object is *this + * @param other Second indexed expression + * @return sum of self and other + * @see ex::simplify_indexed() */ +ex basic::add_indexed(const ex & self, const ex & other) const +{ + return self + other; +} + /** Try to contract two indexed expressions that appear in the same product. * If a contraction exists, the function overwrites one or both of the * expressions and returns true. Otherwise it returns false. It is * guaranteed that both expressions are of class indexed (or a subclass) - * and that at least one dummy index has been found. + * and that at least one dummy index has been found. This functions is + * used internally by simplify_indexed(). * * @param self Pointer to first indexed expression; it's base object is *this * @param other Pointer to second indexed expression * @param v The complete vector of factors - * @return true if the contraction was successful, false otherwise */ + * @return true if the contraction was successful, false otherwise + * @see ex::simplify_indexed() */ bool basic::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { // Do nothing diff --git a/ginac/basic.h b/ginac/basic.h index ee2bef23..4277d949 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -130,6 +130,7 @@ public: // only const functions please (may break reference counting) virtual exvector get_free_indices(void) const; virtual ex simplify_ncmul(const exvector & v) const; virtual ex eval_indexed(const basic & i) const; + virtual ex add_indexed(const ex & self, const ex & other) const; virtual bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; protected: // non-const functions should be called from class ex only virtual ex derivative(const symbol & s) const; diff --git a/ginac/indexed.cpp b/ginac/indexed.cpp index 4fa4256a..e23709d5 100644 --- a/ginac/indexed.cpp +++ b/ginac/indexed.cpp @@ -705,15 +705,17 @@ ex simplify_indexed(const ex & e, exvector & free_indices, const scalar_products // Simplification of sum = sum of simplifications, check consistency of // free indices in each term if (is_ex_exactly_of_type(e_expanded, add)) { - ex sum = _ex0(); + ex sum = simplify_indexed(e_expanded.op(0), free_indices, sp); - for (unsigned i=0; iadd_indexed(sum, term); + else + sum += term; } return sum; diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 1a5d6953..f34534af 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -381,6 +381,41 @@ ex matrix::eval_indexed(const basic & i) const return i.hold(); } +/** Sum of two indexed matrices. */ +ex matrix::add_indexed(const ex & self, const ex & other) const +{ + GINAC_ASSERT(is_ex_of_type(self, indexed)); + GINAC_ASSERT(is_ex_of_type(other, indexed)); + GINAC_ASSERT(self.nops() == 2 || self.nops() == 3); + + // Only add two matrices + if (is_ex_of_type(other.op(0), matrix)) { + GINAC_ASSERT(other.nops() == 2 || other.nops() == 3); + + const matrix &self_matrix = ex_to_matrix(self.op(0)); + const matrix &other_matrix = ex_to_matrix(other.op(0)); + + if (self.nops() == 2 && other.nops() == 2) { // vector + vector + + if (self_matrix.row == other_matrix.row) + return indexed(self_matrix.add(other_matrix), self.op(1)); + else if (self_matrix.row == other_matrix.col) + return indexed(self_matrix.add(other_matrix.transpose()), self.op(1)); + + } else if (self.nops() == 3 && other.nops() == 3) { // matrix + matrix + + if (self.op(1).is_equal(other.op(1)) && self.op(2).is_equal(other.op(2))) + return indexed(self_matrix.add(other_matrix), self.op(1), self.op(2)); + else if (self.op(1).is_equal(other.op(2)) && self.op(2).is_equal(other.op(1))) + return indexed(self_matrix.add(other_matrix.transpose()), self.op(1), self.op(2)); + + } + } + + // Don't know what to do, return unevaluated sum + return self + other; +} + /** Contraction of an indexed matrix with something else. */ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const { diff --git a/ginac/matrix.h b/ginac/matrix.h index 8c24dcf9..53afeb2e 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -52,6 +52,7 @@ public: ex evalf(int level=0) const; ex subs(const lst & ls, const lst & lr) const; ex eval_indexed(const basic & i) const; + ex add_indexed(const ex & self, const ex & other) const; bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const; protected: unsigned return_type(void) const { return return_types::noncommutative; }; -- 2.44.0