]> www.ginac.de Git - ginac.git/commitdiff
 author Christian Bauer Thu, 22 Mar 2001 18:52:08 +0000 (18:52 +0000) committer Christian Bauer Thu, 22 Mar 2001 18:52:08 +0000 (18:52 +0000)
 doc/tutorial/ginac.texi patch | blob | history ginac/basic.cpp patch | blob | history ginac/basic.h patch | blob | history ginac/indexed.cpp patch | blob | history ginac/matrix.cpp patch | blob | history ginac/matrix.h patch | blob | history

index fe09793d0f8b91e24b5f867d004be03c8870dc98..ce36a1b5e93a630be2657e294b13831cbaf8a505 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
+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 9fc4fa308526fdecb00235ab6f528210034759e0..6415f33fde5a79456ebb97ba0e81f36139110c97 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
@@ -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 4fa4256abfb286b86c95d427efa9a5c8de27c4c8..e23709d5f506647280eca2f54fb528f807a5dc8f 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
-               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))
+                       else
+                               sum += term;
}

return sum;
index 1a5d6953c3d994365d9d6b05a60ebf351119e7bd..f34534afdb1ab3ee39b567bd11c7ed7408010855 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)
+                       else if (self_matrix.row == other_matrix.col)
+
+               } 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)))
+                       else if (self.op(1).is_equal(other.op(2)) && self.op(2).is_equal(other.op(1)))
+
+               }
+       }
+
+       // 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 8c24dcf9043a423c21521841928327eab22586a4..53afeb2e6e05c8ab4448ec4f746c3bd39775db82 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; };