sums of indexed matrices are now possible
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Thu, 22 Mar 2001 18:52:08 +0000 (18:52 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Thu, 22 Mar 2001 18:52:08 +0000 (18:52 +0000)
doc/tutorial/ginac.texi
ginac/basic.cpp
ginac/basic.h
ginac/indexed.cpp
ginac/matrix.cpp
ginac/matrix.h

index fe09793..ce36a1b 100644 (file)
@@ -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
index 9fc4fa3..6415f33 100644 (file)
@@ -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
index ee2bef2..4277d94 100644 (file)
@@ -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;
index 4fa4256..e23709d 100644 (file)
@@ -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; i<e_expanded.nops(); i++) {
+               for (unsigned i=1; i<e_expanded.nops(); i++) {
                        exvector free_indices_of_term;
-                       sum += simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
-                       if (i == 0)
-                               free_indices = free_indices_of_term;
-                       else if (!indices_consistent(free_indices, free_indices_of_term))
+                       ex term = simplify_indexed(e_expanded.op(i), free_indices_of_term, sp);
+                       if (!indices_consistent(free_indices, free_indices_of_term))
                                throw (std::runtime_error("simplify_indexed: inconsistent indices in sum"));
+                       if (is_ex_of_type(sum, indexed) && is_ex_of_type(term, indexed))
+                               sum = sum.op(0).bp->add_indexed(sum, term);
+                       else
+                               sum += term;
                }
 
                return sum;
index 1a5d695..f34534a 100644 (file)
@@ -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
 {
index 8c24dcf..53afeb2 100644 (file)
@@ -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; };